p-phung/flood hazard logic

data/.gitignore

@@ -55,3 +55,6 @@ share/python-wheels/
 .installed.cfg
 *.egg
 MANIFEST
+
+# temp pipeline source data
+pipelines/flood/bronze

data/pipelines/flood/README.md

@@ -0,0 +1,91 @@
+# Flood Pipeline Logic (WIP)
+
+This folder contains the flood-specific forecast logic used by the pipeline framework.
+
+## Main script
+
+- `flood/forecast.py`
+  - Entry point for flood hazard logic via `calculate_flood_forecasts(...)`:
+  - Loads station points and admin areas through `DataProvider`, then combines them with local bronze flood inputs.
+  - Builds alerts, severity time series, admin-area exposure, and raster exposure through `DataSubmitter`.
+
+## Accompanying scripts in this folder
+
+- `extract_forecast.py`
+  - Samples GloFAS discharge values from (sliced) NetCDF rasters at station coordinates.
+  - Produces per-station, per-lead-time ensemble discharge series.
+
+- `compute_alert_extent.py`
+  - Resolves the flood extent raster to use for an alert from the available flood extent files.
+  - Selects the highest matched return period, falls back to the closest lower return period, then to `*_empty.tif`.
+
+- `determine_exposure.py`
+  - Reads the station-to-admin-area mapping and filters to place codes present in the loaded admin areas.
+  - Clips the selected flood extent raster to affected admin areas for raster exposure output.
+  - Computes an exposed-population raster and aggregates exposed population per place code.
+
+- `utils_raster.py`
+  - Utility functions for geospatial preprocessing:
+    - derive country bounding box from admin geometries,
+    - slice NetCDF to country bounds,
+    - clip rasters to bounding boxes,
+    - get raster extent for output metadata.
+
+## Bronze input data
+
+`bronze/` holds preprocessed flood inputs used by `forecast.py`:
+
+- `bronze/glofas/`
+  - GloFAS discharge NetCDF files (ensemble forecast source).
+  - Current code points to a local file such as `dis_00_YYYYMMDDHH.nc` and creates `_sliced.nc` files per country run.
+
+- `bronze/thresholds/`
+  - Country-specific threshold JSON files (`*_{ISO3}.json`).
+  - Used to convert station discharge values into return-period-based alert severities.
+
+- `bronze/station-district/`
+  - Station-to-admin mapping JSON (`{ISO3}_station_district_mapping.json`).
+  - Links station IDs to impacted place codes before flood extent clipping and exposure aggregation.
+
+- `bronze/population/`
+  - Country population raster (`{ISO3}.tif`).
+  - Used to create `{ISO3}_exposed.tif`, which is then aggregated over affected admin areas.
+
+- `bronze/flood_extents/`
+  - Flood extent rasters by return period (`flood_map_{ISO3}_rp*.tif`) plus empty fallback.
+  - Used to select, clip, and attach alert extent rasters to events.
+
+## `forecast.py` flow (read -> output)
+
+1. Load core inputs:
+   - Load GloFAS station metadata and target admin areas through `DataProvider`.
+   - Load threshold JSON, station-district mapping JSON, population raster path, and flood extent raster paths from local `bronze/` files.
+   - Stop early and record an error if stations or admin areas are missing.
+
+2. Build country spatial extent
+   - Compute country bounding box from target admin areas.
+   - Slice NetCDF files once to this bounding box.
+
+3. Process discharge per station
+   - Iterate through stations and currently limit processing to the first two station entries.
+   - Extract discharge ensemble values per lead time.
+   - Derive lead-time severities from thresholds.
+   - Skip stations with no threshold exceedance.
+
+4. Build alert payload
+   - Select the alert extent raster based on the matched return periods.
+   - Clip the flood extent to mapped admin areas and collect exposed place codes.
+   - Compute the exposed-population raster and aggregate exposed population per place code.
+
+5. Compute exposure
+   - Create one alert event per alerting station.
+   - Add severity time-series data for ensemble runs and the median discharge.
+   - Add admin-area population exposure per place code.
+   - Add raster exposure metadata for the generated `alert_extent_{station_code}.tif`.
+
+6. Write final output to local forecast folder
+   - `forecast.py` fills `DataSubmitter`.
+   - `pipelines/infra/run_forecasts.py` finalizes and writes `forecast.json`.
+   - Default local base path is `pipelines/output`, resulting in paths like:
+     - `pipelines/output/floods/{ISO3}/{timestamp}/forecast.json`
+   - In this repository this appears under `data/pipelines/output/floods/...`.

data/pipelines/flood/compute_alert_extent.py

