From 2bcf5ae8cd8fc5058d0e039b84c086718d47bc1b Mon Sep 17 00:00:00 2001
From: aaTman <mandelbaum.taylor@gmail.com>
Date: Tue, 7 Apr 2026 21:25:33 +0000
Subject: [PATCH 01/12] update heat and cold bounds, update docstrings and
 methodology. add base_temp_events.yaml for reproducibility

---
 data_prep/heat_cold_bounds_case.py            | 808 +++++++++++++++---
 data_prep/heat_cold_bounds_global.py          | 344 ++++++--
 data_prep/plot_temperature_events.py          | 205 ++++-
 .../data/base_temp_events.yaml                | 661 ++++++++++++++
 src/extremeweatherbench/data/events.yaml      | 456 +++++-----
 5 files changed, 2028 insertions(+), 446 deletions(-)
 create mode 100644 src/extremeweatherbench/data/base_temp_events.yaml

diff --git a/data_prep/heat_cold_bounds_case.py b/data_prep/heat_cold_bounds_case.py
index 3f398331..74153e82 100644
--- a/data_prep/heat_cold_bounds_case.py
+++ b/data_prep/heat_cold_bounds_case.py
@@ -1,10 +1,18 @@
 #!/usr/bin/env python3
-"""Validate and expand heat wave / cold snap bounding boxes from events.yaml.
-
-For each heat_wave or cold snap event in events.yaml, iteratively
-grows the bounding box by 2 degrees in each direction until fewer
-than 50% of grid points on an edge exceed the climatological
-threshold (or 10 iterations).
+"""Validate and expand heat wave / freeze bounding boxes.
+
+Reads heat_wave and freeze cases from base_temp_events.yaml (which uses
+centered_region),
+then:
+  1. Uses the centered_region center as the initial box.
+  2. Grows the time window forward from start_date-3 one day at a
+     time until the single-day exceedance fraction drops below 50 %
+     of its peak.
+  3. Iteratively grows the spatial box by 2 degrees in each direction
+     until fewer than 50 % of grid points on each edge exceed the
+     climatological threshold (or 10 iterations).
+  4. Writes the final bounds as bounded_region entries into
+     events.yaml.
 
 Usage:
     python heat_cold_bounds_case.py \\
