Evaluating Observed Flood Extent from the GFM against Forecast Flood Extent from the EFAS or GloFAS Rapid Flood Mapping layer

With this notebook, we show how to use the new GFM STAC service to derive the maximum flood extent for a flood event, extract the corresponding forecast product from the EFAS or GloFAS Rapid Flood Mapping (RFM) layer and perform a simple evaluation.

This example notebook is structured as follows:

  1. (Import and install the necessary python packages)

  2. Extract GFM Data using the GFM STAC service (ensemble flood extent, reference water mask, exclusion mask)

  3. Download the data as a Zarr store for later use (recommended, but optional)

  4. Mask and save data as GeoTiff

  5. Download Forecasted Flood Extent from the EFAS or GloFAS Rapid Flood Mapping layer

  6. Compare the Observed and Forecast Flood Extent

We present an example analysis of the flooding that occurred in southern Mozambique between the 7th - 28th January 2026.

Another use case in Morocco between 1st - 15th February 2026 is also shown.

Please note you must have an EFAS account to access EFAS data. GloFAS data is publicly accessible and does not require an account.

0. Install and import necessary python modules

!pip install pyproj xarray shapely pystac_client odc-stac matplotlib rioxarray geopandas ipykernel

# Some necessary imports

import shutil
from pystac_client import Client
from datetime import datetime
from odc import stac as odc_stac
from pathlib import Path
import pyproj
import rioxarray # noqa
import xarray as xr
from shapely.geometry import box
import geopandas as gpd
import rasterio
from rasterio.features import rasterize
from rasterio.transform import from_bounds

import numpy as np
import requests
from zipfile import ZipFile

1. Extract GFM Data using the GFM STAC service

Firstly define the GFM product you wish to download.

Then define a time range of interest and the coordinates of the area of interest (AOI). The coordinates of the AOI can be retrieved either via the GFM Portal (portal.gfm.eodc.eu) or a tool like bboxfinder.com.

These parameters are passed to the GFM STAC service to find the available data.

Use Case 1: Southern Mozambique January 2026

use_case_name = "Southern_Mozambique_Jan_2026"

# Define bounding box (West, South, East, North) - Southern Mozambique
aoi = box(32.048307386, -25.272982457, 34.109853235, -23.221383919)

# Define time range
time_range = (datetime(2026, 1, 7), datetime(2026, 1, 28))

# Flood forecast date
forecast_date = datetime(2026, 1, 7)

Use Case 2: Morocco February 2026

use_case_name = "Morocco_Feb_2026"

# Define bounding box (West, South, East, North) - Southern Mozambique
aoi = box(32.048307386, -25.272982457, 34.109853235, -23.221383919)

# Define time range
time_range = (datetime(2026, 2, 1), datetime(2026, 2, 15))

# Flood forecast date
forecast_date = datetime(2026, 2, 1)

Extract GFM Data

# Define the output path, where you want to save temporary data and final results
output_path = Path("./gfm_rfm_evaluation") / use_case_name
output_path.mkdir(parents=True, exist_ok=True)

# eodc STAC API URL
eodc_catalog = Client.open("https://stac.eodc.eu/api/v1")

# Define search query using pystac_client
search = eodc_catalog.search(
    max_items=1000,
    collections=["GFM"],
    bbox=aoi.bounds,
    datetime=time_range
)

# Get all items
items = search.item_collection()
if not items:
    raise RuntimeError("No items found — check your AOI, time range, or collection name.")

print(f"Found {len(items)} items matching filter criteria.")

# Retrieve the coordinate reference system (CRS) from one of the found items
crs_equi7 = pyproj.CRS.from_wkt(items[0].properties["proj:wkt2"])

Lazy-load the returned data

Now the returned data are lazily loaded into an xarray dataset in the Equi7Grid projection.

We define the asset names we want to use in this notebook. In our case, that is ‘ensemble_flood_extent’, ‘reference water mask’ and ‘exclusion_mask’.

# Define the asset names you want to use
asset_names = ["ensemble_flood_extent", "reference_water_mask", "exclusion_mask"]

# Retrieve the resolution of the data
resolution = items[0].properties.get("gsd", 20)  # fallback to 20m if gsd missing

# By setting chunks, the odc library "lazy loads" the data. -1 for time means
# that all time steps are included in one chunk. Reduce the chunk size for x
# and/or y if your kernel crashes due to out of memory issues
raw_data = odc_stac.load(
    items,
    bbox=aoi.bounds,   # Define the bounding box for the area of interest
    crs=crs_equi7,   # Set the coordinate reference system
    bands=asset_names,   # Specify the bands to load
    resolution=resolution,   # Set the resolution of the data
    dtype='uint8',   # Define the data type
    groupby="solar_day",
    chunks={"x": 4096, "y": 4096, "time": -1},
)

