Skip to content

food_security

Vegetation and climate stress assessment for food-security risk. Entry point: FoodSecurityUseCase.

from climate_change import FoodSecurityUseCase

Use case

use_case

FoodSecurityUseCase

FoodSecurityUseCase(dask_engine: DaskEngine)

Bases: BaseUseCase

Entry point for the food security assessment domain.

Minimal config (all optional — defaults match Marsabit County, Kenya 2018–2023):

{ "aoi_geojson": {"type": "Polygon", "coordinates": [...]}, "gee_project": "your-gee-project-id", "start_date": "2018-01-01", "end_date": "2023-12-31", "lt_baseline_start": "2001-01-01", # long-term baseline for VCI/TCI "s2_start": "2020-01-01", # Sentinel-2 MNDWI start "model_type": "rf", # "rf" | "xgboost" | "ensemble" "n_pixels": 3000, "scale": 1000, # metres — MODIS / CHIRPS native "output_dir": "outputs", "prefix": "food_security", }

Source code in food_security/use_case.py
def __init__(self, dask_engine: DaskEngine) -> None:
    super().__init__(dask_engine)

fetch_data

fetch_data(config: AnalysisConfig) -> dict

Authenticate GEE and parse the AOI geometry (see _fetch_from_dict).

Source code in food_security/use_case.py
def fetch_data(self, config: AnalysisConfig) -> dict:
    """Authenticate GEE and parse the AOI geometry (see `_fetch_from_dict`)."""
    flat: dict = {
        "aoi_geojson": config.aoi_geojson,
        "start_date": config.start_date,
        "end_date": config.end_date,
        "country": config.country,
        **config.extra_params,
    }
    return self._fetch_from_dict(flat)

preprocess

preprocess(raw_data: dict, config: AnalysisConfig) -> dict

Build the feature stack and sample training pixels (see _preprocess_raw).

Source code in food_security/use_case.py
def preprocess(self, raw_data: dict, config: AnalysisConfig) -> dict:
    """Build the feature stack and sample training pixels (see `_preprocess_raw`)."""
    return self._preprocess_raw(raw_data)

run_model

run_model(features: dict, config: AnalysisConfig) -> AnalysisOutput

Train models, export the risk COG, and package an AnalysisOutput.

Source code in food_security/use_case.py
def run_model(self, features: dict, config: AnalysisConfig) -> AnalysisOutput:
    """Train models, export the risk COG, and package an `AnalysisOutput`."""
    model_type = config.extra_params.get("model_type", "rf")
    dict_config = {"model_type": model_type, **config.extra_params}
    result = self._run_model_dict(features, dict_config)
    total_area_ha = _aoi_area_ha(config.aoi_geojson)
    spatial_grid = predict_food_security_grid(
        rf_model=features["_rf"],
        xgb_model=features["_xgb"],
        scaler=features["_scaler"],
        datasets=features["datasets"],
        model_type=model_type,
        aoi_geojson=config.aoi_geojson,
    )
    self._apply_spatial_risk_summary(result, spatial_grid, total_area_ha)
    if total_area_ha:
        result["stats"]["total_area_ha"] = round(total_area_ha, 1)
        _attach_risk_area(result.get("charts", {}), total_area_ha)
    raster_paths: dict[str, str] | None = None
    raster_error: str | None = None
    try:
        raster_paths = export_food_security_cog(
            rf_model=features["_rf"],
            xgb_model=features["_xgb"],
            scaler=features["_scaler"],
            datasets=features["datasets"],
            output_dir=dict_config.get("output_dir", "outputs"),
            prefix=dict_config.get("prefix", "food_security"),
            model_type=model_type,
            aoi_geojson=config.aoi_geojson,
        )
    except Exception as exc:
        raster_error = str(exc)

    geojson_features = self._risk_grid_to_geojson(spatial_grid)
    for p in result.pop("_sample_points", []):
        geojson_features.append(
            {
                "type": "Feature",
                "geometry": {"type": "Point", "coordinates": [p["lon"], p["lat"]]},
                "properties": {"risk_class": p["risk_class"]},
            }
        )
    metadata = {
        "model": model_type,
        "country": config.country,
        "start_date": config.start_date,
        "end_date": config.end_date,
        "raster": raster_paths or {},
        "spatial_resolution_m": dict_config.get("scale"),
        "n_pixels_sampled": dict_config.get("n_pixels"),
    }
    if raster_error:
        metadata["raster_error"] = raster_error

    return AnalysisOutput(
        module="food_security",
        geojson={"type": "FeatureCollection", "features": geojson_features},
        raster_path=raster_paths,
        stats={**result["stats"], "country": config.country},
        shap=result.get("charts", {}).get("shap"),
        charts=result.get("charts", {}),
        metadata=metadata,
    )

run

run(config: dict) -> dict

Single AOI analysis — full pipeline in one call.

Source code in food_security/use_case.py
def run(self, config: dict) -> dict:
    """Single AOI analysis — full pipeline in one call."""
    raw_data = self._fetch_from_dict(config)
    features = self._preprocess_raw(raw_data)
    result = self._run_model_dict(features, config)
    total_area_ha = _aoi_area_ha(config.get("aoi_geojson"))
    spatial_grid = predict_food_security_grid(
        rf_model=features["_rf"],
        xgb_model=features["_xgb"],
        scaler=features["_scaler"],
        datasets=features["datasets"],
        model_type=config.get("model_type", "rf"),
        aoi_geojson=config.get("aoi_geojson"),
    )
    self._apply_spatial_risk_summary(result, spatial_grid, total_area_ha)
    if total_area_ha:
        result["stats"]["total_area_ha"] = round(total_area_ha, 1)
        _attach_risk_area(result.get("charts", {}), total_area_ha)

    cog_paths: dict[str, str] | None = None
    try:
        cog_paths = export_food_security_cog(
            rf_model=features["_rf"],
            xgb_model=features["_xgb"],
            scaler=features["_scaler"],
            datasets=features["datasets"],
            output_dir=config.get("output_dir", "outputs"),
            prefix=config.get("prefix", "food_security"),
            model_type=config.get("model_type", "rf"),
            aoi_geojson=config.get("aoi_geojson"),
        )
    except Exception as exc:
        result["raster_error"] = str(exc)
    result["raster"] = cog_paths or {}
    result["geojson"] = {
        "type": "FeatureCollection",
        "features": self._risk_grid_to_geojson(spatial_grid),
    }
    result["stats"].update(
        {
            "bbox": features["bbox"],
            "prefix": config.get("prefix", "food_security"),
        }
    )
    return result

run_date_ranges

run_date_ranges(config: dict, date_ranges: list[dict]) -> list[dict]

Run the same AOI over multiple date-range configurations in parallel.

Source code in food_security/use_case.py
def run_date_ranges(self, config: dict, date_ranges: list[dict]) -> list[dict]:
    """Run the same AOI over multiple date-range configurations in parallel."""
    merged = [{**config, **dr} for dr in date_ranges]
    return self._run_parallel(merged)

run_multi_regions

run_multi_regions(configs: list[dict]) -> list[dict]

Run multiple AOI configs (same timeframe, different regions) in parallel.

Source code in food_security/use_case.py
def run_multi_regions(self, configs: list[dict]) -> list[dict]:
    """Run multiple AOI configs (same timeframe, different regions) in parallel."""
    return self._run_parallel(configs)

Features

features

fetch_vci_tci

fetch_vci_tci(aoi: Geometry, lt_start: str, study_start: str, study_end: str, scale: int = 1000) -> xr.Dataset

Compute VCI and TCI from MODIS NDVI and LST relative to the long-term baseline.

VCI = (NDVI_current − NDVI_min) / (NDVI_max − NDVI_min) × 100 TCI = (LST_max − LST_current) / (LST_max − LST_min) × 100

Returns Dataset with variables 'vci' and 'tci'.

Source code in food_security/features.py
def fetch_vci_tci(
    aoi: ee.Geometry,
    lt_start: str,
    study_start: str,
    study_end: str,
    scale: int = 1000,
) -> xr.Dataset:
    """
    Compute VCI and TCI from MODIS NDVI and LST relative to the long-term baseline.

    VCI = (NDVI_current − NDVI_min) / (NDVI_max − NDVI_min) × 100
    TCI = (LST_max − LST_current) / (LST_max − LST_min) × 100

    Returns Dataset with variables 'vci' and 'tci'.
    """

    def _scale_ndvi(img: ee.Image) -> ee.Image:
        return img.multiply(0.0001).copyProperties(img, ["system:time_start"])

    def _scale_lst(img: ee.Image) -> ee.Image:
        return img.multiply(0.02).subtract(273.15).copyProperties(img, ["system:time_start"])

    # Long-term NDVI baseline
    modis_lt = (
        ee.ImageCollection("MODIS/061/MOD13A3")
        .filterBounds(aoi)
        .filterDate(lt_start, study_start)
        .select("NDVI")
        .map(_scale_ndvi)
    )
    ndvi_min = modis_lt.min().rename("ndvi_min")
    ndvi_max = modis_lt.max().rename("ndvi_max")

    # Study-period NDVI current
    ndvi_current = (
        ee.ImageCollection("MODIS/061/MOD13A3")
        .filterBounds(aoi)
        .filterDate(study_start, study_end)
        .select("NDVI")
        .map(_scale_ndvi)
        .mean()
        .rename("ndvi_current")
    )
    vci_img = (
        ndvi_current.subtract(ndvi_min)
        .divide(ndvi_max.subtract(ndvi_min).add(1e-6))
        .multiply(100)
        .rename("vci")
        .clip(aoi)
    )

    # Long-term LST baseline (10th and 90th percentiles)
    lst_lt = (
        ee.ImageCollection("MODIS/061/MOD11A2")
        .filterBounds(aoi)
        .filterDate(lt_start, study_end)
        .select("LST_Day_1km")
        .map(_scale_lst)
    )
    lst_min = lst_lt.reduce(ee.Reducer.percentile([10])).rename("lst_min")
    lst_max = lst_lt.reduce(ee.Reducer.percentile([90])).rename("lst_max")

    # Study-period LST current
    lst_current = (
        ee.ImageCollection("MODIS/061/MOD11A2")
        .filterBounds(aoi)
        .filterDate(study_start, study_end)
        .select("LST_Day_1km")
        .map(_scale_lst)
        .mean()
        .rename("lst_current")
    )
    tci_img = (
        lst_max.subtract(lst_current)
        .divide(lst_max.subtract(lst_min).add(1e-6))
        .multiply(100)
        .rename("tci")
        .clip(aoi)
    )

    stack = vci_img.addBands(tci_img)
    url = stack.getDownloadURL(
        {"region": aoi, "scale": scale, "crs": "EPSG:4326", "format": "GEO_TIFF"}
    )
    raw = _download_band(url)
    raw = raw.assign_coords(band=["vci", "tci"])
    return raw.to_dataset(dim="band").rename({"x": "lon", "y": "lat"})

fetch_rainfall_anomaly_pct

fetch_rainfall_anomaly_pct(aoi: Geometry, lt_start: str, study_start: str, study_end: str, scale: int = 1000) -> xr.Dataset

Pixel-wise CHIRPS rainfall anomaly as percentage vs long-term daily mean. anomaly_pct = (study_mean − LT_mean) / (LT_mean + ε) × 100

Returns Dataset with variable 'rainfall_anom_pct'.

Source code in food_security/features.py
def fetch_rainfall_anomaly_pct(
    aoi: ee.Geometry,
    lt_start: str,
    study_start: str,
    study_end: str,
    scale: int = 1000,
) -> xr.Dataset:
    """
    Pixel-wise CHIRPS rainfall anomaly as percentage vs long-term daily mean.
    anomaly_pct = (study_mean − LT_mean) / (LT_mean + ε) × 100

    Returns Dataset with variable 'rainfall_anom_pct'.
    """
    chirps = ee.ImageCollection("UCSB-CHG/CHIRPS/DAILY").filterBounds(aoi)
    lt_mean = (
        chirps.filterDate(lt_start, study_start)
        .select("precipitation")
        .mean()
        .multiply(365)
        .rename("lt_mean")
    )
    study_mean = (
        chirps.filterDate(study_start, study_end)
        .select("precipitation")
        .mean()
        .multiply(365)
        .rename("study_mean")
    )
    rain_anom = (
        study_mean.subtract(lt_mean)
        .divide(lt_mean.add(1e-6))
        .multiply(100)
        .rename("rainfall_anom_pct")
        .clip(aoi)
    )
    url = rain_anom.getDownloadURL(
        {"region": aoi, "scale": scale, "crs": "EPSG:4326", "format": "GEO_TIFF"}
    )
    da = _download_band(url).squeeze()
    return xr.Dataset({"rainfall_anom_pct": da.rename({"x": "lon", "y": "lat"})})

fetch_ndvi_slope_img

fetch_ndvi_slope_img(aoi: Geometry, lt_start: str, study_end: str, scale: int = 1000) -> xr.Dataset

Pixel-wise MODIS NDVI linear trend (NDVI per year) from lt_start to study_end. Returns Dataset with variable 'ndvi_slope'.

Source code in food_security/features.py
def fetch_ndvi_slope_img(
    aoi: ee.Geometry,
    lt_start: str,
    study_end: str,
    scale: int = 1000,
) -> xr.Dataset:
    """
    Pixel-wise MODIS NDVI linear trend (NDVI per year) from lt_start to study_end.
    Returns Dataset with variable 'ndvi_slope'.
    """

    def _scale_ndvi(img: ee.Image) -> ee.Image:
        return img.multiply(0.0001).copyProperties(img, ["system:time_start"])

    def _add_time(img: ee.Image) -> ee.Image:
        t = img.date().difference(ee.Date(lt_start), "year")
        return img.addBands(ee.Image(t).rename("time").float())

    modis = (
        ee.ImageCollection("MODIS/061/MOD13A3")
        .filterBounds(aoi)
        .filterDate(lt_start, study_end)
        .select("NDVI")
        .map(_scale_ndvi)
    )
    trend = modis.map(_add_time).select(["time", "NDVI"]).reduce(ee.Reducer.linearFit())
    ndvi_slope = trend.select("scale").rename("ndvi_slope").clip(aoi)

    url = ndvi_slope.getDownloadURL(
        {"region": aoi, "scale": scale, "crs": "EPSG:4326", "format": "GEO_TIFF"}
    )
    da = _download_band(url).squeeze()
    return xr.Dataset({"ndvi_slope": da.rename({"x": "lon", "y": "lat"})})

fetch_mndwi

fetch_mndwi(aoi: Geometry, start: str, end: str, scale: int = 1000) -> xr.Dataset

Sentinel-2 MNDWI = (Green − SWIR1) / (Green + SWIR1) median composite. Returns Dataset with variable 'mndwi'.

Source code in food_security/features.py
def fetch_mndwi(
    aoi: ee.Geometry,
    start: str,
    end: str,
    scale: int = 1000,
) -> xr.Dataset:
    """
    Sentinel-2 MNDWI = (Green − SWIR1) / (Green + SWIR1) median composite.
    Returns Dataset with variable 'mndwi'.
    """

    def _add_mndwi(img: ee.Image) -> ee.Image:
        return (
            img.normalizedDifference(["B3", "B11"])
            .rename("mndwi")
            .copyProperties(img, ["system:time_start"])
        )

    s2 = (
        ee.ImageCollection("COPERNICUS/S2_SR_HARMONIZED")
        .filterBounds(aoi)
        .filterDate(start, end)
        .filter(ee.Filter.lt("CLOUDY_PIXEL_PERCENTAGE", 30))
        .map(_add_mndwi)
    )
    mndwi_img = s2.median().clip(aoi)
    url = mndwi_img.getDownloadURL(
        {"region": aoi, "scale": scale, "crs": "EPSG:4326", "format": "GEO_TIFF"}
    )
    da = _download_band(url).squeeze()
    return xr.Dataset({"mndwi": da.rename({"x": "lon", "y": "lat"})})

fetch_terrain_slope

fetch_terrain_slope(aoi: Geometry, scale: int = 1000) -> xr.Dataset

SRTM-derived terrain slope in degrees. Returns Dataset with variable 'slope_terrain'.

Source code in food_security/features.py
def fetch_terrain_slope(aoi: ee.Geometry, scale: int = 1000) -> xr.Dataset:
    """SRTM-derived terrain slope in degrees. Returns Dataset with variable 'slope_terrain'."""
    slope = ee.Terrain.slope(ee.Image("USGS/SRTMGL1_003")).clip(aoi).rename("slope_terrain")
    url = slope.getDownloadURL(
        {"region": aoi, "scale": scale, "crs": "EPSG:4326", "format": "GEO_TIFF"}
    )
    da = _download_band(url).squeeze()
    return xr.Dataset({"slope_terrain": da.rename({"x": "lon", "y": "lat"})})

fetch_landcover

fetch_landcover(aoi: Geometry, scale: int = 1000) -> xr.Dataset

ESA WorldCover v200 normalised to [0, 1] (class value / 100). Returns Dataset with variable 'land_cover'.

Source code in food_security/features.py
def fetch_landcover(aoi: ee.Geometry, scale: int = 1000) -> xr.Dataset:
    """
    ESA WorldCover v200 normalised to [0, 1] (class value / 100).
    Returns Dataset with variable 'land_cover'.
    """
    lc = (
        ee.ImageCollection("ESA/WorldCover/v200")
        .first()
        .select("Map")
        .divide(100)
        .rename("land_cover")
        .clip(aoi)
    )
    url = lc.getDownloadURL(
        {"region": aoi, "scale": scale, "crs": "EPSG:4326", "format": "GEO_TIFF"}
    )
    da = _download_band(url).squeeze()
    return xr.Dataset({"land_cover": da.rename({"x": "lon", "y": "lat"})})

build_feature_datasets

build_feature_datasets(aoi: Geometry, config: dict) -> dict[str, xr.Dataset]

Download all seven food-security feature bands from GEE in parallel. The six fetch calls are independent HTTP requests; they run concurrently via DaskEngine.run_io_parallel (ThreadPoolExecutor) sharing the GEE session. Returns a dict keyed by band group name, consumed by cog_export.

Source code in food_security/features.py
def build_feature_datasets(aoi: ee.Geometry, config: dict) -> dict[str, xr.Dataset]:
    """
    Download all seven food-security feature bands from GEE in parallel.
    The six fetch calls are independent HTTP requests; they run concurrently
    via DaskEngine.run_io_parallel (ThreadPoolExecutor) sharing the GEE session.
    Returns a dict keyed by band group name, consumed by cog_export.
    """
    from climate_change.core.dask_engine import DaskEngine
    from climate_change.core.population import fetch_population_count_safe

    scale = config.get("scale", 1000)
    study_start = config.get("start_date", "2018-01-01")
    study_end = config.get("end_date", "2023-12-31")
    lt_start = config.get("lt_baseline_start", LT_BASELINE_START)
    s2_start = config.get("s2_start", "2020-01-01")

    return DaskEngine.run_io_parallel(
        {
            "vci_tci": lambda: fetch_vci_tci(aoi, lt_start, study_start, study_end, scale=scale),
            "rainfall": lambda: fetch_rainfall_anomaly_pct(
                aoi, lt_start, study_start, study_end, scale=scale
            ),
            "ndvi_slope": lambda: fetch_ndvi_slope_img(aoi, lt_start, study_end, scale=scale),
            "mndwi": lambda: fetch_mndwi(aoi, s2_start, study_end, scale=scale),
            "terrain": lambda: fetch_terrain_slope(aoi, scale=scale),
            "landcover": lambda: fetch_landcover(aoi, scale=scale),
            # Raw population COUNT for exposure reporting (population within
            # medium/high food-insecurity risk zones). Best-effort — see
            # core.population.fetch_population_count_safe.
            "population_count": lambda: fetch_population_count_safe(
                aoi, scale=scale, country=config.get("country")
            ),
        }
    )

build_gee_feature_stack

build_gee_feature_stack(aoi: Geometry, config: dict) -> ee.Image

Assemble the 7-band GEE image used for stratified pixel sampling. Band order matches FEATURE_COLS.

Source code in food_security/features.py
def build_gee_feature_stack(aoi: ee.Geometry, config: dict) -> ee.Image:
    """
    Assemble the 7-band GEE image used for stratified pixel sampling.
    Band order matches FEATURE_COLS.
    """
    lt_start = config.get("lt_baseline_start", LT_BASELINE_START)
    study_start = config.get("start_date", "2018-01-01")
    study_end = config.get("end_date", "2023-12-31")
    s2_start = config.get("s2_start", "2020-01-01")
    scale = config.get("scale", 1000)

    def _scale_ndvi(img: ee.Image) -> ee.Image:
        return img.multiply(0.0001).copyProperties(img, ["system:time_start"])

    def _scale_lst(img: ee.Image) -> ee.Image:
        return img.multiply(0.02).subtract(273.15).copyProperties(img, ["system:time_start"])

    def _add_time(img: ee.Image) -> ee.Image:
        t = img.date().difference(ee.Date(lt_start), "year")
        return img.addBands(ee.Image(t).rename("time").float())

    # NDVI long-term min/max for VCI
    modis_lt = (
        ee.ImageCollection("MODIS/061/MOD13A3")
        .filterBounds(aoi)
        .filterDate(lt_start, study_start)
        .select("NDVI")
        .map(_scale_ndvi)
    )
    ndvi_min = modis_lt.min()
    ndvi_max = modis_lt.max()
    ndvi_current = (
        ee.ImageCollection("MODIS/061/MOD13A3")
        .filterBounds(aoi)
        .filterDate(study_start, study_end)
        .select("NDVI")
        .map(_scale_ndvi)
        .mean()
    )
    vci = (
        ndvi_current.subtract(ndvi_min)
        .divide(ndvi_max.subtract(ndvi_min).add(1e-6))
        .multiply(100)
        .rename("vci")
        .clip(aoi)
    )

    # LST long-term baseline for TCI
    lst_lt = (
        ee.ImageCollection("MODIS/061/MOD11A2")
        .filterBounds(aoi)
        .filterDate(lt_start, study_end)
        .select("LST_Day_1km")
        .map(_scale_lst)
    )
    lst_min = lst_lt.reduce(ee.Reducer.percentile([10]))
    lst_max = lst_lt.reduce(ee.Reducer.percentile([90]))
    lst_current = (
        ee.ImageCollection("MODIS/061/MOD11A2")
        .filterBounds(aoi)
        .filterDate(study_start, study_end)
        .select("LST_Day_1km")
        .map(_scale_lst)
        .mean()
    )
    tci = (
        lst_max.subtract(lst_current)
        .divide(lst_max.subtract(lst_min).add(1e-6))
        .multiply(100)
        .rename("tci")
        .clip(aoi)
    )

    # CHIRPS rainfall anomaly %
    chirps = ee.ImageCollection("UCSB-CHG/CHIRPS/DAILY").filterBounds(aoi)
    lt_rain_mean = (
        chirps.filterDate(lt_start, study_start).select("precipitation").mean().multiply(365)
    )
    study_rain_mean = (
        chirps.filterDate(study_start, study_end).select("precipitation").mean().multiply(365)
    )
    rain_anom = (
        study_rain_mean.subtract(lt_rain_mean)
        .divide(lt_rain_mean.add(1e-6))
        .multiply(100)
        .rename("rainfall_anom_pct")
        .clip(aoi)
    )

    # MODIS NDVI linearFit slope
    modis_full = (
        ee.ImageCollection("MODIS/061/MOD13A3")
        .filterBounds(aoi)
        .filterDate(lt_start, study_end)
        .select("NDVI")
        .map(_scale_ndvi)
    )
    trend = modis_full.map(_add_time).select(["time", "NDVI"]).reduce(ee.Reducer.linearFit())
    ndvi_slope = trend.select("scale").rename("ndvi_slope").clip(aoi)

    # Sentinel-2 MNDWI
    def _add_mndwi(img: ee.Image) -> ee.Image:
        return (
            img.normalizedDifference(["B3", "B11"])
            .rename("mndwi")
            .copyProperties(img, ["system:time_start"])
        )

    mndwi = (
        ee.ImageCollection("COPERNICUS/S2_SR_HARMONIZED")
        .filterBounds(aoi)
        .filterDate(s2_start, study_end)
        .filter(ee.Filter.lt("CLOUDY_PIXEL_PERCENTAGE", 30))
        .map(_add_mndwi)
        .median()
        .clip(aoi)
    )

    # SRTM slope
    slope = (
        ee.Terrain.slope(ee.Image("USGS/SRTMGL1_003"))
        .clip(aoi)
        .reproject("EPSG:4326", None, scale)
        .rename("slope_terrain")
    )

    # ESA WorldCover normalised
    land_cover = (
        ee.ImageCollection("ESA/WorldCover/v200")
        .first()
        .select("Map")
        .divide(100)
        .rename("land_cover")
        .clip(aoi)
        .reproject("EPSG:4326", None, scale)
    )

    stack = ee.Image.cat([vci, tci, rain_anom, ndvi_slope, mndwi, slope, land_cover])
    # Exclude water bodies and built-up areas from feature preparation/sampling —
    # food-insecurity risk only applies to agricultural/vegetated land.
    valid_land = exclusion_mask(aoi, [WATER_CLASS, BUILTUP_CLASS])
    return stack.updateMask(valid_land)

fetch_ndvi_monthly_timeseries

fetch_ndvi_monthly_timeseries(aoi: Geometry, start: str, end: str) -> pd.DataFrame

Monthly area-mean MODIS NDVI over [start, end]. Returns DataFrame with DatetimeIndex and column 'ndvi'.

Source code in food_security/features.py
def fetch_ndvi_monthly_timeseries(
    aoi: ee.Geometry,
    start: str,
    end: str,
) -> pd.DataFrame:
    """
    Monthly area-mean MODIS NDVI over [start, end].
    Returns DataFrame with DatetimeIndex and column 'ndvi'.
    """

    def _scale_ndvi(img: ee.Image) -> ee.Image:
        return img.multiply(0.0001).copyProperties(img, ["system:time_start"])

    def _mean_feat(img: ee.Image) -> ee.Feature:
        v = img.reduceRegion(ee.Reducer.mean(), aoi, 1000, maxPixels=int(1e9)).get("NDVI")
        return ee.Feature(None, {"date": img.date().format("YYYY-MM"), "ndvi": v})

    modis = (
        ee.ImageCollection("MODIS/061/MOD13A3")
        .filterBounds(aoi)
        .filterDate(start, end)
        .select("NDVI")
        .map(_scale_ndvi)
    )
    records = cast(dict, ee.FeatureCollection(modis.map(_mean_feat)).getInfo())["features"]
    return (
        pd.DataFrame([f["properties"] for f in records])
        .dropna()
        .groupby("date")[["ndvi"]]
        .mean()
        .sort_index()
    )

fetch_monthly_rainfall

fetch_monthly_rainfall(aoi: Geometry, start: str, end: str) -> pd.DataFrame

Monthly area-mean CHIRPS rainfall sum over [start, end]. Returns DataFrame with DatetimeIndex and column 'rain_mm'.

Source code in food_security/features.py
def fetch_monthly_rainfall(
    aoi: ee.Geometry,
    start: str,
    end: str,
) -> pd.DataFrame:
    """
    Monthly area-mean CHIRPS rainfall sum over [start, end].
    Returns DataFrame with DatetimeIndex and column 'rain_mm'.
    """
    n_months = int(cast(int, ee.Date(end).difference(ee.Date(start), "month").round().getInfo()))
    months = ee.List.sequence(0, n_months - 1)

    def _monthly_sum(offset: ee.Number) -> ee.Feature:
        s = ee.Date(start).advance(offset, "month")
        e = s.advance(1, "month")
        v = (
            ee.ImageCollection("UCSB-CHG/CHIRPS/DAILY")
            .filterBounds(aoi)
            .filterDate(s, e)
            .select("precipitation")
            .sum()
            .reduceRegion(ee.Reducer.mean(), aoi, 1000, maxPixels=int(1e9))
            .get("precipitation", 0)
        )
        return ee.Feature(None, {"date": s.format("YYYY-MM"), "rain_mm": v})

    records = cast(dict, ee.FeatureCollection(months.map(_monthly_sum)).getInfo())["features"]
    return (
        pd.DataFrame([f["properties"] for f in records])
        .dropna()
        .groupby("date")[["rain_mm"]]
        .mean()
        .sort_index()
    )

sample_training_data

sample_training_data(feature_stack: Image, aoi: Geometry, n_pixels: int = 3000, scale: int = 1000, seed: int = 42) -> pd.DataFrame

Sample n_pixels from the GEE feature stack and assign 3-class food insecurity labels. Labels: fixed absolute thresholds on the composite food stress score (0-100 scale). 0 = Low Risk (score < 33) 1 = Medium Risk (33 <= score < 66) 2 = High Risk (score >= 66) Returns DataFrame with FEATURE_COLS + ['food_score', 'label'].

Source code in food_security/features.py
def sample_training_data(
    feature_stack: ee.Image,
    aoi: ee.Geometry,
    n_pixels: int = 3000,
    scale: int = 1000,
    seed: int = 42,
) -> pd.DataFrame:
    """
    Sample n_pixels from the GEE feature stack and assign 3-class food insecurity labels.
    Labels: fixed absolute thresholds on the composite food stress score (0-100 scale).
      0 = Low Risk    (score < 33)
      1 = Medium Risk (33 <= score < 66)
      2 = High Risk   (score >= 66)
    Returns DataFrame with FEATURE_COLS + ['food_score', 'label'].
    """
    samples = feature_stack.sample(
        region=aoi,
        scale=scale,
        numPixels=n_pixels,
        seed=seed,
        geometries=False,
        dropNulls=False,
    )
    records = cast(dict, samples.getInfo())["features"]
    df = (
        pd.DataFrame([f["properties"] for f in records])
        .dropna(subset=FEATURE_COLS)[FEATURE_COLS]
        .reset_index(drop=True)
    )

    scores = _compute_food_score(df)
    labels = np.zeros(len(df), dtype=np.intp)
    labels[scores >= RISK_SCORE_THRESHOLDS[0]] = 1
    labels[scores >= RISK_SCORE_THRESHOLDS[1]] = 2

    df["food_score"] = scores
    df["label"] = labels
    return df

align_datasets

align_datasets(datasets: dict[str, Dataset], ref_key: str = 'vci_tci', method_continuous: InterpOptions = 'linear', method_categorical: InterpOptions = 'nearest') -> dict[str, xr.Dataset]

Interpolate all datasets onto the VCI/TCI reference grid. Land cover is treated as categorical (nearest-neighbour). population_count is passed through entirely unchanged, at its own native ~100 m resolution — it is never interpolated onto this (coarser) reference grid. Summing population within a risk zone requires the opposite: upsampling the (categorical) classification onto population's native grid — see core.population.population_exposure, which cog_export.py uses for this instead of relying on align_datasets. Interpolating population here would silently corrupt totals (verified against live GEE data — see core.population.fetch_population_count's docstring).

