Flood risk assessment in Himalayan urban landscape¶
Pithoragarh sub-watershed¶
This notebook attemots to bridge the doctoal reasearch work on biophydical modeling of flood regulation ecosystem services in rapidly urbanising Himalayan landscape. It adpots the quantitative flood risk framework central to global climate adapation science.
Phyiscal vulnerability is derived here for the coldpsots emeprically derived from Getis-Ord Gi* spatial statistics on the Soil Conservation Service-Curve Number layer -- producing statistically validated, locally-calibrated vulnerability evidence.
RISK = HAZARD × EXPOSURE × VULNERABILITY
| 📍Study area | Pithoragarh sub-watershed, Uttarakhand, India (EPSG:32644) |
| 🌧️ Hazard | SCS-CN runoff model · AMC III · Extreme rainfall event 232 mm/24hr (Aug 2011; over past two decades) · Threshold: P75 of runoff distribution |
| 👥 Exposure | Sentinel-2 LULC built-up class (10m) aggregated to 100m fraction × WorldPop 2020 (100m) |
| 📊 Vulnerability | Getis-Ord Gi* coldspots of FRC — statistically significant high-runoff clusters (p ≤ 0.01) from doctoral work |
Translating biophysically modelled flood regulation ecosystem service to disaster adaptation¶
To overcome generic depth-damage curves (e.g., Huizinga et al. 2017 for Asian residential buildings) that translate hydraulic inundation depth into a building damage fraction using lookup tables calibrated globally, this notebook takes into account local factors.
In this analysis, physical vulnerability is derived empirically using the Getis-Ord Gi* spatial statistic applied to the SCS-CN flood regulation capacity (FRC) raster : an output which undertakes local factors relative to slope, soil-type, land cover and rainfall events; produced during the doctoral research. Pixels identified as statistically significant FRC coldspots represent spatial clusters where runoff is very high relative to the neighbourhood : i.e., locations where the landscape structurally fails to buffer/permeate surface water. This is, by definition, the spatial signature of physical vulnerability.
This approach produces vulnerability evidence that is statistically validated, landscape-specific, and internally consistent with the hazard data : no external calibration dataset is required.
| § | Section | Key output |
|---|---|---|
| 1 | Setup | Libraries, paths, colour schemes |
| 2 | Data inputs | Raster metadata, CRS verification |
| 3 | Hazard | P75 flood threshold, sensitivity analysis |
| 4 | Exposure | Built-up fraction × WorldPop |
| 5 | Vulnerability | Gi* FRC coldspot map + LULC crosstab |
| 6 | Risk | Integrated risk map + summary statistics |
| 7 | Outputs | GeoTIFFs + summary CSV |
1. Setup¶
In [1]:
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors
import matplotlib.patches as mpatches
from matplotlib.colors import TwoSlopeNorm
import pandas as pd
from pathlib import Path
from scipy import stats
import warnings
warnings.filterwarnings('ignore')
import rasterio
from rasterio.enums import Resampling
from rasterio.warp import reproject, Resampling as WarpResampling
print(f'✅ All libraries ready')
✅ All libraries ready
2. Data inputs¶
In [2]:
# ── Input/Outpout file paths ───────────────────────────────────────────
RUNOFF_PATH = r"E:\Portfolio\Flood exposure\Pithoragarh\Run_off_2021_AMC_III.tif"
LULC_PATH = r"E:\Portfolio\Flood exposure\Pithoragarh\LULC_2021_PG.tif"
WORLDPOP_PATH = r"E:\Portfolio\Flood exposure\Pithoragarh\worldpop_pithoragarh_2020.tif"
HOTCOLDSPOT_PATH = r"E:\Portfolio\Flood exposure\Pithoragarh\FRC_2021_GiStar_class.tif"
OUTPUT_DIR = r"E:\Portfolio\Flood exposure\Pithoragarh\outputs"
Path(OUTPUT_DIR).mkdir(parents=True, exist_ok=True)
# ── Vulnerability weights per Gi* FRC class ───────────────────────────────
# Logic: FRC coldspot = high runoff cluster = landscape fails to buffer floods
# = structurally high physical vulnerability for built assets
# Weights are calibrated to the significance level of the spatial cluster.
# Applied ONLY within the flood-prone zone (P75 threshold).
VULN_WEIGHTS = {
'coldspot_99' : 0.80, # class −3 · p < 0.01 · Very high vulnerability
'coldspot_95' : 0.60, # class −2 · p < 0.05 · High vulnerability
'coldspot_90' : 0.40, # class −1 · p < 0.10 · Moderate vulnerability
'not_sig' : 0.10, # class 0 · Background vulnerability
'hotspot_90' : 0.05, # class +1 · High FRC = low vulnerability
'hotspot_95' : 0.03, # class +2
'hotspot_99' : 0.01, # class +3 · Very high FRC = minimal vulnerability
}
# ── Gi* FRC colour scheme (Blue = coldspot = high vulnerability) ──────────
GISTAR_COLORS = {
-3: '#053061', # Coldspot** 99%
-2: '#2166ac', # Coldspot* 95%
-1: '#92c5de', # Coldspot 90%
0: '#f7f7f7', # Not significant
1: '#f4a582', # Hotspot 90%
2: '#d6604d', # Hotspot* 95%
3: '#67001f', # Hotspot** 99%
}
GISTAR_LABELS = {
-3: 'Coldspot** (99%) — Very High Vulnerability',
-2: 'Coldspot* (95%) — High Vulnerability',
-1: 'Coldspot (90%) — Moderate Vulnerability',
0: 'Not Significant',
1: 'Hotspot (90%) — Low Vulnerability',
2: 'Hotspot* (95%) — Very Low Vulnerability',
3: 'Hotspot** (99%) — Minimal Vulnerability',
}
# ── LULC legend (class value → label, colour) ─────────────────────────
LULC_LEGEND = {
1: ('Oak Forest', '#2d6a4f'),
2: ('Pine Forest', '#52b788'),
3: ('Vegetation', '#74c69d'),
4: ('Waterbody', '#4895ef'),
5: ('Barren Land', '#adb5bd'),
6: ('Dry River Bed', '#00FFFF'),
7: ('Cropland', '#f9c74f'),
8: ('Open Area', '#e9ecef'),
9: ('Built-up', '#e63946'),
}
BUILTUP_CLASS = 9
Data checks!!¶
All inputs rasters must share CRS (EPSG:32644, UTM Zone 44N) and spatial resolution before any analysis.
In [3]:
def load_raster(path, name, is_lulc=False, is_gi=False):
with rasterio.open(path) as src:
data = src.read(1).astype(np.float32)
nodata = src.nodata
profile = src.profile
if nodata is not None:
data[data == nodata] = np.nan
print(f' {name}')
print(f' CRS: {src.crs} | Res: {src.res[0]:.0f}m × {src.res[1]:.0f}m | Shape: {src.height}×{src.width}')
print(f' Range: {np.nanmin(data):.2f} – {np.nanmax(data):.2f}')
if is_lulc or is_gi:
print(f' Classes: {np.unique(data[~np.isnan(data)]).astype(int)}')
print()
return data, profile
runoff_data, runoff_profile = load_raster(RUNOFF_PATH, 'Runoff — SCS-CN AMC III (mm)')
lulc_data, lulc_profile = load_raster(LULC_PATH, 'LULC 2021 — Sentinel 10m', is_lulc=True)
gistar_data, gistar_profile = load_raster(HOTCOLDSPOT_PATH, 'Gi* FRC Class — pre-computed 10m', is_gi=True)
Runoff — SCS-CN AMC III (mm) CRS: EPSG:32644 | Res: 10m × 10m | Shape: 986×1301 Range: 113.41 – 229.94 LULC 2021 — Sentinel 10m CRS: EPSG:32644 | Res: 10m × 10m | Shape: 995×1308 Range: 1.00 – 9.00 Classes: [1 2 3 4 5 6 7 8 9] Gi* FRC Class — pre-computed 10m CRS: EPSG:32644 | Res: 10m × 10m | Shape: 986×1301 Range: -3.00 – 3.00 Classes: [-3 -2 -1 0 1 2 3]
In [4]:
# ── Figure 1: Input data overview ───────────────────────────────────────────────
fig, axes = plt.subplots(1, 3, figsize=(20, 7))
fig.suptitle(f'Figure 1 — Input Data Overview',
fontsize=14, fontweight='bold', y=1.01)
# Panel A: LULC
cids = list(LULC_LEGEND.keys())
cmap_l = mcolors.ListedColormap([LULC_LEGEND[c][1] for c in cids])
norm_l = mcolors.BoundaryNorm([c-0.5 for c in cids]+[cids[-1]+0.5], cmap_l.N)
axes[0].imshow(lulc_data, cmap=cmap_l, norm=norm_l, interpolation='none')
axes[0].set_title('A) LULC 2021 — Sentinel-2 10m', fontweight='bold')
axes[0].axis('off')
axes[0].legend(
handles=[mpatches.Patch(color=LULC_LEGEND[c][1], label=LULC_LEGEND[c][0]) for c in cids],
loc='upper right', fontsize=7.5, framealpha=0.9)
# Panel B: Runoff
im0 = axes[1].imshow(runoff_data, cmap='Blues', interpolation='none')
axes[1].set_title(
'B) SCS-CN Runoff Depth (mm)\nAMC III · Extreme event: 232 mm/24hr',
fontweight='bold'
)
axes[1].axis('off')
# Panel C: Gi* FRC
gi_classes = sorted(GISTAR_COLORS.keys())
cmap_gi = mcolors.ListedColormap([GISTAR_COLORS[c] for c in gi_classes])
norm_gi = mcolors.BoundaryNorm([c-0.5 for c in gi_classes]+[gi_classes[-1]+0.5], cmap_gi.N)
axes[2].imshow(gistar_data, cmap=cmap_gi, norm=norm_gi, interpolation='none')
axes[2].set_title('C) Flood regulation capacity\nBlue = Coldspot (high vulnerability)',
fontweight='bold')
axes[2].axis('off')
axes[2].legend(
handles=[mpatches.Patch(color=GISTAR_COLORS[c], label=GISTAR_LABELS[c]) for c in gi_classes],
loc='upper right', fontsize=6.5, framealpha=0.9, title='Gi* FRC class')
cax = fig.add_axes([0.39, 0.10, 0.22, 0.025])
# [left, bottom, width, height]
cbar = fig.colorbar(im0, cax=cax, orientation='horizontal')
cbar.set_label('Runoff (mm)', fontsize=9)
plt.tight_layout(rect=[0, 0.08, 1, 1])
plt.savefig(f'{OUTPUT_DIR}/Fig1_inputs.png', dpi=300, bbox_inches='tight')
plt.show()
print('✅ Figure 1 saved')
✅ Figure 1 saved
3. Hazard¶
How the flood threshold matters?¶
The SCS-CN model was run under AMC III (Antecedent Moisture Condition III: saturated soils, worst-case monsoon conditions) for Pithoragarh's maximum observed 24-hour rainfall event (232mm, August 16, 2011).
Under such conditions, even forested mountain slopes generate substantial runoff, compressing the entire landscape into a narrow high-runoff range. An arbitrary threshold like 50mm/100mm would classify 100% of the watershed as flood-prone — physically meaningless for risk assessment.
Approach adopted: Pixels exceeding the 75th percentile (P75 = 205mm) of the runoff distribution are classified as flood-prone. This identifies the upper quartile where infiltration capacity is most severely overwhelmed — consistent with the 'Very Low' and 'Low' flood regulation classes in Pithoragarh subwatershed.
In [5]:
# Identify Threshold
runoff_flat = runoff_data[~np.isnan(runoff_data)].flatten()
pcts = {p: float(np.nanpercentile(runoff_flat, p)) for p in [25,50,75,90,95,99]}
FLOOD_THRESHOLD_MM = pcts[75]
print('Runoff distribution — Pithoragarh 2021')
for p, v in pcts.items():
tag = ' ◀ FLOOD THRESHOLD' if p == 75 else ''
print(f' P{p:2d} : {v:7.2f} mm{tag}')
print(f' Min : {runoff_flat.min():.2f} mm')
print(f' Max : {runoff_flat.max():.2f} mm')
print(f'\n Threshold: {FLOOD_THRESHOLD_MM:.1f} mm (P75)')
Runoff distribution — Pithoragarh 2021 P25 : 182.66 mm P50 : 195.20 mm P75 : 205.09 mm ◀ FLOOD THRESHOLD P90 : 210.34 mm P95 : 218.91 mm P99 : 223.48 mm Min : 113.41 mm Max : 229.94 mm Threshold: 205.1 mm (P75)
In [6]:
## Now since we have World population layer at 100x100 m resolution, we need to
## resample inputs raster at same resolution
# ── Build 100m common reference grid from LULC ────────────────────────────
scale = 10 # 100m
with rasterio.open(LULC_PATH) as src:
base_h = src.height
base_w = src.width
base_t = src.transform
base_crs = src.crs
new_h = base_h // scale
new_w = base_w // scale
new_t = base_t * base_t.scale(base_w / new_w, base_h / new_h)
profile_100m = lulc_profile.copy()
profile_100m.update(height=new_h, width=new_w, transform=new_t,
dtype=rasterio.float32, nodata=-9999)
# ── Resample runoff to 100m ─────────────────────────────────────
with rasterio.open(RUNOFF_PATH) as src:
runoff_100m = np.zeros((new_h, new_w), dtype=np.float32)
reproject(source=rasterio.band(src, 1), destination=runoff_100m,
src_transform=src.transform, src_crs=src.crs,
dst_transform=new_t, dst_crs=base_crs,
resampling=WarpResampling.average)
runoff_100m = np.nan_to_num(runoff_100m, nan=0.0)
flood_prone = (runoff_100m >= FLOOD_THRESHOLD_MM).astype(np.float32)
flood_area = float(flood_prone.sum() * 1e4 / 1e6) # 100m² → km²
flood_pct = float(flood_prone.mean() * 100)
print(f'Flood threshold : {FLOOD_THRESHOLD_MM:.1f} mm (P75)')
print(f'Flood-prone area : {flood_area:.2f} km² ({flood_pct:.1f}% of landscape)')
print('✅ Hazard layer is prepared')
Flood threshold : 205.1 mm (P75) Flood-prone area : 31.47 km² (24.5% of landscape) ✅ Hazard layer is prepared
4. Exposure¶
What/Who are at risk?¶
Beyond natural and built infrastructure, it is the population living in built-up areas.
Method¶
- Extract built-up pixels from the 10m Sentinel LULC (class 9)
- Aggregate to 100m via
averageresampling → built-up fraction per pixel (0–1) - Align WorldPop 2020 (100m) to the same grid
Exposure = built-up fraction × population
In [7]:
# ── Built-up binary → resample at 100m ───────────────────────────────────
builtup_bin = (lulc_data == BUILTUP_CLASS).astype(np.float32)
builtup_km2 = builtup_bin.sum() * 100 / 1e6
print(f'Built-up area: {builtup_km2:.2f} km² ')
# Write 10m built-up to temp file for resampling
tmp_bu = str(Path(OUTPUT_DIR) / 'tmp_builtup_10m.tif')
with rasterio.open(tmp_bu, 'w', driver='GTiff',
height=lulc_profile['height'], width=lulc_profile['width'],
count=1, dtype=rasterio.float32,
crs=lulc_profile['crs'], transform=lulc_profile['transform'],
nodata=-9999) as dst:
dst.write(builtup_bin, 1)
with rasterio.open(tmp_bu) as src:
builtup_frac = src.read(1, out_shape=(new_h, new_w),
resampling=Resampling.average).astype(np.float32)
builtup_frac = np.clip(builtup_frac, 0, 1)
# ── WorldPop data alignment ────────────────────────────────────────────────────
with rasterio.open(WORLDPOP_PATH) as wp:
population = np.zeros((new_h, new_w), dtype=np.float32)
reproject(source=rasterio.band(wp, 1), destination=population,
src_transform=wp.transform, src_crs=wp.crs,
dst_transform=new_t, dst_crs=base_crs,
resampling=WarpResampling.average)
# Remove Na values
population = np.where(np.isnan(population), 0, population)
population = np.where(population < 0, 0, population)
# ── Exposure ──────────────────────────────────────────────────────────────
exposure = builtup_frac * population
total_pop = float(population.sum())
total_exposed = float(exposure.sum())
print(f'\nTotal population : {total_pop:,.0f} persons')
print(f'In built-up areas : {total_exposed:,.0f} persons ({total_exposed/total_pop*100:.1f}%)')
print('✅ Exposure layer is ready')
Built-up area: 10.37 km² Total population : 68,054 persons In built-up areas : 14,051 persons (20.6%) ✅ Exposure layer is ready
In [8]:
# ── Figure 3: Exposure layer inputs & components ────────────────────────────────────────────────────
fig, axes = plt.subplots(1, 3, figsize=(17, 5.5))
fig.suptitle('Figure 3 — Exposure: Built-up Fraction × WorldPop Population',
fontsize=13, fontweight='bold')
im0 = axes[0].imshow(builtup_frac, cmap='Reds', vmin=0, vmax=1, interpolation='none')
axes[0].set_title('A) Built-up fraction (100m)\nAggregated from 10m Sentinel LULC',
fontweight='bold')
axes[0].axis('off')
plt.colorbar(im0, ax=axes[0], label='Fraction (0–1)', shrink=0.85)
im1 = axes[1].imshow(population, cmap='YlOrRd', interpolation='none')
axes[1].set_title('B) WorldPop population (100m)\nPersons per pixel · 2020', fontweight='bold')
axes[1].axis('off')
plt.colorbar(im1, ax=axes[1], label='Persons / pixel', shrink=0.85)
im2 = axes[2].imshow(np.where(exposure==0, np.nan, exposure), cmap='YlOrRd', interpolation='none')
axes[2].set_title(
f'C) EXPOSURE = A × B\nPopulation in built-up areas: {total_exposed:,.0f} persons',
fontweight='bold')
axes[2].axis('off')
plt.colorbar(im2, ax=axes[2], label='Exposed persons / pixel', shrink=0.85)
plt.tight_layout()
plt.savefig(f'{OUTPUT_DIR}/Fig3_exposure.png', dpi=300, bbox_inches='tight')
plt.show()
print('✅ Figure 3 saved')
✅ Figure 3 saved
5. Vulnerability¶
As stated earlier, here we utilise the following layer as vulnerability layer.
The Getis-Ord Gi* classification of Flood Regulation Capacity (FRC) was produced as part of the doctoral work using the SCS-CN runoff layer.
These FRC coldspots are precisely the zones of highest physical vulnerability: built assets located in areas where the landscape cannot moderate flood intensity are more severely damaged for a given rainfall event. The Gi* significance level (90%, 95%, 99%) is used to weight the vulnerability score — stronger statistical evidence of low FRC corresponds to higher vulnerability.
In [9]:
# ── Inspect Gi* FRC class distribution ───────────────────────────────────
gi_class = gistar_data.copy() # integer classes: −3 to +3
# Summary: total hot/coldspot area
coldspot_area = sum(np.nansum(gi_class == c) * 100 / 1e6 for c in [-3,-2,-1])
hotspot_area = sum(np.nansum(gi_class == c) * 100 / 1e6 for c in [1,2,3])
print(f'\n Total coldspot area : {coldspot_area:.2f} km²')
print(f' Total hotspot area : {hotspot_area:.2f} km²')
Total coldspot area : 38.39 km² Total hotspot area : 19.24 km²
In [10]:
# ── Figure 4: Gi* FRC map ─────────────────────────────────────────────────
gi_classes = sorted(GISTAR_COLORS.keys())
cmap_gi = mcolors.ListedColormap([GISTAR_COLORS[c] for c in gi_classes])
norm_gi = mcolors.BoundaryNorm([c-0.5 for c in gi_classes]+[gi_classes[-1]+0.5], cmap_gi.N)
fig, axes = plt.subplots(1, 2, figsize=(16, 10))
fig.suptitle(
f'Figure 4 — Flood Regulation Capacity (FRC) Gi* Classification\n'
f'Vulnerability Source',
fontsize=13, fontweight='bold'
)
# Panel A: Spatial map
axes[0].imshow(gi_class, cmap=cmap_gi, norm=norm_gi, interpolation='none')
axes[0].set_title(
'A) FRC Classification Map\n'
'Blue = Coldspot (low FRC = high vulnerability) | '
'Red = Hotspot (high FRC = low vulnerability)',
fontweight='bold')
axes[0].axis('off')
axes[0].legend(
handles=[mpatches.Patch(color=GISTAR_COLORS[c], label=GISTAR_LABELS[c]) for c in gi_classes],
loc='upper right', fontsize=7.5, framealpha=0.92, title='Gi* FRC class', title_fontsize=8)
# Panel B: Class area bar chart
class_labels = [GISTAR_LABELS[c].split('—')[0].strip() for c in gi_classes]
class_areas = [np.nansum(gi_class == c) * 100 / 1e6 for c in gi_classes]
bar_colors = [GISTAR_COLORS[c] for c in gi_classes]
bars = axes[1].barh(class_labels, class_areas, color=bar_colors, edgecolor='grey', linewidth=0.4)
axes[1].set_xlabel('Area (km²)')
axes[1].set_title('B) Area per Gi* FRC Class (km²)', fontweight='bold')
axes[1].invert_yaxis()
for bar, area in zip(bars, class_areas):
axes[1].text(bar.get_width() + 0.05, bar.get_y() + bar.get_height()/2,
f'{area:.1f} km²', va='center', fontsize=8)
axes[1].grid(True, axis='x', alpha=0.3)
axes[1].spines[['top','right']].set_visible(False)
plt.tight_layout()
plt.savefig(f'{OUTPUT_DIR}/Fig4_gistar_frc_map.png', dpi=300, bbox_inches='tight')
plt.show()
print('✅ Figure 4 saved')
✅ Figure 4 saved
In [11]:
# ── Convert Gi* FRC classes → vulnerability weights, resample to 100m ─────
vuln_10m = np.full_like(gi_class, np.nan, dtype=np.float32)
vuln_10m[gi_class == -3] = VULN_WEIGHTS['coldspot_99']
vuln_10m[gi_class == -2] = VULN_WEIGHTS['coldspot_95']
vuln_10m[gi_class == -1] = VULN_WEIGHTS['coldspot_90']
vuln_10m[gi_class == 0] = VULN_WEIGHTS['not_sig']
vuln_10m[gi_class == 1] = VULN_WEIGHTS['hotspot_90']
vuln_10m[gi_class == 2] = VULN_WEIGHTS['hotspot_95']
vuln_10m[gi_class == 3] = VULN_WEIGHTS['hotspot_99']
tmp_vu = str(Path(OUTPUT_DIR) / 'tmp_vuln_10m.tif')
with rasterio.open(tmp_vu, 'w', driver='GTiff',
height=gistar_profile['height'], width=gistar_profile['width'],
count=1, dtype=rasterio.float32,
crs=gistar_profile['crs'], transform=gistar_profile['transform'],
nodata=-9999) as dst:
dst.write(np.where(np.isnan(vuln_10m), -9999, vuln_10m), 1)
with rasterio.open(tmp_vu) as src:
vulnerability = src.read(1, out_shape=(new_h, new_w),
resampling=Resampling.average).astype(np.float32)
vulnerability = np.where(vulnerability == -9999, 0, vulnerability)
vulnerability = np.clip(vulnerability, 0, 1)
print('Vulnerability weights at 100m (Gi* FRC coldspot-derived):')
print(f' Full landscape mean : {vulnerability.mean():.3f}')
print(f' In flood zones only : {vulnerability[flood_prone==1].mean():.3f}')
print(f' Range : {vulnerability.min():.3f} – {vulnerability.max():.3f}')
print('\n✅ Vulnerability layer is ready')
Vulnerability weights at 100m (Gi* FRC coldspot-derived): Full landscape mean : 0.196 In flood zones only : 0.550 Range : 0.000 – 0.800 ✅ Vulnerability layer is ready
In [12]:
# ── Figure 5: Vulnerability weight map ───────────────────────────────────
fig, axes = plt.subplots(1, 2, figsize=(14, 6))
fig.suptitle(
f'Figure 5 — Physical Vulnerability: Gi* FRC Coldspot Weights \n'
'Coldspot confidence level → damage weight · applied within flood-prone zone only',
fontsize=12, fontweight='bold')
# ── Panel A: Full landscape vulnerability ─────────────────────────────────
im0 = axes[0].imshow(vulnerability, cmap='Blues_r', vmin=0, vmax=0.85,
interpolation='none')
axes[0].set_title(
'A) Vulnerability weights — full landscape\n'
'Derived from Gi* FRC coldspot significance level',
fontweight='bold')
axes[0].axis('off')
cb0 = plt.colorbar(im0, ax=axes[0], shrink=0.85, pad=0.02)
cb0.set_label('Damage weight (0 = resilient, 0.80 = very high)', fontsize=8)
# Add weight reference ticks on colorbar
for cls, w in VULN_WEIGHTS.items():
cb0.ax.axhline(w, color='black', lw=0.6, linestyle=':')
cb0.ax.text(1.15, w, f'{cls.replace("_"," ")} → {w}',
fontsize=6.5, va='center', color='#333')
mean_vuln = float(vulnerability[flood_prone == 1].mean())
# ── Panel B: Vulnerability within flood-prone zone only ───────────────────
vuln_flood = np.where(flood_prone == 1, vulnerability, np.nan)
im1 = axes[1].imshow(vuln_flood, cmap='YlOrRd', vmin=0, vmax=0.85,
interpolation='none')
# Overlay flood zone boundary
axes[1].contour(flood_prone, levels=[0.5], colors='#1d3557', linewidths=1.2)
axes[1].set_title(
f'B) Vulnerability within flood-prone zone (P75)\n'
f'Blue outline = flood boundary · Mean weight: {mean_vuln:.2f}',
fontweight='bold')
axes[1].axis('off')
cb1 = plt.colorbar(im1, ax=axes[1], shrink=0.85, pad=0.02)
cb1.set_label('Damage weight (0–1)', fontsize=8)
# Weight legend as text box
weight_lines = '\n'.join([
f' Coldspot** (99%) → {VULN_WEIGHTS["coldspot_99"]}',
f' Coldspot* (95%) → {VULN_WEIGHTS["coldspot_95"]}',
f' Coldspot (90%) → {VULN_WEIGHTS["coldspot_90"]}',
f' Not significant → {VULN_WEIGHTS["not_sig"]}',
f' Hotspot → {VULN_WEIGHTS["hotspot_90"]}–{VULN_WEIGHTS["hotspot_99"]}',
])
plt.tight_layout()
plt.savefig(f'{OUTPUT_DIR}/Fig5_vulnerability_weights.png', dpi=300, bbox_inches='tight')
plt.show()
print('✅ Figure 5 saved')
✅ Figure 5 saved
6. Flood Risk Integration¶
$$\text{Risk}_{i} = \text{Flood-prone}_{i} \times \text{Exposure}_{i} \times \text{Vulnerability}^{\text{Gi*}}_{i}$$
Risk is computed only within the flood-prone zone (P75 threshold). The output is in person-damage equivalents : the population in built-up areas within statistically significant low-FRC clusters, weighted by vulnerability confidence.
In [13]:
risk = flood_prone * exposure * vulnerability
exposed_flood = float(exposure[flood_prone == 1].sum())
total_risk = float(risk.sum())
mean_vuln = float(vulnerability[flood_prone == 1].mean())
high_risk_area = float((risk > np.nanpercentile(risk[risk>0], 75)).sum() * 1e4 / 1e6)
print('=' * 62)
print(f' FLOOD RISK SUMMARY')
print('=' * 62)
print(f' Hazard model : SCS-CN AMC III · 232 mm/24hr')
print(f' Flood threshold : {FLOOD_THRESHOLD_MM:.0f} mm (P75)')
print(f' Flood-prone area : {flood_area:.2f} km² ({flood_pct:.1f}% of landscape)')
print(f' Vulnerability source : Gi* FRC coldspots · doctoral thesis pipeline')
print(f' ─────────────────────────────────────────────────────────────')
print(f' Total population : {total_pop:,.0f} persons (WorldPop 2020)')
print(f' Exposed in flood zone: {exposed_flood:,.0f} persons')
print(f' Mean vulnerability : {mean_vuln:.3f} (Gi* FRC weighted)')
print(f' Total risk score : {total_risk:,.1f} person-damage units')
print(f' High-risk area (Q75) : {high_risk_area:.2f} km²')
print('=' * 62)
============================================================== FLOOD RISK SUMMARY ============================================================== Hazard model : SCS-CN AMC III · 232 mm/24hr Flood threshold : 205 mm (P75) Flood-prone area : 31.47 km² (24.5% of landscape) Vulnerability source : Gi* FRC coldspots · doctoral thesis pipeline ───────────────────────────────────────────────────────────── Total population : 68,054 persons (WorldPop 2020) Exposed in flood zone: 2,291 persons Mean vulnerability : 0.550 (Gi* FRC weighted) Total risk score : 1,470.4 person-damage units High-risk area (Q75) : 3.51 km² ==============================================================
In [14]:
# ── Figure 6: Main framework — 4-panel output ────────────────────────────
fig, axes = plt.subplots(2, 2, figsize=(14, 12))
fig.suptitle(
f'Figure 6 — Flood Risk Assessment: {"Pithoragarh, 2021"}\n'
'Risk = Hazard × Exposure × Vulnerability (Gi* FRC coldspots)',
fontsize=14, fontweight='bold', y=1.02
)
# A: Hazard
im0 = axes[0,0].imshow(runoff_100m, cmap='Blues', interpolation='none')
axes[0,0].contour(flood_prone, levels=[0.5], colors='#e63946', linewidths=1.8)
axes[0,0].set_title(
f'A) HAZARD — SCS-CN Runoff Depth (mm)\n'
f'Red outline = flood-prone zone (P75 · {flood_area:.1f} km² · {flood_pct:.0f}%)',
fontweight='bold')
axes[0,0].axis('off')
plt.colorbar(im0, ax=axes[0,0], label='Runoff depth (mm)', shrink=0.85)
# B: Exposure
im1 = axes[0,1].imshow(np.where(exposure==0, np.nan, exposure), cmap='YlOrRd', interpolation='none')
axes[0,1].set_title(
f'B) EXPOSURE — Population in Built-up Areas\n'
f'Total: {total_pop:,.0f} · Urban footprint: {total_exposed:,.0f} persons',
fontweight='bold')
axes[0,1].axis('off')
plt.colorbar(im1, ax=axes[0,1], label='Persons per 100m pixel', shrink=0.85)
# C: Vulnerability
im2 = axes[1,0].imshow(
np.where(flood_prone==1, vulnerability, np.nan),
cmap='YlOrRd', vmin=0, vmax=0.85, interpolation='none')
axes[1,0].set_title(
f'C) VULNERABILITY — Gi* FRC Coldspot Weights\n'
f'p ≤ 0.01 clusters of low flood regulation · Mean: {mean_vuln:.2f}',
fontweight='bold')
axes[1,0].axis('off')
plt.colorbar(im2, ax=axes[1,0], label='Vulnerability weight (0–1)', shrink=0.85)
# D: Risk
risk_plot = np.where(risk == 0, np.nan, risk)
vmax_r = np.nanpercentile(risk[risk > 0], 95)
im3 = axes[1,1].imshow(risk_plot, cmap='YlOrRd', vmin=0, vmax=vmax_r, interpolation='none')
axes[1,1].set_title(
f'D) FLOOD RISK = A × B × C\n'
f'Exposed: {exposed_flood:,.0f} persons · Score: {total_risk:,.0f}',
fontweight='bold')
axes[1,1].axis('off')
plt.colorbar(im3, ax=axes[1,1], label='Risk (person-damage units)', shrink=0.85)
fig.text(0.01, -0.01,
'Vulnerability: Gi* FRC coldspots (Sharma, doctoral thesis Ch. 5) · '
'Exposure: Sentinel-2 LULC (10m) × WorldPop 2020 (100m) · '
'Hazard: SCS-CN AMC III, 232 mm/24hr · Flood zone: P75 runoff threshold',
fontsize=7.5, color='grey', style='italic')
plt.tight_layout()
plt.savefig(f'{OUTPUT_DIR}/Fig6_flood_risk_main.png', dpi=300, bbox_inches='tight')
plt.show()
print('✅ Figure 6 (main output) saved')
✅ Figure 6 (main output) saved
In [15]:
def save_tif_100m(array, filename):
"""Write a numpy array to a GeoTIFF at 100m resolution."""
p = profile_100m.copy()
p.update(dtype=rasterio.float32, count=1, nodata=-9999)
arr = np.where(np.isnan(array.astype(np.float32)), -9999, array.astype(np.float32))
with rasterio.open(f'{OUTPUT_DIR}/{filename}', 'w', **p) as dst:
dst.write(arr, 1)
print(f' ✅ {filename}')
print('Saving GeoTIFFs (100m) ...')
save_tif_100m(builtup_frac, 'exposure_builtup_fraction_100m.tif')
save_tif_100m(population, 'exposure_population_100m.tif')
save_tif_100m(exposure, 'exposure_combined_100m.tif')
save_tif_100m(flood_prone, 'hazard_flood_prone_P75_100m.tif')
save_tif_100m(vulnerability, 'vulnerability_gistar_frc_100m.tif')
save_tif_100m(risk, 'FLOOD_RISK_gistar_frc_2021_100m.tif')
# ── Summary CSV ───────────────────────────────────────────────────────────
summary = pd.DataFrame([{
'site' : "Pithoragarh",
'year' : "2021",
'hazard_model' : 'SCS-CN AMC-III · 232 mm/24hr',
'flood_threshold_mm' : round(FLOOD_THRESHOLD_MM, 1),
'flood_threshold_method' : 'P75 runoff percentile',
'exposure_lulc' : f'Sentinel-2 10m class {BUILTUP_CLASS}',
'exposure_population' : 'WorldPop 100m 2020',
'vulnerability_method' : 'Gi* FRC coldspots',
'flood_area_km2' : round(flood_area, 2),
'flood_pct_landscape' : round(flood_pct, 1),
'total_population' : int(round(total_pop, 0)),
'exposed_in_flood_zones' : int(round(exposed_flood, 0)),
'exposed_pct_total' : round(exposed_flood / total_pop * 100, 1),
'mean_vulnerability_gi' : round(mean_vuln, 3),
'total_risk_score' : round(total_risk, 1),
'high_risk_area_km2_q75' : round(high_risk_area, 2),
}])
csv_path = f'{OUTPUT_DIR}/flood_risk_summary.csv'
summary.to_csv(csv_path, index=False)
print(f'\n ✅ flood_risk_summary.csv')
print()
print(summary.T.to_string(header=False))
Saving GeoTIFFs (100m) ... ✅ exposure_builtup_fraction_100m.tif ✅ exposure_population_100m.tif ✅ exposure_combined_100m.tif ✅ hazard_flood_prone_P75_100m.tif ✅ vulnerability_gistar_frc_100m.tif ✅ FLOOD_RISK_gistar_frc_2021_100m.tif ✅ flood_risk_summary.csv site Pithoragarh year 2021 hazard_model SCS-CN AMC-III · 232 mm/24hr flood_threshold_mm 205.1 flood_threshold_method P75 runoff percentile exposure_lulc Sentinel-2 10m class 9 exposure_population WorldPop 100m 2020 vulnerability_method Gi* FRC coldspots flood_area_km2 31.47 flood_pct_landscape 24.5 total_population 68054 exposed_in_flood_zones 2291 exposed_pct_total 3.4 mean_vulnerability_gi 0.55 total_risk_score 1470.4 high_risk_area_km2_q75 3.51
In [ ]:
!jupyter nbconvert --to html Himalayan_flood_risk_assessment.ipynb --output-dir "E:\Portfolio\notebooks"