3. Mask and save data as GeoTiff

We mask the data to exclude invalid values such as nodata (=255) and non-flood (=0).

After masking, we save a GeoTiff file for each GFM layer for easy visualization in e.g. QGIS.

# Only include data which is not nodata and not 0
nodata = 255
valid_data = raw_data.where((raw_data != nodata) & (raw_data != 0))

# Create the sum over the time dimension
# This is done for all specified assets
data = valid_data.sum(dim="time").astype("uint8")

# Binarize all layers: any positive value → 1, zero stays 0
for var in asset_names:
    data[var] = data[var].clip(0, 1)

rio_options = {
    "compress": "LZW",
    "tiled": True,
    "blockxsize": 512,
    "blockysize": 512
}

# Save as Raster in Equi7Grid
data["ensemble_flood_extent"].rio.to_raster(output_path / "ensemble_flood_extent_equi7.tif", **rio_options)
data["reference_water_mask"].rio.to_raster(output_path / "reference_water_mask_equi7.tif", **rio_options)
data["exclusion_mask"].rio.to_raster(output_path / "exclusion_mask_equi7.tif", **rio_options)

4. Download Forecasted Flood Extent from the EFAS or GloFAS Rapid Flood Mapping layer

Specify the data source ("efas" or "glofas") and the date of the forecast.

For EFAS, you will also need a web token, which you can generate by logging in to https://european-flood.emergency.copernicus.eu/en/efas-web-services and copying the web token shown in the table. Note: the token needs to be regenerated after every use.

GloFAS data is publicly accessible and does not require a token.

source = "glofas"  # Set to "efas" or "glofas"
fcast_date = forecast_date.strftime("%Y%m%d00")  # date of the forecast as string in YYYYMMDDHH format

# EFAS requires a token; GloFAS does not
# only needed when source == "efas", token must be regenerated after every use
efas_token = 'your_efas_token'

date_str = f'{fcast_date[0:4]}-{fcast_date[4:6]}-{fcast_date[6:8]}T{fcast_date[8:10]}:00Z'
base_url = f'https://european-flood.emergency.copernicus.eu/api/fms/download/{source}/RapidFloodMapping/{date_str}'
url_rfm = f'{base_url}?token={efas_token}' if source == "efas" else base_url

Download the forecast flood extent file, load it in then clip to AOI

This will save a zip file into the same folder where you have saved this notebook. An ESRI shapefile containing the forecast flood extent will be extracted from this zip file into the same folder.

response = requests.get(url_rfm)
response.raise_for_status()

out_zip = output_path / f'RFM_{fcast_date}.zip'
with open(out_zip, 'wb') as file:
  file.write(response.content)

with ZipFile(out_zip, 'r') as zip_ref:
  zip_ref.extractall(output_path)

Load the Forecasted Flood Extent file

The forecast flood extent shapefile is loaded as a GeoPandas GeoDataFrame. It is then clipped to the extent of the aoi which you defined previously

rfm_file = output_path / f'FloodMaskMerged{fcast_date}.shp'
rfm_data = gpd.read_file(rfm_file)

# Clip the Rapid Flood Mapping to your drawn area
rfm_intersect = gpd.clip(rfm_data, gpd.GeoDataFrame(index=[0], crs='epsg:4326', geometry=[aoi]))

if rfm_intersect.empty:
    raise RuntimeError("No RFM data found within the AOI — check your bounding box or forecast date.")

# Get the bounding coordinates of the GFM data
dx = data.x.values[1]-data.x.values[0]
dy = data.y.values[1]-data.y.values[0]

xmin = data.x.values.min() - (dx/2)
xmax = data.x.values.max() + (dx/2)

ymin = data.y.values.min() - abs(dy/2)
ymax = data.y.values.max() + abs(dy/2)

gfm_flood = data["ensemble_flood_extent"].values

transform = from_bounds(xmin, ymin, xmax, ymax, gfm_flood.shape[1], gfm_flood.shape[0])

rfm_intersect_equi7 = rfm_intersect.to_crs(crs_equi7)

# Get geometry and corresponding value from the GeoDataFrame
shapes = ((geom, 1) for geom in rfm_intersect_equi7.geometry)

