The quantile calculations were refactored due to the conceptual change in logic between realtime monitoring and historical analysis. Additionally, there are different data sets being used for historical analysis and realtime monitoring. Therefore, this script is meant to sanity check the calculations used for real time monitoring.
During the development of the hurricane monitoring system, two different approaches emerged for calculating rainfall quantiles:
Historical Analysis: Uses pre-computed quantile values stored in zma_stats_imerg_quantiles.parquet, calculated from comprehensive historical storm data (2000-2024) using IBTrACS tracks.
Real-time Monitoring: Uses on-the-fly calculations via IMERGRasterProcessor with 0.8 quantile and 2-day rolling sums, applied to IMERG rainfall data during active storm monitoring.
This validation compares the two calculation methods to ensure consistency between historical baselines and real-time assessments. The test:
IMERGRasterProcessor to compute rainfall metrics for qualifying stormsThe validation framework includes several optimizations and debugging features:
This ensures that real-time monitoring calculations produce results consistent with historical analysis, maintaining accuracy and reliability in operational hurricane monitoring.
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, THRESHSprint("🔍 Rainfall Aggregation Validation")
print("=" * 50)
# Load pre-computed aggregations from exploration
blob_name = (
f"{PROJECT_PREFIX}/processed/storm_stats/zma_stats_imerg_quantiles.parquet"
)
print(f"📊 Loading pre-computed aggregations from: {blob_name}")
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)
print(f"\n🌀 Loading historical observational tracks from IBTrACS...")
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)
print("📋 Loading SID to ATCF mapping...")
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"\n🎯 Applying monitoring criteria to identify qualifying storms...")
print(f" Wind threshold: {THRESHS['obsv']['s']} kt")
print(f" Must intersect ZMA")
# Create monitor instance to use existing methods
monitor = CubaHurricaneMonitor()
qualifying_storms = []
raster_processor = IMERGRasterProcessor(quantile=0.8)
# monitor._filter_by_zma
# monitor.filter
# Check each storm against monitoring criteria
print(f"\n🔄 Processing storms to find qualifying ones...")
for i, (sid, group) in enumerate(
obsv_tracks.groupby("sid")
): # Group by SID instead of atcf_id
if i >= 10: # Limit to first 10 storms for testing
print(f" Limiting to first 10 storms for validation...")
break
print(f" Checking {sid}...")
try:
# First, apply wind threshold filter (cheap operation)
max_wind = group["intensity"].max()
if max_wind < THRESHS["obsv"]["s"]:
print(
f" ❌ {sid}: Wind {max_wind:.1f} kt < {THRESHS['obsv']['s']} kt (pre-filter)"
)
continue
print(f" ✅ {sid}: Passes wind threshold ({max_wind:.1f} kt)")
# Quick boolean check: does this storm have ANY track points in ZMA? (cheap operation)
group_clean = monitor._remove_track_duplicates(
group.sort_values("lastUpdate"), "lastUpdate"
)
# Create minimal GeoDataFrame just for ZMA intersection check
gdf_simple = monitor._create_track_geodataframe(group_clean)
has_zma_intersection = not monitor._filter_by_zma(gdf_simple).empty
if not has_zma_intersection:
print(f" ❌ {sid}: No ZMA intersection (pre-filter)")
continue
print(f" ✅ {sid}: Has ZMA intersection, processing full track...")
# Now do the expensive interpolation only for storms that pass both filters
df_interp = monitor._interpolate_track(
group_clean,
"lastUpdate",
["latitude", "longitude", "intensity", "pressure"],
)
# Create full GeoDataFrame and get ZMA-intersecting segments
gdf = monitor._create_track_geodataframe(df_interp)
gdf_zma = monitor._filter_by_zma(gdf)
print(
f" ✅ {sid}: QUALIFIES (wind={max_wind:.1f} kt, ZMA intersection)"
)
# Store storm info for comparison - use original group data for dates
storm_start = group["lastUpdate"].min()
storm_end = group["lastUpdate"].max()
# Debug output to check date values
print(
f" Debug: storm_start = {storm_start}, storm_end = {storm_end}"
)
print(
f" Debug: date range = {storm_start.date() if pd.notna(storm_start) else 'NaT'} to {storm_end.date() if pd.notna(storm_end) else 'NaT'}"
)
qualifying_storms.append(
{
"sid": sid,
"atcf_id": (
group["atcf_id"].iloc[0]
if "atcf_id" in group.columns
else None
), # Keep ATCF ID for reference
"storm_start": storm_start,
"storm_end": storm_end,
"max_wind": max_wind,
"min_dist": gdf["distance"].min(),
}
)
except Exception as e:
print(f" ⚠️ {sid}: Error processing - {e}")
print(f"\n📊 Summary:")
print(
f" Total storms processed: {min(10, obsv_tracks['atcf_id'].nunique())}"
)
print(f" Qualifying storms: {len(qualifying_storms)}")
if len(qualifying_storms) > 0:
print(
f"\n🔄 Running on-the-fly rainfall calculations for qualifying storms..."
)
validation_results = []
for storm_info in qualifying_storms[:3]: # Limit to first 3 for speed
sid = storm_info["sid"]
atcf_id = storm_info["atcf_id"]
storm_start = storm_info["storm_start"]
storm_end = storm_info["storm_end"]
print(f" Processing {sid} (ATCF: {atcf_id})...")
# 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)
print(f" Storm period: {storm_start.date()} to {storm_end.date()}")
print(
f" Analysis period: {analysis_start.date()} to {analysis_end.date()}"
)
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()
print(f" On-the-fly calculation: {max_rainfall:.2f} mm")
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"\n🔍 Comparing 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"\n📊 Validation 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}")
print(f"\n✅ Validation complete!")
print(
f"📋 Note: This is a basic validation. Full comparison would require:"
)
print(f" - Exact SID matching between datasets")
print(f" - Same date ranges and quantile calculations")
print(f" - Processing all storms, not just a subset")
else:
print(f"\n⚠️ No 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
==================================================
📊 Loading pre-computed aggregations from: ds-aa-cub-hurricanes/processed/storm_stats/zma_stats_imerg_quantiles.parquet
✅ Loaded 62 storms from pre-computed data
Date range: 2000-08-24 to 2024-11-07
🌀 Loading historical observational tracks from IBTrACS...
📋 Loading SID to ATCF mapping...
✅ Loaded 347 track points from 62 storms
Date range: 2000-08-24 to 2024-11-07
Storms with ATCF mapping: 347/347
🎯 Applying monitoring criteria to identify qualifying storms...
Wind threshold: 105 kt
Must intersect ZMA
🔄 Processing storms to find qualifying ones...
Checking 2000233N12316...
❌ 2000233N12316: Wind 40.0 kt < 105 kt (pre-filter)
Checking 2000260N15308...
❌ 2000260N15308: Wind 30.0 kt < 105 kt (pre-filter)
Checking 2001303N13276...
✅ 2001303N13276: Passes wind threshold (120.0 kt)
✅ 2001303N13276: Has ZMA intersection, processing full track...
✅ 2001303N13276: QUALIFIES (wind=120.0 kt, ZMA intersection)
Debug: storm_start = 2001-11-04 06:00:00.000040, storm_end = 2001-11-05 06:00:00.000040
Debug: date range = 2001-11-04 to 2001-11-05
Checking 2002258N10300...
✅ 2002258N10300: Passes wind threshold (110.0 kt)
✅ 2002258N10300: Has ZMA intersection, processing full track.../Users/zackarno/.pyenv/versions/3.11.4/envs/ds-aa-cub-hurricanes/lib/python3.11/site-packages/shapely/measurement.py:81: RuntimeWarning:
invalid value encountered in distance
/Users/zackarno/.pyenv/versions/3.11.4/envs/ds-aa-cub-hurricanes/lib/python3.11/site-packages/shapely/measurement.py:81: RuntimeWarning:
invalid value encountered in distance
/Users/zackarno/.pyenv/versions/3.11.4/envs/ds-aa-cub-hurricanes/lib/python3.11/site-packages/shapely/measurement.py:81: RuntimeWarning:
invalid value encountered in distance
✅ 2002258N10300: QUALIFIES (wind=110.0 kt, ZMA intersection)
Debug: storm_start = 2002-09-19 18:00:00.000040, storm_end = 2002-09-21 18:00:00.000040
Debug: date range = 2002-09-19 to 2002-09-21
Checking 2002265N10315...
❌ 2002265N10315: Wind 90.0 kt < 105 kt (pre-filter)
Checking 2002288N17277...
❌ 2002288N17277: Wind 30.0 kt < 105 kt (pre-filter)
Checking 2004217N13306...
❌ 2004217N13306: Wind 25.0 kt < 105 kt (pre-filter)
Checking 2004223N11301...
✅ 2004223N11301: Passes wind threshold (105.0 kt)
✅ 2004223N11301: Has ZMA intersection, processing full track...
✅ 2004223N11301: QUALIFIES (wind=105.0 kt, ZMA intersection)
Debug: storm_start = 2004-08-12 18:00:00.000040, storm_end = 2004-08-13 06:00:00.000040
Debug: date range = 2004-08-12 to 2004-08-13
Checking 2004247N10332...
✅ 2004247N10332: Passes wind threshold (140.0 kt)
✅ 2004247N10332: Has ZMA intersection, processing full track...
✅ 2004247N10332: QUALIFIES (wind=140.0 kt, ZMA intersection)
Debug: storm_start = 2004-09-13 12:00:00.000040, storm_end = 2004-09-14 06:00:00.000040
Debug: date range = 2004-09-13 to 2004-09-14
Checking 2004258N16300...
❌ 2004258N16300: Wind 45.0 kt < 105 kt (pre-filter)
Limiting to first 10 storms for validation...
📊 Summary:
Total storms processed: 10
Qualifying storms: 4
🔄 Running on-the-fly rainfall calculations for qualifying storms...
Processing 2001303N13276 (ATCF: AL152001)...
Storm period: 2001-11-04 to 2001-11-05
Analysis period: 2001-11-03 to 2001-11-06
2025-08-06 15:09:29 [INFO] src.monitoring.monitoring_utils: Loading IMERG raster data for 4 dates/Users/zackarno/.pyenv/versions/3.11.4/envs/ds-aa-cub-hurricanes/lib/python3.11/site-packages/shapely/measurement.py:81: RuntimeWarning:
invalid value encountered in distance
2025-08-06 15:09:55 [INFO] src.monitoring.monitoring_utils: Successfully loaded IMERG data for 4 dates
On-the-fly calculation: 196.17 mm
Processing 2002258N10300 (ATCF: AL102002)...
Storm period: 2002-09-19 to 2002-09-21
Analysis period: 2002-09-18 to 2002-09-22
2025-08-06 15:09:56 [INFO] src.monitoring.monitoring_utils: Loading IMERG raster data for 5 dates
2025-08-06 15:10:18 [INFO] src.monitoring.monitoring_utils: Successfully loaded IMERG data for 5 dates
On-the-fly calculation: 96.10 mm
Processing 2004223N11301 (ATCF: AL032004)...
Storm period: 2004-08-12 to 2004-08-13
Analysis period: 2004-08-11 to 2004-08-14
2025-08-06 15:10:18 [INFO] src.monitoring.monitoring_utils: Loading IMERG raster data for 4 dates
2025-08-06 15:10:29 [INFO] src.monitoring.monitoring_utils: Successfully loaded IMERG data for 4 dates
On-the-fly calculation: 44.66 mm
🔍 Comparing with pre-computed values...
On-the-fly results: 3 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
✅ Validation complete!
📋 Note: This is a basic validation. Full comparison would require:
- Exact SID matching between datasets
- Same date ranges and quantile calculations
- Processing all storms, not just a subset