Each dataset is chunked before interpolation so that xarray produces Dask-backed lazy arrays; dask.compute() then materialises them all in a single parallel scheduler pass.

Source code in food_security/features.py
def align_datasets(
    datasets: dict[str, xr.Dataset],
    ref_key: str = "vci_tci",
    method_continuous: InterpOptions = "linear",
    method_categorical: InterpOptions = "nearest",
) -> dict[str, xr.Dataset]:
    """
    Interpolate all datasets onto the VCI/TCI reference grid. Land cover is
    treated as categorical (nearest-neighbour). population_count is passed
    through entirely unchanged, at its own native ~100 m resolution — it is
    never interpolated onto this (coarser) reference grid. Summing
    population within a risk zone requires the opposite: upsampling the
    (categorical) classification onto population's native grid — see
    core.population.population_exposure, which cog_export.py uses for this
    instead of relying on align_datasets. Interpolating population here
    would silently corrupt totals (verified against live GEE data — see
    core.population.fetch_population_count's docstring).

    Each dataset is chunked before interpolation so that xarray produces
    Dask-backed lazy arrays; dask.compute() then materialises them all
    in a single parallel scheduler pass.
    """
    from dask.base import compute

    _CHUNK = {"lat": 256, "lon": 256}
    ref = datasets[ref_key]

    lazy: dict[str, xr.Dataset] = {ref_key: ref}
    for key, ds in datasets.items():
        if key == ref_key or ds is None or key == "population_count":
            continue
        method = method_categorical if key == "landcover" else method_continuous
        lazy[key] = ds.chunk(_CHUNK).interp(lat=ref.lat, lon=ref.lon, method=method)

    keys = list(lazy)
    computed = compute(*[lazy[k] for k in keys])
    result = dict(zip(keys, computed))
    pop_ds = datasets.get("population_count")
    if pop_ds is not None:
        result["population_count"] = pop_ds
    return result

