How to combine hazard, exposure and vulnerability into a risk score

Problem statement

South of Valencia, the European Flood Awareness System publishes three depth grids — the 1-in-10, 1-in-100 and 1-in-500 year floods — on a 100 m grid that does not quite line up between files. Turning them into one number per place that can be compared means putting all three on one grid, converting each depth into a damaged share of the built-up land with a depth–damage curve, and integrating the damage over annual probability. The result is expected annual damage, 36.47 hectare-equivalents a year here, three quarters of it from floods more frequent than one in ten years.

This guide computes expected annual damage over a window south of Valencia: three flood depth grids, built-up area as exposure, and a depth–damage curve. It also shows how much of the answer comes from the rarest scenario, which is the part hardest to estimate.

Quick answer

Damage per scenario, then integrate over annual probability:

import numpy as np, rasterio
from rasterio.enums import Resampling
from rasterio.vrt import WarpedVRT
from rasterio.windows import from_bounds
from rasterio.warp import transform_bounds

BBOX = [-0.46, 39.25, -0.33, 39.45]

# One grid for every scenario: the three files do not share a pixel grid exactly.
with rasterio.open("flood/floodmap_EFAS_RP100_C.tif") as src:
    win = from_bounds(*transform_bounds(4326, src.crs, *BBOX), src.transform)
    GRID = {"crs": src.crs, "transform": src.window_transform(win),
            "width": int(win.width), "height": int(win.height)}

def depth(rp):
    with rasterio.open(f"flood/floodmap_EFAS_RP{rp}_C.tif") as src, WarpedVRT(
            src, crs=GRID["crs"], transform=GRID["transform"], width=GRID["width"],
            height=GRID["height"], resampling=Resampling.bilinear) as vrt:
        d = vrt.read(1, masked=True).filled(np.nan)
    return np.where(np.isfinite(d) & (d > 0), d, 0.0)

def damage_fraction(d):
    """Illustrative depth-damage curve for buildings: 0 at 0 m, 0.5 at 2 m, 0.8 at 4 m.
    Published curves (JRC, FEMA HAZUS, national guidance) differ by country and
    building type; substitute the right one for your study area."""
    return np.clip(np.interp(d, [0, 0.5, 1, 2, 4, 6], [0, 0.20, 0.35, 0.50, 0.80, 0.90]), 0, 1)

scenarios = {10: depth("010"), 100: depth("100"), 500: depth("500")}
for rp, d in scenarios.items():
    wet = d > 0
    print(f"1-in-{rp:<3}: {wet.sum() * 0.01:5.1f} km² wet, mean damage fraction where wet "
          f"{damage_fraction(d[wet]).mean():.2f}")

1-in-10 :  61.7 km² wet, mean damage fraction where wet 0.23
1-in-100:  77.4 km² wet, mean damage fraction where wet 0.32
1-in-500:  84.3 km² wet, mean damage fraction where wet 0.35

The damage fraction turns depth into a proportion of value lost. It is the vulnerability term, and it is the one that has to come from somewhere other than the GIS.

Scenarios with probabilities, exposure, vulnerability, integration

1. Assemble hazard scenarios with their probabilities

Each map has an annual exceedance probability: 1/10, 1/100, 1/500. Those are the weights in the integration.

2. Put exposure on the same grid

Built-up area, building footprints, population or value. Here it is built-up land from a 10 m land-cover map, resampled to the 100 m flood grid as a fraction.

3. Apply the vulnerability curve

Damage fraction from depth, then multiply by exposure and by value per unit if you have it.

4. Integrate over probability

Expected annual damage is the area under the loss–probability curve: sum over scenario pairs of the average loss times the difference in exceedance probability.

5. Report the contribution of each scenario

If most of the expected damage comes from the rarest scenario, the answer rests on the least reliable map.

6. Keep the assumptions with the number

Curve, exposure layer, scenarios and grid. Expected annual damage in "units of built-up hectare-equivalents" is honest; a currency figure without a value per hectare is not.

Expected annual flood damage, and its sensitivity

The data used on this page

