---
jupyter: ds-flood-gfm
execute:
eval: false
---
## Hurricane Melissa.
Here we will build real time snap shots of GFM detected flooding related to Hurricane Melissa for Jamaica, Haiti, and Cuba
```{python}
import pystac_client
import stackstac
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap, BoundaryNorm
import numpy as np
import geopandas as gpd
from fsspec.implementations.http import HTTPFileSystem
from ds_flood_gfm.geo_utils import load_adm0_lowres
import time
import pandas as pd
import ocha_stratus as stratus
from datetime import datetime, timedelta
from exactextract import exact_extract
# Configure matplotlib for better display
%matplotlib inline
plt.rcParams['figure.dpi'] = 100
```
```{python}
# Analysis configuration
notebook_config = {
'combined_range': {
'start': '2025-10-01',
'end': '2025-11-30'
},
'pre_storm': {
'name': 'Pre-Storm Conditions',
'start': '2025-10-01',
'end': '2025-10-14',
'snapshot_dates': ['2024-10-18', '2024-10-25'],
'reference_date': '2024-10-18'
},
'post_storm': {
'name': 'Post-Storm Conditoins',
'start': '2025-10-18',
'end': '2025-11-30',
'snapshot_dates': ['2025-11-25', '2024-11-27'],
'reference_date': '2024-11-18'
}
}
```
## Area of Interest
Analysis focuses on Admin 1 region CM004 in Cameroon, with three custom AOI polygons defined by Copernicus Emergency Management Service (EMSR772 & EMSR779).
```{python}
GLOBAL_ADM1 = (
"https://data.fieldmaps.io/edge-matched/humanitarian/intl/adm1_polygons.parquet"
)
ISO3 = "JAM"
filesystem = HTTPFileSystem()
filters = [("iso_3", "=", ISO3)]
gdf = gpd.read_parquet(GLOBAL_ADM1, filesystem=filesystem, filters=filters)
gdf_aoi = gdf
bbox = gdf_aoi.total_bounds
```
```{python}
#| eval: false
gdf_world_adm0 = load_adm0_lowres()
gdf_subset_admo = gdf_world_adm0[gdf_world_adm0.name == "Cameroon"]
```
## STAC Query
Query GFM flood extent data from EODC STAC API for the combined date range (Oct 14 - Nov 30, 2024).
```{python}
# | cache: true
# | cache-comments: "Cache STAC search results - change dates to invalidate"
stac_api = "https://stac.eodc.eu/api/v1"
client = pystac_client.Client.open(stac_api)
# Query for combined date range covering both October and November events
start_date = notebook_config['combined_range']['start']
end_date = notebook_config['combined_range']['end']
# Step 2: Search for GFM items
datetime_range = f"{start_date}/{end_date}"
search = client.search(
collections=["GFM"],
bbox=gdf_aoi.total_bounds,
datetime=datetime_range,
)
item_collection = search.item_collection()
```
```{python}
#| cache: true
#| cache-lazy: true
stack = stackstac.stack(item_collection, epsg=4326)
```
## Pre-Processing
### GFM Data (STAC) Lazy
Extract flood extent band and create daily composites by taking the maximum flood value per day across overlapping tiles.
```{python}
# | cache: true
# | cache-lazy: true
stack_flood = stack.sel(band="ensemble_flood_extent")
stack_flood_clipped = stack_flood.sel(
x=slice(bbox[0], bbox[2]), y=slice(bbox[3], bbox[1]) # y is reversealri
)
# Group by day and take the maximum flood value for each day
# Use groupby instead of resample to only get days that actually have data
stack_flood_max = stack_flood_clipped.groupby("time.date").max()
# Rename the dimension back to 'time' and convert to datetime
stack_flood_max = stack_flood_max.rename({"date": "time"})
stack_flood_max["time"] = stack_flood_max.time.astype("datetime64[ns]")
```
Split the combined dataset into October and November events before applying forward-fill. This is more efficient since forward-fill is a sequential operation that processes all timesteps.
```{python}
# Split into individual events using config dates (before ffill for performance)
ic_pre = stack_flood_max.sel(
time=slice(notebook_config['pre_storm']['start'], notebook_config['pre_storm']['end'])
)
ic_post = stack_flood_max.sel(
time=slice(notebook_config['post_storm']['start'], notebook_config['post_storm']['end'])
)
# Apply forward-fill separately to each event (much more efficient)
ic_pre_latest = ic_pre.ffill(dim="time")
ic_post_latest = ic_post.ffill(dim="time")
snapshot_dates_post = ic_post_latest.time.values[-2:]
```
# tmp example
```{python}
ic_post_oct27 = ic_post.sel(time= "2025-10-27").compute()
```
```{python}
# Extract tile footprints directly from the stackstac object for October 27
from shapely.geometry import shape, box
# Get the specific time slice from ic_post that corresponds to Oct 27
ic_post_oct27_stack = ic_post.sel(time="2025-10-27")
# Extract unique tile footprints from the stackstac coordinates
# stackstac stores the original item IDs and geometries in the data array attributes
tiles_from_stack = []
# Get the items that were actually used in this time slice
if hasattr(ic_post_oct27_stack, 'item') or 'id' in ic_post_oct27_stack.coords:
# stackstac stores item IDs - we need to match back to item_collection
item_ids_in_stack = ic_post_oct27_stack.coords.get('id', ic_post_oct27_stack.coords.get('item', None))
if item_ids_in_stack is not None:
for item_id in item_ids_in_stack.values:
# Find matching item in collection
matching_items = [item for item in item_collection if item.id == item_id]
if matching_items:
item = matching_items[0]
tile_id = item.properties.get("Equi7Tile", item.id)
tiles_from_stack.append({
'geometry': shape(item.geometry),
'tile_id': tile_id,
'date': str(pd.Timestamp(item.datetime).normalize())[:10],
'item_id': item.id
})
# If we couldn't extract from stackstac coords, fall back to filtering item_collection
if len(tiles_from_stack) == 0:
print("Extracting footprints from item_collection (stackstac coords not available)")
target_date = pd.Timestamp("2025-10-27").normalize()
for item in item_collection:
item_date = pd.Timestamp(item.datetime).normalize().tz_localize(None)
if item_date == target_date:
tile_id = item.properties.get("Equi7Tile", item.id)
tiles_from_stack.append({
'geometry': shape(item.geometry),
'tile_id': tile_id,
'date': str(item_date)[:10],
'item_id': item.id
})
gdf_tiles = gpd.GeoDataFrame(tiles_from_stack, crs='EPSG:4326')
# Deduplicate by unique tile_id (spatial tiles, not temporal duplicates)
gdf_tiles_unique = gdf_tiles.drop_duplicates(subset='tile_id', keep='first').copy()
print(f"Found {len(gdf_tiles)} total items, {len(gdf_tiles_unique)} unique spatial tiles for 2025-10-27")
gdf_tiles_unique
```
```{python}
# Create a map with Jamaica boundary, flood extent, and tile footprints
fig, ax = plt.subplots(1, 1, figsize=(14, 10))
# Define colors for flood data visualization including NaN
colors = ["lightgrey", "darkblue", "red"] # 0=no flood, 1=flood, NaN=nodata
cmap = ListedColormap(colors)
bounds = [0, 0.5, 1.5, 255.5]
norm = BoundaryNorm(bounds, cmap.N)
# Set the color for NaN values (bad data)
cmap.set_bad(color="red", alpha=0.5)
# Plot the flood data - NaN values will now show as red
im = ic_post_oct27.plot(ax=ax, cmap=cmap, norm=norm, add_colorbar=False)
# Add Jamaica boundary
gdf_aoi.boundary.plot(ax=ax, color="black", linewidth=3, alpha=0.8)
# Add tile footprints as polygons
if len(gdf_tiles_unique) > 0:
gdf_tiles_unique.boundary.plot(
ax=ax, color="purple", linewidth=2, linestyle="--", alpha=0.8
)
# Add labels with tile ID and acquisition date inside each polygon
for idx, row in gdf_tiles_unique.iterrows():
bounds = row.geometry.bounds
# Position label at center of tile
label_x = (bounds[0] + bounds[2]) / 2
label_y = (bounds[1] + bounds[3]) / 2
# Create label text with tile ID and date
label_text = f"{row['tile_id']}\n{row['date']}"
ax.text(
label_x,
label_y,
label_text,
fontsize=9,
fontweight="bold",
ha="center",
va="center",
bbox=dict(
boxstyle="round,pad=0.4",
facecolor="white",
edgecolor="purple",
alpha=0.85,
linewidth=1.5,
),
)
# Add colorbar
cbar = plt.colorbar(im, ax=ax, label="Flood Status")
cbar.set_ticks([0, 1])
cbar.set_ticklabels(["No Flood", "Flood"])
# Zoom out to show full data extent with buffer
data_bounds = ic_post_oct27.rio.bounds()
x_buffer = (data_bounds[2] - data_bounds[0]) * 0.15
y_buffer = (data_bounds[3] - data_bounds[1]) * 0.15
ax.set_xlim(data_bounds[0] - x_buffer, data_bounds[2] + x_buffer)
ax.set_ylim(data_bounds[1] - y_buffer, data_bounds[3] + y_buffer)
# Set title and labels
ax.set_title(
"Flood Extent - October 27, 2025\nJamaica\n(Red = No Data/NaN, Purple boxes = Tile footprints)",
fontsize=14,
fontweight="bold",
)
ax.set_xlabel("Longitude")
ax.set_ylabel("Latitude")
ax.grid(True, alpha=0.3, linestyle=":")
plt.tight_layout()
plt.show()
```
```{python}
# | eval: false
# | echo: false
def print_mb(da):
print(f"Size in memory: {da.nbytes / 1024**2:.1f} MB")
```
### Population Data (GHSL)
Load GHSL population data to overlay. Here we are using the estimated 2025 pop produced in 2023.
- We load it just for the AOI bbox, then clip it down further.
```{python}
blob_name = "ghsl/pop/GHS_POP_E2025_GLOBE_R2023A_4326_3ss_V1_0.tif"
da_global = stratus.open_blob_cog(blob_name, container_name="raster").squeeze(drop=True)
# Ensure the global data has the correct CRS
# da_global = da_global.rio.write_crs("EPSG:4326")
da_global.rio.crs
# clip to box (need to do this first, otherwise Python crashes on normal .rio.clip)
min_x, min_y, max_x, max_y = gdf_aoi.total_bounds
da_clip_box = da_global.rio.clip_box(minx=min_x, miny=min_y, maxx=max_x, maxy=max_y)
# clip to admin and preserve CRS
da_clip = da_clip_box.rio.clip(gdf_aoi.geometry)
da_clip = da_clip.rio.write_crs("EPSG:4326")
da_clip = da_clip.where(da_clip != da_clip.rio.nodata)
```
```{python}
#| eval: false
#| echo: false
# Load local GHSL population data
import rioxarray as rxr
# Load the local GeoTIFF file
da_local = rxr.open_rasterio(
"data/ghsl_e2020/GHS_POP_E2020_GLOBE_R2023A_4326_3ss_V1_0_R8_C20.tif"
).squeeze(drop=True)
# Ensure the local data has the correct CRS
da_local = da_local.rio.write_crs("EPSG:4326")
# clip to box (need to do this first, otherwise Python crashes on normal .rio.clip)
min_x, min_y, max_x, max_y = gdf_aoi.total_bounds
da_local_clip_box = da_local.rio.clip_box(
minx=min_x, miny=min_y, maxx=max_x, maxy=max_y
)
# clip to admin and preserve CRS
da_local_clip = da_local_clip_box.rio.clip(gdf_aoi.geometry)
da_local_clip = da_local_clip.rio.write_crs("EPSG:4326")
# Create a map to visualize overlap between local population data and AOI
fig, ax = plt.subplots(1, 1, figsize=(12, 8))
# Check coordinate ranges
print(f"da_local_clip bounds: {da_local_clip.rio.bounds()}")
print(f"gdf_aoi bounds: {gdf_aoi.total_bounds}")
print(f"da_local_clip CRS: {da_local_clip.rio.crs}")
print(f"gdf_aoi CRS: {gdf_aoi.crs}")
# Plot the local population data with manual extent setting
im = da_local_clip.plot(ax=ax, cmap="YlOrRd", vmin=0, vmax=200, add_colorbar=False)
# Add colorbar
cbar = plt.colorbar(im, ax=ax, label="Population Count")
# Add AOI boundary
gdf_aoi.boundary.plot(ax=ax, color="red", linewidth=3, alpha=1.0)
# Set axis limits to the AOI bounds to ensure visibility
ax.set_xlim(gdf_aoi.total_bounds[0], gdf_aoi.total_bounds[2])
ax.set_ylim(gdf_aoi.total_bounds[1], gdf_aoi.total_bounds[3])
# Set title and labels
ax.set_title(
"Local Population Data 2020 - Overlap with AOI", fontsize=12, fontweight="bold"
)
ax.set_xlabel("Longitude")
ax.set_ylabel("Latitude")
plt.tight_layout()
plt.show()
```
```{python}
# | eval: false
# | echo: false
print(f"da_clip fill value: {da_clip.rio.nodata}")
print(f"da_clip has NaNs: {da_clip.isnull().any().values}")
print(f"da_clip min (should be 0): {da_clip.min().values}")
da_clip = da_local_clip
# Or more robust - use the nodata value
# Verify
print(
f"da_clip min after masking: {da_clip.min().values}"
) # Should be 0 or small positive
print(f"da_clip has NaNs: {da_clip.isnull().any().values}") # Should be True
```
Here we show the GHSL population raster for the admin 1 of interest (grey outline) and the specific custom AOI's from CEMS (black outline).
```{python}
# Create a map of the clipped population data
fig, ax = plt.subplots(1, 1, figsize=(12, 8))
# Plot the population data
im = da_clip.plot(
ax=ax,
cmap='YlOrRd',
vmin=0,
vmax=200,
add_colorbar=False
)
# Add colorbar
cbar = plt.colorbar(im, ax=ax, label='Population Count')
# Add admin boundary in grey
gdf_aoi.boundary.plot(ax=ax, color='grey', linewidth=2, alpha=0.8)
# Add Copernicus AOI boundaries in black
gdf_aoi_copernicus.boundary.plot(ax=ax, color='black', linewidth=2, alpha=1.0)
# Set title and labels
ax.set_title('Population Density 2025 - Clipped to AOI', fontsize=12, fontweight='bold')
ax.set_xlabel('Longitude')
ax.set_ylabel('Latitude')
plt.tight_layout()
plt.show()
```
```{python}
#| eval: false
#| echo: false
print(f"da_clip dtype: {da_clip.dtype}")
print(f"da_clip nodata: {da_clip.rio.nodata}")
print(f"da_clip min: {da_clip.min().values}")
print(f"da_clip max: {da_clip.max().values}")
print(f"da_clip sum: {da_clip.sum(skipna=True).values:,.0f}")
print(f"da_clip has NaN: {da_clip.isnull().any().values}")
print(f"da_clip has negative: {(da_clip < 0).any().values}")
```
## Overlay Flood Extent & Population Data
Here we follow the GFM methodology as closely as possible. We need to get the data on comparable grids:
- We take the 100m population data and resample it to the 20m flood data and rescale the population numbers to reflect this.
- We then run zonal statistics (using `exact_extract`)
```{python}
def gfm_exact_workflow_zonal(flood_20m, population_100m, gdf_zones, zone_id_col=None):
"""
Calculate affected population using GFM methodology with zonal statistics.
Parameters
----------
flood_20m : xr.DataArray
Binary flood mask at 20m resolution (0 or 1)
population_100m : xr.DataArray
Population density at 100m resolution
gdf_zones : gpd.GeoDataFrame
Polygons for zonal statistics (e.g., multiple admin areas)
zone_id_col : str, optional
Column name to use as zone identifier
Returns
-------
pd.DataFrame
Results with affected_population per zone
"""
print("=" * 80)
print("GFM Exact Workflow with Zonal Statistics (exactextract)")
print("=" * 80)
print(f"\nNumber of zones: {len(gdf_zones)}")
# Ensure CRS match
if gdf_zones.crs != flood_20m.rio.crs:
print(f"Reprojecting zones from {gdf_zones.crs} to {flood_20m.rio.crs}")
gdf_zones = gdf_zones.to_crs(flood_20m.rio.crs)
# Clean population nodata
if population_100m.rio.nodata is not None:
print(f"Clearing nodata metadata (was: {population_100m.rio.nodata})")
population_100m = population_100m.rio.write_nodata(None)
# Resample population to 20m
print("Upsampling population from 100m to 20m...")
population_20m = population_100m.rio.reproject_match(
flood_20m,
resampling=0 # nearest neighbor
)
# Adjust for pixel area (divide by 25)
print("Adjusting population density for resolution change...")
pixel_area_ratio = (100 / 20) ** 2 # = 25
population_20m_adjusted = population_20m / pixel_area_ratio
print(f" Population adjustment factor: {pixel_area_ratio}")
print(f" Original sum: {population_100m.sum(skipna=True).values:,.0f}")
print(f" After adjustment: {population_20m_adjusted.sum(skipna=True).values:,.0f}")
# Calculate affected population
print("Calculating affected population at 20m resolution...")
affected_pop_20m = flood_20m * population_20m_adjusted
# Zonal statistics using exactextract
print(f"Computing zonal statistics for {len(gdf_zones)} zones...")
# Prepare zone identifier
if zone_id_col and zone_id_col in gdf_zones.columns:
zone_ids = gdf_zones[zone_id_col].tolist()
else:
zone_ids = gdf_zones.index.tolist()
# Extract zonal sums for affected population
results = exact_extract(
affected_pop_20m,
gdf_zones,
'sum',
output='pandas'
)
# Extract flood extent statistics
flood_results = exact_extract(
flood_20m,
gdf_zones,
['sum', 'count'],
output='pandas'
)
# Combine results
results_df = pd.DataFrame({
'zone_id': zone_ids,
'affected_population': results['sum'].values,
'flooded_pixels': flood_results['sum'].values,
'total_pixels': flood_results['count'].values,
})
# Calculate flooded area
pixel_area_m2 = 20 * 20 # 20m resolution
results_df['flooded_area_m2'] = results_df['flooded_pixels'] * pixel_area_m2
results_df['flooded_area_km2'] = results_df['flooded_area_m2'] / 1_000_000
# Add zone names if available
if zone_id_col and zone_id_col in gdf_zones.columns:
results_df['zone_name'] = gdf_zones[zone_id_col].values
# Clean up values
results_df['affected_population'] = results_df['affected_population'].fillna(0).astype(int)
results_df['flooded_pixels'] = results_df['flooded_pixels'].fillna(0).astype(int)
print("\n" + "=" * 80)
print("RESULTS BY ZONE")
print("=" * 80)
print(results_df.to_string(index=False))
return results_df
```
```{python}
def map_flood_extent_snapshots(flood_data_computed, snapshot_dates, gdf_boundaries, event_name):
"""
Create side-by-side flood extent maps for snapshot dates.
Parameters
----------
flood_data_computed : xr.DataArray
Computed flood data array (already filtered and computed)
snapshot_dates : list
List of snapshot dates being visualized
gdf_boundaries : gpd.GeoDataFrame
Boundaries to overlay on maps
event_name : str
Event name for plot context
Returns
-------
None
Displays matplotlib figure
"""
# Define colors for flood data visualization
colors = ["lightgrey", "darkblue", "white"] # 0=no flood, 1=flood, 255=nodata
cmap = ListedColormap(colors)
bounds = [0, 0.5, 1.5, 255.5]
norm = BoundaryNorm(bounds, cmap.N)
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(20, 8))
# Process each time step
for i, ax in enumerate([ax1, ax2]):
time_step = flood_data_computed.isel(time=i)
date_str = str(time_step.time.values)[:10]
# Process data (handle NaN and downsample for visualization)
data_clean = np.where(np.isnan(time_step), 255, time_step)
data_clean = data_clean[::4, ::4] # Downsample for speed
# Get extent from the xarray data
bounds_extent = time_step.rio.bounds()
extent = [bounds_extent[0], bounds_extent[2], bounds_extent[1], bounds_extent[3]]
# Plot the flood data
im = ax.imshow(
data_clean,
cmap=cmap,
norm=norm,
interpolation="nearest",
origin="upper",
extent=extent,
)
# Overlay boundaries
gdf_boundaries.boundary.plot(ax=ax, color='red', linewidth=3, alpha=0.8)
# Set title and labels
ax.set_title(f'{event_name}\nFlood Extent - {date_str}',
fontsize=12, fontweight='bold')
ax.set_xlabel('Longitude')
ax.set_ylabel('Latitude')
plt.tight_layout()
plt.show()
```
```{python}
def create_comparison_barplot(gfm_results, other_systems_data, event_name,
systems=['CEMS', 'FloodScan', 'GFM (STAC)']):
"""
Create comparison bar plot of affected population across systems.
Parameters
----------
gfm_results : pd.DataFrame
GFM zonal statistics results with columns: zone_id, affected_population
other_systems_data : dict
Dict of system data, e.g., {'FloodScan': {'AOI 1': 7580, ...}, 'CEMS': {...}}
event_name : str
Event name for plot title
systems : list
List of systems to include in plot (e.g., ['CEMS', 'GFM (STAC)'])
Returns
-------
None
Displays matplotlib figure and prints comparison stats
"""
# Convert other systems data to DataFrame
comparison_data = []
for system, values in other_systems_data.items():
for aoi, pop in values.items():
comparison_data.append({
'AOI': aoi,
'System': system,
'affected_population': pop
})
df_other_systems = pd.DataFrame(comparison_data)
# Prepare GFM data
gfm_comparison = gfm_results[['zone_id', 'affected_population']].copy()
gfm_comparison['System'] = 'GFM (STAC)'
# Map zone_id to standard AOI format
zone_id_to_aoi = {}
for idx, zone_id in enumerate(gfm_comparison['zone_id']):
if isinstance(zone_id, (int, np.integer)):
zone_id_to_aoi[zone_id] = f"AOI {zone_id + 1}"
elif 'AOI' in str(zone_id).upper():
import re
match = re.search(r'AOI(\d+)', str(zone_id), re.IGNORECASE)
if match:
aoi_num = int(match.group(1))
zone_id_to_aoi[zone_id] = f"AOI {aoi_num}"
else:
zone_id_to_aoi[zone_id] = f"AOI {idx + 1}"
else:
zone_id_to_aoi[zone_id] = f"AOI {idx + 1}"
gfm_comparison['AOI'] = gfm_comparison['zone_id'].map(zone_id_to_aoi)
gfm_comparison = gfm_comparison[['AOI', 'System', 'affected_population']]
# Combine all data
df_combined = pd.concat([df_other_systems, gfm_comparison], ignore_index=True)
# Filter to requested systems
df_combined = df_combined[df_combined['System'].isin(systems)]
# Display comparison table
print("\n" + "=" * 80)
print(f"AFFECTED POPULATION COMPARISON - {event_name}")
print("=" * 80)
print(df_combined.pivot(index='AOI', columns='System', values='affected_population'))
# Create comparison barplot
fig, ax = plt.subplots(figsize=(12, 6))
# Get unique AOIs and systems
aois = sorted(df_combined['AOI'].unique())
systems_plot = [s for s in systems if s in df_combined['System'].values]
# Set up bar positions
x = np.arange(len(aois))
width = 0.8 / len(systems_plot) # Adjust width based on number of systems
# Define consistent colors
color_map = {
'GFM (STAC)': '#3cb371', # Minty green
'FloodScan': '#ff6347', # Tomato red
'CEMS': '#0f52ba' # Sapphire blue
}
# Plot bars for each system
for i, system in enumerate(systems_plot):
system_data = df_combined[df_combined['System'] == system]
values = [system_data[system_data['AOI'] == aoi]['affected_population'].values[0]
if len(system_data[system_data['AOI'] == aoi]) > 0 else 0
for aoi in aois]
bars = ax.bar(x + i * width, values, width, label=system, color=color_map.get(system, 'gray'))
# Add value labels on bars
for bar, val in zip(bars, values):
height = bar.get_height()
ax.text(bar.get_x() + bar.get_width()/2., height,
f'{int(val):,}',
ha='center', va='bottom', fontsize=9)
# Customize plot
ax.set_xlabel('Area of Interest (AOI)', fontsize=12, fontweight='bold')
ax.set_ylabel('Affected Population', fontsize=12, fontweight='bold')
ax.set_title(f'Affected Population Comparison\n{event_name}', fontsize=14, fontweight='bold')
ax.set_xticks(x + width * (len(systems_plot) - 1) / 2)
ax.set_xticklabels(aois)
ax.legend(title='System', loc='upper right', fontsize=11)
ax.grid(axis='y', alpha=0.3, linestyle='--')
# Format y-axis
ax.yaxis.set_major_formatter(plt.FuncFormatter(lambda x, p: f'{int(x):,}'))
plt.tight_layout()
plt.show()
# Print comparison ratios
if 'GFM (STAC)' in systems_plot:
print("\n" + "=" * 80)
print("SYSTEM COMPARISON RATIOS (relative to GFM)")
print("=" * 80)
for aoi in aois:
gfm_val = df_combined[(df_combined['AOI'] == aoi) & (df_combined['System'] == 'GFM (STAC)')]['affected_population'].values
if len(gfm_val) > 0:
gfm_val = gfm_val[0]
print(f"\n{aoi}:")
print(f" GFM (baseline): {gfm_val:,}")
for system in systems_plot:
if system != 'GFM (STAC)':
sys_val = df_combined[(df_combined['AOI'] == aoi) & (df_combined['System'] == system)]['affected_population'].values
if len(sys_val) > 0:
sys_val = sys_val[0]
ratio = sys_val / gfm_val if gfm_val > 0 else np.inf
print(f" {system}: {sys_val:,} ({ratio:.2f}x GFM)")
```
### Flood Extent Visualizations
```{python}
# pull arrays into memory for easier building iteration\
ic_post_snapshot = ic_post_latest.sel(time=snapshot_dates_post)
ic_post_snapshot_computed = ic_post_snapshot.compute()
```
```{python}
#| cache: true
# Visualize flood extent for snapshot dates (with forward-fill compositing)
map_flood_extent_snapshots(
flood_data_computed=ic_post_snapshot_computed,
snapshot_dates=snapshot_dates_post,
gdf_boundaries=gdf_aoi,
event_name=notebook_config['post_storm']['name']
)
```
```{python}
# claude here only
# Extract acquisition footprints showing ACTUAL acquisition dates (not forward-filled dates)
# This shows users which parts of the map are more/less up-to-date
from shapely.geometry import shape
# Convert snapshot dates to datetime for comparison (timezone-naive)
snapshot_dates_dt = [pd.Timestamp(str(date)[:10]) for date in snapshot_dates_post]
print(f"Snapshot dates: {[str(d)[:10] for d in snapshot_dates_dt]}")
# For each snapshot date, find the most recent actual acquisition per tile
footprints_by_snapshot = {}
for snapshot_date in snapshot_dates_dt:
snapshot_str = str(snapshot_date)[:10]
print(f"\nProcessing snapshot: {snapshot_str}")
# Get all items acquired on or before this snapshot date
# Convert item datetime to date-only (timezone-naive) for comparison
items_before = []
for item in item_collection:
item_date = pd.Timestamp(item.datetime).normalize().tz_localize(None)
if item_date <= snapshot_date:
items_before.append(item)
print(f" Found {len(items_before)} items acquired on or before {snapshot_str}")
# Group by tile_id and keep only the most recent acquisition per tile
tile_latest = {}
for item in items_before:
tile_id = item.properties.get("Equi7Tile", item.id)
item_dt = pd.Timestamp(item.datetime).normalize().tz_localize(None)
if tile_id not in tile_latest or item_dt > tile_latest[tile_id]['datetime']:
tile_latest[tile_id] = {
'geometry': shape(item.geometry),
'datetime': item_dt,
'date': str(item_dt)[:10],
'tile_id': tile_id,
'id': item.id,
'snapshot_date': snapshot_str
}
footprints_by_snapshot[snapshot_str] = list(tile_latest.values())
print(f"\n{snapshot_str} snapshot uses {len(tile_latest)} tiles:")
for tile_info in tile_latest.values():
print(f" Tile {tile_info['tile_id']}: acquired {tile_info['date']}")
# Create visualization showing footprints with ACTUAL acquisition dates
fig, axes = plt.subplots(1, 2, figsize=(20, 8))
for i, snapshot_date_dt in enumerate(snapshot_dates_dt):
ax = axes[i]
snapshot_str = str(snapshot_date_dt)[:10]
# Plot Jamaica boundary first
gdf_aoi.boundary.plot(ax=ax, color='black', linewidth=2, alpha=0.8)
gdf_aoi.plot(ax=ax, facecolor='lightgrey', alpha=0.2)
# Get footprints for this snapshot
footprints = footprints_by_snapshot[snapshot_str]
gdf_footprints = gpd.GeoDataFrame(footprints, crs='EPSG:4326')
# Plot footprints
if len(gdf_footprints) > 0:
gdf_footprints.boundary.plot(ax=ax, color='darkblue', linewidth=1.5, alpha=0.7)
gdf_footprints.plot(ax=ax, facecolor='blue', alpha=0.15)
# Add ACTUAL acquisition date labels inside each footprint
for idx, row in gdf_footprints.iterrows():
bounds = row.geometry.bounds
# Position label at top-center of footprint
label_x = (bounds[0] + bounds[2]) / 2
label_y = bounds[3] - (bounds[3] - bounds[1]) * 0.1
# Show actual acquisition date (not snapshot date)
ax.text(label_x, label_y, row['date'],
fontsize=10, fontweight='bold',
ha='center', va='top',
bbox=dict(boxstyle='round,pad=0.3', facecolor='white',
edgecolor='darkblue', alpha=0.8))
# Set title
ax.set_title(f'{notebook_config["post_storm"]["name"]}\nSnapshot: {snapshot_str}\n(Labels show actual acquisition dates, {len(gdf_footprints)} tiles)',
fontsize=12, fontweight='bold')
ax.set_xlabel('Longitude')
ax.set_ylabel('Latitude')
ax.grid(True, alpha=0.3, linestyle='--')
plt.tight_layout()
plt.show()
```
```{python}
# Convert flooded pixels to points for better visualization
from shapely.geometry import Point
# Function to convert flood raster to points
def flood_pixels_to_points(flood_da, time_idx=0):
"""
Convert flooded pixels (value=1) to point geometries.
Parameters
----------
flood_da : xr.DataArray
Flood extent data array
time_idx : int
Time index to extract (default: 0 for first snapshot)
Returns
-------
gpd.GeoDataFrame
Points representing flooded pixels
"""
# Select the time step
flood_time = flood_da.isel(time=time_idx)
# Get flooded pixels only (value == 1)
flooded = flood_time.where(flood_time == 1, drop=True)
# Extract coordinates
coords_list = []
for y_val in flooded.y.values:
for x_val in flooded.x.values:
pixel_val = flooded.sel(x=x_val, y=y_val).values
if pixel_val == 1:
coords_list.append({'x': x_val, 'y': y_val})
# Create GeoDataFrame
if len(coords_list) > 0:
df = pd.DataFrame(coords_list)
geometry = [Point(xy) for xy in zip(df['x'], df['y'])]
gdf_points = gpd.GeoDataFrame(df, geometry=geometry, crs='EPSG:4326')
return gdf_points
else:
return gpd.GeoDataFrame(columns=['x', 'y', 'geometry'], crs='EPSG:4326')
# Convert flood pixels to points for both snapshot dates
gdf_flood_points_0 = flood_pixels_to_points(ic_post_snapshot_computed, time_idx=0)
gdf_flood_points_1 = flood_pixels_to_points(ic_post_snapshot_computed, time_idx=1)
print(f"Snapshot 0: {len(gdf_flood_points_0)} flooded points")
print(f"Snapshot 1: {len(gdf_flood_points_1)} flooded points")
# Create side-by-side visualization
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(20, 8))
for i, (ax, gdf_points) in enumerate([(ax1, gdf_flood_points_0), (ax2, gdf_flood_points_1)]):
# Plot Jamaica boundary
gdf_aoi.boundary.plot(ax=ax, color='grey', linewidth=2, alpha=0.8)
gdf_aoi.plot(ax=ax, facecolor='lightgrey', alpha=0.2)
# Plot flooded points
if len(gdf_points) > 0:
gdf_points.plot(ax=ax, color='darkblue', markersize=5, alpha=0.6)
# Get date for title
date_str = str(ic_post_snapshot_computed.isel(time=i).time.values)[:10]
# Set title and labels
ax.set_title(f'{notebook_config["post_storm"]["name"]}\nFlooded Pixels (as Points) - {date_str}',
fontsize=12, fontweight='bold')
ax.set_xlabel('Longitude')
ax.set_ylabel('Latitude')
# Add grid
ax.grid(True, alpha=0.3, linestyle='--')
plt.tight_layout()
plt.show()
```
## Resulting Flood Extents & Affected Population Zonal Stats
Calculate affected population for the November event using the forward-filled flood stack.
```{python}
#| cache: true
# Run zonal stats for November event using filled/forward-filled stack
flood_day_max = nov_stack_flood_max.sel(time=notebook_config['november']['snapshot_dates'])
flood_day_max_computed= flood_day_max.compute()
# Process each time step in flood_data_latest_computed
list_zonal = []
for i in range(len(flood_day_max_computed.time)):
time_step = flood_day_max_computed.isel(time=i)
date_str = str(time_step.time.values)[:10]
print(f"\nProcessing {date_str}...")
# Apply GFM threshold
flood_day_max_binary = (time_step == 1).astype(int)
flood_day_max_binary = flood_day_max_binary.rio.write_crs(time_step.rio.crs)
# Run zonal analysis
results_time = gfm_exact_workflow_zonal(
flood_20m=flood_day_max_binary,
population_100m=da_clip,
gdf_zones=gdf_aoi_copernicus,
zone_id_col='aoi_code' # Adjust to match your column name
)
# Add date column to results
results_time['date'] = date_str
list_zonal.append(results_time)
# Display results for this time step
print(f"Total affected population for {date_str}: {results_time['affected_population'].sum():,}")
# Combine all results
day_max_zonal = pd.concat(list_zonal, ignore_index=True)
print(f"\nOverall results shape: {day_max_zonal.shape}")
day_max_zonal
```
```{python}
#| cache: true
# Zonal statistics using forward-filled composites (November event)
results_filled = []
for i in range(len(flood_data_latest_computed.time)):
time_step = flood_data_latest_computed.isel(time=i)
date_str = str(time_step.time.values)[:10]
print(f"\nProcessing {date_str}...")
# Apply GFM threshold
flood_binary = (time_step == 1).astype(int)
flood_binary = flood_binary.rio.write_crs(time_step.rio.crs)
# Run zonal analysis
results_time = gfm_exact_workflow_zonal(
flood_20m=flood_binary,
population_100m=da_clip,
gdf_zones=gdf_aoi_copernicus,
zone_id_col='aoi_code' # Adjust to match your column name
)
# Add date column to results
results_time['date'] = date_str
results_filled.append(results_time)
# Display results for this time step
print(f"Total affected population for {date_str}: {results_time['affected_population'].sum():,}")
# Combine all results
nov_zonal_results = pd.concat(results_filled, ignore_index=True)
# Store the latest date result for comparison
nov_latest_img_zonal_date = nov_zonal_results[nov_zonal_results["date"] == notebook_config['november']['snapshot_dates'][0]]
```
### System Comparisons
```{python}
# Hardcoded November 2024 reference data from other systems
nov_other_systems_data = {
'FloodScan': {
'AOI 1': 7580,
'AOI 2': 110655,
'AOI 3': 86461
},
'CEMS': {
'AOI 1': 1500,
'AOI 2': 13000,
'AOI 3': 8200
}
}
# Create comparison plot with all three systems
create_comparison_barplot(
gfm_results=nov_latest_img_zonal_date,
other_systems_data=nov_other_systems_data,
event_name=notebook_config['november']['name'],
systems=['CEMS', 'FloodScan', 'GFM (STAC)']
)
```
## October 2024 Event Analysis
### Zonal Statistics
```{python}
#| cache: true
# Run zonal stats for October event using filled/forward-filled stack
oct_flood_data_latest = oct_ic_latest.sel(time=notebook_config['october']['snapshot_dates'])
oct_flood_data_latest_computed = oct_flood_data_latest.compute()
oct_results_filled = []
for i in range(len(oct_flood_data_latest_computed.time)):
time_step = oct_flood_data_latest_computed.isel(time=i)
date_str = str(time_step.time.values)[:10]
print(f"\nProcessing {date_str}...")
# Apply GFM threshold
flood_binary = (time_step == 1).astype(int)
flood_binary = flood_binary.rio.write_crs(time_step.rio.crs)
# Run zonal analysis
results_time = gfm_exact_workflow_zonal(
flood_20m=flood_binary,
population_100m=da_clip,
gdf_zones=gdf_aoi_copernicus,
zone_id_col='aoi_code'
)
# Add date column to results
results_time['date'] = date_str
oct_results_filled.append(results_time)
# Display results for this time step
print(f"Total affected population for {date_str}: {results_time['affected_population'].sum():,}")
# Combine all results
oct_zonal_results = pd.concat(oct_results_filled, ignore_index=True)
# Store the latest date result for comparison
oct_latest_img_zonal_date = oct_zonal_results[oct_zonal_results["date"] == notebook_config['october']['snapshot_dates'][0]]
```
### Flood Extent Visualizations
```{python}
#| cache: true
# Visualize flood extent for October snapshot dates
map_flood_extent_snapshots(
flood_data_computed=oct_flood_data_latest_computed,
snapshot_dates=notebook_config['october']['snapshot_dates'],
gdf_boundaries=gdf_aoi_copernicus,
event_name=notebook_config['october']['name']
)
```
### System Comparisons
FloodScan data not available for October event.
```{python}
# Hardcoded October 2024 CEMS reference data
# TODO: Replace x, y, z with actual CEMS values for October
oct_other_systems_data = {
'CEMS': {
'AOI 1': 8300,
'AOI 2': 6400, # TODO: Add October CEMS value for AOI 2
'AOI 3': 26000 # TODO: Add October CEMS value for AOI 3
}
}
# Create comparison plot with CEMS and GFM only (no FloodScan)
create_comparison_barplot(
gfm_results=oct_latest_img_zonal_date,
other_systems_data=oct_other_systems_data,
event_name=notebook_config['october']['name'],
systems=['CEMS', 'GFM (STAC)'] # Only CEMS and GFM
)
```
## Appendix
```{python}
#| eval: false
#| echo: false
# All dates in ic_latest are actual observation dates (no filling needed)
# Create binary flood masks and calculate population impacts for each time step
ic_ready = stack_flood_max
population_impacts = []
dates = []
for i, time_step in enumerate(ic_ready.time):
print(f"\nProcessing date: {str(time_step.values)[:10]}")
try:
# Get flood data for this time step
flood_data = ic_ready.sel(time=time_step)
# Compute the data first to handle network errors
print(" Computing flood data from EODC...")
flood_data_computed = flood_data.compute()
# Create binary flood mask (1 = flooded, 0 = not flooded)
# Use GFM's official binary classification: exactly 1 = flood
# (ensemble requires 2 out of 3 algorithms to agree)
flood_binary = (flood_data_computed == 1).astype(int)
flood_binary = flood_binary.rio.write_crs(flood_data_computed.rio.crs)
# DEBUG: Check flood mask BEFORE resampling
print(
f" DEBUG - Unique values in flood_binary: {np.unique(flood_binary.values[~np.isnan(flood_binary.values)])}"
)
print(
f" DEBUG - Flooded cells BEFORE resample: {(flood_binary == 1).sum().values:,}"
)
# Resample flood data to match population grid resolution
# Population data is at higher resolution, so we need to align grids
# Use nearest neighbor (resampling=0) per GFM methodology
# https://rasterio.readthedocs.io/en/stable/api/rasterio.enums.html#rasterio.enums.Resampling
print(" Reprojecting to match population grid...")
flood_resampled = flood_binary.rio.reproject_match(
da_clip, resampling=0
) # 0 = Resampling.nearest (GFM spec)
except Exception as e:
print(f" ⚠️ ERROR: Failed to process {str(time_step.values)[:10]}: {e}")
print(f" Skipping this date...")
continue
# Calculate flooded population (multiply binary flood by population)
print(" Calculating affected population...")
# DEBUG: Check spatial extents and masking
print(f" DEBUG - flood_resampled extent: {flood_resampled.rio.bounds()}")
print(f" DEBUG - da_clip extent: {da_clip.rio.bounds()}")
print(f" DEBUG - da_clip has NaNs: {da_clip.isnull().any().values}")
print(f" DEBUG - da_clip non-null cells: {(~da_clip.isnull()).sum().values:,}")
print(f" DEBUG - flooded cells (>0): {(flood_resampled > 0).sum().values:,}")
print(f" DEBUG - flood_resampled sum: {flood_resampled.sum().values:,.0f}")
print(
f" DEBUG - flood_resampled unique values: {np.unique(flood_resampled.values[~np.isnan(flood_resampled.values)])}"
)
pop_flooded = flood_resampled * da_clip
# DEBUG: Check population calculation
print(f" DEBUG - pop_flooded > 0 cells: {(pop_flooded > 0).sum().values:,}")
print(
f" DEBUG - pop_flooded sum (skipna=True): {pop_flooded.sum(skipna=True).values:,.0f}"
)
print(f" DEBUG - pop_flooded max: {pop_flooded.max(skipna=True).values:.2f}")
# Store results
population_impacts.append(pop_flooded)
date_str = str(time_step.values)[:10]
dates.append(date_str)
print(f" ✅ {date_str}: {pop_flooded.sum().values:,.0f} people affected")
# Calculate zonal statistics (total flooded population per admin area)
zonal_stats = []
for i, pop_impact in enumerate(population_impacts):
# Sum population within the admin boundary
total_pop_flooded = pop_impact.sum().values
zonal_stats.append(
{
"date": dates[i],
"total_population_flooded": total_pop_flooded,
"admin_code": gdf_aoi.iloc[0]["adm1_src"],
}
)
# Display results
import pandas as pd
df_results = pd.DataFrame(zonal_stats)
print("\n=== Zonal Statistics Summary ===")
print(df_results)
```
```{python}
#| eval: false
#| echo: false
def gfm_exact_workflow(flood_20m, population_100m, aoi_geometry=None):
"""
Replicate GFM's exact workflow for affected population calculation.
Parameters
----------
flood_20m : xr.DataArray
Binary flood mask at 20m resolution (0 or 1)
population_100m : xr.DataArray
Population density at 100m resolution
aoi_geometry : gpd.GeoSeries, optional
Geometry to clip final results
Returns
-------
tuple
(affected_pop_raster_20m, total_affected_count)
"""
print("GFM Exact Workflow:")
print("1. Flood data at 20m resolution")
print("2. Upsampling population from 100m to 20m...")
# Step 1: Resample population TO match flood grid (100m → 20m)
# Use nearest neighbor per GFM spec
population_20m = population_100m.rio.reproject_match(
flood_20m,
resampling=0 # nearest neighbor
)
# Step 2: Adjust population density for pixel area change
# When going from 100m to 20m, each 100m cell becomes 25 cells of 20m
# Population density needs to be divided by 25
print("3. Adjusting population density for resolution change...")
pixel_area_ratio = (100 / 20) ** 2 # = 25
population_20m_adjusted = population_20m / pixel_area_ratio
print(f" Population adjustment factor: {pixel_area_ratio}")
print(f" Original sum: {population_100m.sum(skipna=True).values:,.0f}")
print(f" After resampling: {population_20m.sum(skipna=True).values:,.0f}")
print(f" After adjustment: {population_20m_adjusted.sum(skipna=True).values:,.0f}")
# Step 3: Multiply flood mask (20m) × population (20m)
print("4. Calculating affected population at 20m resolution...")
affected_pop_20m = flood_20m * population_20m_adjusted
# Step 4: Clip to AOI if provided
if aoi_geometry is not None:
print("5. Clipping to AOI boundary...")
affected_pop_20m = affected_pop_20m.rio.clip(
aoi_geometry.values,
affected_pop_20m.rio.crs,
drop=False
)
# Step 5: Calculate total
total_affected = float(affected_pop_20m.sum(skipna=True).values)
print(f"\n✅ Total affected population: {int(total_affected):,}")
return affected_pop_20m, total_affected
```
```{python}
#| eval: false
#| echo: false
# Process the latest time step as an example
time_step = ic_ready.time[1] # Get the last time step
print(f"\nProcessing date: {str(time_step.values)[:10]}")
da_clip_nd = da_clip.rio.write_nodata(None)
try:
# Get flood data for this time step
flood_data = ic_ready.sel(time=time_step)
# Compute the data first to handle network errors
print(" Computing flood data from EODC...")
flood_data_computed = flood_data.compute()
# Create binary flood mask (1 = flooded, 0 = not flooded)
# Use GFM's official binary classification: exactly 1 = flood
# (ensemble requires 2 out of 3 algorithms to agree)
flood_binary = (flood_data_computed == 1).astype(int)
flood_binary = flood_binary.rio.write_crs(flood_data_computed.rio.crs)
# Run the GFM exact workflow
affected_pop_raster, total_affected = gfm_exact_workflow(
flood_20m=flood_binary,
population_100m=da_clip_nd,
aoi_geometry=gdf_aoi.geometry,
)
print(f"Analysis complete for {str(time_step.values)[:10]}")
except Exception as e:
print(f"Error processing data: {e}")
```
### Google Earth Engine GHSL stuff
```{python}
#| eval: false
#| echo: false
import ee
import geemap
# Initialize Earth Engine
ee.Initialize()
# | cache: true
# | cache-lazy: true
# Download GHSL 2025 population data from Earth Engine
# NOTE: EE's computePixels API has a 50MB limit. This AOI at 100m = 67MB (too large)
# Using 200m resolution = 17MB (works fine). For 100m you'd need to use Task API + Drive export.
pop_2025_ee = ee.Image('JRC/GHSL/P2023A/GHS_POP/2025').select('population_count')
# Use AOI bounds for region
region = ee.Geometry.Rectangle([bbox[0], bbox[1], bbox[2], bbox[3]])
# Download as numpy array at 200m (still good resolution for population analysis)
pop_array = geemap.ee_to_numpy(
pop_2025_ee,
region=region,
scale=200 # 200m to stay within 50MB API limit
)
# Squeeze out the band dimension
pop_2025 = np.squeeze(pop_array)
print(f"Population data shape: {pop_2025.shape}")
print(f"Resolution: 200m")
print(f"Min: {pop_2025.min():.2f}, Max: {pop_2025.max():.2f}")
print(f"Non-zero pixels: {np.count_nonzero(pop_2025 > 0):,}")
```
```{python}
#| eval: false
#| echo: false
#
# Create a static visualization combining population and flood data
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(20, 8))
# Left panel: Population data
im1 = ax1.imshow(
pop_2025,
cmap='YlOrRd',
interpolation='nearest',
origin='upper',
extent=[bbox[0], bbox[2], bbox[1], bbox[3]],
vmin=0,
vmax=200
)
gdf_aoi.boundary.plot(ax=ax1, color='black', linewidth=2, alpha=1.0)
ax1.set_title('Population Density 2025\n(people per grid cell)', fontsize=12, fontweight='bold')
ax1.set_xlabel('Longitude')
ax1.set_ylabel('Latitude')
plt.colorbar(im1, ax=ax1, label='Population Count')
# Right panel: Latest flood extent
data_flood = np.where(np.isnan(img_tlater), 255, img_tlater)
data_flood = data_flood[::8, ::8]
extent_flood = [
img_tlater.x.min().item(),
img_tlater.x.max().item(),
img_tlater.y.min().item(),
img_tlater.y.max().item(),
]
im2 = ax2.imshow(
data_flood,
cmap=cmap,
norm=norm,
interpolation='nearest',
origin='upper',
extent=extent_flood,
)
gdf_aoi.boundary.plot(ax=ax2, color='black', linewidth=2, alpha=1.0)
date_str = str(img_tlater.time.values)[:10]
ax2.set_title(f'Latest Flood Extent\n{date_str}', fontsize=12, fontweight='bold')
ax2.set_xlabel('Longitude')
ax2.set_ylabel('Latitude')
plt.tight_layout()
plt.show()
```