The live monitoring code computes storm rainfall differently from the historical analysis. The historical analysis uses pre-computed values in zma_stats_imerg_quantiles.parquet, calculated from IBTrACS tracks for 2000-2024. Monitoring computes them on the fly with IMERGRasterProcessor, as the 80th percentile of two-day rolling sums of IMERG rainfall over Cuba. This appendix checks that the two give the same numbers.
The check takes the historical storms that reach the observational wind threshold (105 kt) and pass through the ZMA, recomputes their rainfall with the monitoring code, and compares it with the pre-computed q80_roll2 value for the same IBTrACS storm ID. To keep the run short it processes the first 10 of the 11 qualifying storms.
All 10 match the pre-computed values, with a difference of 0.00 mm.
import pandas as pd
import numpy as np
import ocha_stratus as stratus
from src.datasources import nhc, ibtracs
from src.monitoring.monitoring_utils import (
CubaHurricaneMonitor,
IMERGRasterProcessor,
)
from src.constants import PROJECT_PREFIX, THRESHS
import logging
logging.getLogger("src.monitoring.monitoring_utils").setLevel(logging.WARNING)print("Rainfall Aggregation Validation")
print("=" * 50)
# Load pre-computed aggregations from exploration
blob_name = (
f"{PROJECT_PREFIX}/processed/storm_stats/zma_stats_imerg_quantiles.parquet"
)
try:
df_precomputed = stratus.load_parquet_from_blob(blob_name)
print(f"Loaded {len(df_precomputed)} storms from pre-computed data")
print(
f" Date range: {df_precomputed['valid_time_min'].min().date()} to {df_precomputed['valid_time_max'].max().date()}"
)
except Exception as e:
print(f"Error loading pre-computed data: {e}")
df_precomputed = None
# Load historical observational tracks from IBTrACS (same as used in imerg_aggregation.md)
try:
# Load the same historical data used to create the pre-computed aggregations
ibtracs_blob_name = (
f"{PROJECT_PREFIX}/processed/ibtracs/zma_tracks_2000-2024.parquet"
)
obsv_tracks = stratus.load_parquet_from_blob(ibtracs_blob_name)
# Load storm metadata to get SID to ATCF mapping (following your approach)
df_storms = ibtracs.load_storms()
# Filter to relevant storms and get SID-ATCF mapping
storm_mapping = df_storms[["sid", "atcf_id", "name"]].drop_duplicates()
storm_mapping = storm_mapping[
storm_mapping["sid"].notna() & storm_mapping["atcf_id"].notna()
]
# Merge tracks with storm mapping to get ATCF IDs
obsv_tracks = obsv_tracks.merge(storm_mapping, on="sid", how="left")
print(
f"Loaded {len(obsv_tracks)} track points from {obsv_tracks['sid'].nunique()} storms"
)
print(
f" Date range: {obsv_tracks['valid_time'].min().date()} to {obsv_tracks['valid_time'].max().date()}"
)
print(
f" Storms with ATCF mapping: {obsv_tracks['atcf_id'].notna().sum()}/{len(obsv_tracks)}"
)
# Rename columns to match what monitoring_utils expects
obsv_tracks = obsv_tracks.rename(
columns={"valid_time": "lastUpdate", "wind_speed": "intensity"}
)
# Keep original SID for exact matching with pre-computed data
obsv_tracks["original_sid"] = obsv_tracks["sid"].copy()
except Exception as e:
print(f"Error loading IBTrACS historical data: {e}")
print(" Falling back to recent NHC data...")
obsv_tracks = nhc.load_recent_glb_obsv()
obsv_tracks = obsv_tracks[obsv_tracks["basin"] == "al"]
obsv_tracks = obsv_tracks.rename(columns={"id": "atcf_id"})
obsv_tracks["original_sid"] = None # No SID for recent NHC data
print(
f"Loaded {len(obsv_tracks)} recent track points from {obsv_tracks['atcf_id'].nunique()} storms"
)
print(f"\nPre-filtering storms by monitoring criteria...")
print(f" Wind threshold: {THRESHS['obsv']['s']} kt")
print(f" Must intersect ZMA")
# Create monitor instance to use existing methods
monitor = CubaHurricaneMonitor()
raster_processor = IMERGRasterProcessor(quantile=0.8)
# Step 1: Pre-filter to find ALL qualifying storms first
qualifying_storms = []
total_storms = obsv_tracks["sid"].nunique()
# Debug: Let's analyze the filtering process step by step
wind_qualified = []
zma_qualified = []
both_qualified = []
for sid, group in obsv_tracks.groupby("sid"):
try:
# First, apply wind threshold filter (cheap operation)
max_wind = group["intensity"].max()
if max_wind >= THRESHS["obsv"]["s"]:
wind_qualified.append(sid)
# Let's also check distance to Cuba for all storms
group_clean = monitor._remove_track_duplicates(
group.sort_values("lastUpdate"), "lastUpdate"
)
gdf_simple = monitor._create_track_geodataframe(group_clean)
min_distance = gdf_simple["distance"].min()
# Check ZMA intersection
has_zma_intersection = not monitor._filter_by_zma(gdf_simple).empty
if has_zma_intersection:
zma_qualified.append(sid)
# Only proceed with storms that meet BOTH criteria
if max_wind >= THRESHS["obsv"]["s"] and has_zma_intersection:
both_qualified.append(sid)
# Store qualifying storm info
storm_start = group["lastUpdate"].min()
storm_end = group["lastUpdate"].max()
qualifying_storms.append(
{
"sid": sid,
"atcf_id": (
group["atcf_id"].iloc[0]
if "atcf_id" in group.columns
else None
),
"storm_start": storm_start,
"storm_end": storm_end,
"max_wind": max_wind,
"min_distance": min_distance,
"group_data": group, # Store for later processing
}
)
except Exception as e:
print(f" {sid}: Error in pre-filtering - {e}")
print(f"\nFiltering Analysis:")
print(f" Total storms: {total_storms}")
print(f" Storms with wind ≥{THRESHS['obsv']['s']} kt: {len(wind_qualified)}")
print(f" Storms intersecting ZMA: {len(zma_qualified)}")
print(f" Storms meeting BOTH criteria: {len(both_qualified)}")
print(f" Final qualifying storms: {len(qualifying_storms)}")
# If we have too few qualifying storms, let's try a distance-based approach
if len(qualifying_storms) < 10:
print(
f"\nToo few qualifying storms. Let's try distance-based filtering..."
)
print(f" Looking for storms within 200km of Cuba with wind ≥80kt...")
distance_qualified = []
for sid, group in obsv_tracks.groupby("sid"):
try:
max_wind = group["intensity"].max()
if max_wind < 80: # Lower wind threshold
continue
group_clean = monitor._remove_track_duplicates(
group.sort_values("lastUpdate"), "lastUpdate"
)
gdf_simple = monitor._create_track_geodataframe(group_clean)
min_distance = gdf_simple["distance"].min()
if min_distance < 200: # 200km distance threshold
storm_start = group["lastUpdate"].min()
storm_end = group["lastUpdate"].max()
distance_qualified.append(
{
"sid": sid,
"atcf_id": (
group["atcf_id"].iloc[0]
if "atcf_id" in group.columns
else None
),
"storm_start": storm_start,
"storm_end": storm_end,
"max_wind": max_wind,
"min_distance": min_distance,
"group_data": group,
}
)
except Exception as e:
continue
print(f" Distance-qualified storms: {len(distance_qualified)}")
if len(distance_qualified) > len(qualifying_storms):
print(f" Using distance-qualified storms for analysis...")
qualifying_storms = distance_qualified
print(f"\nFinal dataset: {len(qualifying_storms)} storms for analysis")
# Step 2: Now sample from qualifying storms for detailed processing
sample_size = min(10, len(qualifying_storms))
for i, storm_info in enumerate(qualifying_storms[:sample_size]):
sid = storm_info["sid"]
atcf_id = storm_info["atcf_id"]
max_wind = storm_info["max_wind"]
group = storm_info["group_data"]
try:
# Now do the expensive interpolation for detailed analysis
group_clean = monitor._remove_track_duplicates(
group.sort_values("lastUpdate"), "lastUpdate"
)
df_interp = monitor._interpolate_track(
group_clean,
"lastUpdate",
["latitude", "longitude", "intensity", "pressure"],
)
# Create full GeoDataFrame and get detailed stats
gdf = monitor._create_track_geodataframe(df_interp)
gdf_zma = monitor._filter_by_zma(gdf)
# Store storm info for comparison - use original group data for dates
storm_start = storm_info["storm_start"]
storm_end = storm_info["storm_end"]
# Update the storm info with detailed stats
storm_info.update(
{
"min_dist": gdf["distance"].min(),
}
)
except Exception as e:
print(f" {sid}: Error in detailed processing - {e}")
print(f"\nFinal Summary:")
print(f" Total storms in dataset: {total_storms}")
print(f" Storms available for analysis: {len(qualifying_storms)}")
print(f" Storms processed for detailed analysis: {sample_size}")
if len(qualifying_storms) > 0:
validation_results = []
# Process more storms for better validation - use at least 5-10 storms
rainfall_sample_size = min(10, len(qualifying_storms))
for storm_info in qualifying_storms[:rainfall_sample_size]:
sid = storm_info["sid"]
atcf_id = storm_info["atcf_id"]
storm_start = storm_info["storm_start"]
storm_end = storm_info["storm_end"]
# Extend period by 1 day on each side (same as monitoring_utils)
analysis_start = storm_start - pd.Timedelta(days=1)
analysis_end = storm_end + pd.Timedelta(days=1)
try:
# Calculate rainfall using monitoring_utils method
storm_rainfall_df = (
raster_processor.calculate_rainfall_for_storm_period(
analysis_start, analysis_end
)
)
if storm_rainfall_df.empty:
print(f" No rainfall data for {sid}")
continue
# Get max rainfall for the period (same as monitoring_utils)
max_rainfall = storm_rainfall_df["roll2_sum"].max()
validation_results.append(
{
"sid": sid,
"atcf_id": atcf_id,
"storm_start": storm_start,
"storm_end": storm_end,
"onthefly_max_roll2": max_rainfall,
"analysis_start": analysis_start,
"analysis_end": analysis_end,
}
)
except Exception as e:
print(f" Error calculating rainfall for {sid}: {e}")
import traceback
traceback.print_exc()
# Convert to DataFrame for easier comparison
df_validation = pd.DataFrame(validation_results)
if len(df_validation) > 0 and df_precomputed is not None:
print(f"\nComparing with pre-computed values...")
print(f" On-the-fly results: {len(df_validation)} storms")
print(f" Pre-computed results: {len(df_precomputed)} storms")
# Note: This might need adjustment based on the actual SID format in precomputed data
print(f"\nValidation results:")
for _, row in df_validation.iterrows():
print(
f" {row['sid']} (ATCF: {row['atcf_id']}): On-the-fly max roll2 = {row['onthefly_max_roll2']:.2f} mm"
)
# Use exact SID matching with pre-computed data
original_sid = row["sid"]
potential_match = df_precomputed[
df_precomputed["sid"] == original_sid
]
if len(potential_match) > 0:
precomputed_val = potential_match.iloc[0]["q80_roll2"]
diff = abs(row["onthefly_max_roll2"] - precomputed_val)
print(
f" Pre-computed q80_roll2 = {precomputed_val:.2f} mm"
)
print(f" Difference = {diff:.2f} mm")
print(f" Exact SID match: {original_sid}")
else:
print(f" No exact SID match found for: {original_sid}")
else:
print(f"\nNo qualifying storms found in the sample. Try:")
print(f" - Increasing the storm limit (currently 10)")
print(f" - Using more recent data")
print(f" - Checking the wind/ZMA thresholds")Rainfall Aggregation Validation
==================================================
Loaded 62 storms from pre-computed data
Date range: 2000-08-24 to 2024-11-07
Loaded 347 track points from 62 storms
Date range: 2000-08-24 to 2024-11-07
Storms with ATCF mapping: 347/347
Pre-filtering storms by monitoring criteria...
Wind threshold: 105 kt
Must intersect ZMA
Filtering Analysis:
Total storms: 62
Storms with wind ≥105 kt: 11
Storms intersecting ZMA: 62
Storms meeting BOTH criteria: 11
Final qualifying storms: 11
Final dataset: 11 storms for analysisFinal Summary:
Total storms in dataset: 62
Storms available for analysis: 11
Storms processed for detailed analysis: 10Comparing with pre-computed values...
On-the-fly results: 10 storms
Pre-computed results: 62 storms
Validation results:
2001303N13276 (ATCF: AL152001): On-the-fly max roll2 = 196.17 mm
Pre-computed q80_roll2 = 196.17 mm
Difference = 0.00 mm
Exact SID match: 2001303N13276
2002258N10300 (ATCF: AL102002): On-the-fly max roll2 = 96.10 mm
Pre-computed q80_roll2 = 96.10 mm
Difference = 0.00 mm
Exact SID match: 2002258N10300
2004223N11301 (ATCF: AL032004): On-the-fly max roll2 = 44.66 mm
Pre-computed q80_roll2 = 44.66 mm
Difference = 0.00 mm
Exact SID match: 2004223N11301
2004247N10332 (ATCF: AL092004): On-the-fly max roll2 = 65.58 mm
Pre-computed q80_roll2 = 65.58 mm
Difference = 0.00 mm
Exact SID match: 2004247N10332
2005186N12299 (ATCF: AL042005): On-the-fly max roll2 = 153.27 mm
Pre-computed q80_roll2 = 153.27 mm
Difference = 0.00 mm
Exact SID match: 2005186N12299
2008238N13293 (ATCF: AL072008): On-the-fly max roll2 = 71.27 mm
Pre-computed q80_roll2 = 71.27 mm
Difference = 0.00 mm
Exact SID match: 2008238N13293
2008245N17323 (ATCF: AL092008): On-the-fly max roll2 = 169.01 mm
Pre-computed q80_roll2 = 169.01 mm
Difference = 0.00 mm
Exact SID match: 2008245N17323
2008311N14278 (ATCF: AL172008): On-the-fly max roll2 = 103.72 mm
Pre-computed q80_roll2 = 103.72 mm
Difference = 0.00 mm
Exact SID match: 2008311N14278
2016273N13300 (ATCF: AL142016): On-the-fly max roll2 = 96.22 mm
Pre-computed q80_roll2 = 96.22 mm
Difference = 0.00 mm
Exact SID match: 2016273N13300
2017242N16333 (ATCF: AL112017): On-the-fly max roll2 = 145.04 mm
Pre-computed q80_roll2 = 145.04 mm
Difference = 0.00 mm
Exact SID match: 2017242N16333