The European Flood Awareness System's pan-European flood hazard maps from the JRC (data catalogue): water depth in metres on a 100 m grid for the 1-in-10, 1-in-100, 1-in-500-year floods, one zipped GeoTIFF per return period, saved under flood/.

import os
import zipfile

import requests

EFAS = "https://jeodpp.jrc.ec.europa.eu/ftp/jrc-opendata/FLOODS/EuropeanMaps/floodMap_RP{}.zip"
os.makedirs("flood", exist_ok=True)
for rp in ("010", "100", "500"):
    tif = f"floodmap_EFAS_RP{rp}_C.tif"
    if not os.path.exists(f"flood/{tif}"):                      # about 180 MB per return period
        with requests.get(EFAS.format(rp), stream=True, timeout=600) as r:
            r.raise_for_status()
            with open("efas.zip", "wb") as f:
                for chunk in r.iter_content(1 << 20):
                    f.write(chunk)
        zipfile.ZipFile("efas.zip").extract(tif, "flood")
        os.remove("efas.zip")
print(sorted(os.listdir("flood")))

['floodmap_EFAS_RP010_C.tif', 'floodmap_EFAS_RP100_C.tif', 'floodmap_EFAS_RP500_C.tif']

Example 1 — exposure on the hazard grid

import pystac_client, planetary_computer
from rasterio.enums import Resampling
from rasterio.vrt import WarpedVRT

shape = (GRID["height"], GRID["width"])
cat = pystac_client.Client.open("https://planetarycomputer.microsoft.com/api/stac/v1",
                                modifier=planetary_computer.sign_inplace)
wc_item = next(cat.search(collections=["esa-worldcover"], bbox=BBOX).items())
with rasterio.open(wc_item.assets["map"].href) as src, WarpedVRT(
        src, crs=GRID["crs"], transform=GRID["transform"], width=GRID["width"],
        height=GRID["height"], resampling=Resampling.nearest) as vrt:
    wc = vrt.read(1)
built = (wc == 50).astype("float32")                    # 1 where the 100 m cell is built-up

print(f"grid {shape}, built-up cells {int(built.sum()):,} = {built.sum():,.0f} ha")
for rp, d in scenarios.items():
    print(f"  1-in-{rp:<3}: built-up hectares wet {float(((d > 0) * built).sum()):6,.0f}")

grid (234, 138), built-up cells 5,643 = 5,643 ha
  1-in-10 : built-up hectares wet    186
  1-in-100: built-up hectares wet    302
  1-in-500: built-up hectares wet    328

Example 2 — expected annual damage

probs = {10: 1 / 10, 100: 1 / 100, 500: 1 / 500}
loss = {rp: float((damage_fraction(d) * built).sum()) for rp, d in scenarios.items()}
for rp in sorted(loss):
    print(f"1-in-{rp:<3}: damage {loss[rp]:7,.1f} built-up hectare-equivalents "
          f"(annual probability {probs[rp]:.3f})")

pairs = sorted(probs.items(), key=lambda kv: -kv[1])     # most frequent first
ead = 0.0
prev_p, prev_loss = 1.0, 0.0                             # below the 1-in-10, assume no damage
contrib = {}
for rp, p in pairs:
    band = (prev_p - p) * (prev_loss + loss[rp]) / 2
    contrib[rp] = band
    ead += band
    prev_p, prev_loss = p, loss[rp]
ead += probs[500] * loss[500]                            # tail beyond the rarest scenario
print(f"expected annual damage: {ead:,.2f} hectare-equivalents per year")
for rp, c in contrib.items():
    print(f"  band down to 1-in-{rp:<3}: {c:6,.2f} ({c / ead:5.1%})")

1-in-10 : damage    60.9 built-up hectare-equivalents (annual probability 0.100)
1-in-100: damage   113.4 built-up hectare-equivalents (annual probability 0.010)
1-in-500: damage   125.3 built-up hectare-equivalents (annual probability 0.002)
expected annual damage: 36.47 hectare-equivalents per year
  band down to 1-in-10 :  27.42 (75.2%)
  band down to 1-in-100:   7.85 (21.5%)
  band down to 1-in-500:   0.95 ( 2.6%)