@@ -12,12 +20,18 @@
 """
 
 import argparse
+import importlib
 import logging
 import pathlib
 import time as time_module
 from typing import Dict, List, Optional
 
+import cartopy.crs as ccrs
+import cartopy.feature as cfeature
 import joblib
+import matplotlib.colors as mcolors
+import matplotlib.patches as mpatches
+import matplotlib.pyplot as plt
 import numpy as np
 import pandas as pd
 import regionmask
@@ -25,12 +39,15 @@
 import xarray as xr
 from dask.distributed import Client, LocalCluster
 from plot_temperature_events import (
+    VALID_QUANTILES,
     detect_time_dim,
     max_consecutive_days,
     plot_consecutive_map,
+    resolve_op,
 )
+from ruamel.yaml import YAML
 
-from extremeweatherbench import cases, defaults, inputs
+from extremeweatherbench import cases, defaults, inputs, regions
 
 logging.basicConfig(
     level=logging.INFO,
@@ -39,62 +56,245 @@
 logger = logging.getLogger(__name__)
 
 MIN_CONSECUTIVE_DAYS = 3
-EXPANSION_DEGREES = 2.0
-MAX_ITERATIONS = 10
+EXPANSION_DEGREES = 1
+MAX_ITERATIONS = 20
 EDGE_VALIDITY_THRESHOLD = 0.5
 MIN_GRIDPOINTS = 500
+TEMPORAL_LOAD_BUFFER_DAYS = 14
+MAX_TEMPORAL_DAYS = 21
 
 
 def _apply_consecutive_filter(
     mask: np.ndarray,
     min_days: int = MIN_CONSECUTIVE_DAYS,
+    max_grace_days: int = 1,
 ) -> np.ndarray:
-    """Keep only runs of ``min_days``+ True days (axis 0)."""
+    """Keep runs of ``min_days``+ True days along axis 0.
+
+    After ``min_days`` strict consecutive True days are established,
+    gaps of up to ``max_grace_days`` are bridged so the event can
+    continue. Runs that never reach ``min_days`` strict consecutive
+    True days are discarded.
+
+    Args:
+        mask: Boolean array of shape (time, lat, lon).
+        min_days: Minimum run length required to qualify as an event.
+        max_grace_days: Maximum gap length to bridge after the
+            minimum run is established.
+
+    Returns:
+        Boolean array of the same shape with only qualifying runs
+        retained.
+    """
     struct = np.zeros((min_days, 1, 1), dtype=bool)
     struct[:, 0, 0] = True
-    eroded = ndimage.binary_erosion(
-        mask,
-        structure=struct,
-        border_value=False,
+
+    # Strict 3-day runs (no grace)
+    strict = (
+        ndimage.binary_dilation(
+            ndimage.binary_erosion(
+                mask, structure=struct, border_value=False,
+            ),
+            structure=struct, border_value=False,
+        )
+        & mask
     )
-    dilated = ndimage.binary_dilation(
-        eroded,
-        structure=struct,
-        border_value=False,
+
+    if max_grace_days <= 0:
+        return strict
+
+    # Fill gaps of ≤ max_grace_days
+    close_k = np.zeros((2 * max_grace_days + 1, 1, 1), dtype=bool)
+    close_k[:, 0, 0] = True
+    filled = ndimage.binary_closing(
+        mask, structure=close_k, border_value=False,
+    )
+
+    # Runs of min_days+ in the gap-filled mask
+    filled_runs = (
+        ndimage.binary_dilation(
+            ndimage.binary_erosion(
+                filled, structure=struct, border_value=False,
+            ),
+            structure=struct, border_value=False,
+        )
+        & filled
     )
-    return dilated & mask
+
+    # Label each temporal run independently per grid point
+    lbl_struct = np.zeros((3, 3, 3), dtype=int)
+    lbl_struct[0, 1, 1] = 1
+    lbl_struct[1, 1, 1] = 1
+    lbl_struct[2, 1, 1] = 1
+    labels, _ = ndimage.label(filled_runs, structure=lbl_struct)
+
+    # Keep only runs that contain a strict 3-day block
+    valid = np.unique(labels[strict & (labels > 0)])
+    return np.isin(labels, valid) & filled_runs
 
 
 def _edge_valid_fraction(
-    filtered: np.ndarray,
+    mask_2d: np.ndarray,
+    land_2d: np.ndarray,
     edge: str,
     band_pts: int,
 ) -> float:
-    """Fraction of grid points on *edge* with any valid day."""
+    """Return the fraction of land grid points on an edge that are active.
+
+    Ocean/masked points are excluded from both numerator and denominator
+    so coastal edges aren't penalised.
+
+    Args:
+        mask_2d: 2-D boolean activity array (lat, lon).
+        land_2d: 2-D boolean land mask (lat, lon); True = land.
+        edge: One of "north", "south", "east", "west".
+        band_pts: Width of the edge strip in grid points.
+
+    Returns:
+        Fraction in [0, 1] of land points in the strip that are active,
+        or 0.0 if the strip contains no land points.
+
+    Raises:
+        ValueError: If edge is not one of the recognised values.
+    """
     if edge == "north":
-        strip = filtered[:, -band_pts:, :]
+        strip = mask_2d[-band_pts:, :]
+        land_strip = land_2d[-band_pts:, :]
     elif edge == "south":
-        strip = filtered[:, :band_pts, :]
+        strip = mask_2d[:band_pts, :]
+        land_strip = land_2d[:band_pts, :]
     elif edge == "west":
-        strip = filtered[:, :, :band_pts]
+        strip = mask_2d[:, :band_pts]
+        land_strip = land_2d[:, :band_pts]
     elif edge == "east":
-        strip = filtered[:, :, -band_pts:]
+        strip = mask_2d[:, -band_pts:]
+        land_strip = land_2d[:, -band_pts:]
     else:
         raise ValueError(f"Unknown edge: {edge}")
 
-    has_event = strip.any(axis=0)
-    if has_event.size == 0:
+    n_land = int(land_strip.sum())
+    if n_land == 0:
         return 0.0
-    return float(has_event.mean())
+    return float((strip & land_strip).sum()) / n_land
+
+
+def _to_plot_lon(lon: float) -> float:
+    """Convert 0-360 longitude to -180..180 for plotting.
+
+    Args:
+        lon: Longitude in 0-360 degrees.
+
+    Returns:
+        Longitude in -180..180 degrees.
+    """
+    return lon - 360 if lon > 180 else lon
+
+
+def plot_peak_day_with_bounds(
+    peak_mask: np.ndarray,
+    all_lats: np.ndarray,
+    all_lons: np.ndarray,
+    initial_box: tuple,
+    final_box: tuple,
+    event_type: str,
+    title: str,
+    output_path: str,
+) -> None:
+    """Plot the peak-footprint day with initial and final bounding boxes.
+
+    Args:
+        peak_mask: 2-D boolean array (lat, lon) for the peak day.
+        all_lats: 1-D latitude array matching peak_mask rows.
+        all_lons: 1-D longitude array (0-360) matching peak_mask cols.
+        initial_box: (lat_min, lat_max, lon_min, lon_max) of the
+            initial centered region.
+        final_box: (lat_min, lat_max, lon_min, lon_max) after spatial
+            expansion.
+        event_type: ``"heat_wave"`` or ``"freeze"``.
+        title: Figure title.
+        output_path: Destination PNG file path.
+    """
+    plot_lons = np.array([_to_plot_lon(x) for x in all_lons])
+    lon_min = float(plot_lons.min()) - 1
+    lon_max = float(plot_lons.max()) + 1
+    lat_min = float(all_lats.min()) - 1
+    lat_max = float(all_lats.max()) + 1
+
+    cmap = mcolors.ListedColormap(
+        ["white", "firebrick" if event_type == "heat_wave" else "steelblue"]
+    )
+    norm = mcolors.BoundaryNorm([0, 0.5, 1], cmap.N)
+
+    fig, ax = plt.subplots(
+        subplot_kw={"projection": ccrs.PlateCarree()},
+        figsize=(10, 8),
+    )
+    ax.set_extent(
+        [lon_min, lon_max, lat_min, lat_max],
+        crs=ccrs.PlateCarree(),
+    )
+    ax.add_feature(cfeature.OCEAN, facecolor="lightblue", zorder=0)
+
+    ax.pcolormesh(
+        plot_lons,
+        all_lats,
+        peak_mask.astype(float),
+        cmap=cmap,
+        norm=norm,
+        transform=ccrs.PlateCarree(),
+        shading="auto",
+    )
+    ax.add_feature(cfeature.COASTLINE, linewidth=0.5)
+    ax.add_feature(cfeature.BORDERS, linewidth=0.3)
+
+    def _add_rect(box, color, label, ls="-"):
+        blat_min, blat_max, blon_min, blon_max = box
+        blon_min = _to_plot_lon(blon_min)
+        blon_max = _to_plot_lon(blon_max)
+        rect = mpatches.Rectangle(
+            (blon_min, blat_min),
+            blon_max - blon_min,
+            blat_max - blat_min,
+            linewidth=2,
+            edgecolor=color,
+            facecolor="none",
+            linestyle=ls,
+            transform=ccrs.PlateCarree(),
+            label=label,
+        )
+        ax.add_patch(rect)
+
+    _add_rect(initial_box, "blue", "Initial box", ls="--")
+    _add_rect(final_box, "green", "Final bounds", ls="-")
+
+    ax.legend(loc="lower left", fontsize=9)
+    ax.set_title(title, loc="left", fontsize=12)
+    fig.savefig(output_path, dpi=200, bbox_inches="tight")
+    plt.close(fig)
+    logger.info("  Saved peak-day plot: %s", output_path)
 
 
 def process_event(
     single_case: cases.IndividualCase,
+    out_dir: pathlib.Path = pathlib.Path("."),
+    quantile: float | None = None,
+    op_str: str | None = None,
 ) -> Optional[Dict]:
-    """Process one event: compute mask and iteratively expand.
-
-    Opens ERA5 and climatology inside the worker to avoid
-    serialisation of lazy zarr handles across joblib processes.
+    """Process one event: find the time window then expand bounds.
+
+    Args:
+        single_case: The individual case to process; must have a
+            CenteredRegion location.
+        out_dir: Directory in which plots are saved.
+        quantile: Climatology quantile. Defaults to 0.85 for
+            heat_wave, 0.15 for freeze.
+        op_str: Comparison operator string (e.g. ">", ">=", "<",
+            "<="). Defaults to ">" for heat_wave, "<" for freeze.
+
+    Returns:
+        A dict with keys event_type, start_date, end_date,
+        latitude_min/max, longitude_min/max, case_id, title, and
+        _peak_gridpoints, or None if the case was skipped.
     """
     is_heatwave = single_case.event_type == "heat_wave"
     logger.info(
@@ -104,22 +304,60 @@ def process_event(
         single_case.event_type,
     )
 
-    bounds = single_case.location.as_geopandas().total_bounds
-    lon_min, lat_min, lon_max, lat_max = bounds
+    loc = single_case.location
+    if not isinstance(loc, regions.CenteredRegion):
+        logger.warning(
+            "  Case %d location is %s, not CenteredRegion — skipping",
+            single_case.case_id_number,
+            type(loc).__name__,
+        )
+        return None
+
+    center_lat = loc.latitude
+    center_lon = loc.longitude
+    if isinstance(loc.bounding_box_degrees, tuple):
+        half_lat = loc.bounding_box_degrees[0] / 2.0
+        half_lon = loc.bounding_box_degrees[1] / 2.0
+    else:
+        half_lat = loc.bounding_box_degrees / 2.0
+        half_lon = loc.bounding_box_degrees / 2.0
 
+    box_lat_min = center_lat - half_lat
+    box_lat_max = center_lat + half_lat
+    box_lon_min = center_lon - half_lon
+    box_lon_max = center_lon + half_lon
+
+    logger.info(
+        "  Initial box (%.1f°): center (%.2f, %.2f)"
+        "  lat [%.2f, %.2f], lon [%.2f, %.2f]",
+        loc.bounding_box_degrees
+        if not isinstance(loc.bounding_box_degrees, tuple)
+        else loc.bounding_box_degrees[0],
+        center_lat,
+        center_lon,
+        box_lat_min,
+        box_lat_max,
+        box_lon_min,
+        box_lon_max,
+    )
+
+    # Time range and spatial pre-fetch extent
     start_date = pd.Timestamp(single_case.start_date) - pd.Timedelta(days=3)
-    end_date = pd.Timestamp(single_case.end_date) + pd.Timedelta(days=3)
+    end_date = (
+        pd.Timestamp(single_case.end_date)
+        + pd.Timedelta(days=TEMPORAL_LOAD_BUFFER_DAYS)
+    )
 
     pot_lat_min = max(
         -90,
-        lat_min - MAX_ITERATIONS * EXPANSION_DEGREES,
+        box_lat_min - MAX_ITERATIONS * EXPANSION_DEGREES,
     )
     pot_lat_max = min(
         90,
-        lat_max + MAX_ITERATIONS * EXPANSION_DEGREES,
+        box_lat_max + MAX_ITERATIONS * EXPANSION_DEGREES,
     )
-    pot_lon_min = lon_min - MAX_ITERATIONS * EXPANSION_DEGREES
-    pot_lon_max = lon_max + MAX_ITERATIONS * EXPANSION_DEGREES
+    pot_lon_min = box_lon_min - MAX_ITERATIONS * EXPANSION_DEGREES
+    pot_lon_max = box_lon_max + MAX_ITERATIONS * EXPANSION_DEGREES
 
     ds = xr.open_zarr(
         inputs.ARCO_ERA5_FULL_URI,
@@ -133,85 +371,83 @@ def process_event(
     six_hourly = t2m[tdim].dt.hour.isin([0, 6, 12, 18])
     t2m = t2m.sel({tdim: six_hourly}).sortby("latitude")
 
-    # Normalise ERA5 longitudes from 0-360 to -180/180 so that the
-    # case bounding boxes (always in -180/180 from geopandas) can be
-    # used directly for slicing without wrap-around issues.
     if float(t2m.longitude.max()) > 180:
         t2m = t2m.assign_coords(
             longitude=(t2m.longitude.values + 180) % 360 - 180,
         ).sortby("longitude")
 
+    # Convert center lon to match ERA5 convention for slicing
+    pot_lon_min_era = pot_lon_min
+    pot_lon_max_era = pot_lon_max
+    if pot_lon_min_era > 180:
+        pot_lon_min_era -= 360
+    if pot_lon_max_era > 180:
+        pot_lon_max_era -= 360
+
     t2m = t2m.sel(
         latitude=slice(pot_lat_min, pot_lat_max),
-        longitude=slice(pot_lon_min, pot_lon_max),
+        longitude=slice(
+            min(pot_lon_min_era, pot_lon_max_era),
+            max(pot_lon_min_era, pot_lon_max_era),
+        ),
     )
 
     if t2m.latitude.size == 0 or t2m.longitude.size == 0:
         logger.warning(
-            "  Case %d: empty spatial selection"
-            " (lon=[%.2f, %.2f], lat=[%.2f, %.2f]) — skipping",
+            "  Case %d: empty spatial selection — skipping",
             single_case.case_id_number,
-            pot_lon_min,
-            pot_lon_max,
-            pot_lat_min,
-            pot_lat_max,
         )
         return None
 
-    if is_heatwave:
-        daily = t2m.resample({tdim: "1D"}).max()
-        clim = (
-            defaults.get_climatology(0.85)
-            .max(
-                dim="hour",
-            )
-            .sortby("latitude")
-        )
-    else:
-        daily = t2m.resample({tdim: "1D"}).min()
-        clim = (
-            defaults.get_climatology(0.15)
-            .min(
-                dim="hour",
-            )
-            .sortby("latitude")
-        )
+    # 6-hourly exceedance then daily all-pass
+    if quantile is None:
+        quantile = 0.85 if is_heatwave else 0.15
+    if op_str is None:
+        op_str = ">" if is_heatwave else "<"
+    cmp = resolve_op(op_str)
+
+    clim = defaults.get_climatology(quantile).sortby("latitude")
 
-    # Match climatology longitude convention to the (possibly
-    # remapped) ERA5 data so reindex_like aligns correctly.
     if float(clim.longitude.max()) > 180:
         clim = clim.assign_coords(
             longitude=(clim.longitude.values + 180) % 360 - 180,
         ).sortby("longitude")
 
-    doy = daily[tdim].dt.dayofyear
+    doy = t2m[tdim].dt.dayofyear
+    hour = t2m[tdim].dt.hour
     max_clim_doy = int(clim.dayofyear.max())
     doy_capped = doy.clip(max=max_clim_doy)
+
     clim_aligned = clim.sel(
         dayofyear=doy_capped,
-    ).reindex_like(daily, method="nearest")
+        hour=hour,
+    ).reindex_like(t2m, method="nearest")
+
+    exc_6h = cmp(t2m, clim_aligned)
+
+    daily_all_pass = (
+        exc_6h.resample({tdim: "1D"}).min().astype(bool)
+    )
 
     land_reg = regionmask.defined_regions.natural_earth_v5_0_0.land_110
     land_mask = (
         land_reg.mask(
-            daily.longitude,
-            daily.latitude,
+            daily_all_pass.longitude,
+            daily_all_pass.latitude,
         )
         == 0
     )
 
-    if is_heatwave:
-        exc_mask = (daily > clim_aligned) & land_mask
-    else:
-        exc_mask = (daily < clim_aligned) & land_mask
+    exc_mask = daily_all_pass & land_mask
 
     logger.info("  Computing exceedance mask...")
     mask_np = exc_mask.compute().values.astype(bool)
-    filtered = _apply_consecutive_filter(mask_np)
 
-    all_lats = daily.latitude.values
-    all_lons = daily.longitude.values
-    grid_res = np.abs(np.diff(all_lats[:2]))[0] if len(all_lats) > 1 else 0.25
+    all_lats = daily_all_pass.latitude.values
+    all_lons = daily_all_pass.longitude.values
+    grid_res = (
+        np.abs(np.diff(all_lats[:2]))[0] if len(all_lats) > 1 else 0.25
+    )
     band_pts = max(1, int(round(EXPANSION_DEGREES / grid_res)))
 
     def _lat_idx(val: float) -> int:
@@ -220,22 +456,159 @@ def _lat_idx(val: float) -> int:
     def _lon_idx(val: float) -> int:
         return int(np.argmin(np.abs(all_lons - val)))
 
-    idx_s = _lat_idx(lat_min)
-    idx_n = _lat_idx(lat_max)
-    idx_w = _lon_idx(lon_min)
-    idx_e = _lon_idx(lon_max)
+    # Convert box bounds to ERA5 lon convention for index lookup
+    box_lon_min_era = box_lon_min
+    box_lon_max_era = box_lon_max
+    if box_lon_min_era > 180:
+        box_lon_min_era -= 360
+    if box_lon_max_era > 180:
+        box_lon_max_era -= 360
+
+    idx_s0 = _lat_idx(box_lat_min)
+    idx_n0 = _lat_idx(box_lat_max)
+    idx_w0 = _lon_idx(box_lon_min_era)
+    idx_e0 = _lon_idx(box_lon_max_era)
+
+    # Temporal iteration
+    # Land mask within the initial box (exclude ocean from denominator)
+    land_mask_np = land_mask.compute().values.astype(bool)
+    box_land = land_mask_np[idx_s0 : idx_n0 + 1, idx_w0 : idx_e0 + 1]
+    n_land_pts = int(box_land.sum())
+    if n_land_pts == 0:
+        logger.warning(
+            "  Case %d: no land points in initial box — skipping",
+            single_case.case_id_number,
+        )
+        return None
+
+    n_days_total = min(mask_np.shape[0], MAX_TEMPORAL_DAYS)
+    final_t = n_days_total
+    established = False
 
-    edges_active = {
-        "north": True,
-        "south": True,
-        "east": True,
-        "west": True,
-    }
+    for t in range(n_days_total):
+        if t < 2:
+            continue
+
+        # "Currently active": last 3 days all exceed threshold
+        box_sl = (
+            slice(idx_s0, idx_n0 + 1),
+            slice(idx_w0, idx_e0 + 1),
+        )
+        currently_active = (
+            mask_np[t][box_sl]
+            & mask_np[t - 1][box_sl]
+            & mask_np[t - 2][box_sl]
+        )
+        active_land = int((currently_active & box_land).sum())
+        frac = active_land / n_land_pts
+
+        if t == 2 and frac < EDGE_VALIDITY_THRESHOLD:
+            logger.warning(
+                "  Day 3: only %.1f%% of land points have"
+                " 3 consecutive days (< 50%%)",
+                frac * 100,
+            )
+
+        if not established:
+            if frac >= EDGE_VALIDITY_THRESHOLD:
+                established = True
+                logger.info(
+                    "  Event established at day %d"
+                    " (%.1f%% of land points active)",
+                    t,
+                    frac * 100,
+                )
+        else:
+            if frac < EDGE_VALIDITY_THRESHOLD:
+                final_t = t + 1
+                logger.info(
+                    "  Temporal stop at day %d"
+                    " (%.1f%% < 50%% of land points)",
+                    t,
+                    frac * 100,
+                )
+                break
+
+        if t == 9:
+            logger.warning(
+                "  Exceeded 10 days (frac %.1f%%), continuing...",
+                frac * 100,
+            )
+
+    if not established:
+        logger.warning(
+            "  Case %d: event never reached 50%% of land"
+            " points — using all %d days",
+            single_case.case_id_number,
+            final_t,
+        )
+
+    logger.info(
+        "  Using %d of %d available days",
+        final_t,
+        mask_np.shape[0],
+    )
+
+    mask_np = mask_np[:final_t]
+    filtered = _apply_consecutive_filter(mask_np)
+
+    # Spatial expansion on peak-footprint day
+    # Find the timestep(s) with the most active grid points.
+    # If tied: 2nd of 2, middle of odd count, first-middle of even.
+    daily_counts = filtered.sum(axis=(1, 2))
+    max_count = daily_counts.max()
+    (tied_days,) = np.where(daily_counts == max_count)
+    n_tied = len(tied_days)
+    if n_tied <= 2:
+        peak_day = int(tied_days[-1]) if n_tied == 2 else int(tied_days[0])
+    else:
+        peak_day = int(tied_days[(n_tied - 1) // 2])
+    peak_mask = filtered[peak_day]
+    logger.info(
+        "  Peak footprint on day %d (%d active grid points,"
+        " %d tied days)",
+        peak_day,
+        int(max_count),
+        n_tied,
+    )
+
+    idx_s = idx_s0
+    idx_n = idx_n0
+    idx_w = idx_w0
+    idx_e = idx_e0
+
+    # Pre-check: disable edges that are >= 95% ocean
+    # in the initial box (prevents expansion through water).
+    init_region_land = land_mask_np[
+        idx_s0 : idx_n0 + 1, idx_w0 : idx_e0 + 1
+    ]
+    edges_active = {}
+    for edge in ("north", "south", "east", "west"):
+        if edge == "north":
+            strip = init_region_land[-band_pts:, :]
+        elif edge == "south":
+            strip = init_region_land[:band_pts, :]
+        elif edge == "west":
+            strip = init_region_land[:, :band_pts]
+        else:
+            strip = init_region_land[:, -band_pts:]
+        land_frac = strip.sum() / max(strip.size, 1)
+        edges_active[edge] = land_frac >= 0.25
+        if not edges_active[edge]:
+            logger.info(
+                "  %s edge disabled (%.1f%% land"
+                " in initial box)",
+                edge.capitalize(),
+                land_frac * 100,
+            )
     n_iter = 0
 
     for iteration in range(MAX_ITERATIONS):
         n_iter = iteration + 1
-        region = filtered[:, idx_s : idx_n + 1, idx_w : idx_e + 1]
+        region = peak_mask[idx_s : idx_n + 1, idx_w : idx_e + 1]
+        land_region = land_mask_np[
+            idx_s : idx_n + 1, idx_w : idx_e + 1
+        ]
 
         all_done = True
         for edge in list(edges_active.keys()):
@@ -244,6 +617,7 @@ def _lon_idx(val: float) -> int:
 
             frac = _edge_valid_fraction(
                 region,
+                land_region,
                 edge,
                 band_pts,
             )
@@ -252,23 +626,27 @@ def _lon_idx(val: float) -> int:
             else:
                 all_done = False
                 if edge == "north":
-                    idx_n = min(len(all_lats) - 1, idx_n + band_pts)
+                    idx_n = min(
+                        len(all_lats) - 1, idx_n + band_pts,
+                    )
                 elif edge == "south":
                     idx_s = max(0, idx_s - band_pts)
                 elif edge == "east":
-                    idx_e = min(len(all_lons) - 1, idx_e + band_pts)
+                    idx_e = min(
+                        len(all_lons) - 1, idx_e + band_pts,
+                    )
                 elif edge == "west":
                     idx_w = max(0, idx_w - band_pts)
 
         if all_done:
             logger.info(
-                "  Converged at iteration %d",
+                "  Spatial expansion converged at iteration %d",
                 n_iter,
             )
             break
     else:
         logger.info(
-            "  Reached max iterations (%d)",
+            "  Reached max spatial iterations (%d)",
             MAX_ITERATIONS,
         )
 
@@ -276,7 +654,9 @@ def _lon_idx(val: float) -> int:
     fin_lons = all_lons[idx_w : idx_e + 1]
     fin_filtered = filtered[:, idx_s : idx_n + 1, idx_w : idx_e + 1]
     consec = max_consecutive_days(fin_filtered)
-    peak_gridpoints = int(fin_filtered.any(axis=0).sum())
+    peak_gridpoints = int(
+        (consec >= MIN_CONSECUTIVE_DAYS).sum()
+    )
 
     result = {
         "case_id": single_case.case_id_number,
@@ -289,25 +669,90 @@ def _lon_idx(val: float) -> int:
         "longitude_min": float(all_lons[idx_w]),
         "longitude_max": float(all_lons[idx_e]),
         "n_iterations": n_iter,
+        "final_days": final_t,
         "_consec": consec,
         "_lats": fin_lats,
         "_lons": fin_lons,
         "_peak_gridpoints": peak_gridpoints,
     }
     logger.info(
-        "  Final bounds: lat [%.2f, %.2f], lon [%.2f, %.2f]",
+        "  Final bounds: lat [%.2f, %.2f], lon [%.2f, %.2f]"
+        " (%d days, %d iterations)",
+        result["latitude_min"],
+        result["latitude_max"],
+        result["longitude_min"],
+        result["longitude_max"],
+        final_t,
+        n_iter,
+    )
+
+    kind = "heatwave" if is_heatwave else "freeze"
+    start = result["start_date"][:10]
+    end = result["end_date"][:10]
+    out_png = str(
+        out_dir / f"case_{result['case_id']}_consecutive_{kind}_days.png"
+    )
+    plot_consecutive_map(
+        consec,
+        fin_lats,
+        fin_lons,
+        single_case.event_type,
+        title=(
+            f"Consecutive {kind.capitalize()} Days"
+            f" — {result['title']}\n{start} to {end}"
+        ),
+        output_path=out_png,
+    )
+    logger.info("  Saved plot: %s", out_png)
+
+    initial_box = (box_lat_min, box_lat_max, box_lon_min, box_lon_max)
+    final_box = (
         result["latitude_min"],
         result["latitude_max"],
         result["longitude_min"],
         result["longitude_max"],
     )
+    peak_png = str(
+        out_dir
+        / f"case_{result['case_id']}_peak_day_{kind}.png"
+    )
+    plot_peak_day_with_bounds(
+        peak_mask,
+        all_lats,
+        all_lons,
+        initial_box,
+        final_box,
+        single_case.event_type,
+        title=(
+            f"Peak Footprint (day {peak_day})"
+            f" — {result['title']}\n{start} to {end}"
+        ),
+        output_path=peak_png,
+    )
+
     return result
 
 
 def results_to_dataframe(results: List[Dict]) -> pd.DataFrame:
-    """Convert result dicts to a labelled DataFrame."""
+    """Convert result dicts to a labelled DataFrame.
+
+    Args:
+        results: List of dicts returned by process_event (may
+            contain None entries which are ignored). Each dict must
+            have keys event_type, start_date, end_date,
+            latitude_min/max, longitude_min/max, case_id, title, and
+            _peak_gridpoints.
+
+    Returns:
+        DataFrame with columns label, case_id, title, event_type,
+        start_date, end_date, latitude_min, latitude_max,
+        longitude_min, longitude_max, sorted by start_date.
+        Events below MIN_GRIDPOINTS are excluded.
+    """
     columns = [
         "label",
+        "case_id",
+        "title",
         "event_type",
         "start_date",
         "end_date",
@@ -346,10 +791,81 @@ def results_to_dataframe(results: List[Dict]) -> pd.DataFrame:
     return df
 
 
+def write_bounds_to_yaml(
+    results: List[Dict],
+    yaml_path: pathlib.Path,
+) -> None:
+    """Update events.yaml with computed bounded_region bounds.
+
+    Uses ruamel.yaml round-trip mode so indentation, quoting style,
+    key order, blank lines, and comments are preserved exactly.
+    Only heat_wave / freeze cases with a valid result are modified.
+
+    Args:
+        results: List of dicts returned by process_event (None
+            entries are ignored).
+        yaml_path: Path to the events.yaml file to update in-place.
+    """
+    result_map = {
+        r["case_id"]: r for r in results if r is not None
+    }
+    if not result_map:
+        logger.warning(
+            "write_bounds_to_yaml: no valid results, skipping."
+        )
+        return
+
+    yaml = YAML(typ="rt")
+    yaml.preserve_quotes = True
+
+    with yaml_path.open("r") as fh:
+        data = yaml.load(fh)
+
+    updated = 0
+    for entry in data:
+        cid = entry.get("case_id_number")
+        if cid not in result_map:
+            continue
+        r = result_map[cid]
+        params = entry["location"]["parameters"]
+        params["latitude_min"] = round(r["latitude_min"], 2)
+        params["latitude_max"] = round(r["latitude_max"], 2)
+        params["longitude_min"] = round(r["longitude_min"], 2)
+        params["longitude_max"] = round(r["longitude_max"], 2)
+        updated += 1
+
+    with yaml_path.open("w") as fh:
+        yaml.dump(data, fh)
+
+    logger.info(
+        "write_bounds_to_yaml: updated %d cases in %s",
+        updated,
+        yaml_path,
+    )
+
+
+def _load_base_temp_events() -> list[cases.IndividualCase]:
+    """Load all cases from base_temp_events.yaml.
+
+    Returns:
+        List of IndividualCase objects parsed from the bundled
+        base_temp_events.yaml resource file.
+    """
+    import extremeweatherbench.data
+
+    old_yaml = importlib.resources.files(
+        extremeweatherbench.data,
+    ).joinpath("base_temp_events.yaml")
+    with importlib.resources.as_file(old_yaml) as f:
+        raw = cases.read_incoming_yaml(f)
+    return cases.load_individual_cases(raw)
+
+
 def main():
     parser = argparse.ArgumentParser(
         description=(
-            "Validate / expand heat wave and cold snap bounds from events.yaml."
+            "Validate / expand heat wave and freeze bounds"
+            " using base_temp_events.yaml centered_region data."
         ),
     )
     parser.add_argument(
@@ -363,6 +879,37 @@ def main():
         default=4,
         help="Number of parallel workers (joblib)",
     )
+    parser.add_argument(
+        "--case-min",
+        type=int,
+        default=None,
+        help="Minimum case_id_number to process (inclusive)",
+    )
+    parser.add_argument(
+        "--case-max",
+        type=int,
+        default=None,
+        help="Maximum case_id_number to process (inclusive)",
+    )
+    parser.add_argument(
+        "--quantile",
+        type=float,
+        default=None,
+        help=(
+            "Climatology quantile "
+            f"({VALID_QUANTILES}; "
+            "default: 0.85 for heat_wave, 0.15 for freeze)"
+        ),
+    )
+    parser.add_argument(
+        "--operator",
+        default=None,
+        help=(
+            "Comparison operator "
+            "(>, >=, <, <=, ==; "
+            "default: > for heat_wave, < for freeze)"
+        ),
+    )
     args = parser.parse_args()
 
     wall_start = time_module.time()
@@ -370,43 +917,46 @@ def main():
     client = Client(LocalCluster(n_workers=args.n_workers))
     logger.info("Dask dashboard: %s", client.dashboard_link)
 
-    events_yaml = cases.load_ewb_events_yaml_into_case_list()
-    hw_fz = [e for e in events_yaml if e.event_type in ("heat_wave", "cold_snap")]
+    old_cases = _load_base_temp_events()
+    hw_fz = [
+        e for e in old_cases
+        if e.event_type in ("heat_wave", "freeze")
+        and (args.case_min is None
+             or e.case_id_number >= args.case_min)
+        and (args.case_max is None
+             or e.case_id_number <= args.case_max)
+    ]
     logger.info(
-        "Found %d heat_wave / cold snap events",
+        "Found %d heat_wave / freeze events in base_temp_events.yaml",
         len(hw_fz),
     )
 
+    out_dir = pathlib.Path(args.output).parent
+
     results = joblib.Parallel(n_jobs=args.n_workers)(
-        joblib.delayed(process_event)(c) for c in hw_fz
+        joblib.delayed(process_event)(
+            c, out_dir,
+            quantile=args.quantile,
+            op_str=args.operator,
+        )
+        for c in hw_fz
     )
 
     df = results_to_dataframe(results)
-    df.to_csv(args.output, index=False)
-
-    out_dir = pathlib.Path(args.output).parent
-    for r in results:
-        if r is None:
-            continue
-        consec = r["_consec"]
-        lats = r["_lats"]
-        lons = r["_lons"]
-        case_id = r["case_id"]
-        event_type = r["event_type"]
-        kind = "heatwave" if event_type == "heat_wave" else "cold_snap"
-        start = r["start_date"][:10]
-        end = r["end_date"][:10]
-        out_png = str(out_dir / f"case_{case_id}_consecutive_{kind}_days.png")
-        plot_consecutive_map(
-            consec,
-            lats,
-            lons,
-            event_type,
-            title=(
-                f"Consecutive {kind.capitalize()} Days — {r['title']}\n{start} to {end}"
-            ),
-            output_path=out_png,
-        )
+    out_path = pathlib.Path(args.output)
+    if out_path.exists():
+        df.to_csv(out_path, index=False, mode="a", header=False)
+    else:
+        df.to_csv(out_path, index=False)
+
+    yaml_path = (
+        pathlib.Path(__file__).parent.parent
+        / "src"
+        / "extremeweatherbench"
+        / "data"
+        / "events.yaml"
+    )
+    write_bounds_to_yaml(results, yaml_path)
 
     elapsed = time_module.time() - wall_start
     logger.info(
diff --git a/data_prep/heat_cold_bounds_global.py b/data_prep/heat_cold_bounds_global.py
index 1e81845a..c5f54d8f 100644
--- a/data_prep/heat_cold_bounds_global.py
+++ b/data_prep/heat_cold_bounds_global.py
@@ -29,11 +29,12 @@
 import xarray as xr
 from dask.distributed import Client, LocalCluster
 from plot_temperature_events import (
-    _align_climatology,
+    VALID_QUANTILES,
     detect_time_dim,
     max_consecutive_days,
     open_era5_t2m,
     plot_consecutive_map,
+    resolve_op,
 )
 
 from extremeweatherbench import defaults
@@ -49,26 +50,62 @@
 MIN_GRIDPOINTS = 500
 
 
-def get_daily_climatology_thresholds() -> Tuple[xr.DataArray, xr.DataArray]:
-    """Derive daily thresholds from 6-hourly percentile climatology.
+def get_climatology_thresholds(
+    q_hw: float = 0.85,
+    q_fz: float = 0.15,
+    q_hw_upper: Optional[float] = None,
+    q_fz_lower: Optional[float] = None,
+) -> Tuple[
+    xr.DataArray,
+    xr.DataArray,
+    Optional[xr.DataArray],
+    Optional[xr.DataArray],
+]:
+    """Return percentile climatology DataArrays for heat/freeze detection.
 
-    Returns the max-over-hours of the 85th percentile (heat wave
-    threshold) and min-over-hours of the 15th percentile (cold snap
-    threshold), each indexed by dayofyear.
+    Args:
+        q_hw: Lower-bound quantile for heat wave detection.
+        q_fz: Upper-bound quantile for freeze detection.
+        q_hw_upper: Upper-bound quantile for heat waves. When set,
+            only days where temp > q_hw AND temp < q_hw_upper are
+            flagged.
+        q_fz_lower: Lower-bound quantile for freezes. When set, only
+            days where temp < q_fz AND temp > q_fz_lower are flagged.
+
+    Returns:
+        A tuple of (clim_hw, clim_fz, clim_hw_upper, clim_fz_lower).
+        The last two elements are None when the corresponding optional
+        quantile argument is not supplied.
     """
-    clim_85 = defaults.get_climatology(0.85)
-    clim_15 = defaults.get_climatology(0.15)
-    return (
-        clim_85.max(dim="hour").sortby("latitude"),
-        clim_15.min(dim="hour").sortby("latitude"),
+    clim_hw = defaults.get_climatology(q_hw).sortby("latitude")
+    clim_fz = defaults.get_climatology(q_fz).sortby("latitude")
+    clim_hw_upper = (
+        defaults.get_climatology(q_hw_upper).sortby("latitude")
+        if q_hw_upper is not None
+        else None
+    )
+    clim_fz_lower = (
+        defaults.get_climatology(q_fz_lower).sortby("latitude")
+        if q_fz_lower is not None
+        else None
     )
+    return clim_hw, clim_fz, clim_hw_upper, clim_fz_lower
 
 
 def build_land_mask(
     lons: xr.DataArray,
     lats: xr.DataArray,
 ) -> xr.DataArray:
-    """Boolean land mask (True = land) for the given grid."""
+    """Build a boolean land mask (True = land) for the given grid.
+
+    Args:
+        lons: 1-D DataArray of longitudes.
+        lats: 1-D DataArray of latitudes.
+
+    Returns:
+        Boolean DataArray of the same shape as the meshgrid of lons
+        and lats, where True indicates a land grid point.
+    """
     land = regionmask.defined_regions.natural_earth_v5_0_0.land_110
     mask = land.mask(lons, lats)
     return mask == 0
@@ -76,61 +113,141 @@ def build_land_mask(
 
 def build_exceedance_masks(
     t2m: xr.DataArray,
-    clim_daily_max: xr.DataArray,
-    clim_daily_min: xr.DataArray,
+    clim_hw: xr.DataArray,
+    clim_fz: xr.DataArray,
     land_mask: xr.DataArray,
+    op_hw: str = ">",
+    op_fz: str = "<",
+    clim_hw_upper: Optional[xr.DataArray] = None,
+    clim_fz_lower: Optional[xr.DataArray] = None,
 ) -> Tuple[xr.DataArray, xr.DataArray]:
-    """Build lazy boolean masks for heat wave / cold snap exceedance.
+    """Build daily exceedance masks from 6-hourly temperature data.
 
-    Daily max/min aggregation stays in the dask graph. The
-    climatology is loaded into memory once (small) and then
-    broadcast against the lazy daily arrays to avoid chunk
-    multiplication warnings.
-    """
-    tdim = detect_time_dim(t2m)
+    Each 6-hourly timestep is compared to its matching
+    (dayofyear, hour) climatology. A day passes only if all four
+    6-hourly timesteps exceed the threshold.
 
-    daily_max = t2m.resample({tdim: "1D"}).max()
-    daily_min = t2m.resample({tdim: "1D"}).min()
+    When clim_hw_upper is provided, heat-wave days must also satisfy
+    temp < upper bound on every 6-hourly step. When clim_fz_lower is
+    provided, freeze days must also satisfy temp > lower bound on
+    every 6-hourly step.
 
-    clim_max_aligned = _align_climatology(
-        clim_daily_max,
-        daily_max,
-        tdim,
-    )
-    clim_min_aligned = _align_climatology(
-        clim_daily_min,
-        daily_min,
-        tdim,
-    )
+    Args:
+        t2m: 6-hourly 2m temperature DataArray.
+        clim_hw: Heat-wave lower-bound climatology indexed by
+            (dayofyear, hour).
+        clim_fz: Freeze upper-bound climatology indexed by
+            (dayofyear, hour).
+        land_mask: Boolean DataArray (True = land) matching t2m grid.
+        op_hw: Comparison operator string for heat waves.
+        op_fz: Comparison operator string for freezes.
+        clim_hw_upper: Optional upper-bound climatology for heat
+            waves (exclusive cap).
+        clim_fz_lower: Optional lower-bound climatology for freezes
+            (exclusive floor).
+
+    Returns:
+        A tuple (hw, fz) of daily boolean DataArrays masked to land,
+        where True indicates an exceedance day.
+    """
+    cmp_hw = resolve_op(op_hw)
+    cmp_fz = resolve_op(op_fz)
+    tdim = detect_time_dim(t2m)
 
-    hw = (daily_max > clim_max_aligned) & land_mask
-    fz = (daily_min < clim_min_aligned) & land_mask
+    doy = t2m[tdim].dt.dayofyear
+    hour = t2m[tdim].dt.hour
+    max_clim_doy = int(clim_hw.dayofyear.max())
+    doy_capped = doy.clip(max=max_clim_doy)
+
+    clim_hw_aligned = clim_hw.sel(
+        dayofyear=doy_capped, hour=hour,
+    ).reindex_like(t2m, method="nearest")
+    clim_fz_aligned = clim_fz.sel(
+        dayofyear=doy_capped, hour=hour,
+    ).reindex_like(t2m, method="nearest")
+
+    hw_6h = cmp_hw(t2m, clim_hw_aligned)
+    fz_6h = cmp_fz(t2m, clim_fz_aligned)
+
+    if clim_hw_upper is not None:
+        clim_hw_upper_aligned = clim_hw_upper.sel(
+            dayofyear=doy_capped, hour=hour,
+        ).reindex_like(t2m, method="nearest")
+        hw_6h = hw_6h & (t2m < clim_hw_upper_aligned)
+
+    if clim_fz_lower is not None:
+        clim_fz_lower_aligned = clim_fz_lower.sel(
+            dayofyear=doy_capped, hour=hour,
+        ).reindex_like(t2m, method="nearest")
+        fz_6h = fz_6h & (t2m > clim_fz_lower_aligned)
+
+    hw = hw_6h.resample({tdim: "1D"}).min().astype(bool) & land_mask
+    fz = fz_6h.resample({tdim: "1D"}).min().astype(bool) & land_mask
     return hw, fz
 
 
 def apply_consecutive_filter(
     mask: np.ndarray,
     min_days: int = MIN_CONSECUTIVE_DAYS,
+    max_grace_days: int = 1,
 ) -> np.ndarray:
-    """Keep only grid points in runs of ``min_days``+ True days.
+    """Keep runs of ``min_days``+ True days along axis 0.
+
+    After ``min_days`` strict consecutive True days are established,
+    gaps of up to ``max_grace_days`` are bridged so the event can
+    continue. Runs that never reach ``min_days`` strict consecutive
+    True days are discarded.
 
-    Binary erosion removes runs shorter than ``min_days``;
-    dilation restores surviving runs to their original extent.
-    Operates only along axis 0 (time).
+    Args:
+        mask: Boolean array of shape (time, lat, lon).
+        min_days: Minimum run length required to qualify as an event.
+        max_grace_days: Maximum gap length to bridge after the
+            minimum run is established.
+
+    Returns:
+        Boolean array of the same shape with only qualifying runs
+        retained.
     """
     struct = np.zeros((min_days, 1, 1), dtype=bool)
     struct[:, 0, 0] = True
-    eroded = ndimage.binary_erosion(
-        mask,
-        structure=struct,
-        border_value=False,
+
+    strict = (
+        ndimage.binary_dilation(
+            ndimage.binary_erosion(
+                mask, structure=struct, border_value=False,
+            ),
+            structure=struct, border_value=False,
+        )
+        & mask
     )
-    dilated = ndimage.binary_dilation(
-        eroded,
-        structure=struct,
-        border_value=False,
+
+    if max_grace_days <= 0:
+        return strict
+
+    close_k = np.zeros((2 * max_grace_days + 1, 1, 1), dtype=bool)
+    close_k[:, 0, 0] = True
+    filled = ndimage.binary_closing(
+        mask, structure=close_k, border_value=False,
+    )
+
+    filled_runs = (
+        ndimage.binary_dilation(
+            ndimage.binary_erosion(
+                filled, structure=struct, border_value=False,
+            ),
+            structure=struct, border_value=False,
+        )
+        & filled
     )
-    return dilated & mask
+
+    lbl_struct = np.zeros((3, 3, 3), dtype=int)
+    lbl_struct[0, 1, 1] = 1
+    lbl_struct[1, 1, 1] = 1
+    lbl_struct[2, 1, 1] = 1
+    labels, _ = ndimage.label(filled_runs, structure=lbl_struct)
+
+    valid = np.unique(labels[strict & (labels > 0)])
+    return np.isin(labels, valid) & filled_runs
 
 
 @nb.njit(cache=True)
@@ -140,11 +257,20 @@ def _count_overlaps_nb(
     n_a: int,
     n_b: int,
 ) -> np.ndarray:
-    """Pixel-count overlap matrix between two consecutive-day label grids.
+    """Compute a pixel-count overlap matrix between two label grids.
 
-    Returns int32 (n_a, n_b) where result[i, j] is the number of
-    pixels where labels_a == i+1 and labels_b == j+1. Single-pass
-    JIT loop avoids the per-blob numpy scan in the Python tracker.
+    A single-pass JIT loop avoids the per-blob numpy scan in the
+    Python tracker.
+
+    Args:
+        labels_a: Integer label array for day t (shape lat x lon).
+        labels_b: Integer label array for day t+1 (shape lat x lon).
+        n_a: Number of blobs in labels_a.
+        n_b: Number of blobs in labels_b.
+
+    Returns:
+        Int32 array of shape (n_a, n_b) where element [i, j] is the
+        number of pixels where labels_a == i+1 and labels_b == j+1.
     """
     mat = np.zeros((n_a, n_b), dtype=np.int32)
     rows, cols = labels_a.shape
@@ -164,10 +290,23 @@ def _resolve_event(
     events: Dict[int, Dict],
     cur_map: Dict[int, int],
 ) -> Optional[int]:
-    """Find and merge events from the previous day for blob oid.
+    """Find and merge prior-day events overlapping with blob oid.
 
-    Uses a precomputed overlap matrix row instead of scanning the
+    Uses a precomputed overlap matrix column instead of scanning the
     full prev_labels array per blob.
+
+    Args:
+        oid: Current-day blob label (1-indexed).
+        overlap_mat: Pixel-count overlap matrix from
+            _count_overlaps_nb, or None if no previous day.
+        prev_map: Mapping from previous-day blob label to event ID.
+        events: Mutable dict of all live events keyed by event ID.
+        cur_map: Mutable mapping from current-day blob label to event
+            ID; updated in-place when events are merged.
+
+    Returns:
+        The surviving event ID after merging, or None if no overlap
+        with a prior-day event was found.
     """
     if overlap_mat is None:
         return None
@@ -198,7 +337,13 @@ def _terminate_declined_events(
     events: Dict[int, Dict],
     cur_map: Dict[int, int],
 ) -> None:
-    """Terminate events absent today or below 50% of peak."""
+    """Mark events as done if absent today or below 50% of peak area.
+
+    Args:
+        events: Mutable dict of all live events keyed by event ID.
+        cur_map: Mapping from current-day blob label to event ID;
+            used to determine which events are still active today.
+    """
     active = set(cur_map.values())
     for eid, ev in events.items():
         if ev["done"]:
@@ -225,6 +370,19 @@ def detect_events(
 
     An event's bounding box is the union of all its daily extents.
     Events terminate when their active area drops below 50% of peak.
+
+    Args:
+        filtered_mask: Boolean array of shape (time, lat, lon) with
+            consecutive-day filtering already applied.
+        dates: 1-D array of date labels aligned with axis 0.
+        lats: 1-D latitude array aligned with axis 1.
+        lons: 1-D longitude array aligned with axis 2.
+        event_type: Label string stored in each returned event dict
+            (e.g. "heat_wave" or "cold_snap").
+
+    Returns:
+        List of event dicts with keys type, start, end, lat_min,
+        lat_max, lon_min, lon_max, peak, area, done.
     """
     n_days = filtered_mask.shape[0]
     events: Dict[int, Dict] = {}
@@ -329,9 +487,14 @@ def events_to_dataframe(
 
     Args:
         events: Raw event dicts from ``detect_events``.
-        min_gridpoints: Drop events whose peak spatial extent
-            (in grid points) is below this threshold. Defaults
-            to the module-level ``MIN_GRIDPOINTS`` constant.
+        min_gridpoints: Drop events whose peak spatial extent (in
+            grid points) is below this threshold.
+
+    Returns:
+        DataFrame with columns label, event_type, start_date,
+        end_date, latitude_min, latitude_max, longitude_min,
+        longitude_max, sorted by start_date. Events below
+        min_gridpoints are excluded.
     """
     columns = [
         "label",
@@ -400,6 +563,54 @@ def main():
         default=4,
         help="Number of dask workers",
     )
+    parser.add_argument(
+        "--quantile-hw",
+        type=float,
+        default=0.85,
+        help=(
+            "Climatology quantile for heat waves "
+            f"({VALID_QUANTILES}; default: 0.85)"
+        ),
+    )
+    parser.add_argument(
+        "--quantile-fz",
+        type=float,
+        default=0.15,
+        help=(
+            "Climatology quantile for freezes "
+            f"({VALID_QUANTILES}; default: 0.15)"
+        ),
+    )
+    parser.add_argument(
+        "--quantile-hw-upper",
+        type=float,
+        default=None,
+        help=(
+            "Upper-bound quantile for heat waves; days must be "
+            "> --quantile-hw AND < this value "
+            f"({VALID_QUANTILES}; default: None)"
+        ),
+    )
+    parser.add_argument(
+        "--quantile-fz-lower",
+        type=float,
+        default=None,
+        help=(
+            "Lower-bound quantile for freezes; days must be "
+            "< --quantile-fz AND > this value "
+            f"({VALID_QUANTILES}; default: None)"
+        ),
+    )
+    parser.add_argument(
+        "--operator-hw",
+        default=">",
+        help="Comparison operator for heat waves (default: >)",
+    )
+    parser.add_argument(
+        "--operator-fz",
+        default="<",
+        help="Comparison operator for freezes (default: <)",
+    )
     args = parser.parse_args()
 
     wall_start = time_module.time()
@@ -420,7 +631,14 @@ def main():
     logger.info("  sizes=%s", dict(t2m.sizes))
 
     logger.info("Loading climatology thresholds...")
-    clim_max, clim_min = get_daily_climatology_thresholds()
+    clim_hw, clim_fz, clim_hw_upper, clim_fz_lower = (
+        get_climatology_thresholds(
+            q_hw=args.quantile_hw,
+            q_fz=args.quantile_fz,
+            q_hw_upper=args.quantile_hw_upper,
+            q_fz_lower=args.quantile_fz_lower,
+        )
+    )
 
     logger.info("Building land mask...")
     land_mask = build_land_mask(t2m.longitude, t2m.latitude)
@@ -428,9 +646,13 @@ def main():
     logger.info("Building exceedance masks (lazy)...")
     hw_lazy, fz_lazy = build_exceedance_masks(
         t2m,
-        clim_max,
-        clim_min,
+        clim_hw,
+        clim_fz,
         land_mask,
+        op_hw=args.operator_hw,
+        op_fz=args.operator_fz,
+        clim_hw_upper=clim_hw_upper,
+        clim_fz_lower=clim_fz_lower,
     )
 
     tdim = detect_time_dim(hw_lazy)
diff --git a/data_prep/plot_temperature_events.py b/data_prep/plot_temperature_events.py
index 6517fdd9..ed33ac39 100644
--- a/data_prep/plot_temperature_events.py
+++ b/data_prep/plot_temperature_events.py
@@ -16,9 +16,10 @@
 
 import argparse
 import logging
+import operator as op_module
 import pathlib
 import time as time_module
-from typing import Literal, cast
+from typing import Callable, Literal, cast
 
 import cartopy.crs as ccrs
 import cartopy.feature as cfeature
@@ -41,11 +42,52 @@
 
 MIN_CONSECUTIVE_DAYS = 3
 
+_OP_MAP: dict[str, Callable] = {
+    ">": op_module.gt,
+    ">=": op_module.ge,
+    "<": op_module.lt,
+    "<=": op_module.le,
+    "==": op_module.eq,
+}
+
+VALID_QUANTILES = [0.10, 0.15, 0.25, 0.50, 0.75, 0.85, 0.90]
+
+
+def resolve_op(op_str: str) -> Callable:
+    """Convert a string operator to a callable.
+
+    Args:
+        op_str: One of ">", ">=", "<", "<=", "==".
+
+    Returns:
+        The corresponding operator callable.
+
+    Raises:
+        ValueError: If op_str is not a recognised operator.
+    """
+    if op_str not in _OP_MAP:
+        raise ValueError(
+            f"Unknown operator {op_str!r}; "
+            f"choose from {list(_OP_MAP)}"
+        )
+    return _OP_MAP[op_str]
+
+
 # ── shared ERA5 utilities ─────────────────────────────────────────────
 
 
 def detect_time_dim(obj: xr.Dataset | xr.DataArray) -> str:
-    """Return the name of the time dimension."""
+    """Return the name of the time dimension.
+
+    Args:
+        obj: An xarray Dataset or DataArray.
+
+    Returns:
+        The name of the first matching time dimension.
+
+    Raises:
+        ValueError: If no recognised time dimension is found.
+    """
     for name in ("valid_time", "time"):
         if name in obj.dims:
             return name
@@ -57,6 +99,13 @@ def open_era5_t2m(start_date: str, end_date: str) -> xr.DataArray:
 
     Selects 6-hourly timesteps (0/6/12/18 UTC) to match the
     climatology base and sorts latitude to ascending order.
+
+    Args:
+        start_date: Inclusive start date string (YYYY-MM-DD).
+        end_date: Inclusive end date string (YYYY-MM-DD).
+
+    Returns:
+        Lazy DataArray of 2m temperature with ascending latitude.
     """
     ds = xr.open_zarr(
         inputs.ARCO_ERA5_FULL_URI,
@@ -80,6 +129,15 @@ def _align_climatology(
 
     Computes the climatology into memory (366 x lat x lon) to avoid
     dask chunk multiplication, then indexes by dayofyear via numpy.
+
+    Args:
+        clim: Climatology DataArray indexed by dayofyear.
+        daily: Daily ERA5 DataArray to align against.
+        tdim: Name of the time dimension in daily.
+
+    Returns:
+        DataArray with the same shape as daily, containing the
+        climatology value for each timestep's day-of-year.
     """
     clim_vals = clim.compute().values
     clim_doy = clim.dayofyear.values
@@ -119,12 +177,27 @@ def _align_climatology(
 
 
 def _to_plot_lon(lon: float) -> float:
-    """Wrap 0-360 longitude to -180..180 for PlateCarree."""
+    """Wrap 0-360 longitude to -180..180 for PlateCarree.
+
+    Args:
+        lon: Longitude in 0-360 degrees.
+
+    Returns:
+        Longitude in -180..180 degrees.
+    """
     return lon - 360.0 if lon > 180.0 else lon
 
 
 def max_consecutive_days(mask_3d: np.ndarray) -> np.ndarray:
-    """Max consecutive True days per grid point along axis 0."""
+    """Compute max consecutive True days per grid point.
+
+    Args:
+        mask_3d: Boolean array of shape (time, lat, lon).
+
+    Returns:
+        Int32 array of shape (lat, lon) with the maximum number of
+        consecutive True values along axis 0 for each grid point.
+    """
     nt, nlat, nlon = mask_3d.shape
     flat = mask_3d.reshape(nt, -1).astype(np.int8)
     result = np.zeros(flat.shape[1], dtype=np.int32)
@@ -140,7 +213,11 @@ def max_consecutive_days(mask_3d: np.ndarray) -> np.ndarray:
 
 
 def _add_map_features(ax) -> None:
-    """Add standard cartopy features to an axis."""
+    """Add standard cartopy features to an axis.
+
+    Args:
+        ax: A cartopy GeoAxes instance.
+    """
     ax.coastlines(linewidth=0.5, zorder=10)
     ax.add_feature(
         cfeature.BORDERS,
@@ -187,7 +264,7 @@ def plot_consecutive_map(
         title: Two-line figure title (left-aligned).
         output_path: Destination PNG file path.
         extent: (lon_min, lon_max, lat_min, lat_max) in -180..180.
-            Derived from lats/lons when None.
+            Derived from lats/lons when None. Default is None.
     """
     plot_data = consec.astype(float)
     plot_data[consec < MIN_CONSECUTIVE_DAYS] = np.nan
@@ -258,6 +335,14 @@ def plot_event_bounds(
     """Draw bounding boxes for all events on a global Robinson map.
 
     Saves a PNG alongside the CSV with the same stem.
+
+    Args:
+        df: DataFrame with columns event_type, longitude_min,
+            longitude_max, latitude_min, latitude_max.
+        csv_path: Path to the output CSV; the PNG is saved with the
+            same stem.
+        title: Figure title. Default is "Detected Heat Wave and
+            Cold Snap Events".
     """
     if df.empty:
         logger.warning("No events to plot -- skipping bounds plot.")
@@ -352,7 +437,17 @@ def plot_event_bounds(
 
 
 def load_case(case_id_number: int) -> cases.IndividualCase:
-    """Load a single case from events.yaml by case_id_number."""
+    """Load a single case from events.yaml by case_id_number.
+
+    Args:
+        case_id_number: Integer identifier for the case.
+
+    Returns:
+        The matching IndividualCase object.
+
+    Raises:
+        ValueError: If no case with the given ID exists.
+    """
     all_cases = cases.load_ewb_events_yaml_into_case_list()
     for c in all_cases:
         if c.case_id_number == case_id_number:
@@ -372,13 +467,33 @@ def load_case(case_id_number: int) -> cases.IndividualCase:
 
 def compute_consecutive_field(
     single_case: cases.IndividualCase,
+    quantile: float | None = None,
+    op_str: str | None = None,
 ) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
     """Return (consecutive_days, lats, lons) for a case.
 
-    Works for both heat_wave (daily max > 85th pct) and
-    cold_snap (daily min < 15th pct) event types.
+    Args:
+        single_case: The case to compute the consecutive field for.
+        quantile: Climatology quantile. Default is None, which
+            resolves to 0.85 for heat_wave or 0.15 for
+            freeze/cold_snap.
+        op_str: Comparison operator string (e.g. ">", ">=",
+            "<", "<="). Default is None, which resolves to ">"
+            for heat_wave or "<" for freeze.
+
+    Returns:
+        A tuple of (consec, lats, lons) where consec is an int32
+        array of shape (lat, lon) containing max consecutive event
+        days, lats is the 1-D latitude array, and lons is the 1-D
+        longitude array (0-360).
     """
     is_heatwave = single_case.event_type == "heat_wave"
+    if quantile is None:
+        quantile = 0.85 if is_heatwave else 0.15
+    if op_str is None:
+        op_str = ">" if is_heatwave else "<"
+    cmp = resolve_op(op_str)
+
     start = str(single_case.start_date.date())
     end = str(single_case.end_date.date())
     bounds = single_case.location.as_geopandas().total_bounds
@@ -399,36 +514,47 @@ def compute_consecutive_field(
 
     tdim = detect_time_dim(t2m)
 
-    if is_heatwave:
-        logger.info("Loading 85th-percentile climatology...")
-        daily = t2m.resample({tdim: "1D"}).max()
-        clim = defaults.get_climatology(0.85).max(dim="hour").sortby("latitude")
-    else:
-        logger.info("Loading 15th-percentile climatology...")
-        daily = t2m.resample({tdim: "1D"}).min()
-        clim = defaults.get_climatology(0.15).min(dim="hour").sortby("latitude")
+    logger.info(
+        "Loading q=%.2f climatology (op=%s)...",
+        quantile, op_str,
+    )
+    clim = defaults.get_climatology(quantile).sortby("latitude")
 
     clim = clim.sel(
-        latitude=daily.latitude,
-        longitude=daily.longitude,
+        latitude=t2m.latitude,
+        longitude=t2m.longitude,
         method="nearest",
     )
-    clim_aligned = _align_climatology(clim, daily, tdim)
+
+    doy = t2m[tdim].dt.dayofyear
+    hour = t2m[tdim].dt.hour
+    max_clim_doy = int(clim.dayofyear.max())
+    doy_capped = doy.clip(max=max_clim_doy)
+
+    clim_aligned = clim.sel(
+        dayofyear=doy_capped,
+        hour=hour,
+    )
+
+    logger.info("Computing 6-hourly exceedance...")
+    exc_6h = cmp(t2m, clim_aligned)
+
+    daily_all_pass = (
+        exc_6h.resample({tdim: "1D"}).min().astype(bool)
+    )
 
     logger.info("Building land mask...")
     land = regionmask.defined_regions.natural_earth_v5_0_0.land_110
-    land_mask = land.mask(daily.longitude, daily.latitude) == 0
+    land_mask = land.mask(
+        daily_all_pass.longitude, daily_all_pass.latitude
+    ) == 0
 
-    logger.info("Computing exceedance mask...")
-    if is_heatwave:
-        exc = (daily > clim_aligned) & land_mask
-    else:
-        exc = (daily < clim_aligned) & land_mask
+    exc = daily_all_pass & land_mask
     mask_np = exc.compute().values.astype(bool)
 
     logger.info("Computing max consecutive days...")
     consec = max_consecutive_days(mask_np)
-    return consec, daily.latitude.values, daily.longitude.values
+    return consec, daily_all_pass.latitude.values, daily_all_pass.longitude.values
 
 
 # ── CLI ───────────────────────────────────────────────────────────────
@@ -459,6 +585,25 @@ def main() -> None:
             "(default: case_N_consecutive_{heatwave|cold_snap}_days.png)"
         ),
     )
+    parser.add_argument(
+        "--quantile",
+        type=float,
+        default=None,
+        help=(
+            "Climatology quantile "
+            f"(one of {VALID_QUANTILES}; "
+            "default: 0.85 for heat_wave, 0.15 for freeze)"
+        ),
+    )
+    parser.add_argument(
+        "--operator",
+        default=None,
+        help=(
+            "Comparison operator string "
+            "(>, >=, <, <=, ==; "
+            "default: > for heat_wave, < for freeze)"
+        ),
+    )
     args = parser.parse_args()
 
     t0 = time_module.time()
@@ -476,7 +621,11 @@ def main() -> None:
             single_case.event_type,
         )
 
-    consec, lats, lons = compute_consecutive_field(single_case)
+    consec, lats, lons = compute_consecutive_field(
+        single_case,
+        quantile=args.quantile,
+        op_str=args.operator,
+    )
 
     kind = "Heat Wave" if single_case.event_type == "heat_wave" else "Cold Snap"
     start_str = str(single_case.start_date.date())
diff --git a/src/extremeweatherbench/data/base_temp_events.yaml b/src/extremeweatherbench/data/base_temp_events.yaml
new file mode 100644
index 00000000..5cf4ff1f
--- /dev/null
+++ b/src/extremeweatherbench/data/base_temp_events.yaml
@@ -0,0 +1,661 @@
+- case_id_number: 1
+  title: 2021 Pacific Northwest
+  start_date: 2021-06-20 00:00:00
+  end_date: 2021-07-03 00:00:00
+  location:
+    type: centered_region
+    parameters:
+      latitude: 47.6062
+      longitude: 237.6679
+      bounding_box_degrees: 5
+  event_type: heat_wave
+- case_id_number: 2
+  title: 2022 Upper Midwest
+  start_date: 2022-05-07 00:00:00
+  end_date: 2022-05-17 00:00:00
+  location:
+    type: centered_region
+    parameters:
+      latitude: 41.8781
+      longitude: 272.3702
+      bounding_box_degrees: 5
+  event_type: heat_wave
+- case_id_number: 3
+  title: 2022 California
+  start_date: 2022-06-07 00:00:00
+  end_date: 2022-06-15 00:00:00
+  location:
+    type: centered_region
+    parameters:
+      latitude: 34.0522
+      longitude: 241.7563
+      bounding_box_degrees: 5
+  event_type: heat_wave
+- case_id_number: 4
+  title: 2022 Texas
+  start_date: 2022-06-30 00:00:00
+  end_date: 2022-07-18 00:00:00
+  location:
+    type: centered_region
+    parameters:
+      latitude: 32.7767
+      longitude: 263.203
+      bounding_box_degrees: 5
+  event_type: heat_wave
+- case_id_number: 5
+  title: 2023 Pacific Northwest
+  start_date: 2023-05-10 00:00:00
+  end_date: 2023-05-23 00:00:00
+  location:
+    type: centered_region
+    parameters:
+      latitude: 47.6062
+      longitude: 237.6679
+      bounding_box_degrees: 5
+  event_type: heat_wave
+- case_id_number: 6
+  title: 2022 Mid-Atlantic
+  start_date: 2022-05-17 00:00:00
+  end_date: 2022-05-24 00:00:00
+  location:
+    type: centered_region
+    parameters:
+      latitude: 39.2904
+      longitude: 283.38779999999997
+      bounding_box_degrees: 5
+  event_type: heat_wave
+- case_id_number: 7
+  title: 2023 Australia
+  start_date: 2023-11-18 00:00:00
+  end_date: 2023-11-28 00:00:00
+  location:
+    type: centered_region
+    parameters:
+      latitude: -31.9505
+      longitude: 115.8605
+      bounding_box_degrees: 5
+  event_type: heat_wave
+- case_id_number: 8
+  title: 2023 Ireland
+  start_date: 2023-09-02 00:00:00
+  end_date: 2023-09-13 00:00:00
+  location:
+    type: centered_region
+    parameters:
+      latitude: 53.1424
+      longitude: 352.3079
+      bounding_box_degrees: 5
+  event_type: heat_wave
+- case_id_number: 9
+  title: 2023 Italy
+  start_date: 2023-07-07 00:00:00
+  end_date: 2023-07-27 00:00:00
+  location:
+    type: centered_region
+    parameters:
+      latitude: 41.9028
+      longitude: 12.4964
+      bounding_box_degrees: 5
+  event_type: heat_wave
+- case_id_number: 10
+  title: 2023 SW Europe
+  start_date: 2023-08-17 00:00:00
+  end_date: 2023-08-28 00:00:00
+  location:
+    type: centered_region
+    parameters:
+      latitude: 40.4637
+      longitude: 356.2508
+      bounding_box_degrees: 5
+  event_type: heat_wave
+- case_id_number: 11
+  title: 2023 South America
+  start_date: 2023-07-29 00:00:00
+  end_date: 2023-08-04 00:00:00
+  location:
+    type: centered_region
+    parameters:
+      latitude: -27.0
+      longitude: 294.5
+      bounding_box_degrees: 5
+  event_type: heat_wave
+- case_id_number: 12
+  title: 2023 China (Shanghai)
+  start_date: 2023-05-24 00:00:00
+  end_date: 2023-06-01 00:00:00
+  location:
+    type: centered_region
+    parameters:
+      latitude: 31.2304
+      longitude: 121.4737
+      bounding_box_degrees: 5
+  event_type: heat_wave
+- case_id_number: 13
+  title: 2023 China (Yunnan Province)
+  start_date: 2023-04-14 00:00:00
+  end_date: 2023-04-23 00:00:00
+  location:
+    type: centered_region
+    parameters:
+      latitude: 25.0458
+      longitude: 102.71
+      bounding_box_degrees: 5
+  event_type: heat_wave
+- case_id_number: 14
+  title: 2023 Algeria
+  start_date: 2023-07-05 00:00:00
+  end_date: 2023-07-27 00:00:00
+  location:
+    type: centered_region
+    parameters:
+      latitude: 36.7372
+      longitude: 3.0869
+      bounding_box_degrees: 5
+  event_type: heat_wave
+- case_id_number: 15
+  title: 2023 Iberian Peninsula
+  start_date: 2023-04-22 00:00:00
+  end_date: 2023-05-01 00:00:00
+  location:
+    type: centered_region
+    parameters:
+      latitude: 41.1496
+      longitude: 352.389
+      bounding_box_degrees: 5
+  event_type: heat_wave
+- case_id_number: 16
+  title: 2023 Southeast Asia
+  start_date: 2023-04-16 00:00:00
+  end_date: 2023-04-22 00:00:00
+  location:
+    type: bounded_region
+    parameters:
+      latitude_min: 9.0
+      latitude_max: 31.5
+      longitude_min: 61.5
+      longitude_max: 113.75
+  event_type: heat_wave
+- case_id_number: 17
+  title: 2023 India
+  start_date: 2023-02-15 00:00:00
+  end_date: 2023-03-01 00:00:00
+  location:
+    type: centered_region
+    parameters:
+      latitude: 28.6139
+      longitude: 77.209
+      bounding_box_degrees: 5
+  event_type: heat_wave
+- case_id_number: 18
+  title: 2021 Western Russia
+  start_date: 2021-06-18 00:00:00
+  end_date: 2021-06-30 00:00:00
+  location:
+    type: centered_region
+    parameters:
+      latitude: 55.7558
+      longitude: 37.6173
+      bounding_box_degrees: 5
+  event_type: heat_wave
+- case_id_number: 19
+  title: 2022 Germany
+  start_date: 2022-12-23 00:00:00
+  end_date: 2022-12-31 00:00:00
+  location:
+    type: centered_region
+    parameters:
+      latitude: 52.52
+      longitude: 13.405
+      bounding_box_degrees: 5
+  event_type: heat_wave
+- case_id_number: 20
+  title: 2022 UK (August)
+  start_date: 2022-08-08 00:00:00
+  end_date: 2022-08-16 00:00:00
+  location:
+    type: centered_region
+    parameters:
+      latitude: 51.5074
+      longitude: 359.8722
+      bounding_box_degrees: 5
+  event_type: heat_wave
+- case_id_number: 21
+  title: 2022 UK (July)
+  start_date: 2022-07-15 00:00:00
+  end_date: 2022-07-23 00:00:00
+  location:
+    type: centered_region
+    parameters:
+      latitude: 51.5074
+      longitude: 359.8722
+      bounding_box_degrees: 5
+  event_type: heat_wave
+- case_id_number: 22
+  title: 2022 France
+  start_date: 2022-06-09 00:00:00
+  end_date: 2022-06-21 00:00:00
+  location:
+    type: centered_region
+    parameters:
+      latitude: 43.4832
+      longitude: 358.4414
+      bounding_box_degrees: 5
+  event_type: heat_wave
+- case_id_number: 23
+  title: 2022 Japan
+  start_date: 2022-06-20 00:00:00
+  end_date: 2022-07-05 00:00:00
+  location:
+    type: centered_region
+    parameters:
+      latitude: 35.6895
+      longitude: 139.6917
+      bounding_box_degrees: 5
+  event_type: heat_wave
+- case_id_number: 24
+  title: 2022 India
+  start_date: 2022-04-24 00:00:00
+  end_date: 2022-05-04 00:00:00
+  location:
+    type: centered_region
+    parameters:
+      latitude: 28.2769
+      longitude: 68.4376
+      bounding_box_degrees: 5
+  event_type: heat_wave
+- case_id_number: 25
+  title: 2022 East Antarctica
+  start_date: 2022-03-12 00:00:00
+  end_date: 2022-03-26 00:00:00
+  location:
+    type: centered_region
+    parameters:
+      latitude: -75.1
+      longitude: 123.35
+      bounding_box_degrees: 5
+  event_type: heat_wave
+- case_id_number: 26
+  title: 2022 West Australia
+  start_date: 2022-01-15 00:00:00
+  end_date: 2022-01-25 00:00:00
+  location:
+    type: centered_region
+    parameters:
+      latitude: -31.9505
+      longitude: 115.8605
+      bounding_box_degrees: 5
+  event_type: heat_wave
+- case_id_number: 27
+  title: 2021 Canada Plains
+  start_date: 2021-05-30 00:00:00
+  end_date: 2021-06-09 00:00:00
+  location:
+    type: centered_region
+    parameters:
+      latitude: 49.0
+      longitude: 262.43330000000003
+      bounding_box_degrees: 5
+  event_type: heat_wave
+- case_id_number: 28
+  title: 2021 New Zealand
+  start_date: 2021-01-12 18:00:00
+  end_date: 2021-01-18 18:00:00
+  location:
+    type: centered_region
+    parameters:
+      latitude: -43.8983
+      longitude: 171.731
+      bounding_box_degrees: 5
+  event_type: heat_wave
+- case_id_number: 29
+  title: 2020 Australia
+  start_date: 2020-11-25 00:00:00
+  end_date: 2020-12-01 00:00:00
+  location:
+    type: centered_region
+    parameters:
+      latitude: -33.8245
+      longitude: 150.9448
+      bounding_box_degrees: 5
+  event_type: heat_wave
+- case_id_number: 30
+  title: 2021 Texas
+  start_date: 2021-02-10 12:00:00
+  end_date: 2021-02-22 00:00:00
+  location:
+    type: centered_region
+    parameters:
+      latitude: 30.2672
+      longitude: 262.2569
+      bounding_box_degrees: 5
+  event_type: freeze
+- case_id_number: 31
+  title: 2022 Arkansas
+  start_date: 2022-02-17 18:00:00
+  end_date: 2022-03-01 06:00:00
+  location:
+    type: centered_region
+    parameters:
+      latitude: 34.7445
+      longitude: 267.7104
+      bounding_box_degrees: 5
+  event_type: freeze
+- case_id_number: 32
+  title: 2023 Germany
+  start_date: 2023-12-02 06:00:00
+  end_date: 2023-12-08 06:00:00
+  location:
+    type: centered_region
+    parameters:
+      latitude: 48.1351
+      longitude: 11.582
+      bounding_box_degrees: 5
+  event_type: freeze
+- case_id_number: 33
+  title: 2023 China
+  start_date: 2023-12-15 06:00:00
+  end_date: 2023-12-26 18:00:00
+  location:
+    type: centered_region
+    parameters:
+      latitude: 31.2304
+      longitude: 121.4737
+      bounding_box_degrees: 5
+  event_type: freeze
+- case_id_number: 34
+  title: 2023 Afghanistan
+  start_date: 2023-01-11 06:00:00
+  end_date: 2023-01-28 00:00:00
+  location:
+    type: centered_region
+    parameters:
+      latitude: 34.5553
+      longitude: 69.2075
+      bounding_box_degrees: 5
+  event_type: freeze
+- case_id_number: 35
+  title: 2022 Colorado
+  start_date: 2022-04-06 00:00:00
+  end_date: 2022-04-16 00:00:00
+  location:
+    type: centered_region
+    parameters:
+      latitude: 39.7392
+      longitude: 255.0097
+      bounding_box_degrees: 5
+  event_type: freeze
+- case_id_number: 72
+  title: May 2024 Texas
+  start_date: 2024-05-25 00:00:00
+  end_date: 2024-05-31 00:00:00
+  location:
+    type: centered_region
+    parameters:
+      latitude: 25.9017
+      longitude: 262.50260000000003
+      bounding_box_degrees: 5
+  event_type: heat_wave
+- case_id_number: 73
+  title: June 2024 Northeast US
+  start_date: 2024-06-17 00:00:00
+  end_date: 2024-06-23 00:00:00
+  location:
+    type: centered_region
+    parameters:
+      latitude: 41.8781
+      longitude: 286.771
+      bounding_box_degrees: 6
+  event_type: heat_wave
+- case_id_number: 74
+  title: July 2024 Southwest US
+  start_date: 2024-07-04 00:00:00
+  end_date: 2024-07-10 00:00:00
+  location:
+    type: centered_region
+    parameters:
+      latitude: 33.7701
+      longitude: 243.78539999999998
+      bounding_box_degrees: 6
+  event_type: heat_wave
+- case_id_number: 75
+  title: July 2024 Northeast US
+  start_date: 2024-07-07 00:00:00
+  end_date: 2024-07-13 00:00:00
+  location:
+    type: centered_region
+    parameters:
+      latitude: 40.7128
+      longitude: 285.994
+      bounding_box_degrees: 6
+  event_type: heat_wave
+- case_id_number: 76
+  title: July 2024 Mid-Atlantic US
+  start_date: 2024-07-15 00:00:00
+  end_date: 2024-07-21 00:00:00
+  location:
+    type: centered_region
+    parameters:
+      latitude: 39.9526
+      longitude: 284.8348
+      bounding_box_degrees: 6
+  event_type: heat_wave
+- case_id_number: 77
+  title: August 2024 Midwest US
+  start_date: 2024-08-25 00:00:00
+  end_date: 2024-08-31 00:00:00
+  location:
+    type: centered_region
+    parameters:
+      latitude: 40.1106
+      longitude: 271.79269999999997
+      bounding_box_degrees: 5
+  event_type: heat_wave
+- case_id_number: 78
+  title: July 2024 Antarctica
+  start_date: 2024-07-01 00:00:00
+  end_date: 2024-07-31 00:00:00
+  location:
+    type: centered_region
+    parameters:
+      latitude: -75.0
+      longitude: 15.0
+      bounding_box_degrees: 5
+  event_type: heat_wave
+- case_id_number: 79
+  title: August 2024 Canada
+  start_date: 2024-08-09 00:00:00
+  end_date: 2024-08-15 00:00:00
+  location:
+    type: centered_region
+    parameters:
+      latitude: 67.0
+      longitude: 248.0
+      bounding_box_degrees: 5
+  event_type: heat_wave
+- case_id_number: 80
+  title: August 2024 Australia
+  start_date: 2024-08-22 00:00:00
+  end_date: 2024-08-30 00:00:00
+  location:
+    type: centered_region
+    parameters:
+      latitude: -20.0
+      longitude: 120.0
+      bounding_box_degrees: 5
+  event_type: heat_wave
+- case_id_number: 81
+  title: July/August 2024 Japan
+  start_date: 2024-07-25 00:00:00
+  end_date: 2024-08-05 00:00:00
+  location:
+    type: centered_region
+    parameters:
+      latitude: 36.0
+      longitude: 138.0
+      bounding_box_degrees: 6
+  event_type: heat_wave
+- case_id_number: 82
+  title: June 2024 Saudi Arabia
+  start_date: 2024-06-16 00:00:00
+  end_date: 2024-06-22 00:00:00
+  location:
+    type: centered_region
+    parameters:
+      latitude: 24.0
+      longitude: 45.0
+      bounding_box_degrees: 6
+  event_type: heat_wave
+- case_id_number: 83
+  title: August 2024 Europe
+  start_date: 2024-08-10 00:00:00
+  end_date: 2024-08-16 00:00:00
+  location:
+    type: centered_region
+    parameters:
+      latitude: 48.3794
+      longitude: 10.8978
+      bounding_box_degrees: 10
+  event_type: heat_wave
+- case_id_number: 84
+  title: July 2024 Ukraine
+  start_date: 2024-07-12 00:00:00
+  end_date: 2024-07-18 00:00:00
+  location:
+    type: centered_region
+    parameters:
+      latitude: 50.4501
+      longitude: 30.5234
+      bounding_box_degrees: 5
+  event_type: heat_wave
+- case_id_number: 85
+  title: June 2024 Europe
+  start_date: 2024-06-20 00:00:00
+  end_date: 2024-06-30 00:00:00
+  location:
+    type: centered_region
+    parameters:
+      latitude: 52.3676
+      longitude: 4.9041
+      bounding_box_degrees: 6
+  event_type: heat_wave
+- case_id_number: 86
+  title: May 2024 Central Mexico
+  start_date: 2024-05-23 00:00:00
+  end_date: 2024-05-31 00:00:00
+  location:
+    type: centered_region
+    parameters:
+      latitude: 19.4326
+      longitude: 260.8668
+      bounding_box_degrees: 5
+  event_type: heat_wave
+- case_id_number: 87
+  title: May 2024 Pakistan/India
+  start_date: 2024-05-23 00:00:00
+  end_date: 2024-05-31 00:00:00
+  location:
+    type: centered_region
+    parameters:
+      latitude: 34.0
+      longitude: 76.0
+      bounding_box_degrees: 5
+  event_type: heat_wave
+- case_id_number: 88
+  title: August 2023 Chile
+  start_date: 2023-08-01 00:00:00
+  end_date: 2023-08-07 00:00:00
+  location:
+    type: centered_region
+    parameters:
+      latitude: -33.4489
+      longitude: 289.3307
+      bounding_box_degrees: 5
+  event_type: heat_wave
+- case_id_number: 89
+  title: January 2024 Pacific Northwest
+  start_date: 2024-01-11 12:00:00
+  end_date: 2024-01-20 00:00:00
+  location:
+    type: centered_region
+    parameters:
+      latitude: 47.6062
+      longitude: 237.6679
+      bounding_box_degrees: 5
+  event_type: freeze
+- case_id_number: 90
+  title: April 2022 Nevada
+  start_date: 2022-04-11 18:00:00
+  end_date: 2022-04-17 18:00:00
+  location:
+    type: centered_region
+    parameters:
+      latitude: 40.8324
+      longitude: 244.2369
+      bounding_box_degrees: 5
+  event_type: freeze
+- case_id_number: 91
+  title: February/March 2021 New England
+  start_date: 2021-03-04 18:00:00
+  end_date: 2021-03-10 18:00:00
+  location:
+    type: centered_region
+    parameters:
+      latitude: 42.8509
+      longitude: 287.4421
+      bounding_box_degrees: 5
+  event_type: freeze
+- case_id_number: 92
+  title: January 2024 Northern Europe
+  start_date: 2024-01-01 06:00:00
+  end_date: 2024-01-09 12:00:00
+  location:
+    type: centered_region
+    parameters:
+      latitude: 59.3293
+      longitude: 18.0686
+      bounding_box_degrees: 5
+  event_type: freeze
+- case_id_number: 93
+  title: December 2022 Europe
+  start_date: 2022-12-10 06:00:00
+  end_date: 2022-12-20 00:00:00
+  location:
+    type: centered_region
+    parameters:
+      latitude: 52.52
+      longitude: 13.405
+      bounding_box_degrees: 5
+  event_type: freeze
+- case_id_number: 94
+  title: April 2024 Europe
+  start_date: 2024-04-16 12:00:00
+  end_date: 2024-04-27 12:00:00
+  location:
+    type: centered_region
+    parameters:
+      latitude: 48.8566
+      longitude: 2.3522
+      bounding_box_degrees: 5
+  event_type: freeze
+- case_id_number: 95
+  title: January 2023 North China
+  start_date: 2024-01-20 00:00:00
+  end_date: 2024-02-04 18:00:00
+  location:
+    type: centered_region
+    parameters:
+      latitude: 42.5246
+      longitude: 122.3853
+      bounding_box_degrees: 5
+  event_type: freeze
+- case_id_number: 96
+  title: April 2024 Sweden
+  start_date: 2024-04-02 00:00:00
+  end_date: 2024-04-08 00:00:00
+  location:
+    type: centered_region
+    parameters:
+      latitude: 60.1282
+      longitude: 18.6435
+      bounding_box_degrees: 5
+  event_type: freeze
diff --git a/src/extremeweatherbench/data/events.yaml b/src/extremeweatherbench/data/events.yaml
index e1486e95..97a753b6 100644
--- a/src/extremeweatherbench/data/events.yaml
+++ b/src/extremeweatherbench/data/events.yaml
@@ -5,10 +5,10 @@
   location:
     type: bounded_region
     parameters:
-      latitude_min: 33.25
-      latitude_max: 68.0
-      longitude_min: 215.25
-      longitude_max: 256.0
+      latitude_min: 39.0
+      latitude_max: 65.0
+      longitude_min: -136.75
+      longitude_max: -105.75
   event_type: heat_wave
 - case_id_number: 2
   title: 2022 Upper Midwest
@@ -17,10 +17,10 @@
   location:
     type: bounded_region
     parameters:
-      latitude_min: 19.5
-      latitude_max: 48.25
-      longitude_min: 252.0
-      longitude_max: 276.75
+      latitude_min: 31.5
+      latitude_max: 47.5
+      longitude_min: -102.25
+      longitude_max: -82.25
   event_type: heat_wave
 - case_id_number: 3
   title: 2022 California
@@ -29,10 +29,10 @@
   location:
     type: bounded_region
     parameters:
-      latitude_min: 23.75
-      latitude_max: 44.5
-      longitude_min: 235.75
-      longitude_max: 266.25
+      latitude_min: 31.5
+      latitude_max: 42.5
+      longitude_min: -123.75
+      longitude_max: -104.75
   event_type: heat_wave
 - case_id_number: 4
   title: 2022 Texas
@@ -41,10 +41,10 @@
   location:
     type: bounded_region
     parameters:
-      latitude_min: 20.5
-      latitude_max: 43.25
-      longitude_min: 254.75
-      longitude_max: 275.5
+      latitude_min: 25.25
+      latitude_max: 37.25
+      longitude_min: -103.25
+      longitude_max: -90.25
   event_type: heat_wave
 - case_id_number: 5
   title: 2023 Pacific Northwest
@@ -53,10 +53,10 @@
   location:
     type: bounded_region
     parameters:
-      latitude_min: 35.25
-      latitude_max: 72.0
-      longitude_min: 219.25
-      longitude_max: 262.0
+      latitude_min: 42.0
+      latitude_max: 70.0
+      longitude_min: -130.75
+      longitude_max: -111.75
   event_type: heat_wave
 - case_id_number: 6
   title: 2022 Mid-Atlantic
@@ -65,10 +65,10 @@
   location:
     type: bounded_region
     parameters:
-      latitude_min: 31.0
-      latitude_max: 42.25
-      longitude_min: 277.0
-      longitude_max: 287.75
+      latitude_min: 33.75
+      latitude_max: 41.75
+      longitude_min: -83.0
+      longitude_max: -74.0
   event_type: heat_wave
 - case_id_number: 7
   title: 2023 Australia
@@ -77,9 +77,9 @@
   location:
     type: bounded_region
     parameters:
-      latitude_min: -35.0
-      latitude_max: -20.0
-      longitude_min: 113.5
+      latitude_min: -34.5
+      latitude_max: -23.5
+      longitude_min: 113.25
       longitude_max: 120.25
   event_type: heat_wave
 - case_id_number: 8
@@ -89,10 +89,10 @@
   location:
     type: bounded_region
     parameters:
-      latitude_min: 40.75
-      latitude_max: 67.5
-      longitude_min: -9.75
-      longitude_max: 16.75
+      latitude_min: 50.75
+      latitude_max: 55.75
+      longitude_min: -10.25
+      longitude_max: -5.25
   event_type: heat_wave
 - case_id_number: 9
   title: 2023 Italy
@@ -101,10 +101,10 @@
   location:
     type: bounded_region
     parameters:
-      latitude_min: 17.5
-      latitude_max: 52.25
-      longitude_min: -12.0
-      longitude_max: 32.75
+      latitude_min: 39.5
+      latitude_max: 52.5
+      longitude_min: 3.0
+      longitude_max: 17.0
   event_type: heat_wave
 - case_id_number: 10
   title: 2023 SW Europe
@@ -113,10 +113,10 @@
   location:
     type: bounded_region
     parameters:
-      latitude_min: 36.0
-      latitude_max: 56.75
-      longitude_min: -9.25
-      longitude_max: 20.75
+      latitude_min: 35.0
+      latitude_max: 50.0
+      longitude_min: -10.25
+      longitude_max: 18.75
   event_type: heat_wave
 - case_id_number: 11
   title: 2023 South America
@@ -125,10 +125,10 @@
   location:
     type: bounded_region
     parameters:
-      latitude_min: -44.25
-      latitude_max: -20.25
-      longitude_min: 285.75
-      longitude_max: 308.0
+      latitude_min: -34.5
+      latitude_max: -18.5
+      longitude_min: -73.0
+      longitude_max: -55.0
   event_type: heat_wave
 - case_id_number: 12
   title: 2023 China (Shanghai)
@@ -137,10 +137,10 @@
   location:
     type: bounded_region
     parameters:
-      latitude_min: 16.75
-      latitude_max: 35.25
-      longitude_min: 99.0
-      longitude_max: 122.0
+      latitude_min: 23.75
+      latitude_max: 33.75
+      longitude_min: 119.0
+      longitude_max: 124.0
   event_type: heat_wave
 - case_id_number: 13
   title: 2023 China (Yunnan Province)
@@ -149,10 +149,10 @@
   location:
     type: bounded_region
     parameters:
-      latitude_min: 12.75
-      latitude_max: 43.5
-      longitude_min: 92.25
-      longitude_max: 121.0
+      latitude_min: 14.5
+      latitude_max: 40.5
+      longitude_min: 96.25
+      longitude_max: 118.25
   event_type: heat_wave
 - case_id_number: 14
   title: 2023 Algeria
@@ -161,10 +161,10 @@
   location:
     type: bounded_region
     parameters:
-      latitude_min: 16.25
-      latitude_max: 51.0
-      longitude_min: -17.0
-      longitude_max: 27.5
+      latitude_min: 29.25
+      latitude_max: 39.25
+      longitude_min: -11.5
+      longitude_max: 11.5
   event_type: heat_wave
 - case_id_number: 15
   title: 2023 Iberian Peninsula
@@ -174,9 +174,9 @@
     type: bounded_region
     parameters:
       latitude_min: 22.75
-      latitude_max: 46.75
-      longitude_min: -15.75
-      longitude_max: 8.75
+      latitude_max: 44.75
+      longitude_min: -10.0
+      longitude_max: 9.0
   event_type: heat_wave
 - case_id_number: 16
   title: 2023 Southeast Asia
@@ -197,10 +197,10 @@
   location:
     type: bounded_region
     parameters:
-      latitude_min: 22.25
-      latitude_max: 53.0
-      longitude_min: 56.75
-      longitude_max: 101.5
+      latitude_min: 25.0
+      latitude_max: 37.0
+      longitude_min: 63.75
+      longitude_max: 93.75
   event_type: heat_wave
 - case_id_number: 18
   title: 2021 Western Russia
@@ -209,10 +209,10 @@
   location:
     type: bounded_region
     parameters:
-      latitude_min: 35.5
-      latitude_max: 65.75
-      longitude_min: 13.25
-      longitude_max: 60.0
+      latitude_min: 40.25
+      latitude_max: 65.25
+      longitude_min: 15.25
+      longitude_max: 56.0
   event_type: heat_wave
 - case_id_number: 19
   title: 2022 Germany
@@ -221,10 +221,10 @@
   location:
     type: bounded_region
     parameters:
-      latitude_min: 36.25
-      latitude_max: 56.25
-      longitude_min: -3.0
-      longitude_max: 21.75
+      latitude_min: 46.0
+      latitude_max: 55.0
+      longitude_min: -2.0
+      longitude_max: 17.0
   event_type: heat_wave
 - case_id_number: 20
   title: 2022 UK (August)
@@ -233,9 +233,9 @@
   location:
     type: bounded_region
     parameters:
-      latitude_min: 39.25
+      latitude_min: 41.0
       latitude_max: 56.0
-      longitude_min: -9.75
+      longitude_min: -6.75
       longitude_max: 8.25
   event_type: heat_wave
 - case_id_number: 21
@@ -245,10 +245,10 @@
   location:
     type: bounded_region
     parameters:
-      latitude_min: 37.25
-      latitude_max: 58.0
-      longitude_min: -9.75
-      longitude_max: 8.25
+      latitude_min: 37.0
+      latitude_max: 54.0
+      longitude_min: -7.75
+      longitude_max: 2.25
   event_type: heat_wave
 - case_id_number: 22
   title: 2022 France
@@ -257,10 +257,10 @@
   location:
     type: bounded_region
     parameters:
-      latitude_min: 29.0
-      latitude_max: 51.75
-      longitude_min: -9.75
-      longitude_max: 12.75
+      latitude_min: 31.0
+      latitude_max: 51.0
+      longitude_min: -9.0
+      longitude_max: 10.0
   event_type: heat_wave
 - case_id_number: 23
   title: 2022 Japan
@@ -269,10 +269,10 @@
   location:
     type: bounded_region
     parameters:
-      latitude_min: 31.25
-      latitude_max: 46.0
-      longitude_min: 125.25
-      longitude_max: 145.25
+      latitude_min: 33.25
+      latitude_max: 42.25
+      longitude_min: 126.25
+      longitude_max: 142.25
   event_type: heat_wave
 - case_id_number: 24
   title: 2022 India
@@ -281,10 +281,10 @@
   location:
     type: bounded_region
     parameters:
-      latitude_min: 18.0
-      latitude_max: 52.75
+      latitude_min: 20.75
+      latitude_max: 41.75
       longitude_min: 58.0
-      longitude_max: 78.75
+      longitude_max: 76.0
   event_type: heat_wave
 - case_id_number: 25
   title: 2022 East Antarctica
@@ -294,9 +294,9 @@
     type: bounded_region
     parameters:
       latitude_min: -85.5
-      latitude_max: -65.5
-      longitude_min: 99.0
-      longitude_max: 147.75
+      latitude_max: -64.5
+      longitude_min: 101.0
+      longitude_max: 145.75
   event_type: heat_wave
 - case_id_number: 26
   title: 2022 West Australia
@@ -305,10 +305,10 @@
   location:
     type: bounded_region
     parameters:
-      latitude_min: -35.0
-      latitude_max: -27.5
-      longitude_min: 114.25
-      longitude_max: 119.75
+      latitude_min: -35.5
+      latitude_max: -28.5
+      longitude_min: 113.25
+      longitude_max: 119.25
   event_type: heat_wave
 - case_id_number: 27
   title: 2021 Canada Plains
@@ -317,10 +317,10 @@
   location:
     type: bounded_region
     parameters:
-      latitude_min: 38.5
-      latitude_max: 59.5
-      longitude_min: 238.0
-      longitude_max: 272.75
+      latitude_min: 39.5
+      latitude_max: 53.5
+      longitude_min: -120.0
+      longitude_max: -87.0
   event_type: heat_wave
 - case_id_number: 28
   title: 2021 New Zealand
@@ -330,9 +330,9 @@
     type: bounded_region
     parameters:
       latitude_min: -46.5
-      latitude_max: -40.75
-      longitude_min: 167.5
-      longitude_max: 175.25
+      latitude_max: -41.5
+      longitude_min: 149.0
+      longitude_max: 149.0
   event_type: heat_wave
 - case_id_number: 29
   title: 2020 Australia
@@ -341,10 +341,10 @@
   location:
     type: bounded_region
     parameters:
-      latitude_min: -36.75
-      latitude_max: -17.5
-      longitude_min: 126.5
-      longitude_max: 152.5
+      latitude_min: -36.25
+      latitude_max: -22.25
+      longitude_min: 128.5
+      longitude_max: 153.5
   event_type: heat_wave
 - case_id_number: 30
   title: 2021 Texas
@@ -353,10 +353,10 @@
   location:
     type: bounded_region
     parameters:
-      latitude_min: 24.0
-      latitude_max: 54.75
-      longitude_min: 250.0
-      longitude_max: 278.75
+      latitude_min: 23.75
+      latitude_max: 52.75
+      longitude_min: -107.25
+      longitude_max: -82.25
   event_type: freeze
 - case_id_number: 31
   title: 2022 Arkansas
@@ -365,10 +365,10 @@
   location:
     type: bounded_region
     parameters:
-      latitude_min: 24.25
-      latitude_max: 59.0
-      longitude_min: 243.25
-      longitude_max: 272.0
+      latitude_min: 30.25
+      latitude_max: 57.0
+      longitude_min: -114.75
+      longitude_max: -88.75
   event_type: freeze
 - case_id_number: 32
   title: 2023 Germany
@@ -377,10 +377,10 @@
   location:
     type: bounded_region
     parameters:
-      latitude_min: 44.25
-      latitude_max: 68.5
-      longitude_min: 7.25
-      longitude_max: 22.0
+      latitude_min: 45.75
+      latitude_max: 68.75
+      longitude_min: 9.0
+      longitude_max: 14.0
   event_type: freeze
 - case_id_number: 33
   title: 2023 China
@@ -389,10 +389,10 @@
   location:
     type: bounded_region
     parameters:
-      latitude_min: 14.75
-      latitude_max: 55.5
-      longitude_min: 101.0
-      longitude_max: 141.75
+      latitude_min: 26.75
+      latitude_max: 53.5
+      longitude_min: 100.0
+      longitude_max: 124.0
   event_type: freeze
 - case_id_number: 34
   title: 2023 Afghanistan
@@ -401,10 +401,10 @@
   location:
     type: bounded_region
     parameters:
-      latitude_min: 22.25
-      latitude_max: 49.0
+      latitude_min: 32.0
+      latitude_max: 47.0
       longitude_min: 50.75
-      longitude_max: 85.5
+      longitude_max: 71.75
   event_type: freeze
 - case_id_number: 35
   title: 2022 Colorado
@@ -413,10 +413,10 @@
   location:
     type: bounded_region
     parameters:
-      latitude_min: 35.25
-      latitude_max: 64.0
-      longitude_min: 236.75
-      longitude_max: 261.5
+      latitude_min: 37.25
+      latitude_max: 42.25
+      longitude_min: -107.5
+      longitude_max: -102.5
   event_type: freeze
 - case_id_number: 36
   title: July 2024 South Dakota
@@ -857,10 +857,10 @@
   location:
     type: bounded_region
     parameters:
-      latitude_min: 30.0
-      latitude_max: 42.75
-      longitude_min: 235.5
-      longitude_max: 248.75
+      latitude_min: 14.5
+      latitude_max: 31.5
+      longitude_min: -105.0
+      longitude_max: -95.0
   event_type: heat_wave
 - case_id_number: 73
   title: July 2024 Northeast US
@@ -869,10 +869,10 @@
   location:
     type: bounded_region
     parameters:
-      latitude_min: 35.75
-      latitude_max: 47.5
-      longitude_min: 281.0
-      longitude_max: 300.0
+      latitude_min: 39.0
+      latitude_max: 48.0
+      longitude_min: -82.25
+      longitude_max: -69.25
   event_type: heat_wave
 - case_id_number: 74
   title: July 2024 Mid-Atlantic US
@@ -881,10 +881,10 @@
   location:
     type: bounded_region
     parameters:
-      latitude_min: 35.25
-      latitude_max: 44.75
-      longitude_min: 280.0
-      longitude_max: 292.75
+      latitude_min: 30.75
+      latitude_max: 37.75
+      longitude_min: -124.25
+      longitude_max: -113.25
   event_type: heat_wave
 - case_id_number: 75
   title: August 2024 Midwest US
@@ -893,10 +893,10 @@
   location:
     type: bounded_region
     parameters:
-      latitude_min: 36.0
-      latitude_max: 46.25
-      longitude_min: 265.5
-      longitude_max: 278.25
+      latitude_min: 37.75
+      latitude_max: 52.75
+      longitude_min: -78.0
+      longitude_max: -61.0
   event_type: heat_wave
 - case_id_number: 76
   title: July 2024 Antarctica
@@ -905,10 +905,10 @@
   location:
     type: bounded_region
     parameters:
-      latitude_min: -90.0
-      latitude_max: -68.75
-      longitude_min: -9.5
-      longitude_max: 39.5
+      latitude_min: 37.0
+      latitude_max: 43.0
+      longitude_min: -81.25
+      longitude_max: -68.25
   event_type: heat_wave
 - case_id_number: 77
   title: August 2024 Canada
@@ -917,10 +917,10 @@
   location:
     type: bounded_region
     parameters:
-      latitude_min: 54.5
-      latitude_max: 73.5
-      longitude_min: 231.5
-      longitude_max: 272.5
+      latitude_min: 37.5
+      latitude_max: 44.5
+      longitude_min: -90.75
+      longitude_max: -85.75
   event_type: heat_wave
 - case_id_number: 78
   title: August 2024 Australia
@@ -929,10 +929,10 @@
   location:
     type: bounded_region
     parameters:
-      latitude_min: -32.5
-      latitude_max: -9.5
-      longitude_min: 116.25
-      longitude_max: 144.5
+      latitude_min: -83.5
+      latitude_max: -68.5
+      longitude_min: 4.5
+      longitude_max: 37.5
   event_type: heat_wave
 - case_id_number: 79
   title: July/August 2024 Japan
@@ -941,10 +941,10 @@
   location:
     type: bounded_region
     parameters:
-      latitude_min: 31.25
-      latitude_max: 43.0
-      longitude_min: 123.0
-      longitude_max: 142.75
+      latitude_min: 58.5
+      latitude_max: 76.5
+      longitude_min: -133.5
+      longitude_max: -89.5
   event_type: heat_wave
 - case_id_number: 80
   title: June 2024 Saudi Arabia
@@ -953,10 +953,10 @@
   location:
     type: bounded_region
     parameters:
-      latitude_min: 19.0
-      latitude_max: 35.0
-      longitude_min: 30.0
-      longitude_max: 54.0
+      latitude_min: -27.5
+      latitude_max: -17.5
+      longitude_min: 114.5
+      longitude_max: 142.5
   event_type: heat_wave
 - case_id_number: 81
   title: August 2024 Europe
@@ -965,10 +965,10 @@
   location:
     type: bounded_region
     parameters:
-      latitude_min: 35.5
-      latitude_max: 55.25
-      longitude_min: 0.0
-      longitude_max: 25.75
+      latitude_min: 33.0
+      latitude_max: 39.0
+      longitude_min: 132.0
+      longitude_max: 141.0
   event_type: heat_wave
 - case_id_number: 82
   title: July 2024 Ukraine
@@ -977,10 +977,10 @@
   location:
     type: bounded_region
     parameters:
-      latitude_min: 40.0
-      latitude_max: 56.75
-      longitude_min: 14.25
-      longitude_max: 45.0
+      latitude_min: 18.0
+      latitude_max: 36.0
+      longitude_min: 22.0
+      longitude_max: 55.0
   event_type: heat_wave
 - case_id_number: 83
   title: June 2024 Europe
@@ -989,10 +989,10 @@
   location:
     type: bounded_region
     parameters:
-      latitude_min: 45.5
-      latitude_max: 63.25
-      longitude_min: -4.0
-      longitude_max: 25.75
+      latitude_min: 34.5
+      latitude_max: 53.5
+      longitude_min: 3.0
+      longitude_max: 22.0
   event_type: heat_wave
 - case_id_number: 84
   title: May 2024 Central Mexico
@@ -1001,10 +1001,10 @@
   location:
     type: bounded_region
     parameters:
-      latitude_min: 7.25
-      latitude_max: 29.75
-      longitude_min: 252.5
-      longitude_max: 283.25
+      latitude_min: 39.0
+      latitude_max: 55.0
+      longitude_min: 15.0
+      longitude_max: 39.0
   event_type: heat_wave
 - case_id_number: 85
   title: May 2024 Pakistan/India
@@ -1013,10 +1013,10 @@
   location:
     type: bounded_region
     parameters:
-      latitude_min: 23.5
-      latitude_max: 42.5
-      longitude_min: 71.5
-      longitude_max: 98.5
+      latitude_min: 48.25
+      latitude_max: 55.25
+      longitude_min: 2.0
+      longitude_max: 9.0
   event_type: heat_wave
 - case_id_number: 86
   title: August 2023 Chile
@@ -1025,10 +1025,10 @@
   location:
     type: bounded_region
     parameters:
-      latitude_min: -37.25
-      latitude_max: -19.0
-      longitude_min: 288.5
-      longitude_max: 305.75
+      latitude_min: 14.0
+      latitude_max: 31.0
+      longitude_min: -104.75
+      longitude_max: -82.75
   event_type: heat_wave
 - case_id_number: 87
   title: January 2024 Pacific Northwest
@@ -1037,10 +1037,10 @@
   location:
     type: bounded_region
     parameters:
-      latitude_min: 41.25
-      latitude_max: 66.0
-      longitude_min: 221.25
-      longitude_max: 262.0
+      latitude_min: 31.5
+      latitude_max: 36.5
+      longitude_min: 73.5
+      longitude_max: 78.5
   event_type: freeze
 - case_id_number: 88
   title: April 2022 Nevada
@@ -1049,10 +1049,10 @@
   location:
     type: bounded_region
     parameters:
-      latitude_min: 36.5
-      latitude_max: 65.25
-      longitude_min: 225.75
-      longitude_max: 268.5
+      latitude_min: -42.0
+      latitude_max: -20.0
+      longitude_min: -73.25
+      longitude_max: -62.25
   event_type: freeze
 - case_id_number: 89
   title: February/March 2021 New England
@@ -1061,10 +1061,10 @@
   location:
     type: bounded_region
     parameters:
-      latitude_min: 38.5
-      latitude_max: 47.0
-      longitude_min: 283.25
-      longitude_max: 291.5
+      latitude_min: 45.0
+      latitude_max: 69.0
+      longitude_min: -133.75
+      longitude_max: -100.0
   event_type: freeze
 - case_id_number: 90
   title: January 2024 Northern Europe
@@ -1073,10 +1073,10 @@
   location:
     type: bounded_region
     parameters:
-      latitude_min: 51.25
-      latitude_max: 70.5
-      longitude_min: 5.25
-      longitude_max: 42.5
+      latitude_min: 38.25
+      latitude_max: 63.25
+      longitude_min: -118.25
+      longitude_max: -113.25
   event_type: freeze
 - case_id_number: 91
   title: December 2022 Europe
@@ -1085,10 +1085,10 @@
   location:
     type: bounded_region
     parameters:
-      latitude_min: 46.25
-      latitude_max: 70.75
-      longitude_min: -9.75
-      longitude_max: 25.75
+      latitude_min: 40.25
+      latitude_max: 46.25
+      longitude_min: -75.0
+      longitude_max: -68.0
   event_type: freeze
 - case_id_number: 92
   title: April 2024 Europe
@@ -1097,10 +1097,10 @@
   location:
     type: bounded_region
     parameters:
-      latitude_min: 40.5
-      latitude_max: 71.0
-      longitude_min: -4.0
-      longitude_max: 26.75
+      latitude_min: 56.75
+      latitude_max: 71.75
+      longitude_min: 3.5
+      longitude_max: 20.5
   event_type: freeze
 - case_id_number: 93
   title: January 2023 North China
@@ -1109,10 +1109,10 @@
   location:
     type: bounded_region
     parameters:
-      latitude_min: 38.25
-      latitude_max: 47.0
-      longitude_min: 118.0
-      longitude_max: 130.75
+      latitude_min: 48.0
+      latitude_max: 55.0
+      longitude_min: -9.0
+      longitude_max: 16.0
   event_type: freeze
 - case_id_number: 94
   title: April 2024 Sweden
@@ -1121,10 +1121,10 @@
   location:
     type: bounded_region
     parameters:
-      latitude_min: 58.25
-      latitude_max: 71.0
-      longitude_min: 8.25
-      longitude_max: 31.0
+      latitude_min: 42.25
+      latitude_max: 54.25
+      longitude_min: -3.25
+      longitude_max: 23.75
   event_type: freeze
 - case_id_number: 95
   title: November 2024 West Coast US
@@ -1133,10 +1133,10 @@
   location:
     type: bounded_region
     parameters:
-      latitude_min: 27.8
-      latitude_max: 54.2
-      longitude_min: 222.4
-      longitude_max: 245.3
+      latitude_min: 40.0
+      latitude_max: 45.0
+      longitude_min: 120.0
+      longitude_max: 125.0
   event_type: atmospheric_river
 - case_id_number: 96
   title: October 2024 British Columbia
@@ -1145,10 +1145,10 @@
   location:
     type: bounded_region
     parameters:
-      latitude_min: 32.8
-      latitude_max: 56.0
-      longitude_min: 217.3
-      longitude_max: 247.0
+      latitude_min: 57.75
+      latitude_max: 71.75
+      longitude_min: 10.25
+      longitude_max: 21.25
   event_type: atmospheric_river
 - case_id_number: 97
   title: September 2024 SW Alaska and British Columbia

From 8328fe8e58d5e2ab8cf1a43c7af6030288f43bc1 Mon Sep 17 00:00:00 2001
From: aaTman <mandelbaum.taylor@gmail.com>
Date: Tue, 7 Apr 2026 21:50:28 +0000
Subject: [PATCH 02/12] add defaults; remove top python line

---
 data_prep/heat_cold_bounds_case.py   | 16 +++++++++-------
 data_prep/heat_cold_bounds_global.py | 22 ++++++++++++++--------
 2 files changed, 23 insertions(+), 15 deletions(-)

diff --git a/data_prep/heat_cold_bounds_case.py b/data_prep/heat_cold_bounds_case.py
index 74153e82..7de458d3 100644
--- a/data_prep/heat_cold_bounds_case.py
+++ b/data_prep/heat_cold_bounds_case.py
@@ -1,4 +1,3 @@
-#!/usr/bin/env python3
 """Validate and expand heat wave / freeze bounding boxes.
 
 Reads heat_wave and freeze cases from base_temp_events.yaml (which uses
@@ -78,9 +77,10 @@ def _apply_consecutive_filter(
 
     Args:
         mask: Boolean array of shape (time, lat, lon).
-        min_days: Minimum run length required to qualify as an event.
+        min_days: Minimum run length required to qualify as an
+            event. Default is 3 (MIN_CONSECUTIVE_DAYS).
         max_grace_days: Maximum gap length to bridge after the
-            minimum run is established.
+            minimum run is established. Default is 1.
 
     Returns:
         Boolean array of the same shape with only qualifying runs
@@ -285,11 +285,13 @@ def process_event(
     Args:
         single_case: The individual case to process; must have a
             CenteredRegion location.
-        out_dir: Directory in which plots are saved.
-        quantile: Climatology quantile. Defaults to 0.85 for
-            heat_wave, 0.15 for freeze.
+        out_dir: Directory in which plots are saved. Default is
+            pathlib.Path(".").
+        quantile: Climatology quantile. Default is None, which
+            resolves to 0.85 for heat_wave or 0.15 for freeze.
         op_str: Comparison operator string (e.g. ">", ">=", "<",
-            "<="). Defaults to ">" for heat_wave, "<" for freeze.
+            "<="). Default is None, which resolves to ">" for
+            heat_wave or "<" for freeze.
 
     Returns:
         A dict with keys event_type, start_date, end_date,
diff --git a/data_prep/heat_cold_bounds_global.py b/data_prep/heat_cold_bounds_global.py
index c5f54d8f..ffb56be9 100644
--- a/data_prep/heat_cold_bounds_global.py
+++ b/data_prep/heat_cold_bounds_global.py
@@ -1,4 +1,3 @@
-#!/usr/bin/env python3
 """Detect heat waves and cold snaps globally from ERA5 reanalysis.
 
 Scans ERA5 2m temperature over an input date range and identifies
@@ -65,12 +64,15 @@ def get_climatology_thresholds(
 
     Args:
         q_hw: Lower-bound quantile for heat wave detection.
+            Default is 0.85.
         q_fz: Upper-bound quantile for freeze detection.
+            Default is 0.15.
         q_hw_upper: Upper-bound quantile for heat waves. When set,
             only days where temp > q_hw AND temp < q_hw_upper are
-            flagged.
+            flagged. Default is None.
         q_fz_lower: Lower-bound quantile for freezes. When set, only
-            days where temp < q_fz AND temp > q_fz_lower are flagged.
+            days where temp < q_fz AND temp > q_fz_lower are
+            flagged. Default is None.
 
     Returns:
         A tuple of (clim_hw, clim_fz, clim_hw_upper, clim_fz_lower).
@@ -140,11 +142,13 @@ def build_exceedance_masks(
             (dayofyear, hour).
         land_mask: Boolean DataArray (True = land) matching t2m grid.
         op_hw: Comparison operator string for heat waves.
+            Default is ">".
         op_fz: Comparison operator string for freezes.
+            Default is "<".
         clim_hw_upper: Optional upper-bound climatology for heat
-            waves (exclusive cap).
+            waves (exclusive cap). Default is None.
         clim_fz_lower: Optional lower-bound climatology for freezes
-            (exclusive floor).
+            (exclusive floor). Default is None.
 
     Returns:
         A tuple (hw, fz) of daily boolean DataArrays masked to land,
@@ -200,9 +204,10 @@ def apply_consecutive_filter(
 
     Args:
         mask: Boolean array of shape (time, lat, lon).
-        min_days: Minimum run length required to qualify as an event.
+        min_days: Minimum run length required to qualify as an
+            event. Default is 3 (MIN_CONSECUTIVE_DAYS).
         max_grace_days: Maximum gap length to bridge after the
-            minimum run is established.
+            minimum run is established. Default is 1.
 
     Returns:
         Boolean array of the same shape with only qualifying runs
@@ -488,7 +493,8 @@ def events_to_dataframe(
     Args:
         events: Raw event dicts from ``detect_events``.
         min_gridpoints: Drop events whose peak spatial extent (in
-            grid points) is below this threshold.
+            grid points) is below this threshold. Default is 500
+            (MIN_GRIDPOINTS).
 
     Returns:
         DataFrame with columns label, event_type, start_date,

From 5cb96d57c1e7a00f0f91e591508bc503bda54ba9 Mon Sep 17 00:00:00 2001
From: aaTman <mandelbaum.taylor@gmail.com>
Date: Wed, 8 Apr 2026 00:11:58 +0000
Subject: [PATCH 03/12] add lat min and max args

---
 data_prep/heat_cold_bounds_global.py | 25 +++++++++++++++++++++++++
 1 file changed, 25 insertions(+)

diff --git a/data_prep/heat_cold_bounds_global.py b/data_prep/heat_cold_bounds_global.py
index ffb56be9..5df61202 100644
--- a/data_prep/heat_cold_bounds_global.py
+++ b/data_prep/heat_cold_bounds_global.py
@@ -617,6 +617,24 @@ def main():
         default="<",
         help="Comparison operator for freezes (default: <)",
     )
+    parser.add_argument(
+        "--lat-min",
+        type=float,
+        default=-90.0,
+        help=(
+            "Minimum latitude to include in detection. "
+            "Default is -90.0"
+        ),
+    )
+    parser.add_argument(
+        "--lat-max",
+        type=float,
+        default=90.0,
+        help=(
+            "Maximum latitude to include in detection. "
+            "Default is 90.0"
+        ),
+    )
     args = parser.parse_args()
 
     wall_start = time_module.time()
@@ -634,6 +652,13 @@ def main():
 
     logger.info("Opening ERA5 data...")
     t2m = open_era5_t2m(args.start_date, args.end_date)
+    if args.lat_min != -90.0 or args.lat_max != 90.0:
+        t2m = t2m.sel(latitude=slice(args.lat_min, args.lat_max))
+        logger.info(
+            "  Latitude filtered to [%.1f, %.1f]",
+            args.lat_min,
+            args.lat_max,
+        )
     logger.info("  sizes=%s", dict(t2m.sizes))
 
     logger.info("Loading climatology thresholds...")

From 996510327684b0ecfbd2a2b8264b76be0ebd4743 Mon Sep 17 00:00:00 2001
From: aaTman <mandelbaum.taylor@gmail.com>
Date: Wed, 8 Apr 2026 00:14:58 +0000
Subject: [PATCH 04/12] ruff + mypy

---
 data_prep/heat_cold_bounds_case.py   | 110 +++++++++++----------------
 data_prep/heat_cold_bounds_global.py |  60 +++++++--------
 data_prep/plot_temperature_events.py |  16 ++--
 3 files changed, 80 insertions(+), 106 deletions(-)

diff --git a/data_prep/heat_cold_bounds_case.py b/data_prep/heat_cold_bounds_case.py
index 7de458d3..dba02627 100644
--- a/data_prep/heat_cold_bounds_case.py
+++ b/data_prep/heat_cold_bounds_case.py
@@ -23,7 +23,7 @@
 import logging
 import pathlib
 import time as time_module
-from typing import Dict, List, Optional
+from typing import Dict, List, Literal, Optional, cast
 
 import cartopy.crs as ccrs
 import cartopy.feature as cfeature
@@ -93,9 +93,12 @@ def _apply_consecutive_filter(
     strict = (
         ndimage.binary_dilation(
             ndimage.binary_erosion(
-                mask, structure=struct, border_value=False,
+                mask,
+                structure=struct,
+                border_value=False,
             ),
-            structure=struct, border_value=False,
+            structure=struct,
+            border_value=False,
         )
         & mask
     )
@@ -107,16 +110,21 @@ def _apply_consecutive_filter(
     close_k = np.zeros((2 * max_grace_days + 1, 1, 1), dtype=bool)
     close_k[:, 0, 0] = True
     filled = ndimage.binary_closing(
-        mask, structure=close_k, border_value=False,
+        mask,
+        structure=close_k,
+        border_value=False,
     )
 
     # Runs of min_days+ in the gap-filled mask
     filled_runs = (
         ndimage.binary_dilation(
             ndimage.binary_erosion(
-                filled, structure=struct, border_value=False,
+                filled,
+                structure=struct,
+                border_value=False,
             ),
-            structure=struct, border_value=False,
+            structure=struct,
+            border_value=False,
         )
         & filled
     )
@@ -345,9 +353,8 @@ def process_event(
 
     # Time range and spatial pre-fetch extent
     start_date = pd.Timestamp(single_case.start_date) - pd.Timedelta(days=3)
-    end_date = (
-        pd.Timestamp(single_case.end_date)
-        + pd.Timedelta(days=TEMPORAL_LOAD_BUFFER_DAYS)
+    end_date = pd.Timestamp(single_case.end_date) + pd.Timedelta(
+        days=TEMPORAL_LOAD_BUFFER_DAYS
     )
 
     pot_lat_min = max(
@@ -427,9 +434,7 @@ def process_event(
 
     exc_6h = cmp(t2m, clim_aligned)
 
-    daily_all_pass = (
-        exc_6h.resample({tdim: "1D"}).min().astype(bool)
-    )
+    daily_all_pass = exc_6h.resample({tdim: "1D"}).min().astype(bool)
 
     land_reg = regionmask.defined_regions.natural_earth_v5_0_0.land_110
     land_mask = (
@@ -447,9 +452,7 @@ def process_event(
 
     all_lats = daily_all_pass.latitude.values
     all_lons = daily_all_pass.longitude.values
-    grid_res = (
-        np.abs(np.diff(all_lats[:2]))[0] if len(all_lats) > 1 else 0.25
-    )
+    grid_res = np.abs(np.diff(all_lats[:2]))[0] if len(all_lats) > 1 else 0.25
     band_pts = max(1, int(round(EXPANSION_DEGREES / grid_res)))
 
     def _lat_idx(val: float) -> int:
@@ -497,17 +500,14 @@ def _lon_idx(val: float) -> int:
             slice(idx_w0, idx_e0 + 1),
         )
         currently_active = (
-            mask_np[t][box_sl]
-            & mask_np[t - 1][box_sl]
-            & mask_np[t - 2][box_sl]
+            mask_np[t][box_sl] & mask_np[t - 1][box_sl] & mask_np[t - 2][box_sl]
         )
         active_land = int((currently_active & box_land).sum())
         frac = active_land / n_land_pts
 
         if t == 2 and frac < EDGE_VALIDITY_THRESHOLD:
             logger.warning(
-                "  Day 3: only %.1f%% of land points have"
-                " 3 consecutive days (< 50%%)",
+                "  Day 3: only %.1f%% of land points have 3 consecutive days (< 50%%)",
                 frac * 100,
             )
 
@@ -515,8 +515,7 @@ def _lon_idx(val: float) -> int:
             if frac >= EDGE_VALIDITY_THRESHOLD:
                 established = True
                 logger.info(
-                    "  Event established at day %d"
-                    " (%.1f%% of land points active)",
+                    "  Event established at day %d (%.1f%% of land points active)",
                     t,
                     frac * 100,
                 )
@@ -524,8 +523,7 @@ def _lon_idx(val: float) -> int:
             if frac < EDGE_VALIDITY_THRESHOLD:
                 final_t = t + 1
                 logger.info(
-                    "  Temporal stop at day %d"
-                    " (%.1f%% < 50%% of land points)",
+                    "  Temporal stop at day %d (%.1f%% < 50%% of land points)",
                     t,
                     frac * 100,
                 )
@@ -539,8 +537,7 @@ def _lon_idx(val: float) -> int:
 
     if not established:
         logger.warning(
-            "  Case %d: event never reached 50%% of land"
-            " points — using all %d days",
+            "  Case %d: event never reached 50%% of land points — using all %d days",
             single_case.case_id_number,
             final_t,
         )
@@ -567,8 +564,7 @@ def _lon_idx(val: float) -> int:
         peak_day = int(tied_days[(n_tied - 1) // 2])
     peak_mask = filtered[peak_day]
     logger.info(
-        "  Peak footprint on day %d (%d active grid points,"
-        " %d tied days)",
+        "  Peak footprint on day %d (%d active grid points, %d tied days)",
         peak_day,
         int(max_count),
         n_tied,
@@ -581,9 +577,7 @@ def _lon_idx(val: float) -> int:
 
     # Pre-check: disable edges that are >= 95% ocean
     # in the initial box (prevents expansion through water).
-    init_region_land = land_mask_np[
-        idx_s0 : idx_n0 + 1, idx_w0 : idx_e0 + 1
-    ]
+    init_region_land = land_mask_np[idx_s0 : idx_n0 + 1, idx_w0 : idx_e0 + 1]
     edges_active = {}
     for edge in ("north", "south", "east", "west"):
         if edge == "north":
@@ -598,8 +592,7 @@ def _lon_idx(val: float) -> int:
         edges_active[edge] = land_frac >= 0.25
         if not edges_active[edge]:
             logger.info(
-                "  %s edge disabled (%.1f%% land"
-                " in initial box)",
+                "  %s edge disabled (%.1f%% land in initial box)",
                 edge.capitalize(),
                 land_frac * 100,
             )
@@ -608,9 +601,7 @@ def _lon_idx(val: float) -> int:
     for iteration in range(MAX_ITERATIONS):
         n_iter = iteration + 1
         region = peak_mask[idx_s : idx_n + 1, idx_w : idx_e + 1]
-        land_region = land_mask_np[
-            idx_s : idx_n + 1, idx_w : idx_e + 1
-        ]
+        land_region = land_mask_np[idx_s : idx_n + 1, idx_w : idx_e + 1]
 
         all_done = True
         for edge in list(edges_active.keys()):
@@ -629,13 +620,15 @@ def _lon_idx(val: float) -> int:
                 all_done = False
                 if edge == "north":
                     idx_n = min(
-                        len(all_lats) - 1, idx_n + band_pts,
+                        len(all_lats) - 1,
+                        idx_n + band_pts,
                     )
                 elif edge == "south":
                     idx_s = max(0, idx_s - band_pts)
                 elif edge == "east":
                     idx_e = min(
-                        len(all_lons) - 1, idx_e + band_pts,
+                        len(all_lons) - 1,
+                        idx_e + band_pts,
                     )
                 elif edge == "west":
                     idx_w = max(0, idx_w - band_pts)
@@ -656,9 +649,7 @@ def _lon_idx(val: float) -> int:
     fin_lons = all_lons[idx_w : idx_e + 1]
     fin_filtered = filtered[:, idx_s : idx_n + 1, idx_w : idx_e + 1]
     consec = max_consecutive_days(fin_filtered)
-    peak_gridpoints = int(
-        (consec >= MIN_CONSECUTIVE_DAYS).sum()
-    )
+    peak_gridpoints = int((consec >= MIN_CONSECUTIVE_DAYS).sum())
 
     result = {
         "case_id": single_case.case_id_number,
@@ -678,8 +669,7 @@ def _lon_idx(val: float) -> int:
         "_peak_gridpoints": peak_gridpoints,
     }
     logger.info(
-        "  Final bounds: lat [%.2f, %.2f], lon [%.2f, %.2f]"
-        " (%d days, %d iterations)",
+        "  Final bounds: lat [%.2f, %.2f], lon [%.2f, %.2f] (%d days, %d iterations)",
         result["latitude_min"],
         result["latitude_max"],
         result["longitude_min"],
@@ -691,14 +681,12 @@ def _lon_idx(val: float) -> int:
     kind = "heatwave" if is_heatwave else "freeze"
     start = result["start_date"][:10]
     end = result["end_date"][:10]
-    out_png = str(
-        out_dir / f"case_{result['case_id']}_consecutive_{kind}_days.png"
-    )
+    out_png = str(out_dir / f"case_{result['case_id']}_consecutive_{kind}_days.png")
     plot_consecutive_map(
         consec,
         fin_lats,
         fin_lons,
-        single_case.event_type,
+        cast(Literal["heat_wave", "cold_snap"], single_case.event_type),
         title=(
             f"Consecutive {kind.capitalize()} Days"
             f" — {result['title']}\n{start} to {end}"
@@ -714,10 +702,7 @@ def _lon_idx(val: float) -> int:
         result["longitude_min"],
         result["longitude_max"],
     )
-    peak_png = str(
-        out_dir
-        / f"case_{result['case_id']}_peak_day_{kind}.png"
-    )
+    peak_png = str(out_dir / f"case_{result['case_id']}_peak_day_{kind}.png")
     plot_peak_day_with_bounds(
         peak_mask,
         all_lats,
@@ -726,8 +711,7 @@ def _lon_idx(val: float) -> int:
         final_box,
         single_case.event_type,
         title=(
-            f"Peak Footprint (day {peak_day})"
-            f" — {result['title']}\n{start} to {end}"
+            f"Peak Footprint (day {peak_day}) — {result['title']}\n{start} to {end}"
         ),
         output_path=peak_png,
     )
@@ -808,13 +792,9 @@ def write_bounds_to_yaml(
             entries are ignored).
         yaml_path: Path to the events.yaml file to update in-place.
     """
-    result_map = {
-        r["case_id"]: r for r in results if r is not None
-    }
+    result_map = {r["case_id"]: r for r in results if r is not None}
     if not result_map:
-        logger.warning(
-            "write_bounds_to_yaml: no valid results, skipping."
-        )
+        logger.warning("write_bounds_to_yaml: no valid results, skipping.")
         return
 
     yaml = YAML(typ="rt")
@@ -921,12 +901,11 @@ def main():
 
     old_cases = _load_base_temp_events()
     hw_fz = [
-        e for e in old_cases
+        e
+        for e in old_cases
         if e.event_type in ("heat_wave", "freeze")
-        and (args.case_min is None
-             or e.case_id_number >= args.case_min)
-        and (args.case_max is None
-             or e.case_id_number <= args.case_max)
+        and (args.case_min is None or e.case_id_number >= args.case_min)
+        and (args.case_max is None or e.case_id_number <= args.case_max)
     ]
     logger.info(
         "Found %d heat_wave / freeze events in base_temp_events.yaml",
@@ -937,7 +916,8 @@ def main():
 
     results = joblib.Parallel(n_jobs=args.n_workers)(
         joblib.delayed(process_event)(
-            c, out_dir,
+            c,
+            out_dir,
             quantile=args.quantile,
             op_str=args.operator,
         )
diff --git a/data_prep/heat_cold_bounds_global.py b/data_prep/heat_cold_bounds_global.py
index 5df61202..442f66b9 100644
--- a/data_prep/heat_cold_bounds_global.py
+++ b/data_prep/heat_cold_bounds_global.py
@@ -164,10 +164,12 @@ def build_exceedance_masks(
     doy_capped = doy.clip(max=max_clim_doy)
 
     clim_hw_aligned = clim_hw.sel(
-        dayofyear=doy_capped, hour=hour,
+        dayofyear=doy_capped,
+        hour=hour,
     ).reindex_like(t2m, method="nearest")
     clim_fz_aligned = clim_fz.sel(
-        dayofyear=doy_capped, hour=hour,
+        dayofyear=doy_capped,
+        hour=hour,
     ).reindex_like(t2m, method="nearest")
 
     hw_6h = cmp_hw(t2m, clim_hw_aligned)
@@ -175,13 +177,15 @@ def build_exceedance_masks(
 
     if clim_hw_upper is not None:
         clim_hw_upper_aligned = clim_hw_upper.sel(
-            dayofyear=doy_capped, hour=hour,
+            dayofyear=doy_capped,
+            hour=hour,
         ).reindex_like(t2m, method="nearest")
         hw_6h = hw_6h & (t2m < clim_hw_upper_aligned)
 
     if clim_fz_lower is not None:
         clim_fz_lower_aligned = clim_fz_lower.sel(
-            dayofyear=doy_capped, hour=hour,
+            dayofyear=doy_capped,
+            hour=hour,
         ).reindex_like(t2m, method="nearest")
         fz_6h = fz_6h & (t2m > clim_fz_lower_aligned)
 
@@ -219,9 +223,12 @@ def apply_consecutive_filter(
     strict = (
         ndimage.binary_dilation(
             ndimage.binary_erosion(
-                mask, structure=struct, border_value=False,
+                mask,
+                structure=struct,
+                border_value=False,
             ),
-            structure=struct, border_value=False,
+            structure=struct,
+            border_value=False,
         )
         & mask
     )
@@ -232,15 +239,20 @@ def apply_consecutive_filter(
     close_k = np.zeros((2 * max_grace_days + 1, 1, 1), dtype=bool)
     close_k[:, 0, 0] = True
     filled = ndimage.binary_closing(
-        mask, structure=close_k, border_value=False,
+        mask,
+        structure=close_k,
+        border_value=False,
     )
 
     filled_runs = (
         ndimage.binary_dilation(
             ndimage.binary_erosion(
-                filled, structure=struct, border_value=False,
+                filled,
+                structure=struct,
+                border_value=False,
             ),
-            structure=struct, border_value=False,
+            structure=struct,
+            border_value=False,
         )
         & filled
     )
@@ -574,18 +586,14 @@ def main():
         type=float,
         default=0.85,
         help=(
-            "Climatology quantile for heat waves "
-            f"({VALID_QUANTILES}; default: 0.85)"
+            f"Climatology quantile for heat waves ({VALID_QUANTILES}; default: 0.85)"
         ),
     )
     parser.add_argument(
         "--quantile-fz",
         type=float,
         default=0.15,
-        help=(
-            "Climatology quantile for freezes "
-            f"({VALID_QUANTILES}; default: 0.15)"
-        ),
+        help=(f"Climatology quantile for freezes ({VALID_QUANTILES}; default: 0.15)"),
     )
     parser.add_argument(
         "--quantile-hw-upper",
@@ -621,19 +629,13 @@ def main():
         "--lat-min",
         type=float,
         default=-90.0,
-        help=(
-            "Minimum latitude to include in detection. "
-            "Default is -90.0"
-        ),
+        help=("Minimum latitude to include in detection. Default is -90.0"),
     )
     parser.add_argument(
         "--lat-max",
         type=float,
         default=90.0,
-        help=(
-            "Maximum latitude to include in detection. "
-            "Default is 90.0"
-        ),
+        help=("Maximum latitude to include in detection. Default is 90.0"),
     )
     args = parser.parse_args()
 
@@ -662,13 +664,11 @@ def main():
     logger.info("  sizes=%s", dict(t2m.sizes))
 
     logger.info("Loading climatology thresholds...")
-    clim_hw, clim_fz, clim_hw_upper, clim_fz_lower = (
-        get_climatology_thresholds(
-            q_hw=args.quantile_hw,
-            q_fz=args.quantile_fz,
-            q_hw_upper=args.quantile_hw_upper,
-            q_fz_lower=args.quantile_fz_lower,
-        )
+    clim_hw, clim_fz, clim_hw_upper, clim_fz_lower = get_climatology_thresholds(
+        q_hw=args.quantile_hw,
+        q_fz=args.quantile_fz,
+        q_hw_upper=args.quantile_hw_upper,
+        q_fz_lower=args.quantile_fz_lower,
     )
 
     logger.info("Building land mask...")
diff --git a/data_prep/plot_temperature_events.py b/data_prep/plot_temperature_events.py
index ed33ac39..6eb2bb17 100644
--- a/data_prep/plot_temperature_events.py
+++ b/data_prep/plot_temperature_events.py
@@ -66,10 +66,7 @@ def resolve_op(op_str: str) -> Callable:
         ValueError: If op_str is not a recognised operator.
     """
     if op_str not in _OP_MAP:
-        raise ValueError(
-            f"Unknown operator {op_str!r}; "
-            f"choose from {list(_OP_MAP)}"
-        )
+        raise ValueError(f"Unknown operator {op_str!r}; choose from {list(_OP_MAP)}")
     return _OP_MAP[op_str]
 
 
@@ -516,7 +513,8 @@ def compute_consecutive_field(
 
     logger.info(
         "Loading q=%.2f climatology (op=%s)...",
-        quantile, op_str,
+        quantile,
+        op_str,
     )
     clim = defaults.get_climatology(quantile).sortby("latitude")
 
@@ -539,15 +537,11 @@ def compute_consecutive_field(
     logger.info("Computing 6-hourly exceedance...")
     exc_6h = cmp(t2m, clim_aligned)
 
-    daily_all_pass = (
-        exc_6h.resample({tdim: "1D"}).min().astype(bool)
-    )
+    daily_all_pass = exc_6h.resample({tdim: "1D"}).min().astype(bool)
 
     logger.info("Building land mask...")
     land = regionmask.defined_regions.natural_earth_v5_0_0.land_110
-    land_mask = land.mask(
-        daily_all_pass.longitude, daily_all_pass.latitude
-    ) == 0
+    land_mask = land.mask(daily_all_pass.longitude, daily_all_pass.latitude) == 0
 
     exc = daily_all_pass & land_mask
     mask_np = exc.compute().values.astype(bool)

From c22e9d6cedc23b6abb2ef731822abf510f03f171 Mon Sep 17 00:00:00 2001
From: aaTman <mandelbaum.taylor@gmail.com>
Date: Wed, 8 Apr 2026 01:08:44 +0000
Subject: [PATCH 05/12] add restriction criteria and method

---
 data_prep/heat_cold_bounds_global.py | 182 +++++++++++++++++++++++++++
 1 file changed, 182 insertions(+)

diff --git a/data_prep/heat_cold_bounds_global.py b/data_prep/heat_cold_bounds_global.py
index 442f66b9..59f35361 100644
--- a/data_prep/heat_cold_bounds_global.py
+++ b/data_prep/heat_cold_bounds_global.py
@@ -47,6 +47,9 @@
 MIN_CONSECUTIVE_DAYS = 3
 AREA_DECLINE_FRACTION = 0.5
 MIN_GRIDPOINTS = 500
+EXPANSION_DEGREES = 1
+MAX_SPATIAL_ITERATIONS = 20
+EDGE_VALIDITY_THRESHOLD = 0.5
 
 
 def get_climatology_thresholds(
@@ -113,6 +116,167 @@ def build_land_mask(
     return mask == 0
 
 
+def _edge_valid_fraction(
+    mask_2d: np.ndarray,
+    land_2d: np.ndarray,
+    edge: str,
+    band_pts: int,
+) -> float:
+    """Return the fraction of land grid points on an edge that are active.
+
+    Ocean/masked points are excluded from both numerator and denominator
+    so coastal edges are not penalised.
+
+    Args:
+        mask_2d: 2-D boolean activity array (lat, lon).
+        land_2d: 2-D boolean land mask (lat, lon); True = land.
+        edge: One of "north", "south", "east", "west".
+        band_pts: Width of the edge strip in grid points.
+
+    Returns:
+        Fraction in [0, 1] of land points in the strip that are active,
+        or 0.0 if the strip contains no land points.
+
+    Raises:
+        ValueError: If edge is not one of the recognised values.
+    """
+    if edge == "north":
+        strip = mask_2d[-band_pts:, :]
+        land_strip = land_2d[-band_pts:, :]
+    elif edge == "south":
+        strip = mask_2d[:band_pts, :]
+        land_strip = land_2d[:band_pts, :]
+    elif edge == "west":
+        strip = mask_2d[:, :band_pts]
+        land_strip = land_2d[:, :band_pts]
+    elif edge == "east":
+        strip = mask_2d[:, -band_pts:]
+        land_strip = land_2d[:, -band_pts:]
+    else:
+        raise ValueError(f"Unknown edge: {edge}")
+    n_land = int(land_strip.sum())
+    if n_land == 0:
+        return 0.0
+    return float((strip & land_strip).sum()) / n_land
+
+
+def expand_event_bounds(
+    event: Dict,
+    filt_mask: np.ndarray,
+    dates: np.ndarray,
+    lats: np.ndarray,
+    lons: np.ndarray,
+    land_mask_np: np.ndarray,
+) -> Dict:
+    """Expand a detected event's bounding box using case-script logic.
+
+    Starting from the blob-tracked bounding box, expands each edge
+    outward while >= EDGE_VALIDITY_THRESHOLD of land points on that
+    edge are active on the peak-footprint day. Also computes and stores
+    max_consecutive_days within the expanded region.
+
+    The expansion matches the spatial-growth algorithm used in
+    heat_cold_bounds_case.py: 1-degree steps, ocean-heavy edges
+    disabled upfront, convergence when all active edges drop below 50%.
+
+    Args:
+        event: Event dict from detect_events (keys: start, end,
+            lat_min, lat_max, lon_min, lon_max).
+        filt_mask: Boolean array (time, lat, lon) with consecutive-day
+            filtering already applied, aligned with dates/lats/lons.
+        dates: 1-D datetime64 array of daily timestamps (axis 0).
+        lats: 1-D latitude array (axis 1).
+        lons: 1-D longitude array (axis 2).
+        land_mask_np: 2-D boolean array (lat, lon); True = land.
+
+    Returns:
+        The same event dict with updated lat_min/lat_max/lon_min/
+        lon_max and a new max_consecutive_days key.
+    """
+    grid_res = float(np.abs(np.diff(lats[:2]))[0]) if len(lats) > 1 else 0.25
+    band_pts = max(1, int(round(EXPANSION_DEGREES / grid_res)))
+
+    start_date = np.datetime64(event["start"])
+    end_date = np.datetime64(event["end"])
+    t_mask = (dates >= start_date) & (dates <= end_date)
+    ev_filt = filt_mask[t_mask]
+
+    if ev_filt.shape[0] == 0:
+        event["max_consecutive_days"] = 0
+        return event
+
+    def _lat_idx(val: float) -> int:
+        return int(np.argmin(np.abs(lats - val)))
+
+    def _lon_idx(val: float) -> int:
+        return int(np.argmin(np.abs(lons - val)))
+
+    idx_s0 = _lat_idx(event["lat_min"])
+    idx_n0 = _lat_idx(event["lat_max"])
+    idx_w0 = _lon_idx(event["lon_min"])
+    idx_e0 = _lon_idx(event["lon_max"])
+
+    daily_counts = ev_filt.sum(axis=(1, 2))
+    max_count = daily_counts.max()
+    tied_days = np.where(daily_counts == max_count)[0]
+    n_tied = len(tied_days)
+    if n_tied <= 2:
+        peak_day = int(tied_days[-1]) if n_tied == 2 else int(tied_days[0])
+    else:
+        peak_day = int(tied_days[(n_tied - 1) // 2])
+    peak_mask = ev_filt[peak_day]
+
+    idx_s, idx_n, idx_w, idx_e = idx_s0, idx_n0, idx_w0, idx_e0
+
+    init_region_land = land_mask_np[idx_s0 : idx_n0 + 1, idx_w0 : idx_e0 + 1]
+    edges_active: Dict[str, bool] = {}
+    for edge in ("north", "south", "east", "west"):
+        if edge == "north":
+            strip = init_region_land[-band_pts:, :]
+        elif edge == "south":
+            strip = init_region_land[:band_pts, :]
+        elif edge == "west":
+            strip = init_region_land[:, :band_pts]
+        else:
+            strip = init_region_land[:, -band_pts:]
+        land_frac = strip.sum() / max(strip.size, 1)
+        edges_active[edge] = bool(land_frac >= 0.25)
+
+    for _ in range(MAX_SPATIAL_ITERATIONS):
+        region = peak_mask[idx_s : idx_n + 1, idx_w : idx_e + 1]
+        land_region = land_mask_np[idx_s : idx_n + 1, idx_w : idx_e + 1]
+
+        all_done = True
+        for edge in list(edges_active.keys()):
+            if not edges_active[edge]:
+                continue
+            frac = _edge_valid_fraction(region, land_region, edge, band_pts)
+            if frac < EDGE_VALIDITY_THRESHOLD:
+                edges_active[edge] = False
+            else:
+                all_done = False
+                if edge == "north":
+                    idx_n = min(len(lats) - 1, idx_n + band_pts)
+                elif edge == "south":
+                    idx_s = max(0, idx_s - band_pts)
+                elif edge == "east":
+                    idx_e = min(len(lons) - 1, idx_e + band_pts)
+                elif edge == "west":
+                    idx_w = max(0, idx_w - band_pts)
+
+        if all_done:
+            break
+
+    fin_filt = ev_filt[:, idx_s : idx_n + 1, idx_w : idx_e + 1]
+    consec = max_consecutive_days(fin_filt)
+    event["lat_min"] = float(lats[idx_s])
+    event["lat_max"] = float(lats[idx_n])
+    event["lon_min"] = float(lons[idx_w])
+    event["lon_max"] = float(lons[idx_e])
+    event["max_consecutive_days"] = int(consec.max()) if consec.size > 0 else 0
+    return event
+
+
 def build_exceedance_masks(
     t2m: xr.DataArray,
     clim_hw: xr.DataArray,
@@ -523,6 +687,7 @@ def events_to_dataframe(
         "latitude_max",
         "longitude_min",
         "longitude_max",
+        "max_consecutive_days",
     ]
     if not events:
         return pd.DataFrame(columns=columns)
@@ -548,6 +713,7 @@ def events_to_dataframe(
             "latitude_max": e["lat_max"],
             "longitude_min": e["lon_min"],
             "longitude_max": e["lon_max"],
+            "max_consecutive_days": e.get("max_consecutive_days", 0),
         }
         for e in events
     ]
@@ -673,6 +839,7 @@ def main():
 
     logger.info("Building land mask...")
     land_mask = build_land_mask(t2m.longitude, t2m.latitude)
+    land_mask_np = land_mask.values.astype(bool)
 
     logger.info("Building exceedance masks (lazy)...")
     hw_lazy, fz_lazy = build_exceedance_masks(
@@ -743,6 +910,21 @@ def main():
     )
     logger.info("  %d events", len(fz_ev))
 
+    logger.info(
+        "Expanding event bounds (case-script spatial logic, %d-deg steps)...",
+        EXPANSION_DEGREES,
+    )
+    t0 = time_module.time()
+    hw_ev = [
+        expand_event_bounds(ev, hw_filt, dates, lats, lons, land_mask_np)
+        for ev in hw_ev
+    ]
+    fz_ev = [
+        expand_event_bounds(ev, fz_filt, dates, lats, lons, land_mask_np)
+        for ev in fz_ev
+    ]
+    logger.info("  done in %.1f s", time_module.time() - t0)
+
     logger.info("Computing max-consecutive-days fields for plots...")
     stem = str(pathlib.Path(args.output).with_suffix(""))
     hw_consec = max_consecutive_days(hw_filt)

From 51e80892ab020ea753c92488655ca269d5dfbeee Mon Sep 17 00:00:00 2001
From: taylor <mandelbaum.taylor@gmail.com>
Date: Wed, 8 Apr 2026 03:24:20 +0000
Subject: [PATCH 06/12] add animation code

---
 data_prep/animate_heat_cold_global.py | 274 ++++++++++++++++++++++++++
 1 file changed, 274 insertions(+)
 create mode 100644 data_prep/animate_heat_cold_global.py

diff --git a/data_prep/animate_heat_cold_global.py b/data_prep/animate_heat_cold_global.py
new file mode 100644
index 00000000..a8bbea56
--- /dev/null
+++ b/data_prep/animate_heat_cold_global.py
@@ -0,0 +1,274 @@
+"""Animate daily heat-wave and cold-snap exceedance masks from ERA5.
+
+Produces two animated GIFs — one for heat waves, one for cold snaps —
+showing the filtered (3+ consecutive day) exceedance mask day by day on
+a global Robinson projection.
+
+Usage:
+    python animate_heat_cold_global.py \\
+        --start-date 2020-01-01 --end-date 2020-03-01 \\
+        --fps 4 --n-workers 4
+"""
+
+import argparse
+import logging
+import time as time_module
+from typing import Literal
+
+import cartopy.crs as ccrs
+import cartopy.feature as cfeature
+import matplotlib.animation as mpl_anim
+import matplotlib.colors as mcolors
+import matplotlib.pyplot as plt
+import numpy as np
+from dask.distributed import Client, LocalCluster
+
+from heat_cold_bounds_global import (
+    apply_consecutive_filter,
+    build_exceedance_masks,
+    build_land_mask,
+    get_climatology_thresholds,
+)
+from plot_temperature_events import (
+    _add_map_features,
+    detect_time_dim,
+    open_era5_t2m,
+)
+
+logging.basicConfig(
+    level=logging.INFO,
+    format="%(asctime)s %(name)s %(levelname)s %(message)s",
+)
+logger = logging.getLogger(__name__)
+
+_HW_COLOR = "#d73027"
+_FZ_COLOR = "#4575b4"
+
+
+def animate_exceedance(
+    filt_mask: np.ndarray,
+    dates: np.ndarray,
+    lats: np.ndarray,
+    lons: np.ndarray,
+    event_type: Literal["heat_wave", "cold_snap"],
+    output_path: str,
+    fps: int = 4,
+) -> None:
+    """Save an animated GIF of the daily filtered exceedance mask.
+
+    Each frame shows one day's exceedance footprint on a global
+    Robinson projection. Active grid points are colored red (heat
+    wave) or blue (cold snap); inactive land and ocean are rendered
+    in whitesmoke and light blue respectively.
+
+    Args:
+        filt_mask: Boolean array (time, lat, lon) with the
+            consecutive-day filter already applied.
+        dates: 1-D datetime64 array aligned with axis 0 of
+            filt_mask.
+        lats: 1-D latitude array (ascending).
+        lons: 1-D longitude array (0-360 degrees).
+        event_type: ``"heat_wave"`` or ``"cold_snap"``.
+        output_path: Destination ``.gif`` file path.
+        fps: Frames per second for the output GIF. Default 4.
+    """
+    is_hw = event_type == "heat_wave"
+    color = _HW_COLOR if is_hw else _FZ_COLOR
+    kind = "Heat Wave" if is_hw else "Cold Snap"
+    n_days = filt_mask.shape[0]
+
+    rgba = mcolors.to_rgba(color)
+    cmap = mcolors.ListedColormap(["none", rgba])
+
+    fig, ax = plt.subplots(
+        subplot_kw={"projection": ccrs.Robinson()},
+        figsize=(14, 7),
+    )
+    ax.set_global()
+    ax.add_feature(cfeature.OCEAN, facecolor="lightblue", zorder=0)
+    ax.add_feature(cfeature.LAND, facecolor="whitesmoke", zorder=0)
+    _add_map_features(ax)
+
+    data0 = filt_mask[0].astype(float)
+    mesh = ax.pcolormesh(
+        lons,
+        lats,
+        data0,
+        cmap=cmap,
+        vmin=0,
+        vmax=1,
+        transform=ccrs.PlateCarree(),
+        shading="auto",
+        zorder=1,
+    )
+
+    date_str = str(dates[0])[:10]
+    title = ax.set_title(
+        f"{kind} Exceedance \u2014 {date_str}",
+        loc="left",
+        fontsize=13,
+    )
+
+    def _update(di: int):
+        mesh.set_array(filt_mask[di].astype(float).ravel())
+        title.set_text(
+            f"{kind} Exceedance \u2014 {str(dates[di])[:10]}"
+        )
+        if di % 10 == 0:
+            logger.info(
+                "  Rendering frame %d / %d  (%s)",
+                di + 1,
+                n_days,
+                str(dates[di])[:10],
+            )
+        return mesh, title
+
+    anim = mpl_anim.FuncAnimation(
+        fig,
+        _update,
+        frames=n_days,
+        interval=1000 // fps,
+        blit=False,
+    )
+
+    writer = mpl_anim.PillowWriter(fps=fps)
+    logger.info("Saving %s animation to %s ...", kind, output_path)
+    anim.save(output_path, writer=writer, dpi=100)
+    plt.close(fig)
+    logger.info("  Saved %s", output_path)
+
+
+def main() -> None:
+    parser = argparse.ArgumentParser(
+        description=(
+            "Animate daily heat-wave and cold-snap exceedance masks "
+            "from ERA5 reanalysis."
+        ),
+    )
+    parser.add_argument(
+        "--start-date",
+        required=True,
+        help="Start date YYYY-MM-DD",
+    )
+    parser.add_argument(
+        "--end-date",
+        required=True,
+        help="End date YYYY-MM-DD",
+    )
+    parser.add_argument(
+        "--output-heat",
+        default=None,
+        help=(
+            "Output GIF path for heat wave animation "
+            "(default: heat_exceedance_<start>_<end>.gif)"
+        ),
+    )
+    parser.add_argument(
+        "--output-cold",
+        default=None,
+        help=(
+            "Output GIF path for cold snap animation "
+            "(default: cold_exceedance_<start>_<end>.gif)"
+        ),
+    )
+    parser.add_argument(
+        "--fps",
+        type=int,
+        default=4,
+        help="Frames per second for the output GIFs (default: 4)",
+    )
+    parser.add_argument(
+        "--n-workers",
+        type=int,
+        default=4,
+        help="Number of dask workers (default: 4)",
+    )
+    args = parser.parse_args()
+
+    if args.output_heat is None:
+        args.output_heat = (
+            f"heat_exceedance_{args.start_date}_{args.end_date}.gif"
+        )
+    if args.output_cold is None:
+        args.output_cold = (
+            f"cold_exceedance_{args.start_date}_{args.end_date}.gif"
+        )
+
+    wall_start = time_module.time()
+    client = Client(LocalCluster(n_workers=args.n_workers))
+    logger.info("Dask dashboard: %s", client.dashboard_link)
+
+    logger.info("Opening ERA5 data...")
+    t2m = open_era5_t2m(args.start_date, args.end_date)
+    logger.info("  sizes=%s", dict(t2m.sizes))
+
+    logger.info("Loading climatology thresholds...")
+    clim_hw, clim_fz, _, _ = get_climatology_thresholds()
+
+    logger.info("Building land mask...")
+    land_mask = build_land_mask(t2m.longitude, t2m.latitude)
+
+    logger.info("Building exceedance masks (lazy)...")
+    hw_lazy, fz_lazy = build_exceedance_masks(
+        t2m, clim_hw, clim_fz, land_mask
+    )
+
+    tdim = detect_time_dim(hw_lazy)
+
+    logger.info("Computing heat wave mask...")
+    t0 = time_module.time()
+    hw_da = hw_lazy.compute()
+    logger.info("  done in %.1f s", time_module.time() - t0)
+
+    logger.info("Computing cold snap mask...")
+    t0 = time_module.time()
+    fz_da = fz_lazy.compute()
+    logger.info("  done in %.1f s", time_module.time() - t0)
+
+    dates = hw_da[tdim].values
+    lats = hw_da.latitude.values
+    lons = hw_da.longitude.values
+    hw_np = hw_da.values.astype(bool)
+    fz_np = fz_da.values.astype(bool)
+    del hw_da, fz_da
+
+    logger.info("Applying 3-day consecutive filter...")
+    hw_filt = apply_consecutive_filter(hw_np)
+    fz_filt = apply_consecutive_filter(fz_np)
+    logger.info(
+        "  HW active cells: %d -> %d",
+        hw_np.sum(),
+        hw_filt.sum(),
+    )
+    logger.info(
+        "  FZ active cells: %d -> %d",
+        fz_np.sum(),
+        fz_filt.sum(),
+    )
+    del hw_np, fz_np
+
+    animate_exceedance(
+        hw_filt,
+        dates,
+        lats,
+        lons,
+        "heat_wave",
+        args.output_heat,
+        fps=args.fps,
+    )
+    animate_exceedance(
+        fz_filt,
+        dates,
+        lats,
+        lons,
+        "cold_snap",
+        args.output_cold,
+        fps=args.fps,
+    )
+
+    client.close()
+    logger.info("Done in %.1f s", time_module.time() - wall_start)
+
+
+if __name__ == "__main__":
+    main()

From 14caa7124bdb3b9b58c1ac956973d685adc03895 Mon Sep 17 00:00:00 2001
From: aaTman <mandelbaum.taylor@gmail.com>
Date: Wed, 8 Apr 2026 03:40:55 +0000
Subject: [PATCH 07/12] simplify to include climatology bound(s), fix func
 names

---
 data_prep/animate_heat_cold_global.py |  24 +-
 data_prep/heat_cold_bounds_global.py  | 465 +++++++++++++-------------
 2 files changed, 243 insertions(+), 246 deletions(-)

diff --git a/data_prep/animate_heat_cold_global.py b/data_prep/animate_heat_cold_global.py
index a8bbea56..0090bc28 100644
--- a/data_prep/animate_heat_cold_global.py
+++ b/data_prep/animate_heat_cold_global.py
@@ -25,9 +25,9 @@
 
 from heat_cold_bounds_global import (
     apply_consecutive_filter,
-    build_exceedance_masks,
+    build_exceedance_mask,
     build_land_mask,
-    get_climatology_thresholds,
+    get_climatology_bounds,
 )
 from plot_temperature_events import (
     _add_map_features,
@@ -111,9 +111,7 @@ def animate_exceedance(
 
     def _update(di: int):
         mesh.set_array(filt_mask[di].astype(float).ravel())
-        title.set_text(
-            f"{kind} Exceedance \u2014 {str(dates[di])[:10]}"
-        )
+        title.set_text(f"{kind} Exceedance \u2014 {str(dates[di])[:10]}")
         if di % 10 == 0:
             logger.info(
                 "  Rendering frame %d / %d  (%s)",
@@ -186,13 +184,9 @@ def main() -> None:
     args = parser.parse_args()
 
     if args.output_heat is None:
-        args.output_heat = (
-            f"heat_exceedance_{args.start_date}_{args.end_date}.gif"
-        )
+        args.output_heat = f"heat_exceedance_{args.start_date}_{args.end_date}.gif"
     if args.output_cold is None:
-        args.output_cold = (
-            f"cold_exceedance_{args.start_date}_{args.end_date}.gif"
-        )
+        args.output_cold = f"cold_exceedance_{args.start_date}_{args.end_date}.gif"
 
     wall_start = time_module.time()
     client = Client(LocalCluster(n_workers=args.n_workers))
@@ -203,15 +197,15 @@ def main() -> None:
     logger.info("  sizes=%s", dict(t2m.sizes))
 
     logger.info("Loading climatology thresholds...")
-    clim_hw, clim_fz, _, _ = get_climatology_thresholds()
+    clim_hw, _ = get_climatology_bounds(q_lower=0.85)
+    _, clim_fz = get_climatology_bounds(q_upper=0.15)
 
     logger.info("Building land mask...")
     land_mask = build_land_mask(t2m.longitude, t2m.latitude)
 
     logger.info("Building exceedance masks (lazy)...")
-    hw_lazy, fz_lazy = build_exceedance_masks(
-        t2m, clim_hw, clim_fz, land_mask
-    )
+    hw_lazy = build_exceedance_mask(t2m, clim_hw, None, land_mask)
+    fz_lazy = build_exceedance_mask(t2m, None, clim_fz, land_mask)
 
     tdim = detect_time_dim(hw_lazy)
 
diff --git a/data_prep/heat_cold_bounds_global.py b/data_prep/heat_cold_bounds_global.py
index 59f35361..a4d5cb99 100644
--- a/data_prep/heat_cold_bounds_global.py
+++ b/data_prep/heat_cold_bounds_global.py
@@ -1,17 +1,30 @@
-"""Detect heat waves and cold snaps globally from ERA5 reanalysis.
+"""Detect temperature exceedance events globally from ERA5 reanalysis.
 
 Scans ERA5 2m temperature over an input date range and identifies
-heat wave (daily max > 85th percentile for 3+ consecutive days)
-and cold snap (daily min < 15th percentile for 3+ consecutive days)
-events globally over land.
+events where the daily temperature falls within a user-specified
+climatology quantile band for 3+ consecutive days, over land.
 
-Bounding boxes represent the maximum spatial extent of each event.
-Events terminate when area drops below 50% of their peak area.
+At least one of --quantile-lower or --quantile-upper must be given.
+Both can be combined to define a band (e.g. 50th–85th percentile).
+
+Bounding boxes are first derived from blob tracking, then expanded
+using the same edge-validity logic as heat_cold_bounds_case.py:
+each edge grows by 1 degree while >= 50% of its land points are
+active on the peak-footprint day.  Events terminate when their
+active area drops below 50% of peak.
 
 Usage:
+    # Anything above the 85th percentile (heat wave)
     python heat_cold_bounds_global.py \\
         --start-date 2023-06-01 --end-date 2023-09-01 \\
-        --output heat_cold_global.csv --n-workers 4
+        --quantile-lower 0.85 --operator-lower ">=" \\
+        --event-type heat_wave --output heat_cold_global.csv
+
+    # Band between 50th and 85th (moderate heat)
+    python heat_cold_bounds_global.py \\
+        --start-date 2023-06-01 --end-date 2023-09-01 \\
+        --quantile-lower 0.50 --quantile-upper 0.85 \\
+        --event-type heat_wave --output heat_cold_global.csv
 """
 
 import argparse
@@ -20,6 +33,8 @@
 import time as time_module
 from typing import Dict, List, Optional, Tuple
 
+import joblib
+
 import numba as nb
 import numpy as np
 import pandas as pd
@@ -47,54 +62,73 @@
 MIN_CONSECUTIVE_DAYS = 3
 AREA_DECLINE_FRACTION = 0.5
 MIN_GRIDPOINTS = 500
+MIN_AREA_KM2 = 200_000.0
 EXPANSION_DEGREES = 1
 MAX_SPATIAL_ITERATIONS = 20
 EDGE_VALIDITY_THRESHOLD = 0.5
 
 
-def get_climatology_thresholds(
-    q_hw: float = 0.85,
-    q_fz: float = 0.15,
-    q_hw_upper: Optional[float] = None,
-    q_fz_lower: Optional[float] = None,
-) -> Tuple[
-    xr.DataArray,
-    xr.DataArray,
-    Optional[xr.DataArray],
-    Optional[xr.DataArray],
-]:
-    """Return percentile climatology DataArrays for heat/freeze detection.
+def compute_grid_cell_area(
+    lats: np.ndarray,
+    lons: np.ndarray,
+) -> np.ndarray:
+    """Return a 2-D (lat, lon) array of grid-cell areas in km².
+
+    Uses the spherical-Earth approximation:
+        area = R² × Δlat_rad × Δlon_rad × cos(lat)
 
     Args:
-        q_hw: Lower-bound quantile for heat wave detection.
-            Default is 0.85.
-        q_fz: Upper-bound quantile for freeze detection.
-            Default is 0.15.
-        q_hw_upper: Upper-bound quantile for heat waves. When set,
-            only days where temp > q_hw AND temp < q_hw_upper are
-            flagged. Default is None.
-        q_fz_lower: Lower-bound quantile for freezes. When set, only
-            days where temp < q_fz AND temp > q_fz_lower are
-            flagged. Default is None.
+        lats: 1-D latitude array in degrees.
+        lons: 1-D longitude array in degrees; used only for shape.
 
     Returns:
-        A tuple of (clim_hw, clim_fz, clim_hw_upper, clim_fz_lower).
-        The last two elements are None when the corresponding optional
-        quantile argument is not supplied.
+        Float64 array of shape (len(lats), len(lons)) with each
+        cell's surface area in km².
+    """
+    R_KM = 6371.0
+    dlat = float(np.abs(np.diff(lats[:2]))[0]) if len(lats) > 1 else 0.25
+    dlon = float(np.abs(np.diff(lons[:2]))[0]) if len(lons) > 1 else 0.25
+    cell_km2 = R_KM**2 * np.deg2rad(dlat) * np.deg2rad(dlon) * np.cos(np.deg2rad(lats))
+    return np.outer(cell_km2, np.ones(len(lons)))
+
+
+def get_climatology_bounds(
+    q_lower: Optional[float] = None,
+    q_upper: Optional[float] = None,
+) -> Tuple[Optional[xr.DataArray], Optional[xr.DataArray]]:
+    """Return climatology DataArrays for the lower and/or upper bound.
+
+    At least one of q_lower or q_upper must be provided. Each returned
+    DataArray is indexed by (dayofyear, hour) and sorted by latitude.
+
+    Args:
+        q_lower: Quantile for the lower bound (e.g. 0.50 means temp
+            must exceed the 50th-percentile climatology). None skips
+            the lower-bound check.
+        q_upper: Quantile for the upper bound (e.g. 0.85 means temp
+            must not exceed the 85th-percentile climatology). None
+            skips the upper-bound check.
+
+    Returns:
+        Tuple (clim_lower, clim_upper); either element may be None
+        when the corresponding quantile argument is not supplied.
+
+    Raises:
+        ValueError: If both q_lower and q_upper are None.
     """
-    clim_hw = defaults.get_climatology(q_hw).sortby("latitude")
-    clim_fz = defaults.get_climatology(q_fz).sortby("latitude")
-    clim_hw_upper = (
-        defaults.get_climatology(q_hw_upper).sortby("latitude")
-        if q_hw_upper is not None
+    if q_lower is None and q_upper is None:
+        raise ValueError("At least one of q_lower or q_upper must be set.")
+    clim_lower = (
+        defaults.get_climatology(q_lower).sortby("latitude")
+        if q_lower is not None
         else None
     )
-    clim_fz_lower = (
-        defaults.get_climatology(q_fz_lower).sortby("latitude")
-        if q_fz_lower is not None
+    clim_upper = (
+        defaults.get_climatology(q_upper).sortby("latitude")
+        if q_upper is not None
         else None
     )
-    return clim_hw, clim_fz, clim_hw_upper, clim_fz_lower
+    return clim_lower, clim_upper
 
 
 def build_land_mask(
@@ -277,85 +311,63 @@ def _lon_idx(val: float) -> int:
     return event
 
 
-def build_exceedance_masks(
+def build_exceedance_mask(
     t2m: xr.DataArray,
-    clim_hw: xr.DataArray,
-    clim_fz: xr.DataArray,
+    clim_lower: Optional[xr.DataArray],
+    clim_upper: Optional[xr.DataArray],
     land_mask: xr.DataArray,
-    op_hw: str = ">",
-    op_fz: str = "<",
-    clim_hw_upper: Optional[xr.DataArray] = None,
-    clim_fz_lower: Optional[xr.DataArray] = None,
-) -> Tuple[xr.DataArray, xr.DataArray]:
-    """Build daily exceedance masks from 6-hourly temperature data.
-
-    Each 6-hourly timestep is compared to its matching
-    (dayofyear, hour) climatology. A day passes only if all four
-    6-hourly timesteps exceed the threshold.
-
-    When clim_hw_upper is provided, heat-wave days must also satisfy
-    temp < upper bound on every 6-hourly step. When clim_fz_lower is
-    provided, freeze days must also satisfy temp > lower bound on
-    every 6-hourly step.
+    op_lower: str = ">",
+    op_upper: str = "<",
+) -> xr.DataArray:
+    """Build a daily exceedance mask from 6-hourly temperature data.
+
+    Each 6-hourly timestep must satisfy all provided bounds. A day
+    passes only when every 6-hourly step satisfies every active bound.
+    The result is further masked to land points.
+
+    At least one of clim_lower or clim_upper must be provided.
 
     Args:
         t2m: 6-hourly 2m temperature DataArray.
-        clim_hw: Heat-wave lower-bound climatology indexed by
-            (dayofyear, hour).
-        clim_fz: Freeze upper-bound climatology indexed by
-            (dayofyear, hour).
+        clim_lower: Climatology for the lower bound, indexed by
+            (dayofyear, hour). None skips this bound.
+        clim_upper: Climatology for the upper bound, indexed by
+            (dayofyear, hour). None skips this bound.
         land_mask: Boolean DataArray (True = land) matching t2m grid.
-        op_hw: Comparison operator string for heat waves.
+        op_lower: Comparison operator for the lower bound.
             Default is ">".
-        op_fz: Comparison operator string for freezes.
+        op_upper: Comparison operator for the upper bound.
             Default is "<".
-        clim_hw_upper: Optional upper-bound climatology for heat
-            waves (exclusive cap). Default is None.
-        clim_fz_lower: Optional lower-bound climatology for freezes
-            (exclusive floor). Default is None.
 
     Returns:
-        A tuple (hw, fz) of daily boolean DataArrays masked to land,
-        where True indicates an exceedance day.
+        Daily boolean DataArray masked to land where True indicates
+        every 6-hourly step satisfied all active bounds.
     """
-    cmp_hw = resolve_op(op_hw)
-    cmp_fz = resolve_op(op_fz)
     tdim = detect_time_dim(t2m)
-
     doy = t2m[tdim].dt.dayofyear
     hour = t2m[tdim].dt.hour
-    max_clim_doy = int(clim_hw.dayofyear.max())
+    ref = clim_lower if clim_lower is not None else clim_upper
+    assert ref is not None
+    max_clim_doy = int(ref.dayofyear.max())
     doy_capped = doy.clip(max=max_clim_doy)
 
-    clim_hw_aligned = clim_hw.sel(
-        dayofyear=doy_capped,
-        hour=hour,
-    ).reindex_like(t2m, method="nearest")
-    clim_fz_aligned = clim_fz.sel(
-        dayofyear=doy_capped,
-        hour=hour,
-    ).reindex_like(t2m, method="nearest")
-
-    hw_6h = cmp_hw(t2m, clim_hw_aligned)
-    fz_6h = cmp_fz(t2m, clim_fz_aligned)
-
-    if clim_hw_upper is not None:
-        clim_hw_upper_aligned = clim_hw_upper.sel(
-            dayofyear=doy_capped,
-            hour=hour,
-        ).reindex_like(t2m, method="nearest")
-        hw_6h = hw_6h & (t2m < clim_hw_upper_aligned)
-
-    if clim_fz_lower is not None:
-        clim_fz_lower_aligned = clim_fz_lower.sel(
-            dayofyear=doy_capped,
-            hour=hour,
-        ).reindex_like(t2m, method="nearest")
-        fz_6h = fz_6h & (t2m > clim_fz_lower_aligned)
-
-    hw = hw_6h.resample({tdim: "1D"}).min().astype(bool) & land_mask
-    fz = fz_6h.resample({tdim: "1D"}).min().astype(bool) & land_mask
-    return hw, fz
+    mask_6h = xr.ones_like(t2m, dtype=bool)
+
+    if clim_lower is not None:
+        cmp = resolve_op(op_lower)
+        aligned = clim_lower.sel(dayofyear=doy_capped, hour=hour).reindex_like(
+            t2m, method="nearest"
+        )
+        mask_6h = mask_6h & cmp(t2m, aligned)
+
+    if clim_upper is not None:
+        cmp = resolve_op(op_upper)
+        aligned = clim_upper.sel(dayofyear=doy_capped, hour=hour).reindex_like(
+            t2m, method="nearest"
+        )
+        mask_6h = mask_6h & cmp(t2m, aligned)
+
+    return mask_6h.resample({tdim: "1D"}).min().astype(bool) & land_mask
 
 
 def apply_consecutive_filter(
@@ -507,6 +519,7 @@ def _resolve_event(
         tgt["lon_min"] = min(tgt["lon_min"], merged["lon_min"])
         tgt["lon_max"] = max(tgt["lon_max"], merged["lon_max"])
         tgt["peak"] = max(tgt["peak"], merged["peak"])
+        tgt["peak_area_km2"] = max(tgt["peak_area_km2"], merged["peak_area_km2"])
         tgt["start"] = min(tgt["start"], merged["start"])
         for k, v in list(cur_map.items()):
             if v == other:
@@ -541,6 +554,7 @@ def detect_events(
     lats: np.ndarray,
     lons: np.ndarray,
     event_type: str,
+    area_grid: Optional[np.ndarray] = None,
 ) -> List[Dict]:
     """Track spatiotemporal events from a filtered boolean mask.
 
@@ -560,10 +574,12 @@ def detect_events(
         lons: 1-D longitude array aligned with axis 2.
         event_type: Label string stored in each returned event dict
             (e.g. "heat_wave" or "cold_snap").
+        area_grid: Optional 2-D array (lat, lon) of grid-cell areas
+            in km². When provided, peak_area_km2 is tracked per event.
 
     Returns:
         List of event dicts with keys type, start, end, lat_min,
-        lat_max, lon_min, lon_max, peak, area, done.
+        lat_max, lon_min, lon_max, peak, peak_area_km2, area, done.
     """
     n_days = filtered_mask.shape[0]
     events: Dict[int, Dict] = {}
@@ -601,6 +617,9 @@ def detect_events(
             om = labels == oid
             area = int(om.sum())
             li, lo = np.where(om)
+            blob_area_km2 = (
+                float(area_grid[li, lo].sum()) if area_grid is not None else 0.0
+            )
 
             eid = _resolve_event(
                 oid,
@@ -622,12 +641,14 @@ def detect_events(
                     "lon_min": float(lons[lo].min()),
                     "lon_max": float(lons[lo].max()),
                     "peak": area,
+                    "peak_area_km2": blob_area_km2,
                     "area": area,
                     "done": False,
                 }
             else:
                 ev = events[eid]
                 ev["end"] = dates[di]
+                ev["peak_area_km2"] = max(ev["peak_area_km2"], blob_area_km2)
                 ev["lat_min"] = min(
                     ev["lat_min"],
                     float(lats[li].min()),
@@ -663,20 +684,22 @@ def detect_events(
 def events_to_dataframe(
     events: List[Dict],
     min_gridpoints: int = MIN_GRIDPOINTS,
+    min_area_km2: float = MIN_AREA_KM2,
 ) -> pd.DataFrame:
     """Convert event dicts to a labelled DataFrame.
 
     Args:
         events: Raw event dicts from ``detect_events``.
-        min_gridpoints: Drop events whose peak spatial extent (in
-            grid points) is below this threshold. Default is 500
-            (MIN_GRIDPOINTS).
+        min_gridpoints: Drop events whose peak grid-point count is
+            below this threshold. Default is 500 (MIN_GRIDPOINTS).
+        min_area_km2: Drop events whose peak area (km²) is below
+            this threshold. Default is 200 000 (MIN_AREA_KM2).
 
     Returns:
         DataFrame with columns label, event_type, start_date,
         end_date, latitude_min, latitude_max, longitude_min,
-        longitude_max, sorted by start_date. Events below
-        min_gridpoints are excluded.
+        longitude_max, max_consecutive_days, sorted by start_date.
+        Events below either threshold are excluded.
     """
     columns = [
         "label",
@@ -693,11 +716,16 @@ def events_to_dataframe(
         return pd.DataFrame(columns=columns)
 
     n_before = len(events)
-    events = [e for e in events if e["peak"] >= min_gridpoints]
+    events = [
+        e
+        for e in events
+        if e["peak"] >= min_gridpoints and e.get("peak_area_km2", 0.0) >= min_area_km2
+    ]
     logger.info(
-        "  Filtered %d events below %d gridpoints; %d remain",
+        "  Filtered %d events (< %d pts or < %.0f km²); %d remain",
         n_before - len(events),
         min_gridpoints,
+        min_area_km2,
         len(events),
     )
 
@@ -724,7 +752,7 @@ def events_to_dataframe(
 
 def main():
     parser = argparse.ArgumentParser(
-        description=("Detect heat waves and cold snaps globally from ERA5."),
+        description=("Detect temperature exceedance events globally from ERA5."),
     )
     parser.add_argument(
         "--start-date",
@@ -738,7 +766,7 @@ def main():
     )
     parser.add_argument(
         "--output",
-        default="heat_cold_global.csv",
+        default="events_global.csv",
         help="Output CSV path",
     )
     parser.add_argument(
@@ -748,63 +776,74 @@ def main():
         help="Number of dask workers",
     )
     parser.add_argument(
-        "--quantile-hw",
+        "--quantile-lower",
         type=float,
-        default=0.85,
+        default=None,
         help=(
-            f"Climatology quantile for heat waves ({VALID_QUANTILES}; default: 0.85)"
+            f"Lower-bound climatology quantile ({VALID_QUANTILES}). "
+            "Days where temp op_lower clim_lower are candidates. "
+            "At least one of --quantile-lower or --quantile-upper "
+            "is required."
         ),
     )
     parser.add_argument(
-        "--quantile-fz",
-        type=float,
-        default=0.15,
-        help=(f"Climatology quantile for freezes ({VALID_QUANTILES}; default: 0.15)"),
-    )
-    parser.add_argument(
-        "--quantile-hw-upper",
-        type=float,
-        default=None,
-        help=(
-            "Upper-bound quantile for heat waves; days must be "
-            "> --quantile-hw AND < this value "
-            f"({VALID_QUANTILES}; default: None)"
-        ),
+        "--operator-lower",
+        default=">",
+        help=("Comparison operator applied to the lower bound (default: >)"),
     )
     parser.add_argument(
-        "--quantile-fz-lower",
+        "--quantile-upper",
         type=float,
         default=None,
         help=(
-            "Lower-bound quantile for freezes; days must be "
-            "< --quantile-fz AND > this value "
-            f"({VALID_QUANTILES}; default: None)"
+            f"Upper-bound climatology quantile ({VALID_QUANTILES}). "
+            "Days where temp op_upper clim_upper are candidates. "
+            "Combine with --quantile-lower to define a band."
         ),
     )
     parser.add_argument(
-        "--operator-hw",
-        default=">",
-        help="Comparison operator for heat waves (default: >)",
-    )
-    parser.add_argument(
-        "--operator-fz",
+        "--operator-upper",
         default="<",
-        help="Comparison operator for freezes (default: <)",
+        help=("Comparison operator applied to the upper bound (default: <)"),
     )
     parser.add_argument(
         "--lat-min",
         type=float,
         default=-90.0,
-        help=("Minimum latitude to include in detection. Default is -90.0"),
+        help="Minimum latitude to include in detection. Default -90.0",
     )
     parser.add_argument(
         "--lat-max",
         type=float,
         default=90.0,
-        help=("Maximum latitude to include in detection. Default is 90.0"),
+        help="Maximum latitude to include in detection. Default 90.0",
     )
     args = parser.parse_args()
 
+    if args.quantile_lower is None and args.quantile_upper is None:
+        parser.error(
+            "At least one of --quantile-lower or --quantile-upper must be provided."
+        )
+
+    # Derive event_type label and plot colour from the quantiles.
+    # Format: "q{lower}+" / "q{upper}-" / "q{lower}-q{upper}"
+    ql, qu = args.quantile_lower, args.quantile_upper
+    if ql is not None and qu is not None:
+        event_type = f"q{ql:.2f}-q{qu:.2f}"
+    elif ql is not None:
+        event_type = f"q{ql:.2f}+"
+    else:
+        event_type = f"q{qu:.2f}-"
+    # Infer warm vs cold for the plot colourmap:
+    # warm when the lower bound is above the median (high percentile),
+    # cold when the upper bound is below the median (low percentile).
+    if ql is not None and ql >= 0.5:
+        plot_event_type = "heat_wave"
+    elif qu is not None and qu <= 0.5:
+        plot_event_type = "cold_snap"
+    else:
+        plot_event_type = "heat_wave"
+
     wall_start = time_module.time()
     client = Client(
         LocalCluster(n_workers=args.n_workers),
@@ -829,127 +868,91 @@ def main():
         )
     logger.info("  sizes=%s", dict(t2m.sizes))
 
-    logger.info("Loading climatology thresholds...")
-    clim_hw, clim_fz, clim_hw_upper, clim_fz_lower = get_climatology_thresholds(
-        q_hw=args.quantile_hw,
-        q_fz=args.quantile_fz,
-        q_hw_upper=args.quantile_hw_upper,
-        q_fz_lower=args.quantile_fz_lower,
+    logger.info(
+        "Loading climatology bounds (lower=%s, upper=%s)...",
+        args.quantile_lower,
+        args.quantile_upper,
+    )
+    clim_lower, clim_upper = get_climatology_bounds(
+        q_lower=args.quantile_lower,
+        q_upper=args.quantile_upper,
     )
 
-    logger.info("Building land mask...")
+    logger.info("Building land mask and grid-cell area array...")
     land_mask = build_land_mask(t2m.longitude, t2m.latitude)
     land_mask_np = land_mask.values.astype(bool)
+    area_grid = compute_grid_cell_area(t2m.latitude.values, t2m.longitude.values)
 
-    logger.info("Building exceedance masks (lazy)...")
-    hw_lazy, fz_lazy = build_exceedance_masks(
+    logger.info("Building exceedance mask (lazy)...")
+    exc_lazy = build_exceedance_mask(
         t2m,
-        clim_hw,
-        clim_fz,
+        clim_lower,
+        clim_upper,
         land_mask,
-        op_hw=args.operator_hw,
-        op_fz=args.operator_fz,
-        clim_hw_upper=clim_hw_upper,
-        clim_fz_lower=clim_fz_lower,
+        op_lower=args.operator_lower,
+        op_upper=args.operator_upper,
     )
 
-    tdim = detect_time_dim(hw_lazy)
+    tdim = detect_time_dim(exc_lazy)
 
-    logger.info("Computing heat wave mask...")
+    logger.info("Computing exceedance mask...")
     t0 = time_module.time()
-    hw_da = hw_lazy.compute()
+    exc_da = exc_lazy.compute()
     logger.info("  done in %.1f s", time_module.time() - t0)
 
-    logger.info("Computing cold snap mask...")
-    t0 = time_module.time()
-    fz_da = fz_lazy.compute()
-    logger.info("  done in %.1f s", time_module.time() - t0)
-
-    dates = hw_da[tdim].values
-    lats = hw_da.latitude.values
-    lons = hw_da.longitude.values
-    hw_np = hw_da.values.astype(bool)
-    fz_np = fz_da.values.astype(bool)
-    del hw_da, fz_da
+    dates = exc_da[tdim].values
+    lats = exc_da.latitude.values
+    lons = exc_da.longitude.values
+    exc_np = exc_da.values.astype(bool)
+    del exc_da
 
     logger.info(
         "Applying %d-day consecutive filter...",
         MIN_CONSECUTIVE_DAYS,
     )
-    hw_filt = apply_consecutive_filter(hw_np)
-    fz_filt = apply_consecutive_filter(fz_np)
-    logger.info(
-        "  HW: %d -> %d True cells",
-        hw_np.sum(),
-        hw_filt.sum(),
-    )
+    exc_filt = apply_consecutive_filter(exc_np)
     logger.info(
-        "  FZ: %d -> %d True cells",
-        fz_np.sum(),
-        fz_filt.sum(),
+        "  %d -> %d True cells",
+        exc_np.sum(),
+        exc_filt.sum(),
     )
-    del hw_np, fz_np
+    del exc_np
 
-    logger.info("Detecting heat wave events...")
-    hw_ev = detect_events(
-        hw_filt,
-        dates,
-        lats,
-        lons,
-        "heat_wave",
-    )
-    logger.info("  %d events", len(hw_ev))
-
-    logger.info("Detecting cold snap events...")
-    fz_ev = detect_events(
-        fz_filt,
-        dates,
-        lats,
-        lons,
-        "cold_snap",
-    )
-    logger.info("  %d events", len(fz_ev))
+    logger.info("Detecting events...")
+    t0 = time_module.time()
+    events = detect_events(exc_filt, dates, lats, lons, event_type, area_grid=area_grid)
+    logger.info("  %d events (%.1f s)", len(events), time_module.time() - t0)
 
     logger.info(
-        "Expanding event bounds (case-script spatial logic, %d-deg steps)...",
+        "Expanding event bounds (%d-deg steps)...",
         EXPANSION_DEGREES,
     )
     t0 = time_module.time()
-    hw_ev = [
-        expand_event_bounds(ev, hw_filt, dates, lats, lons, land_mask_np)
-        for ev in hw_ev
-    ]
-    fz_ev = [
-        expand_event_bounds(ev, fz_filt, dates, lats, lons, land_mask_np)
-        for ev in fz_ev
-    ]
+    events = joblib.Parallel(n_jobs=-1, prefer="threads")(
+        joblib.delayed(expand_event_bounds)(
+            ev, exc_filt, dates, lats, lons, land_mask_np
+        )
+        for ev in events
+    )
     logger.info("  done in %.1f s", time_module.time() - t0)
 
-    logger.info("Computing max-consecutive-days fields for plots...")
+    logger.info("Computing max-consecutive-days plot...")
     stem = str(pathlib.Path(args.output).with_suffix(""))
-    hw_consec = max_consecutive_days(hw_filt)
+    consec = max_consecutive_days(exc_filt)
     plot_consecutive_map(
-        hw_consec,
+        consec,
         lats,
         lons,
-        "heat_wave",
-        title=(f"Consecutive Heatwave Days\n{args.start_date} to {args.end_date}"),
-        output_path=f"{stem}_heatwave.png",
-    )
-    del hw_consec, hw_filt
-
-    fz_consec = max_consecutive_days(fz_filt)
-    plot_consecutive_map(
-        fz_consec,
-        lats,
-        lons,
-        "cold_snap",
-        title=(f"Consecutive Cold Snap Days\n{args.start_date} to {args.end_date}"),
-        output_path=f"{stem}_cold_snap.png",
+        plot_event_type,
+        title=(
+            f"Consecutive Exceedance Days ({event_type})"
+            f"\n{args.start_date} to {args.end_date}"
+        ),
+        output_path=f"{stem}_consec.png",
     )
-    del fz_consec, fz_filt
+    del consec, exc_filt
 
-    df = events_to_dataframe(hw_ev + fz_ev)
+    df = events_to_dataframe(events)
     df.to_csv(args.output, index=False)
 
     elapsed = time_module.time() - wall_start

From 34f1b16a3ead4659abc2941f49c94914ecd0af25 Mon Sep 17 00:00:00 2001
From: aaTman <mandelbaum.taylor@gmail.com>
Date: Wed, 8 Apr 2026 20:31:34 +0000
Subject: [PATCH 08/12] update animate

---
 data_prep/animate_heat_cold_global.py | 135 +++++++++++---------------
 1 file changed, 55 insertions(+), 80 deletions(-)

diff --git a/data_prep/animate_heat_cold_global.py b/data_prep/animate_heat_cold_global.py
index 0090bc28..3c53c8aa 100644
--- a/data_prep/animate_heat_cold_global.py
+++ b/data_prep/animate_heat_cold_global.py
@@ -7,13 +7,13 @@
 Usage:
     python animate_heat_cold_global.py \\
         --start-date 2020-01-01 --end-date 2020-03-01 \\
-        --fps 4 --n-workers 4
+        --fps 4 --n-workers 4 --quantile-lower 0.85 --quantile-upper 0.15
+        --output events_global.gif
 """
 
 import argparse
 import logging
 import time as time_module
-from typing import Literal
 
 import cartopy.crs as ccrs
 import cartopy.feature as cfeature
@@ -22,7 +22,6 @@
 import matplotlib.pyplot as plt
 import numpy as np
 from dask.distributed import Client, LocalCluster
-
 from heat_cold_bounds_global import (
     apply_consecutive_filter,
     build_exceedance_mask,
@@ -41,8 +40,7 @@
 )
 logger = logging.getLogger(__name__)
 
-_HW_COLOR = "#d73027"
-_FZ_COLOR = "#4575b4"
+_EVENT_COLOR = "#d73027"
 
 
 def animate_exceedance(
@@ -50,15 +48,15 @@ def animate_exceedance(
     dates: np.ndarray,
     lats: np.ndarray,
     lons: np.ndarray,
-    event_type: Literal["heat_wave", "cold_snap"],
+    lower_quantile: float | None,
+    upper_quantile: float | None,
     output_path: str,
     fps: int = 4,
 ) -> None:
     """Save an animated GIF of the daily filtered exceedance mask.
 
     Each frame shows one day's exceedance footprint on a global
-    Robinson projection. Active grid points are colored red (heat
-    wave) or blue (cold snap); inactive land and ocean are rendered
+    Robinson projection. Active grid points are colored red; inactive land and ocean are rendered
     in whitesmoke and light blue respectively.
 
     Args:
@@ -68,22 +66,21 @@ def animate_exceedance(
             filt_mask.
         lats: 1-D latitude array (ascending).
         lons: 1-D longitude array (0-360 degrees).
-        event_type: ``"heat_wave"`` or ``"cold_snap"``.
+        lower_quantile: Lower quantile for the event.
+        upper_quantile: Upper quantile for the event.
         output_path: Destination ``.gif`` file path.
         fps: Frames per second for the output GIF. Default 4.
     """
-    is_hw = event_type == "heat_wave"
-    color = _HW_COLOR if is_hw else _FZ_COLOR
-    kind = "Heat Wave" if is_hw else "Cold Snap"
     n_days = filt_mask.shape[0]
 
-    rgba = mcolors.to_rgba(color)
+    rgba = mcolors.to_rgba(_EVENT_COLOR)
     cmap = mcolors.ListedColormap(["none", rgba])
 
     fig, ax = plt.subplots(
         subplot_kw={"projection": ccrs.Robinson()},
-        figsize=(14, 7),
+        figsize=(12, 5.5),
     )
+    fig.subplots_adjust(left=0.01, right=0.99, top=0.93, bottom=0.01)
     ax.set_global()
     ax.add_feature(cfeature.OCEAN, facecolor="lightblue", zorder=0)
     ax.add_feature(cfeature.LAND, facecolor="whitesmoke", zorder=0)
@@ -101,17 +98,22 @@ def animate_exceedance(
         shading="auto",
         zorder=1,
     )
-
+    if lower_quantile is not None and upper_quantile is not None:
+        kind = f"p{lower_quantile:.2f}-p{upper_quantile:.2f}"
+    elif lower_quantile is not None:
+        kind = f"p{lower_quantile:.2f}+"
+    else:
+        kind = f"p{upper_quantile:.2f}-"
     date_str = str(dates[0])[:10]
     title = ax.set_title(
-        f"{kind} Exceedance \u2014 {date_str}",
+        f"{kind} {date_str}",
         loc="left",
         fontsize=13,
     )
 
     def _update(di: int):
         mesh.set_array(filt_mask[di].astype(float).ravel())
-        title.set_text(f"{kind} Exceedance \u2014 {str(dates[di])[:10]}")
+        title.set_text(f"{kind} {str(dates[di])[:10]}")
         if di % 10 == 0:
             logger.info(
                 "  Rendering frame %d / %d  (%s)",
@@ -139,7 +141,7 @@ def _update(di: int):
 def main() -> None:
     parser = argparse.ArgumentParser(
         description=(
-            "Animate daily heat-wave and cold-snap exceedance masks "
+            "Animate daily exceedance masks "
             "from ERA5 reanalysis."
         ),
     )
@@ -153,22 +155,6 @@ def main() -> None:
         required=True,
         help="End date YYYY-MM-DD",
     )
-    parser.add_argument(
-        "--output-heat",
-        default=None,
-        help=(
-            "Output GIF path for heat wave animation "
-            "(default: heat_exceedance_<start>_<end>.gif)"
-        ),
-    )
-    parser.add_argument(
-        "--output-cold",
-        default=None,
-        help=(
-            "Output GIF path for cold snap animation "
-            "(default: cold_exceedance_<start>_<end>.gif)"
-        ),
-    )
     parser.add_argument(
         "--fps",
         type=int,
@@ -181,13 +167,25 @@ def main() -> None:
         default=4,
         help="Number of dask workers (default: 4)",
     )
+    parser.add_argument(
+        "--quantile-lower",
+        type=float,
+        default=None,
+        help="Lower quantile for the event",
+    )
+    parser.add_argument(
+        "--quantile-upper",
+        type=float,
+        default=None,
+        help="Upper quantile for the event",
+    )
+    parser.add_argument(
+        "--output",
+        default="events_global.gif",
+        help="Output GIF path",
+    )
     args = parser.parse_args()
 
-    if args.output_heat is None:
-        args.output_heat = f"heat_exceedance_{args.start_date}_{args.end_date}.gif"
-    if args.output_cold is None:
-        args.output_cold = f"cold_exceedance_{args.start_date}_{args.end_date}.gif"
-
     wall_start = time_module.time()
     client = Client(LocalCluster(n_workers=args.n_workers))
     logger.info("Dask dashboard: %s", client.dashboard_link)
@@ -197,69 +195,46 @@ def main() -> None:
     logger.info("  sizes=%s", dict(t2m.sizes))
 
     logger.info("Loading climatology thresholds...")
-    clim_hw, _ = get_climatology_bounds(q_lower=0.85)
-    _, clim_fz = get_climatology_bounds(q_upper=0.15)
+    clim_lower, clim_upper = get_climatology_bounds(q_lower=args.quantile_lower, q_upper=args.quantile_upper)
 
     logger.info("Building land mask...")
     land_mask = build_land_mask(t2m.longitude, t2m.latitude)
 
     logger.info("Building exceedance masks (lazy)...")
-    hw_lazy = build_exceedance_mask(t2m, clim_hw, None, land_mask)
-    fz_lazy = build_exceedance_mask(t2m, None, clim_fz, land_mask)
-
-    tdim = detect_time_dim(hw_lazy)
+    event_lazy = build_exceedance_mask(t2m, clim_lower=clim_lower, clim_upper=clim_upper, land_mask=land_mask)
 
-    logger.info("Computing heat wave mask...")
-    t0 = time_module.time()
-    hw_da = hw_lazy.compute()
-    logger.info("  done in %.1f s", time_module.time() - t0)
+    tdim = detect_time_dim(event_lazy)
 
-    logger.info("Computing cold snap mask...")
+    logger.info("Computing event mask...")
     t0 = time_module.time()
-    fz_da = fz_lazy.compute()
+    event_da = event_lazy.compute()
     logger.info("  done in %.1f s", time_module.time() - t0)
 
-    dates = hw_da[tdim].values
-    lats = hw_da.latitude.values
-    lons = hw_da.longitude.values
-    hw_np = hw_da.values.astype(bool)
-    fz_np = fz_da.values.astype(bool)
-    del hw_da, fz_da
+    dates = event_da[tdim].values
+    lats = event_da.latitude.values
+    lons = event_da.longitude.values
+    event_np = event_da.values.astype(bool)
+    del event_da
 
     logger.info("Applying 3-day consecutive filter...")
-    hw_filt = apply_consecutive_filter(hw_np)
-    fz_filt = apply_consecutive_filter(fz_np)
+    event_filt = apply_consecutive_filter(event_np)
     logger.info(
         "  HW active cells: %d -> %d",
-        hw_np.sum(),
-        hw_filt.sum(),
-    )
-    logger.info(
-        "  FZ active cells: %d -> %d",
-        fz_np.sum(),
-        fz_filt.sum(),
+        event_np.sum(),
+        event_filt.sum(),
     )
-    del hw_np, fz_np
+    del event_np
 
     animate_exceedance(
-        hw_filt,
+        event_filt,
         dates,
         lats,
         lons,
-        "heat_wave",
-        args.output_heat,
+        lower_quantile=args.quantile_lower,
+        upper_quantile=args.quantile_upper,
+        output_path=args.output,
         fps=args.fps,
     )
-    animate_exceedance(
-        fz_filt,
-        dates,
-        lats,
-        lons,
-        "cold_snap",
-        args.output_cold,
-        fps=args.fps,
-    )
-
     client.close()
     logger.info("Done in %.1f s", time_module.time() - wall_start)
 

From 8cb22d2597679132c4f45519b6d042ac3065b498 Mon Sep 17 00:00:00 2001
From: Taylor Mandelbaum <mandelbaum.taylor@gmail.com>
Date: Wed, 8 Apr 2026 20:58:48 -0400
Subject: [PATCH 09/12] ruff

---
 data_prep/animate_heat_cold_global.py | 13 +++++++------
 1 file changed, 7 insertions(+), 6 deletions(-)

diff --git a/data_prep/animate_heat_cold_global.py b/data_prep/animate_heat_cold_global.py
index 3c53c8aa..7439e432 100644
--- a/data_prep/animate_heat_cold_global.py
+++ b/data_prep/animate_heat_cold_global.py
@@ -140,10 +140,7 @@ def _update(di: int):
 
 def main() -> None:
     parser = argparse.ArgumentParser(
-        description=(
-            "Animate daily exceedance masks "
-            "from ERA5 reanalysis."
-        ),
+        description=("Animate daily exceedance masks from ERA5 reanalysis."),
     )
     parser.add_argument(
         "--start-date",
@@ -195,13 +192,17 @@ def main() -> None:
     logger.info("  sizes=%s", dict(t2m.sizes))
 
     logger.info("Loading climatology thresholds...")
-    clim_lower, clim_upper = get_climatology_bounds(q_lower=args.quantile_lower, q_upper=args.quantile_upper)
+    clim_lower, clim_upper = get_climatology_bounds(
+        q_lower=args.quantile_lower, q_upper=args.quantile_upper
+    )
 
     logger.info("Building land mask...")
     land_mask = build_land_mask(t2m.longitude, t2m.latitude)
 
     logger.info("Building exceedance masks (lazy)...")
-    event_lazy = build_exceedance_mask(t2m, clim_lower=clim_lower, clim_upper=clim_upper, land_mask=land_mask)
+    event_lazy = build_exceedance_mask(
+        t2m, clim_lower=clim_lower, clim_upper=clim_upper, land_mask=land_mask
+    )
 
     tdim = detect_time_dim(event_lazy)
 

From d09d068679f1e26dc2bc955b739dc683c75ec570 Mon Sep 17 00:00:00 2001
From: Taylor Mandelbaum <mandelbaum.taylor@gmail.com>
Date: Wed, 8 Apr 2026 21:04:30 -0400
Subject: [PATCH 10/12] remove animate script

---
 data_prep/animate_heat_cold_global.py | 244 --------------------------
 1 file changed, 244 deletions(-)
 delete mode 100644 data_prep/animate_heat_cold_global.py

diff --git a/data_prep/animate_heat_cold_global.py b/data_prep/animate_heat_cold_global.py
deleted file mode 100644
index 7439e432..00000000
--- a/data_prep/animate_heat_cold_global.py
+++ /dev/null
@@ -1,244 +0,0 @@
-"""Animate daily heat-wave and cold-snap exceedance masks from ERA5.
-
-Produces two animated GIFs — one for heat waves, one for cold snaps —
-showing the filtered (3+ consecutive day) exceedance mask day by day on
-a global Robinson projection.
-
-Usage:
-    python animate_heat_cold_global.py \\
-        --start-date 2020-01-01 --end-date 2020-03-01 \\
-        --fps 4 --n-workers 4 --quantile-lower 0.85 --quantile-upper 0.15
-        --output events_global.gif
-"""
-
-import argparse
-import logging
-import time as time_module
-
-import cartopy.crs as ccrs
-import cartopy.feature as cfeature
-import matplotlib.animation as mpl_anim
-import matplotlib.colors as mcolors
-import matplotlib.pyplot as plt
-import numpy as np
-from dask.distributed import Client, LocalCluster
-from heat_cold_bounds_global import (
-    apply_consecutive_filter,
-    build_exceedance_mask,
-    build_land_mask,
-    get_climatology_bounds,
-)
-from plot_temperature_events import (
-    _add_map_features,
-    detect_time_dim,
-    open_era5_t2m,
-)
-
-logging.basicConfig(
-    level=logging.INFO,
-    format="%(asctime)s %(name)s %(levelname)s %(message)s",
-)
-logger = logging.getLogger(__name__)
-
-_EVENT_COLOR = "#d73027"
-
-
-def animate_exceedance(
-    filt_mask: np.ndarray,
-    dates: np.ndarray,
-    lats: np.ndarray,
-    lons: np.ndarray,
-    lower_quantile: float | None,
-    upper_quantile: float | None,
-    output_path: str,
-    fps: int = 4,
-) -> None:
-    """Save an animated GIF of the daily filtered exceedance mask.
-
-    Each frame shows one day's exceedance footprint on a global
-    Robinson projection. Active grid points are colored red; inactive land and ocean are rendered
-    in whitesmoke and light blue respectively.
-
-    Args:
-        filt_mask: Boolean array (time, lat, lon) with the
-            consecutive-day filter already applied.
-        dates: 1-D datetime64 array aligned with axis 0 of
-            filt_mask.
-        lats: 1-D latitude array (ascending).
-        lons: 1-D longitude array (0-360 degrees).
-        lower_quantile: Lower quantile for the event.
-        upper_quantile: Upper quantile for the event.
-        output_path: Destination ``.gif`` file path.
-        fps: Frames per second for the output GIF. Default 4.
-    """
-    n_days = filt_mask.shape[0]
-
-    rgba = mcolors.to_rgba(_EVENT_COLOR)
-    cmap = mcolors.ListedColormap(["none", rgba])
-
-    fig, ax = plt.subplots(
-        subplot_kw={"projection": ccrs.Robinson()},
-        figsize=(12, 5.5),
-    )
-    fig.subplots_adjust(left=0.01, right=0.99, top=0.93, bottom=0.01)
-    ax.set_global()
-    ax.add_feature(cfeature.OCEAN, facecolor="lightblue", zorder=0)
-    ax.add_feature(cfeature.LAND, facecolor="whitesmoke", zorder=0)
-    _add_map_features(ax)
-
-    data0 = filt_mask[0].astype(float)
-    mesh = ax.pcolormesh(
-        lons,
-        lats,
-        data0,
-        cmap=cmap,
-        vmin=0,
-        vmax=1,
-        transform=ccrs.PlateCarree(),
-        shading="auto",
-        zorder=1,
-    )
-    if lower_quantile is not None and upper_quantile is not None:
-        kind = f"p{lower_quantile:.2f}-p{upper_quantile:.2f}"
-    elif lower_quantile is not None:
-        kind = f"p{lower_quantile:.2f}+"
-    else:
-        kind = f"p{upper_quantile:.2f}-"
-    date_str = str(dates[0])[:10]
-    title = ax.set_title(
-        f"{kind} {date_str}",
-        loc="left",
-        fontsize=13,
-    )
-
-    def _update(di: int):
-        mesh.set_array(filt_mask[di].astype(float).ravel())
-        title.set_text(f"{kind} {str(dates[di])[:10]}")
-        if di % 10 == 0:
-            logger.info(
-                "  Rendering frame %d / %d  (%s)",
-                di + 1,
-                n_days,
-                str(dates[di])[:10],
-            )
-        return mesh, title
-
-    anim = mpl_anim.FuncAnimation(
-        fig,
-        _update,
-        frames=n_days,
-        interval=1000 // fps,
-        blit=False,
-    )
-
-    writer = mpl_anim.PillowWriter(fps=fps)
-    logger.info("Saving %s animation to %s ...", kind, output_path)
-    anim.save(output_path, writer=writer, dpi=100)
-    plt.close(fig)
-    logger.info("  Saved %s", output_path)
-
-
-def main() -> None:
-    parser = argparse.ArgumentParser(
-        description=("Animate daily exceedance masks from ERA5 reanalysis."),
-    )
-    parser.add_argument(
-        "--start-date",
-        required=True,
-        help="Start date YYYY-MM-DD",
-    )
-    parser.add_argument(
-        "--end-date",
-        required=True,
-        help="End date YYYY-MM-DD",
-    )
-    parser.add_argument(
-        "--fps",
-        type=int,
-        default=4,
-        help="Frames per second for the output GIFs (default: 4)",
-    )
-    parser.add_argument(
-        "--n-workers",
-        type=int,
-        default=4,
-        help="Number of dask workers (default: 4)",
-    )
-    parser.add_argument(
-        "--quantile-lower",
-        type=float,
-        default=None,
-        help="Lower quantile for the event",
-    )
-    parser.add_argument(
-        "--quantile-upper",
-        type=float,
-        default=None,
-        help="Upper quantile for the event",
-    )
-    parser.add_argument(
-        "--output",
-        default="events_global.gif",
-        help="Output GIF path",
-    )
-    args = parser.parse_args()
-
-    wall_start = time_module.time()
-    client = Client(LocalCluster(n_workers=args.n_workers))
-    logger.info("Dask dashboard: %s", client.dashboard_link)
-
-    logger.info("Opening ERA5 data...")
-    t2m = open_era5_t2m(args.start_date, args.end_date)
-    logger.info("  sizes=%s", dict(t2m.sizes))
-
-    logger.info("Loading climatology thresholds...")
-    clim_lower, clim_upper = get_climatology_bounds(
-        q_lower=args.quantile_lower, q_upper=args.quantile_upper
-    )
-
-    logger.info("Building land mask...")
-    land_mask = build_land_mask(t2m.longitude, t2m.latitude)
-
-    logger.info("Building exceedance masks (lazy)...")
-    event_lazy = build_exceedance_mask(
-        t2m, clim_lower=clim_lower, clim_upper=clim_upper, land_mask=land_mask
-    )
-
-    tdim = detect_time_dim(event_lazy)
-
-    logger.info("Computing event mask...")
-    t0 = time_module.time()
-    event_da = event_lazy.compute()
-    logger.info("  done in %.1f s", time_module.time() - t0)
-
-    dates = event_da[tdim].values
-    lats = event_da.latitude.values
-    lons = event_da.longitude.values
-    event_np = event_da.values.astype(bool)
-    del event_da
-
-    logger.info("Applying 3-day consecutive filter...")
-    event_filt = apply_consecutive_filter(event_np)
-    logger.info(
-        "  HW active cells: %d -> %d",
-        event_np.sum(),
-        event_filt.sum(),
-    )
-    del event_np
-
-    animate_exceedance(
-        event_filt,
-        dates,
-        lats,
-        lons,
-        lower_quantile=args.quantile_lower,
-        upper_quantile=args.quantile_upper,
-        output_path=args.output,
-        fps=args.fps,
-    )
-    client.close()
-    logger.info("Done in %.1f s", time_module.time() - wall_start)
-
-
-if __name__ == "__main__":
-    main()

From bb2dd8a20aedfd76f0c130b9f14e43ab14523e69 Mon Sep 17 00:00:00 2001
From: Taylor Mandelbaum <mandelbaum.taylor@gmail.com>
Date: Wed, 8 Apr 2026 21:09:54 -0400
Subject: [PATCH 11/12] rename

---
 ...heat_cold_bounds_global.py => temperature_bounds_global.py} | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)
 rename data_prep/{heat_cold_bounds_global.py => temperature_bounds_global.py} (99%)

diff --git a/data_prep/heat_cold_bounds_global.py b/data_prep/temperature_bounds_global.py
similarity index 99%
rename from data_prep/heat_cold_bounds_global.py
rename to data_prep/temperature_bounds_global.py
index a4d5cb99..d8d10689 100644
--- a/data_prep/heat_cold_bounds_global.py
+++ b/data_prep/temperature_bounds_global.py
@@ -34,7 +34,6 @@
 from typing import Dict, List, Optional, Tuple
 
 import joblib
-
 import numba as nb
 import numpy as np
 import pandas as pd
@@ -62,7 +61,7 @@
 MIN_CONSECUTIVE_DAYS = 3
 AREA_DECLINE_FRACTION = 0.5
 MIN_GRIDPOINTS = 500
-MIN_AREA_KM2 = 200_000.0
+MIN_AREA_KM2 = 200000.0
 EXPANSION_DEGREES = 1
 MAX_SPATIAL_ITERATIONS = 20
 EDGE_VALIDITY_THRESHOLD = 0.5

From 91b5662d9f3e62cbe284ed22dd543152b2bd6225 Mon Sep 17 00:00:00 2001
From: Taylor Mandelbaum <mandelbaum.taylor@gmail.com>
Date: Wed, 8 Apr 2026 21:11:50 -0400
Subject: [PATCH 12/12] update with rename

---
 data_prep/plot_temperature_events.py   | 2 +-
 data_prep/temperature_bounds_global.py | 4 ++--
 2 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/data_prep/plot_temperature_events.py b/data_prep/plot_temperature_events.py
index 6eb2bb17..e8c8708a 100644
--- a/data_prep/plot_temperature_events.py
+++ b/data_prep/plot_temperature_events.py
@@ -4,7 +4,7 @@
 Handles both event types from a single entry point. The event type
 is auto-detected from events.yaml; no flag required.
 
-Exported functions used by heat_cold_bounds_global.py and
+Exported functions used by temperature_bounds_global.py and
 heat_cold_bounds_case.py:
     max_consecutive_days  -- compute field from boolean mask
     plot_consecutive_map  -- pcolormesh map (Reds/Blues, discrete)
diff --git a/data_prep/temperature_bounds_global.py b/data_prep/temperature_bounds_global.py
index d8d10689..ce7e91d6 100644
--- a/data_prep/temperature_bounds_global.py
+++ b/data_prep/temperature_bounds_global.py
@@ -15,13 +15,13 @@
 
 Usage:
     # Anything above the 85th percentile (heat wave)
-    python heat_cold_bounds_global.py \\
+    python temperature_bounds_global.py \\
         --start-date 2023-06-01 --end-date 2023-09-01 \\
         --quantile-lower 0.85 --operator-lower ">=" \\
         --event-type heat_wave --output heat_cold_global.csv
 
     # Band between 50th and 85th (moderate heat)
-    python heat_cold_bounds_global.py \\
+    python temperature_bounds_global.py \\
         --start-date 2023-06-01 --end-date 2023-09-01 \\
         --quantile-lower 0.50 --quantile-upper 0.85 \\
         --event-type heat_wave --output heat_cold_global.csv