# Rasterize the geometries into the numpy array
rfm_raster = rasterize(
    shapes,
    out_shape=gfm_flood.shape,
    transform=transform,
    fill=0,
    dtype=np.uint8
)

5. Compare the Observed and Forecast Flood Extent

Use a contingency table approach to classify each grid cell as one of the following:

  1. Hit - where flooding is present in both the GFM & RFM forecast

  2. False alarm - flooding only present in the RFM forecast

  3. Miss - flooding only present in the GFM observations

  4. Correct negative - no flooding present in either GFM or RFM data

We exclude areas which are either in the reference water mask or the exclusion mask.

From this classification, we compute the following performance scores:

  1. CSI (Critical Success Index), range 0-1 where 1 = perfect agreement between forecast and observations

  2. False Alarm Ratio, range 0-1 where 0 is optimal, shows the fraction of grid cells which were forecasted as being flooded which were incorrect

  3. Hit Rate, range 0-1 where 1 is optimal, shows the fraction of observed flooded grid cells which were correctly forecasted

# compare GFM and RFM data

gfm_em = data["exclusion_mask"].values
gfm_rwm = data["reference_water_mask"].values

# compute the total number of hits, false alarms, misses and correct negatives
tp = np.where((gfm_flood == 1) & (rfm_raster == 1) & (gfm_rwm < 1) & (gfm_em == 0))[0].shape[0] # hits
fp = np.where((gfm_flood == 0) & (rfm_raster == 1) & (gfm_rwm < 1) & (gfm_em == 0))[0].shape[0] # false alarms
fn = np.where((gfm_flood == 1) & (rfm_raster == 0) & (gfm_rwm < 1) & (gfm_em == 0))[0].shape[0] # misses
tn = np.where((gfm_flood == 0) & (rfm_raster == 0) & (gfm_rwm < 1) & (gfm_em == 0))[0].shape[0] # correct negatives

csi = tp / (tp + fp + fn) if (tp + fp + fn) > 0 else float('nan')  # critical success index (CSI)
far = fp / (tp + fp) if (tp + fp) > 0 else float('nan')            # false alarm ratio (FAR)
hr  = tp / (tp + fn) if (tp + fn) > 0 else float('nan')            # hit rate (HR)

print(f'CSI: {"%.2f" % csi}')
print(f'False Alarm Ratio: {"%.2f" % far}')
print(f'Hit Rate: {"%.2f" % hr}')

# Create a grid for the evaluation results to show on a map
# create a new grid which shows whether each grid cell is a hit, false alarm or miss
eval_grid = np.zeros(gfm_flood.shape, dtype=np.uint8)
eval_grid[(gfm_flood == 1) & (rfm_raster == 1) & (gfm_rwm < 1) & (gfm_em == 0)] = 1  # Hits
eval_grid[(gfm_flood == 1) & (rfm_raster == 0) & (gfm_rwm < 1) & (gfm_em == 0)] = 2  # Misses
eval_grid[(gfm_flood == 0) & (rfm_raster == 1) & (gfm_rwm < 1) & (gfm_em == 0)] = 3  # False alarms

# Save RFM data as raster
with rasterio.open(output_path / 'rfm_data_equi7.tif', 'w', driver='GTiff', height=rfm_raster.shape[0],
                   width=rfm_raster.shape[1], count=1, dtype=rfm_raster.dtype,
                   crs=crs_equi7, transform=transform, compress="LZW") as dst:
    dst.write(rfm_raster, 1)

# Save RFM GFM evaluation results as raster
with rasterio.open(output_path / 'rfm_gfm_eval_equi7.tif', 'w', driver='GTiff', height=eval_grid.shape[0],
                   width=eval_grid.shape[1], count=1, dtype=eval_grid.dtype,
                   crs=crs_equi7, transform=transform, compress="LZW") as dst:
    dst.write(eval_grid, 1)

5. Visualize outputs in your favourite GIS application (e.g. QGIS, ArcGIS)

Conclusions

This notebook shows how you can download specific GFM data using the GFM STAC service and the EFAS or GloFAS Rapid Flood Mapping forecast layer and compare the data, with the idea of evaluating the performance of the forecast layer.

In general it could be seen that:

  • The forecasted inundation extent generally over-predicted when compared to GFM

  • Because forecast data is coarser resolution

    • Smoothed representation of floodplain topography, e.g. doesn’t represent levees

  • Therefore you need to be careful regarding the precision for which the forecast information can be used

    • It can highlight general floodplain areas at risk but be careful when identifying specific towns/cities affected