Model

model

FoodSecurityModel

FoodSecurityModel()

Orchestrates the full ML pipeline for a single food security analysis. config['model_type'] selects which model(s) are trained and which predictions drive the primary risk distribution: "rf" — Random Forest only (default) "xgboost" — XGBoost only "ensemble" — both, combined as mean softmax probabilities of RF + XGBoost Only the selected model type is trained. Trained models and scaler are stored on self for use by cog_export (whichever model wasn't needed stays None).

Source code in food_security/model.py
def __init__(self) -> None:
    self.rf: RandomForestClassifier | None = None
    self.xgb: Booster | None = None
    self.scaler: StandardScaler | None = None

predict

predict(df: DataFrame, ndvi_df: DataFrame | None = None, rain_df: DataFrame | None = None, config: dict | None = None) -> dict
Parameters

df : DataFrame with FEATURE_COLS + ['food_score', 'label'] ndvi_df : monthly NDVI DataFrame (date index, 'ndvi' column) rain_df : monthly rainfall DataFrame (date index, 'rain_mm' column) config : optional dict; reads 'model_type' (default 'rf')

Returns

dict with keys: stats, charts

Source code in food_security/model.py
def predict(
    self,
    df: pd.DataFrame,
    ndvi_df: pd.DataFrame | None = None,
    rain_df: pd.DataFrame | None = None,
    config: dict | None = None,
) -> dict:
    """
    Parameters
    ----------
    df       : DataFrame with FEATURE_COLS + ['food_score', 'label']
    ndvi_df  : monthly NDVI DataFrame (date index, 'ndvi' column)
    rain_df  : monthly rainfall DataFrame (date index, 'rain_mm' column)
    config   : optional dict; reads 'model_type' (default 'rf')

    Returns
    -------
    dict with keys: stats, charts
    """
    cfg = config or {}
    model_type = cfg.get("model_type", "rf")
    if model_type not in VALID_MODEL_TYPES:
        raise ValueError(f"model_type must be one of {VALID_MODEL_TYPES}, got '{model_type}'")

    _ndvi = ndvi_df if ndvi_df is not None else pd.DataFrame()
    _rain = rain_df if rain_df is not None else pd.DataFrame()

    X = df[FEATURE_COLS].fillna(df[FEATURE_COLS].median()).to_numpy(dtype=np.float64)
    y = df["label"].to_numpy(dtype=np.intp)

    X_train, X_test, y_train, y_test = train_test_split(
        X, y, test_size=0.20, random_state=42, stratify=_safe_stratify(y)
    )
    self.scaler = StandardScaler()
    X_train_s = self.scaler.fit_transform(X_train)
    X_test_s = self.scaler.transform(X_test)

    # Only train what this selection needs — "ensemble" needs both, "rf"
    # and "xgboost" need only the one selected. When both are needed,
    # train them concurrently on the Dask cluster if it's running,
    # otherwise fall back to sequential execution (e.g. during unit tests).
    from climate_change.core.dask_engine import DaskEngine

    need_rf = model_type in ("rf", "ensemble")
    need_xgb = model_type in ("xgboost", "ensemble")
    rf_meta: dict = {"cv_f1_mean": None, "cv_f1_std": None}
    xgb_meta: dict = {"cv_f1_mean": None, "cv_f1_std": None}

    client = DaskEngine.get_client_if_running() if (need_rf and need_xgb) else None
    if client is not None:
        f_rf = client.submit(train_rf, X_train_s, y_train, pure=False)
        f_xgb = client.submit(train_xgb, X_train_s, y_train, pure=False)
        (self.rf, rf_meta), (self.xgb, xgb_meta) = cast(list, client.gather([f_rf, f_xgb]))
    else:
        if need_rf:
            self.rf, rf_meta = train_rf(X_train_s, y_train)
        if need_xgb:
            self.xgb, xgb_meta = train_xgb(X_train_s, y_train)

    eval_result = evaluate_models(self.rf, self.xgb, X_test_s, y_test)
    if model_type == "rf":
        assert self.rf is not None
        shap_payload = compute_shap_importance(self.rf, X_test_s)
    else:
        assert self.xgb is not None
        shap_payload = compute_shap_importance(self.xgb, X_test_s)

    # Compute VHI summary scalars from the feature DataFrame
    vci_mean = float(cast(float, df["vci"].mean())) if "vci" in df.columns else 0.0
    tci_mean = float(cast(float, df["tci"].mean())) if "tci" in df.columns else 0.0
    vhi_mean = 0.5 * vci_mean + 0.5 * tci_mean

    charts = build_food_security_charts(
        eval_result,
        shap_payload,
        _ndvi,
        _rain,
        vci_mean,
        tci_mean,
        vhi_mean,
        model_type,
    )

    # Risk distribution on all pixels using selected model
    X_all_s = self.scaler.transform(
        df[FEATURE_COLS].fillna(df[FEATURE_COLS].median()).to_numpy(dtype=np.float64)
    )
    _KEY = {"rf": "rf", "xgboost": "xgb", "ensemble": "ensemble"}
    result_key = _KEY.get(model_type, "rf")

    if model_type == "rf":
        assert self.rf is not None
        all_preds = self.rf.predict(X_all_s).astype(int)
    elif model_type == "xgboost":
        assert self.xgb is not None
        all_preds = np.argmax(self.xgb.predict(DMatrix(X_all_s)), axis=1).astype(int)
    else:
        assert self.rf is not None
        assert self.xgb is not None
        proba_rf = pad_rf_proba(self.rf, self.rf.predict_proba(X_all_s))
        proba = (proba_rf + self.xgb.predict(DMatrix(X_all_s))) / 2.0
        all_preds = np.argmax(proba, axis=1).astype(int)

    n_total = len(all_preds)
    counts = np.array([(all_preds == i).sum() for i in range(3)], dtype=np.float64)
    risk_pct = (counts / n_total * 100).round(1)
    high_risk_pct = float(risk_pct[2])

    stats = {
        "model_type": model_type,
        "n_pixels_sampled": int(len(df)),
        "rf_cv_f1": round(rf_meta["cv_f1_mean"], 4)
        if rf_meta["cv_f1_mean"] is not None
        else None,
        "rf_f1": eval_result["rf"]["f1"] if "rf" in eval_result else None,
        "rf_accuracy": eval_result["rf"]["accuracy"] if "rf" in eval_result else None,
        "xgb_cv_f1": round(xgb_meta["cv_f1_mean"], 4)
        if xgb_meta["cv_f1_mean"] is not None
        else None,
        "xgb_f1": eval_result["xgb"]["f1"] if "xgb" in eval_result else None,
        "xgb_accuracy": eval_result["xgb"]["accuracy"] if "xgb" in eval_result else None,
        "ensemble_f1": eval_result["ensemble"]["f1"] if "ensemble" in eval_result else None,
        "selected_f1": eval_result[result_key]["f1"],
        "high_risk_pct": round(high_risk_pct, 1),
        "top_driver": shap_payload["features"][0],
        "vci_mean": round(vci_mean, 1),
        "tci_mean": round(tci_mean, 1),
        "vhi_mean": round(vhi_mean, 1),
    }

    _FOOD_CLASS_NAMES = ["Low Risk", "Medium Risk", "High Risk"]
    if "lon" in df.columns and "lat" in df.columns:
        _sample_points = [
            {
                "lon": round(cast(float, df["lon"].iat[i]), 5),
                "lat": round(cast(float, df["lat"].iat[i]), 5),
                "risk_class": _FOOD_CLASS_NAMES[int(all_preds[i])],
            }
            for i in range(len(df))
        ]
    else:
        _sample_points = []
    return {"stats": stats, "charts": charts, "_sample_points": _sample_points}

train_rf

train_rf(X_train: ndarray, y_train: ndarray, cv_folds: int = 5) -> tuple[RandomForestClassifier, dict]

Fit a balanced Random Forest classifier. Returns (model, metadata).

Source code in food_security/model.py
def train_rf(
    X_train: np.ndarray,
    y_train: np.ndarray,
    cv_folds: int = 5,
) -> tuple[RandomForestClassifier, dict]:
    """Fit a balanced Random Forest classifier. Returns (model, metadata)."""
    clf = RandomForestClassifier(
        n_estimators=RF_N_ESTIMATORS,
        max_depth=RF_MAX_DEPTH,
        min_samples_leaf=RF_MIN_SAMPLES_LEAF,
        class_weight="balanced",
        random_state=42,
        n_jobs=-1,
    )
    clf.fit(X_train, y_train)
    folds = _safe_cv_folds(y_train, cv_folds)
    if folds is None:
        return clf, {"cv_f1_mean": None, "cv_f1_std": None}
    cv = StratifiedKFold(n_splits=folds, shuffle=True, random_state=42)
    cv_f1 = cross_val_score(clf, X_train, y_train, cv=cv, scoring="f1_macro")
    return clf, {"cv_f1_mean": float(cv_f1.mean()), "cv_f1_std": float(cv_f1.std())}

train_xgb

train_xgb(X_train: ndarray, y_train: ndarray, cv_folds: int = 5) -> tuple[Booster, dict]

Fit XGBoost multi-class classifier via the low-level Booster API with sample weights.

Uses xgboost.train() rather than the XGBClassifier sklearn wrapper because the wrapper requires every class in FOOD_CLASSES to appear in y_train (it infers num_class from np.unique(y_train) and rejects gaps). Some AOIs/time windows genuinely have zero samples of a given risk class, which is a legitimate state, not invalid input, so class count is fixed via XGB_PARAMS instead.

Returns (booster, metadata).

Source code in food_security/model.py
def train_xgb(
    X_train: np.ndarray,
    y_train: np.ndarray,
    cv_folds: int = 5,
) -> tuple[Booster, dict]:
    """
    Fit XGBoost multi-class classifier via the low-level Booster API with sample weights.

    Uses xgboost.train() rather than the XGBClassifier sklearn wrapper because the
    wrapper requires every class in FOOD_CLASSES to appear in y_train (it infers
    num_class from np.unique(y_train) and rejects gaps). Some AOIs/time windows
    genuinely have zero samples of a given risk class, which is a legitimate state,
    not invalid input, so class count is fixed via XGB_PARAMS instead.

    Returns (booster, metadata).
    """
    sw = compute_sample_weight("balanced", y_train)
    dtrain = DMatrix(X_train, label=y_train, weight=sw)
    booster = xgb_train(XGB_PARAMS, dtrain, num_boost_round=XGB_N_ESTIMATORS)

    folds = _safe_cv_folds(y_train, cv_folds)
    if folds is None:
        return booster, {"cv_f1_mean": None, "cv_f1_std": None}
    cv = StratifiedKFold(n_splits=folds, shuffle=True, random_state=42)
    cv_f1_scores = []
    for train_idx, val_idx in cv.split(X_train, y_train):
        fold_booster = xgb_train(
            XGB_PARAMS,
            DMatrix(X_train[train_idx], label=y_train[train_idx]),
            num_boost_round=XGB_N_ESTIMATORS,
        )
        fold_pred = np.argmax(fold_booster.predict(DMatrix(X_train[val_idx])), axis=1)
        cv_f1_scores.append(f1_score(y_train[val_idx], fold_pred, average="macro"))

    cv_f1 = np.array(cv_f1_scores)
    return booster, {"cv_f1_mean": float(cv_f1.mean()), "cv_f1_std": float(cv_f1.std())}

pad_rf_proba

pad_rf_proba(rf: RandomForestClassifier, proba: ndarray) -> np.ndarray

Expand RF's predict_proba output to the full FOOD_CLASSES width.

RandomForestClassifier only emits a column per class it actually saw during fit (via rf.classes_), whereas train_xgb's Booster always outputs len(FOOD_CLASSES) columns regardless of what the AOI's data contained. Without padding, an AOI missing one risk class produces mismatched shapes (e.g. (n, 2) vs (n, 3)) the moment RF and XGBoost probabilities are combined for the ensemble.

Source code in food_security/model.py
def pad_rf_proba(rf: RandomForestClassifier, proba: np.ndarray) -> np.ndarray:
    """
    Expand RF's predict_proba output to the full FOOD_CLASSES width.

    RandomForestClassifier only emits a column per class it actually saw
    during fit (via rf.classes_), whereas train_xgb's Booster always outputs
    len(FOOD_CLASSES) columns regardless of what the AOI's data contained.
    Without padding, an AOI missing one risk class produces mismatched shapes
    (e.g. (n, 2) vs (n, 3)) the moment RF and XGBoost probabilities are
    combined for the ensemble.
    """
    n_classes = len(FOOD_CLASSES)
    if proba.shape[1] == n_classes:
        return proba
    padded = np.zeros((proba.shape[0], n_classes), dtype=proba.dtype)
    padded[:, rf.classes_.astype(int)] = proba
    return padded

evaluate_models

evaluate_models(rf: RandomForestClassifier | None, xgb: Booster | None, X_test: ndarray, y_test: ndarray) -> dict

Evaluate whichever of RF/XGBoost were actually trained (both, for an "ensemble" selection; exactly one otherwise) on the held-out test set, plus their mean-proba ensemble when both are present.

Source code in food_security/model.py
def evaluate_models(
    rf: RandomForestClassifier | None,
    xgb: Booster | None,
    X_test: np.ndarray,
    y_test: np.ndarray,
) -> dict:
    """
    Evaluate whichever of RF/XGBoost were actually trained (both, for an
    "ensemble" selection; exactly one otherwise) on the held-out test set,
    plus their mean-proba ensemble when both are present.
    """

    def _metrics(pred: np.ndarray, label: str) -> dict:
        return {
            "label": label,
            "f1": round(float(f1_score(y_test, pred, average="macro")), 4),
            "accuracy": round(float(accuracy_score(y_test, pred)), 4),
            "predictions": pred.tolist(),
        }

    proba_rf = pad_rf_proba(rf, rf.predict_proba(X_test)) if rf is not None else None
    proba_xgb = xgb.predict(DMatrix(X_test)) if xgb is not None else None

    result: dict = {"actuals": y_test.tolist()}
    if proba_rf is not None:
        result["rf"] = _metrics(np.argmax(proba_rf, axis=1).astype(int), "Random Forest")
    if proba_xgb is not None:
        result["xgb"] = _metrics(np.argmax(proba_xgb, axis=1).astype(int), "XGBoost")
    if proba_rf is not None and proba_xgb is not None:
        ens_pred = np.argmax((proba_rf + proba_xgb) / 2.0, axis=1).astype(int)
        result["ensemble"] = _metrics(ens_pred, "Ensemble (mean proba)")
    return result

compute_shap_importance

compute_shap_importance(model: RandomForestClassifier | Booster, X_test: ndarray) -> dict

TreeExplainer SHAP values for whichever model was actually selected/trained — see core.shap_utils for the axis-handling details.

Source code in food_security/model.py
def compute_shap_importance(model: RandomForestClassifier | Booster, X_test: np.ndarray) -> dict:
    """TreeExplainer SHAP values for whichever model was actually selected/trained — see core.shap_utils for the axis-handling details."""
    from climate_change.core.shap_utils import compute_shap_importance as _compute_shap_importance

    X_df = pd.DataFrame(X_test, columns=FEATURE_COLS)
    return _compute_shap_importance(model, X_df, FEATURE_COLS)

build_food_security_charts

build_food_security_charts(eval_result: dict, shap_payload: dict, ndvi_df: DataFrame, rain_df: DataFrame, vci_mean: float, tci_mean: float, vhi_mean: float, model_type: str = 'rf') -> dict

Assemble frontend-ready chart payloads for the food security module.

Source code in food_security/model.py
def build_food_security_charts(
    eval_result: dict,
    shap_payload: dict,
    ndvi_df: pd.DataFrame,
    rain_df: pd.DataFrame,
    vci_mean: float,
    tci_mean: float,
    vhi_mean: float,
    model_type: str = "rf",
) -> dict:
    """Assemble frontend-ready chart payloads for the food security module."""
    _KEY = {"rf": "rf", "xgboost": "xgb", "ensemble": "ensemble"}
    result_key = _KEY.get(model_type, "rf")
    predictions = np.array(eval_result[result_key]["predictions"])

    n_total = len(predictions)
    counts = np.array([(predictions == i).sum() for i in range(3)], dtype=np.float64)
    risk_pct = (counts / n_total * 100).round(1)

    risk_dist = {
        "labels": FOOD_CLASSES,
        "data": risk_pct.tolist(),
        "colors": FOOD_COLORS,
    }

    # Time series: NDVI + monthly rainfall aligned to NDVI index
    ts_labels = ndvi_df.index.tolist() if not ndvi_df.empty else []
    ts_datasets = []
    if not ndvi_df.empty:
        ts_datasets.append(
            {
                "label": "NDVI",
                "data": ndvi_df["ndvi"].round(4).tolist(),
                "color": "#27AE60",
            }
        )
    if not rain_df.empty and ts_labels:
        aligned_rain = rain_df["rain_mm"].reindex(ts_labels).round(1).tolist()
        ts_datasets.append(
            {
                "label": "Monthly rain (mm)",
                "data": aligned_rain,
                "color": "#2980B9",
            }
        )

    return {
        "riskDist": risk_dist,
        "timeSeries": {"labels": ts_labels, "datasets": ts_datasets},
        "shap": shap_payload,
        "indices": {
            "vci_mean": round(vci_mean, 1),
            "tci_mean": round(tci_mean, 1),
            "vhi_mean": round(vhi_mean, 1),
        },
        "model_performance": {
            **{
                key: {"f1": eval_result[key]["f1"], "accuracy": eval_result[key]["accuracy"]}
                for key in ("rf", "xgb", "ensemble")
                if key in eval_result
            },
            "selected": model_type,
        },
    }

COG export

cog_export

predict_food_security_grid

predict_food_security_grid(rf_model: RandomForestClassifier | None, xgb_model: Booster | None, scaler: StandardScaler, datasets: dict[str, Dataset], model_type: str = 'rf', aoi_geojson: dict | None = None) -> dict

Run full-grid inference over the aligned feature datasets and classify every valid pixel into Low/Medium/High food-insecurity risk.

Returns a dict with the classified risk_grid (int-encoded per RISK_INT, 0 = nodata), the raster transform/crs/shape, and per-class pixel counts/percentages. Water bodies and built-up land cover are masked out before classification. Used both by export_food_security_cog (COG write) and the use case's GeoJSON conversion.

Source code in food_security/cog_export.py
def predict_food_security_grid(
    rf_model: RandomForestClassifier | None,
    xgb_model: Booster | None,
    scaler: StandardScaler,
    datasets: dict[str, xr.Dataset],
    model_type: str = "rf",
    aoi_geojson: dict | None = None,
) -> dict:
    """
    Run full-grid inference over the aligned feature datasets and classify
    every valid pixel into Low/Medium/High food-insecurity risk.

    Returns a dict with the classified `risk_grid` (int-encoded per
    `RISK_INT`, 0 = nodata), the raster `transform`/`crs`/shape, and
    per-class pixel counts/percentages.
    Water bodies and built-up land cover are masked out before classification.
    Used both by `export_food_security_cog` (COG write) and the use case's
    GeoJSON conversion.
    """
    aligned = align_datasets(datasets, ref_key="vci_tci")
    vci_tci_ds = aligned["vci_tci"]
    ref_lat = vci_tci_ds.lat
    ref_lon = vci_tci_ds.lon
    n_lat, n_lon = len(ref_lat), len(ref_lon)

    vci = vci_tci_ds["vci"].values.astype(np.float32)
    tci = vci_tci_ds["tci"].values.astype(np.float32)
    rain_anom = aligned["rainfall"]["rainfall_anom_pct"].values.astype(np.float32)
    ndvi_slope = aligned["ndvi_slope"]["ndvi_slope"].values.astype(np.float32)
    mndwi = aligned["mndwi"]["mndwi"].values.astype(np.float32)
    slope = aligned["terrain"]["slope_terrain"].values.astype(np.float32)
    lc = aligned["landcover"]["land_cover"].values.astype(np.float32)

    bands = np.stack([vci, tci, rain_anom, ndvi_slope, mndwi, slope, lc], axis=0)
    X_full = bands.reshape(len(FEATURE_COLS), -1).T
    valid_mask = ~np.isnan(X_full).any(axis=1)
    # Water bodies and built-up areas are outside the food-insecurity domain.
    water_or_builtup = np.isclose(lc.reshape(-1), _WATER_LAND_COVER) | np.isclose(
        lc.reshape(-1), _BUILTUP_LAND_COVER
    )
    valid_mask &= ~water_or_builtup

    transform = vci_tci_ds["vci"].rio.transform()
    crs = vci_tci_ds["vci"].rio.crs or "EPSG:4326"
    geometries = _aoi_geometries(aoi_geojson)
    if geometries:
        inside_aoi = geometry_mask(
            geometries,
            out_shape=(n_lat, n_lon),
            transform=transform,
            invert=True,
        ).reshape(-1)
        valid_mask = valid_mask & inside_aoi

    risk_grid: np.ndarray = np.zeros(n_lat * n_lon, dtype=np.uint8)
    if valid_mask.any():
        X_valid = scaler.transform(X_full[valid_mask])
        if model_type == "rf":
            assert rf_model is not None
            proba = rf_model.predict_proba(X_valid)
        elif model_type == "xgboost":
            assert xgb_model is not None
            proba = xgb_model.predict(DMatrix(X_valid))
        else:
            assert rf_model is not None
            assert xgb_model is not None
            proba_rf = pad_rf_proba(rf_model, rf_model.predict_proba(X_valid))
            proba = (proba_rf + xgb_model.predict(DMatrix(X_valid))) / 2.0
        pred_classes = np.argmax(proba, axis=1).astype(np.uint8) + 1
        risk_grid[valid_mask] = pred_classes

    risk_grid = risk_grid.reshape(n_lat, n_lon)
    valid_values = risk_grid[risk_grid > 0]
    counts = np.array(
        [(valid_values == class_id).sum() for class_id in (1, 2, 3)],
        dtype=np.float64,
    )
    total = int(valid_values.size)
    percentages = (counts / total * 100).round(1).tolist() if total else [0.0, 0.0, 0.0]

    population_by_class: dict[str, float] = {}
    total_population: float | None = None
    pop_ds = aligned.get("population_count")
    if pop_ds is not None:
        by_int = population_exposure(risk_grid, transform, pop_ds, classes=[1, 2, 3])
        population_by_class = {RISK_INT[v]: by_int[str(v)] for v in (1, 2, 3)}
        total_population = round(sum(population_by_class.values()), 1)

    return {
        "risk_grid": risk_grid,
        "transform": transform,
        "crs": crs,
        "height": n_lat,
        "width": n_lon,
        "counts": counts.astype(int).tolist(),
        "percentages": percentages,
        "valid_pixel_count": total,
        "population_by_class": population_by_class,
        "total_population": total_population,
    }

export_food_security_cog

export_food_security_cog(rf_model: RandomForestClassifier | None, xgb_model: Booster | None, scaler: StandardScaler, datasets: dict[str, Dataset], output_dir: str, prefix: str, model_type: str = 'rf', aoi_geojson: dict | None = None) -> dict[str, str]

Build the full-AOI feature matrix from in-memory xarray Datasets, run inference with the selected model, classify into 3 food insecurity risk classes, and write a Cloud Optimised GeoTIFF.

Parameters

rf_model : trained RandomForestClassifier xgb_model : trained XGBoost Booster scaler : fitted StandardScaler (from FoodSecurityModel.scaler) datasets : dict returned by build_feature_datasets output_dir : local directory for COG output prefix : file prefix model_type : "rf" | "xgboost" | "ensemble" aoi_geojson : optional AOI polygon used to mask pixels outside the true AOI

Returns

dict with key 'food_security_risk' → path string

Source code in food_security/cog_export.py
def export_food_security_cog(
    rf_model: RandomForestClassifier | None,
    xgb_model: Booster | None,
    scaler: StandardScaler,
    datasets: dict[str, xr.Dataset],
    output_dir: str,
    prefix: str,
    model_type: str = "rf",
    aoi_geojson: dict | None = None,
) -> dict[str, str]:
    """
    Build the full-AOI feature matrix from in-memory xarray Datasets,
    run inference with the selected model, classify into 3 food insecurity risk classes,
    and write a Cloud Optimised GeoTIFF.

    Parameters
    ----------
    rf_model    : trained RandomForestClassifier
    xgb_model   : trained XGBoost Booster
    scaler      : fitted StandardScaler (from FoodSecurityModel.scaler)
    datasets    : dict returned by build_feature_datasets
    output_dir  : local directory for COG output
    prefix      : file prefix
    model_type  : "rf" | "xgboost" | "ensemble"
    aoi_geojson : optional AOI polygon used to mask pixels outside the true AOI

    Returns
    -------
    dict with key 'food_security_risk' → path string
    """
    output_path = Path(output_dir)
    output_path.mkdir(parents=True, exist_ok=True)
    cog_path = output_path / f"{prefix}_food_security_risk.tif"

    grid = predict_food_security_grid(
        rf_model,
        xgb_model,
        scaler,
        datasets,
        model_type=model_type,
        aoi_geojson=aoi_geojson,
    )
    risk_grid = grid["risk_grid"]

    with rasterio.open(
        cog_path,
        "w",
        driver="COG",
        height=grid["height"],
        width=grid["width"],
        count=1,
        dtype=np.uint8,
        crs=grid["crs"],
        transform=grid["transform"],
        nodata=0,
        compress="lzw",
    ) as dst:
        dst.write(risk_grid, 1)
        dst.update_tags(
            RISK_CLASSES="0=nodata,1=Low Risk,2=Medium Risk,3=High Risk",
            MODEL_TYPE=model_type,
        )

    _log.info("COG exported: %s  (%d×%d px)", cog_path, grid["height"], grid["width"])
    return {"food_security_risk": str(cog_path)}