@@ -0,0 +1,152 @@
+from __future__ import annotations
+
+import logging
+import os
+
+# Supported return period labels used by flood alerts.
+
+RETURN_PERIODS: dict[str, int] = {
+    "5yr": 5,
+    "10yr": 10,
+    "20yr": 20,
+    "25yr": 25,
+    "50yr": 50,
+    "100yr": 100
+}
+
+
+def compute_alert_extent(
+    time_interval_severities: list,
+    flood_extent_paths: list[str],
+) -> str:
+    """
+    Compute the flood extent raster for the alert station by resolving the appropriate return period raster.
+    Returns the path to the computed flood extent raster for the station.
+    """
+
+    return_period = _resolve_requested_return_period_value(time_interval_severities)
+
+    flood_extent_path = _resolve_flood_extent(
+        return_period=return_period,
+        flood_extent_paths=flood_extent_paths,
+    )
+    return flood_extent_path
+
+
+def _extract_return_period_value(path: str) -> int | None:
+    """
+    Extract the return period value from a flood extent raster file name.
+    Expected file name format: flood_map_{country}_rp{value}.tif, e.g. flood_map_ETH_rp20.tif -> 20
+    Returns the return period value as an integer, or None if it cannot be extracted e.g. flood_map_ETH_empty.tif.
+    """
+    filename = os.path.basename(path).lower()
+
+    filename_stem, _ = os.path.splitext(filename)
+    rp_index = filename_stem.rfind("_rp")
+    if rp_index == -1:
+        return None
+
+    value_text = filename_stem[rp_index + len("_rp") :]
+    if value_text.isdigit():
+        return int(value_text)
+
+    return None
+
+def _extract_return_period_label_value(return_period_label: str) -> int | None:
+    """
+    Extract the return period value from a return period label.
+    Expected label format: "Xyr", e.g. "5yr", "20yr".
+    Returns the return period value as an integer, or None if it cannot be extracted.
+    """
+    normalized_label = return_period_label.strip().lower()
+
+    if normalized_label.endswith("yr"):
+        value_text = normalized_label[:-2].strip()
+        if value_text.isdigit():
+            return int(value_text)
+
+    return None
+
+
+def _resolve_empty_flood_extent_path(flood_extent_paths: list[str]) -> str | None:
+    for flood_extent_path in flood_extent_paths:
+        filename = os.path.basename(flood_extent_path).lower()
+
+        if filename.endswith("_empty.tif"):
+            return flood_extent_path
+
+        filename_stem, _ = os.path.splitext(filename)
+        rp_index = filename_stem.rfind("_rp")
+        if rp_index == -1:
+            continue
+
+        base_without_suffix = filename_stem[:rp_index]
+        directory = os.path.dirname(flood_extent_path)
+        empty_tif_path = os.path.join(directory, f"{base_without_suffix}_empty.tif")
+        if os.path.exists(empty_tif_path):
+            return empty_tif_path
+
+    return None
+
+def _resolve_requested_return_period_value(time_interval_severities: list) -> int | None:
+    """
+    Resolve flood extent raster for the highest matched return period
+    """
+    highest_return_period = max(
+        time_interval_severities,
+        key=lambda s: s.median_discharge,
+    ).return_period
+
+    return_period = RETURN_PERIODS.get(highest_return_period)
+    if return_period is not None:
+        return return_period
+
+    parsed_return_period = _extract_return_period_label_value(highest_return_period)
+    if parsed_return_period is None:
+        logging.warning(
+            f"Unknown return period '{highest_return_period}', using empty fallback"
+        )
+
+    return parsed_return_period
+
+
+def _resolve_flood_extent(
+    return_period: int | None,
+    flood_extent_paths: list[str],
+) -> str:
+    """
+    Resolve the flood extent raster using this order:
+    1. Exact return period raster.
+    2. Closest lower-or-equal available return period raster.
+    3. flood_map_{country}_empty.tif fallback raster (guaranteed to exist).
+    """
+
+    available_paths_by_value: dict[int, str] = {}
+    for flood_extent_path in flood_extent_paths:
+        if not os.path.exists(flood_extent_path):
+            continue
+        value = _extract_return_period_value(flood_extent_path)
+        if value is not None:
+            available_paths_by_value[value] = flood_extent_path
+
+    if return_period is not None and return_period in available_paths_by_value:
+        return available_paths_by_value[return_period]
+
+    if return_period is not None:
+        fallback_value = max(
+            (value for value in available_paths_by_value.keys() if value <= return_period),
+            default=None,
+        )
+        if fallback_value is not None:
+            return available_paths_by_value[fallback_value]
+
+    # Return empty flood extent
+    empty_path = _resolve_empty_flood_extent_path(flood_extent_paths)
+    if empty_path is None:
+        raise FileNotFoundError(
+            "Could not resolve flood extent raster: no suitable return period raster "
+            "or no empty fallback raster was found among existing files."
+        )
+
+    return empty_path
+