Three quarters of the expected annual damage comes from the band between certainty and the 1-in-10 flood — the frequent, shallow events. The 1-in-500 scenario does the most damage when it happens and contributes 2.6% of the annual expectation. A risk map showing only the rare scenario shows the least of the answer.

Example 3 — how sensitive is it to the curve?

curves = {
    "as above": [0, 0.20, 0.35, 0.50, 0.80, 0.90],
    "shallower damage": [0, 0.35, 0.55, 0.70, 0.90, 0.95],
    "deeper damage": [0, 0.10, 0.20, 0.35, 0.65, 0.85],
}
for name, ys in curves.items():
    f = lambda d, ys=ys: np.clip(np.interp(d, [0, 0.5, 1, 2, 4, 6], ys), 0, 1)
    l = {rp: float((f(d) * built).sum()) for rp, d in scenarios.items()}
    prev_p, prev_loss, e = 1.0, 0.0, 0.0
    for rp, p in pairs:
        e += (prev_p - p) * (prev_loss + l[rp]) / 2
        prev_p, prev_loss = p, l[rp]
    e += probs[500] * l[500]
    print(f"{name:18} EAD {e:6,.2f}  (1-in-100 damage {l[100]:7,.1f})")

as above           EAD  36.47  (1-in-100 damage   113.4)
shallower damage   EAD  53.52  (1-in-100 damage   162.6)
deeper damage      EAD  23.57  (1-in-100 damage    75.5)

Why expected annual damage beats a scenario

Why expected annual damage and not a scenario

A scenario answers "what if the 1-in-100 flood happens". Expected annual damage answers "what does flooding cost per year on average", which is the number that compares places, supports cost–benefit analysis of defences, and can be summed across hazards. It needs the whole curve, not one point on it.

Why the integration is trapezoidal

Plot loss against annual exceedance probability and the expected annual loss is the area under the curve. With three scenarios there are three points, so the area is approximated by trapezia — and the stretch between the most frequent scenario and certainty, where nothing floods, and the tail beyond the rarest, matter as much as the points themselves.

Why the vulnerability curve dominates the uncertainty

The hazard maps are modelled but at least they are published, versioned products. The damage curve is usually borrowed from another country, another building stock and another decade. Changing it within a plausible range moved expected annual damage here from 23.6 to 53.5 hectare-equivalents a year — a factor of 2.3, from one hazard model. That is why any risk figure should be published with the curve it used.

Values, defences and indirect losses

• Values: multiply hectare-equivalents by a value density to get currency, and say where that density came from.
• Contents and structure damage differ; published curves usually separate them.
• Defended areas need a failure probability, not just the defended map.
• Duration and velocity matter for damage as well as depth; depth-only curves are a simplification.
• Indirect losses — displacement, business interruption — are not in these curves.
• Do not sum EAD across hazards without checking the exposure is not double-counted.

Internal links

• Hazard, exposure and vulnerability explained — the three terms
• Flood hazard maps explained — the scenarios used here
• Return periods and annual exceedance probability explained — the probabilities
• How to count buildings and people exposed to a hazard zone — the exposure term
• Exposure counts are double-counted or missed — overlapping zones
• How to run a weighted site suitability analysis in Python — weighted overlays, and how they differ from this

FAQ

What is expected annual damage?

The average damage per year, integrated over all hazard probabilities rather than computed for one scenario.

How many scenarios do I need?

At least three; more return periods make the integration less dependent on the shape assumed between them.

Where do depth–damage curves come from?

Published sources — JRC for Europe, HAZUS for the US, national flood guidance elsewhere. Not from spatial data.

Can I publish a risk map from one scenario?

You can publish that scenario's damage. Calling it risk without the probability integration is the common mistake.

What if I have no value data?

Report damage in exposure units — built-up hectares, dwellings — and let a reader apply their own values.

Which term dominates the uncertainty?

Usually vulnerability, because the curve is borrowed; the hazard maps at least come with versions and documentation.