data/pipelines/flood/determine_alerts.py

@@ -0,0 +1,145 @@
+from __future__ import annotations
+
+import logging
+from dataclasses import dataclass
+from typing import TypedDict
+
+import numpy as np
+
+from pipelines.flood.extract_forecast import TimeIntervalDischarge
+from pipelines.infra.data_types.location_point import LocationPoint
+
+MINIMUM_RETURN_PERIOD = "1.5yr"
+
+
+@dataclass
+class TimeIntervalSeverity:
+    time_interval_start: str
+    time_interval_end: str
+    median_discharge: float
+    ensemble_discharges: list[float]
+    return_period: str
+
+
+@dataclass
+class AlertStation:
+    station_code: str
+    station: LocationPoint
+    time_interval_severities: list[TimeIntervalSeverity]
+
+
+class ReturnPeriodThresholdValue(TypedDict):
+    return_period: float
+    threshold_value: float
+
+
+class ReturnPeriodThresholds(TypedDict):
+    station_code: str
+    thresholds: list[ReturnPeriodThresholdValue]
+
+
+def determine_temporal_extent(
+    station_code: str,
+    time_interval_discharges: list[TimeIntervalDischarge],
+    thresholds: list[ReturnPeriodThresholds],
+    minimum_return_period: str = MINIMUM_RETURN_PERIOD,
+) -> list[TimeIntervalSeverity]:
+    """
+    Compute lead time severities for one station by comparing the median
+    ensemble discharge against return period thresholds.
+    Returns a list of lead time severities for the station.
+    """
+    station_thresholds = _prepare_station_threshold(
+        station_code, thresholds, minimum_return_period
+    )
+    if station_thresholds is None:
+        return []
+
+    time_interval_severities: list[TimeIntervalSeverity] = []
+    for time_interval_discharge in time_interval_discharges:
+        ensemble_array = np.asarray(
+            time_interval_discharge.ensemble_discharges,
+            dtype=float,
+        )
+        if np.isnan(ensemble_array).all():
+            continue
+        median_discharge = float(np.nanmedian(ensemble_array))
+        return_period = _match_return_period(median_discharge, station_thresholds)
+        if return_period is not None:
+            time_interval_severities.append(
+                TimeIntervalSeverity(
+                    time_interval_start=time_interval_discharge.time_interval_start,
+                    time_interval_end=time_interval_discharge.time_interval_end,
+                    median_discharge=median_discharge,
+                    ensemble_discharges=time_interval_discharge.ensemble_discharges,
+                    return_period=return_period,
+                )
+            )
+
+    return time_interval_severities
+
+
+def _format_return_period_label(return_period: float) -> str:
+    return f"{return_period:g}yr"
+
+
+def _get_station_return_period_thresholds(
+    thresholds: list[ReturnPeriodThresholds],
+    station_code: str,
+) -> dict[str, float] | None:
+    for station_threshold in thresholds:
+        if station_threshold["station_code"] != station_code:
+            continue
+
+        return {
+            _format_return_period_label(threshold["return_period"]): threshold[
+                "threshold_value"
+            ]
+            for threshold in station_threshold["thresholds"]
+        }
+
+    return None
+
+
+def _match_return_period(
+    discharge: float,
+    station_thresholds: dict[str, float],
+) -> str | None:
+    """
+    Find the highest return period whose threshold the discharge exceeds.
+    Thresholds are expected as e.g. {"2yr": 100, "5yr": 200, "10yr": 350, ...}.
+    Returns the label of the highest exceeded return period, or None if none exceeded.
+    """
+    # sort thresholds descending by value to find the highest exceeded efficiently
+    for return_period, threshold_value in sorted(
+        station_thresholds.items(), key=lambda item: item[1], reverse=True
+    ):
+        if discharge > threshold_value:
+            return return_period
+    return None
+
+
+def _prepare_station_threshold(
+    station_code: str,
+    thresholds: list[ReturnPeriodThresholds],
+    minimum_return_period: str = MINIMUM_RETURN_PERIOD,
+) -> dict[str, float] | None:
+    """
+    Retrieve and validate station return-period thresholds.
+    Returns the station thresholds dict, or None if validation fails.
+    """
+    station_thresholds = _get_station_return_period_thresholds(thresholds, station_code)
+    if station_thresholds is None:
+        logging.warning(
+            f"No return period thresholds for station {station_code}, skipping"
+        )
+        return None
+
+    if minimum_return_period not in station_thresholds:
+        logging.warning(
+            f"Return period '{minimum_return_period}' not found for station {station_code}, skipping"
+        )
+        return None
+
+    return station_thresholds
+