Precipitation Assessment API Reference

For land-system classification, see the Classification API Reference.

savana.rainfall.pipeline

RainfallAssessment: the chainable orchestrator tying every savana.rainfall stage together, mirroring :class:savana.pipeline.SavanaClassifier's builder pattern.

Every constructor argument is a default that can be overridden — a different product catalogue, a different gauge network, a different zone scheme, or different application weights all work the same way they do calling the individual stage functions directly. This class is a convenience, not a new capability.

RainfallAssessment

Chainable orchestrator for a precipitation product assessment.

Example (defaults — reproduces the WA study)::

ra = (
    RainfallAssessment()
    .get_observations(source="ee_asset")
    .ingest(start="2001-01-01", end="2020-12-31")
    .extract()
    .assign_zones()
    .validate()
    .score()
)
print(ra.summarize())
ra.export_workbook("decision_tool.xlsx")

Example (a different station network, subset of products)::

ra = (
    RainfallAssessment(
        stations=[(-1.5, 12.4), (2.1, 6.5)],  # or a DataFrame, a
                                                # .geojson/.csv path,
                                                # or a single (lon, lat)
        products={"CHIRPS": config.DEFAULT_PRODUCTS["CHIRPS"],
                  "GPM_IMERG": config.DEFAULT_PRODUCTS["GPM_IMERG"]},
    )
    .get_observations(source="download")
    .ingest(start="2015-01-01", end="2023-12-31")
    .extract()
    .validate()   # no assign_zones() call -> pooled validation
    .score()
)

Source code in savana/rainfall/pipeline.py

class RainfallAssessment:
    """Chainable orchestrator for a precipitation product assessment.

    Example (defaults — reproduces the WA study)::

        ra = (
            RainfallAssessment()
            .get_observations(source="ee_asset")
            .ingest(start="2001-01-01", end="2020-12-31")
            .extract()
            .assign_zones()
            .validate()
            .score()
        )
        print(ra.summarize())
        ra.export_workbook("decision_tool.xlsx")

    Example (a different station network, subset of products)::

        ra = (
            RainfallAssessment(
                stations=[(-1.5, 12.4), (2.1, 6.5)],  # or a DataFrame, a
                                                        # .geojson/.csv path,
                                                        # or a single (lon, lat)
                products={"CHIRPS": config.DEFAULT_PRODUCTS["CHIRPS"],
                          "GPM_IMERG": config.DEFAULT_PRODUCTS["GPM_IMERG"]},
            )
            .get_observations(source="download")
            .ingest(start="2015-01-01", end="2023-12-31")
            .extract()
            .validate()   # no assign_zones() call -> pooled validation
            .score()
        )
    """

    def __init__(
        self,
        products: dict | None = None,
        stations=None,
        zones_gdf=None,
        zones_fc=None,
        app_weights: dict | None = None,
        zone_notes: dict | None = None,
        rain_threshold: float | None = None,
        cache_dir=None,
        ee_project: str | None = None,
    ):
        from . import stations as _stations

        self.products = products if products is not None else config.DEFAULT_PRODUCTS
        # Accepts anything load_stations_any() accepts: a DataFrame, a
        # .geojson/.csv path, a list of (lon, lat)/(id, lon, lat)/dicts,
        # or a single (lon, lat) tuple/list. Always resolves to a real
        # stations_df immediately (never left as None), defaulting to
        # the WA 16 stations if stations=None.
        self.stations_df = _stations.load_stations_any(stations)
        self.zones_gdf = zones_gdf
        self.zones_fc = zones_fc
        self.app_weights = (
            app_weights if app_weights is not None else config.DEFAULT_APP_WEIGHTS
        )
        self.zone_notes = (
            zone_notes if zone_notes is not None else config.DEFAULT_ZONE_NOTES
        )
        self.rain_threshold = (
            rain_threshold
            if rain_threshold is not None
            else config.DEFAULT_RAIN_THRESHOLD_MM_DAY
        )
        self.cache_dir = cache_dir
        self.ee_project = ee_project

        # populated as stages run
        self.start = None
        self.end = None
        self.obs_df = None
        self.products_ic = None
        self.sim_df = None
        self.merged_df = None
        self.validation_by_zone_df = None
        self.validation_overall_df = None
        self.ranking_df = None
        self.threshold_df = None
        self.scores_df = None
        self._facts = None

    def _ensure_ee(self):
        """Initialize Earth Engine with this instance's ``ee_project``,
        called lazily right before any EE-touching operation — not at
        construction time, so a purely offline run (csv observations,
        no ``.ingest()``/``.preview_*`` calls) never
        prompts for EE auth at all. Safe to call repeatedly; a no-op
        after the first successful initialize() in this process (Earth
        Engine's session state is global, not per-object).
        """
        from .. import ee_init

        ee_init.initialize(project=self.ee_project)

    # ────────────────────────────────────────────────────
    # Stages
    # ────────────────────────────────────────────────────

    def get_observations(self, source: str = "download", **kwargs):
        from . import stations as _stations

        if source == "ee_asset":
            self._ensure_ee()

        stations_df = (
            self.stations_df
            if self.stations_df is not None
            else config.default_stations_wa()
        )
        self.stations_df = stations_df

        # If .ingest() already ran and set a date range, and the caller
        # didn't explicitly pass their own start_year/end_year, reuse
        # it -- otherwise "download" silently defaults to the
        # full 2001-2020 range regardless of what .ingest() was told,
        # which is surprising and easy to miss. Only applies to sources
        # that actually take a year range ("csv"/"ee_asset" don't).
        if (
            source == "download"
            and "start_year" not in kwargs
            and "end_year" not in kwargs
            and self.start is not None
            and self.end is not None
        ):
            kwargs["start_year"] = int(self.start[:4])
            kwargs["end_year"] = int(self.end[:4])
            print(
                f"  Using date range from .ingest(): "
                f"{kwargs['start_year']}-{kwargs['end_year']} "
                f"(pass start_year=/end_year= explicitly to override)"
            )

        self.obs_df = _stations.get_observations(stations_df, source=source, **kwargs)
        return self

    def ingest(self, start: str, end: str, roi=None):
        from . import ingestion

        self._ensure_ee()
        self.start, self.end = start, end
        roi = roi if roi is not None else ingestion.build_roi(self.stations_df)
        self.products_ic = ingestion.load_all_products(
            start, end, roi=roi, products=self.products, stations_df=self.stations_df
        )
        return self

    def extract(self, cache_dir=None):
        from . import extraction

        if self.products_ic is None:
            raise RuntimeError("Call .ingest() before .extract().")
        self._ensure_ee()
        cache_dir = cache_dir if cache_dir is not None else self.cache_dir
        self.sim_df = extraction.extract_all_products(
            self.products_ic, self.stations_df, cache_dir=cache_dir
        )
        return self

    def assign_zones(self, zones_gdf=None, zones_fc=None, use_default_if_none=False):
        from . import zones as _zones

        zones_gdf = zones_gdf if zones_gdf is not None else self.zones_gdf
        zones_fc = zones_fc if zones_fc is not None else self.zones_fc
        if zones_fc is not None or use_default_if_none:
            self._ensure_ee()
        self.stations_df = _zones.assign_zones(
            self.stations_df,
            zones_gdf=zones_gdf,
            zones_fc=zones_fc,
            use_default_if_none=use_default_if_none,
        )
        return self

    def merge(self):
        from . import extraction

        if self.sim_df is None or self.obs_df is None:
            raise RuntimeError(
                "Call .get_observations() and .extract() before .merge()."
            )
        self.merged_df = extraction.merge_with_observations(
            self.sim_df, self.obs_df, stations_df=self.stations_df
        )
        return self

    # ────────────────────────────────────────────────────
    # Preview — look before you validate
    # ────────────────────────────────────────────────────

    def preview_stations(self, m=None, zoom: int = 5):
        """Interactive map of station locations. Works as soon as
        stations are set (before ``.get_observations()`` even) — the
        first sanity check: are these actually where you think they are?
        """
        from . import stations as _stations

        self._ensure_ee()
        stations_df = (
            self.stations_df
            if self.stations_df is not None
            else config.default_stations_wa()
        )
        return _stations.preview_map(stations_df, m=m, zoom=zoom)

    def preview_observations(self, station_id: str | None = None):
        """Quick time-series plot of raw GPCC observations. Requires
        ``.get_observations()`` to have run — no product data needed."""
        from . import viz

        if self.obs_df is None:
            raise RuntimeError(
                "Call .get_observations() before .preview_observations()."
            )
        return viz.preview_observations(self.obs_df, station_id=station_id)

    def preview_comparison(
        self, station_id: str | None = None, product: str | None = None
    ):
        """Quick obs-vs-sim scatter, before running formal validation
        metrics. Requires ``.merge()`` (or ``.validate()``, which calls
        it) to have run."""
        from . import viz

        if self.merged_df is None:
            self.merge()
        return viz.preview_comparison(
            self.merged_df, station_id=station_id, product=product
        )

    def compare_table(self):
        """Obs vs. every product's simulated value, side by side — one
        row per (station, year, month), GPCC in its own column, one
        column per product. The plain "just let me look at the numbers"
        table, underlying every bias/KGE/etc. computed later. Requires
        ``.merge()`` (or ``.validate()``, which calls it) to have run.
        """
        if self.merged_df is None:
            self.merge()

        idx_cols = [
            c for c in ("station_id", "year", "month") if c in self.merged_df.columns
        ]
        pivot = self.merged_df.pivot_table(
            index=idx_cols, columns="product", values="sim_mm_day"
        )
        # obs_mm_day is identical across products for the same
        # station/year/month (it's the same real GPCC observation) —
        # any one row's value is the right one to pull in as the GPCC column.
        obs = self.merged_df.drop_duplicates(idx_cols).set_index(idx_cols)["obs_mm_day"]
        pivot.insert(0, "GPCC", obs)
        return pivot.reset_index()

    def preview_map(
        self,
        product: str,
        kind: str = "daily",
        reference: str | None = None,
        show_gpcc: bool = False,
        region=None,
        m=None,
    ):
        """Interactive map of one product's mean rainfall (``kind=
        "daily"`` or ``"annual"``), or its bias against ANOTHER PRODUCT
        if ``reference`` is given — a gridded-vs-gridded comparison,
        never a GPCC comparison (GPCC has no gridded form here).

        Set ``show_gpcc=True`` to overlay real GPCC station values (not
        a rasterized surface — the true point observations, colored on
        the same scale as the raster) on top of the mean map. Requires
        ``.get_observations()`` to have already run. Ignored when
        ``reference`` is also given (the overlay only applies to the
        single-product mean map).

        Requires ``.ingest()`` to have run.
        """
        from . import spatial

        if self.products_ic is None:
            raise RuntimeError("Call .ingest() before .preview_map().")
        if product not in self.products_ic:
            raise ValueError(
                f"Unknown product {product!r}. Ingested: " f"{sorted(self.products_ic)}"
            )
        if region is None:
            from . import ingestion

            region = ingestion.build_roi(self.stations_df)

        if reference is not None:
            if reference not in self.products_ic:
                raise ValueError(
                    f"Unknown reference {reference!r}. Ingested: "
                    f"{sorted(self.products_ic)}"
                )
            return spatial.preview_bias_map(
                self.products_ic[product],
                self.products_ic[reference],
                product_name=product,
                reference_name=reference,
                region=region,
                m=m,
            )

        obs_df, stations_df = None, None
        if show_gpcc:
            if self.obs_df is None:
                raise RuntimeError(
                    "show_gpcc=True requires .get_observations() to " "have run first."
                )
            obs_df, stations_df = self.obs_df, self.stations_df

        return spatial.preview_mean_map(
            self.products_ic[product],
            product_name=product,
            region=region,
            kind=kind,
            m=m,
            obs_df=obs_df,
            stations_df=stations_df,
        )

    def preview_station_bias(self, product: str, m=None, zoom: int = 5):
        """Interactive map of per-station bias against REAL GPCC
        observations for one product — the actual "does this agree with
        ground truth, and where" spatial check. Requires ``.merge()``
        (or ``.validate()``, which calls it) to have run.
        """
        from . import spatial

        if self.merged_df is None:
            self.merge()
        return spatial.preview_station_bias_map(self.merged_df, product, m=m, zoom=zoom)

    def validate(self):
        from . import validation

        if self.merged_df is None:
            self.merge()
        if "zone" in self.merged_df.columns:
            self.validation_by_zone_df = validation.validate_by_zone(
                self.merged_df, threshold=self.rain_threshold
            )
        self.validation_overall_df = validation.validate_overall(
            self.merged_df, threshold=self.rain_threshold
        )
        self.ranking_df = validation.rank_products(
            self.validation_by_zone_df
            if self.validation_by_zone_df is not None
            else self.validation_overall_df
        )
        return self

    def analyze_thresholds(self, thresholds: list[float] | None = None):
        from . import thresholds as _thresholds

        if self.merged_df is None:
            self.merge()
        self.threshold_df = _thresholds.threshold_sensitivity(
            self.merged_df, thresholds
        )
        return self

    def score(self, normalization: str = "fixed"):
        from . import decision

        validation_df = (
            self.validation_by_zone_df
            if self.validation_by_zone_df is not None
            else self.validation_overall_df
        )
        if validation_df is None:
            raise RuntimeError("Call .validate() before .score().")
        self.scores_df = decision.score_products(
            validation_df, weights=self.app_weights, normalization=normalization
        )
        return self

    def run(self, start: str, end: str, obs_source: str = "download", **obs_kwargs):
        """Run every stage end-to-end with sensible defaults."""
        return (
            self.get_observations(source=obs_source, **obs_kwargs)
            .ingest(start=start, end=end)
            .extract()
            .assign_zones()
            .merge()
            .validate()
            .analyze_thresholds()
            .score()
        )

    # ────────────────────────────────────────────────────
    # Insights
    # ────────────────────────────────────────────────────

    def facts(self):
        from . import insights

        if self._facts is None:
            validation_df = (
                self.validation_by_zone_df
                if self.validation_by_zone_df is not None
                else self.validation_overall_df
            )
            if self.scores_df is None or validation_df is None:
                raise RuntimeError("Call .score() before .facts().")
            self._facts = insights.compute_facts(
                self.scores_df,
                validation_df,
                self.ranking_df,
                self.threshold_df,
                zone_notes=self.zone_notes,
            )
        return self._facts

    def summarize(self) -> str:
        from . import insights

        return insights.summarize(self.facts())

    def answer(self, question: str) -> str:
        from . import insights

        return insights.answer(self.facts(), question)

    # ────────────────────────────────────────────────────
    # Outputs
    # ────────────────────────────────────────────────────

    def export_workbook(self, out_path):
        from . import decision

        if self.scores_df is None:
            self.score()
        return decision.build_workbook(
            out_path,
            self.validation_by_zone_df,
            validation_overall_df=self.validation_overall_df,
            ranking_df=self.ranking_df,
            threshold_df=self.threshold_df,
            scores_df=self.scores_df,
            app_weights=self.app_weights,
            zone_notes=self.zone_notes,
        )

    def show(self, kind: str = "recommendation_heatmap", **kwargs):
        from . import viz

        fn = getattr(viz, kind, None)
        if fn is None:
            raise ValueError(f"Unknown figure kind {kind!r}. See savana.rainfall.viz.")
        target_df = (
            self.scores_df
            if "scores_df" in fn.__code__.co_varnames
            else (
                self.validation_by_zone_df
                if self.validation_by_zone_df is not None
                else self.validation_overall_df
            )
        )
        return fn(target_df, **kwargs)

compare_table()

Obs vs. every product's simulated value, side by side — one row per (station, year, month), GPCC in its own column, one column per product. The plain "just let me look at the numbers" table, underlying every bias/KGE/etc. computed later. Requires .merge() (or .validate(), which calls it) to have run.

Source code in savana/rainfall/pipeline.py

def compare_table(self):
    """Obs vs. every product's simulated value, side by side — one
    row per (station, year, month), GPCC in its own column, one
    column per product. The plain "just let me look at the numbers"
    table, underlying every bias/KGE/etc. computed later. Requires
    ``.merge()`` (or ``.validate()``, which calls it) to have run.
    """
    if self.merged_df is None:
        self.merge()

    idx_cols = [
        c for c in ("station_id", "year", "month") if c in self.merged_df.columns
    ]
    pivot = self.merged_df.pivot_table(
        index=idx_cols, columns="product", values="sim_mm_day"
    )
    # obs_mm_day is identical across products for the same
    # station/year/month (it's the same real GPCC observation) —
    # any one row's value is the right one to pull in as the GPCC column.
    obs = self.merged_df.drop_duplicates(idx_cols).set_index(idx_cols)["obs_mm_day"]
    pivot.insert(0, "GPCC", obs)
    return pivot.reset_index()

preview_comparison(station_id=None, product=None)

Quick obs-vs-sim scatter, before running formal validation metrics. Requires .merge() (or .validate(), which calls it) to have run.

Source code in savana/rainfall/pipeline.py

def preview_comparison(
    self, station_id: str | None = None, product: str | None = None
):
    """Quick obs-vs-sim scatter, before running formal validation
    metrics. Requires ``.merge()`` (or ``.validate()``, which calls
    it) to have run."""
    from . import viz

    if self.merged_df is None:
        self.merge()
    return viz.preview_comparison(
        self.merged_df, station_id=station_id, product=product
    )

preview_map(product, kind='daily', reference=None, show_gpcc=False, region=None, m=None)

Interactive map of one product's mean rainfall (kind= "daily" or "annual"), or its bias against ANOTHER PRODUCT if reference is given — a gridded-vs-gridded comparison, never a GPCC comparison (GPCC has no gridded form here).

Set show_gpcc=True to overlay real GPCC station values (not a rasterized surface — the true point observations, colored on the same scale as the raster) on top of the mean map. Requires .get_observations() to have already run. Ignored when reference is also given (the overlay only applies to the single-product mean map).

Requires .ingest() to have run.

Source code in savana/rainfall/pipeline.py

def preview_map(
    self,
    product: str,
    kind: str = "daily",
    reference: str | None = None,
    show_gpcc: bool = False,
    region=None,
    m=None,
):
    """Interactive map of one product's mean rainfall (``kind=
    "daily"`` or ``"annual"``), or its bias against ANOTHER PRODUCT
    if ``reference`` is given — a gridded-vs-gridded comparison,
    never a GPCC comparison (GPCC has no gridded form here).

    Set ``show_gpcc=True`` to overlay real GPCC station values (not
    a rasterized surface — the true point observations, colored on
    the same scale as the raster) on top of the mean map. Requires
    ``.get_observations()`` to have already run. Ignored when
    ``reference`` is also given (the overlay only applies to the
    single-product mean map).

    Requires ``.ingest()`` to have run.
    """
    from . import spatial

    if self.products_ic is None:
        raise RuntimeError("Call .ingest() before .preview_map().")
    if product not in self.products_ic:
        raise ValueError(
            f"Unknown product {product!r}. Ingested: " f"{sorted(self.products_ic)}"
        )
    if region is None:
        from . import ingestion

        region = ingestion.build_roi(self.stations_df)

    if reference is not None:
        if reference not in self.products_ic:
            raise ValueError(
                f"Unknown reference {reference!r}. Ingested: "
                f"{sorted(self.products_ic)}"
            )
        return spatial.preview_bias_map(
            self.products_ic[product],
            self.products_ic[reference],
            product_name=product,
            reference_name=reference,
            region=region,
            m=m,
        )

    obs_df, stations_df = None, None
    if show_gpcc:
        if self.obs_df is None:
            raise RuntimeError(
                "show_gpcc=True requires .get_observations() to " "have run first."
            )
        obs_df, stations_df = self.obs_df, self.stations_df

    return spatial.preview_mean_map(
        self.products_ic[product],
        product_name=product,
        region=region,
        kind=kind,
        m=m,
        obs_df=obs_df,
        stations_df=stations_df,
    )

preview_observations(station_id=None)

Quick time-series plot of raw GPCC observations. Requires .get_observations() to have run — no product data needed.

Source code in savana/rainfall/pipeline.py

def preview_observations(self, station_id: str | None = None):
    """Quick time-series plot of raw GPCC observations. Requires
    ``.get_observations()`` to have run — no product data needed."""
    from . import viz

    if self.obs_df is None:
        raise RuntimeError(
            "Call .get_observations() before .preview_observations()."
        )
    return viz.preview_observations(self.obs_df, station_id=station_id)

preview_station_bias(product, m=None, zoom=5)

Interactive map of per-station bias against REAL GPCC observations for one product — the actual "does this agree with ground truth, and where" spatial check. Requires .merge() (or .validate(), which calls it) to have run.

Source code in savana/rainfall/pipeline.py

def preview_station_bias(self, product: str, m=None, zoom: int = 5):
    """Interactive map of per-station bias against REAL GPCC
    observations for one product — the actual "does this agree with
    ground truth, and where" spatial check. Requires ``.merge()``
    (or ``.validate()``, which calls it) to have run.
    """
    from . import spatial

    if self.merged_df is None:
        self.merge()
    return spatial.preview_station_bias_map(self.merged_df, product, m=m, zoom=zoom)

preview_stations(m=None, zoom=5)

Interactive map of station locations. Works as soon as stations are set (before .get_observations() even) — the first sanity check: are these actually where you think they are?

Source code in savana/rainfall/pipeline.py

def preview_stations(self, m=None, zoom: int = 5):
    """Interactive map of station locations. Works as soon as
    stations are set (before ``.get_observations()`` even) — the
    first sanity check: are these actually where you think they are?
    """
    from . import stations as _stations

    self._ensure_ee()
    stations_df = (
        self.stations_df
        if self.stations_df is not None
        else config.default_stations_wa()
    )
    return _stations.preview_map(stations_df, m=m, zoom=zoom)

run(start, end, obs_source='download', **obs_kwargs)

Run every stage end-to-end with sensible defaults.

Source code in savana/rainfall/pipeline.py

def run(self, start: str, end: str, obs_source: str = "download", **obs_kwargs):
    """Run every stage end-to-end with sensible defaults."""
    return (
        self.get_observations(source=obs_source, **obs_kwargs)
        .ingest(start=start, end=end)
        .extract()
        .assign_zones()
        .merge()
        .validate()
        .analyze_thresholds()
        .score()
    )

validate_against_gpcc(stations=None, products=None, start_year=2001, end_year=2020, obs_source=None, obs_csv=None, cache_dir='savana_rainfall_data', zones_gdf=None, zones_fc=None, rain_threshold=None, ee_project=None)

Validate one or more precipitation products against GPCC gauge observations at one or more stations, over a chosen year range — the one-call version of the whole assessment, matching the original paper's exact logic (16 WA stations, 6 products, 2001-2020) as the default, everything else overridable by simple parameters.

This is the function to reach for first. It does exactly what the original per-station CSV workflow did — extract each requested product at each requested station, save/reuse a per-product CSV (cache_dir/precip_extraction_<PRODUCT>.csv, same as before), merge against GPCC observations (cache_dir/gpcc_obs_<start>_<end>.csv), and write the same result CSVs the original scripts did (validation_by_zone.csv, validation_overall.csv, product_ranking.csv, threshold_sensitivity.csv) — just wrapped in one call instead of six separate scripts.

Parameters:

Name	Type	Description	Default
stations		where to validate. Any of: - None (default): the 16 WA GPCC stations from the paper. - a stations_df, a path to a .geojson/.csv file of station points, a list of (lon, lat) tuples, or a single (lon, lat) tuple — see :func:savana.rainfall.stations.load_stations_any for the full list of accepted shapes. Works the same whether you give it 1 station or 100.	None
products	list[str] | None	which products to check, by name (e.g. ["CHIRPS", "GPM_IMERG"]). None (default) uses all 6 in :data:config.DEFAULT_PRODUCTS. Any subset works.	None
start_year, end_year		inclusive year range (plain ints — the paper used 2001-2020; pick whatever you need).	required
obs_source	str | None	where GPCC observations come from — "ee_asset" (fast, only covers the 16 WA stations), "download" (slower, works for any station anywhere), or "csv" (use obs_csv= — you already have one). Defaults to "ee_asset" when stations is the WA default (fastest path for the paper's own network) and "download" otherwise (since the EE asset only has the WA 16).	None
obs_csv		required if obs_source="csv".	None
cache_dir		where per-product extraction CSVs, the GPCC obs CSV, and the result CSVs are read from / written to. Sits right next to your notebook by default; set to None to skip all file caching and keep everything in memory only.	'savana_rainfall_data'
zones_gdf, zones_fc		optional zone geometry (see :mod:savana.rainfall.zones) for zone-stratified results. Omit for pooled (unzoned) validation.	required
rain_threshold	float | None	mm/day wet/dry threshold for categorical metrics. Defaults to the WMO standard (1.0 mm/day).	None
ee_project	str | None	Google Cloud project registered for Earth Engine use (only needed for obs_source="ee_asset" or the default Earth Engine ingestion — not needed at all if you only use obs_source="csv"). If omitted, uses whatever is already configured for the environment (see savana.ee_init.initialize) — set this explicitly if you have more than one Google Cloud project and the wrong one keeps getting picked up.	None

Returns:

Type	Description
	A fully populated :class:RainfallAssessment — inspect
	.validation_by_zone_df / .validation_overall_df
	directly, or call .summarize(), .answer("..."),
	.show(), .export_workbook(...) on it, same as building
	one by hand.

Example::

from savana.rainfall import validate_against_gpcc

# Reproduce the paper exactly:
result = validate_against_gpcc()

# One station, two products, a shorter recent period:
result = validate_against_gpcc(
    stations=(-1.5, 12.4),
    products=["CHIRPS", "GPM_IMERG"],
    start_year=2018, end_year=2023,
)
print(result.validation_overall_df)

Source code in savana/rainfall/pipeline.py

def validate_against_gpcc(
    stations=None,
    products: list[str] | None = None,
    start_year: int = 2001,
    end_year: int = 2020,
    obs_source: str | None = None,
    obs_csv=None,
    cache_dir="savana_rainfall_data",
    zones_gdf=None,
    zones_fc=None,
    rain_threshold: float | None = None,
    ee_project: str | None = None,
):
    """Validate one or more precipitation products against GPCC gauge
    observations at one or more stations, over a chosen year range —
    the one-call version of the whole assessment, matching the original
    paper's exact logic (16 WA stations, 6 products, 2001-2020) as the
    default, everything else overridable by simple parameters.

    This is the function to reach for first. It does exactly what the
    original per-station CSV workflow did — extract each requested
    product at each requested station, save/reuse a per-product CSV
    (``cache_dir/precip_extraction_<PRODUCT>.csv``, same as before),
    merge against GPCC observations
    (``cache_dir/gpcc_obs_<start>_<end>.csv``), and write the same
    result CSVs the original scripts did (``validation_by_zone.csv``,
    ``validation_overall.csv``, ``product_ranking.csv``,
    ``threshold_sensitivity.csv``) — just wrapped in one call instead of
    six separate scripts.

    Args:
        stations: where to validate. Any of:
            - ``None`` (default): the 16 WA GPCC stations from the paper.
            - a ``stations_df``, a path to a ``.geojson``/``.csv`` file
              of station points, a list of ``(lon, lat)`` tuples, or a
              single ``(lon, lat)`` tuple — see
              :func:`savana.rainfall.stations.load_stations_any` for
              the full list of accepted shapes. Works the same whether
              you give it 1 station or 100.
        products: which products to check, by name (e.g.
            ``["CHIRPS", "GPM_IMERG"]``). ``None`` (default) uses all 6
            in :data:`config.DEFAULT_PRODUCTS`. Any subset works.
        start_year, end_year: inclusive year range (plain ints — the
            paper used 2001-2020; pick whatever you need).
        obs_source: where GPCC observations come from —
            ``"ee_asset"`` (fast, only covers the 16 WA stations),
            ``"download"`` (slower, works for any station anywhere),
            or ``"csv"`` (use ``obs_csv=`` — you already have one).
            Defaults to
            ``"ee_asset"`` when ``stations`` is the WA default (fastest
            path for the paper's own network) and ``"download"``
            otherwise (since the EE asset only has the WA 16).
        obs_csv: required if ``obs_source="csv"``.
        cache_dir: where per-product extraction CSVs, the GPCC obs CSV,
            and the result CSVs are read from / written to. Sits right
            next to your notebook by default; set to ``None`` to skip
            all file caching and keep everything in memory only.
        zones_gdf, zones_fc: optional zone geometry (see
            :mod:`savana.rainfall.zones`) for zone-stratified results.
            Omit for pooled (unzoned) validation.
        rain_threshold: mm/day wet/dry threshold for categorical
            metrics. Defaults to the WMO standard (1.0 mm/day).
        ee_project: Google Cloud project registered for Earth Engine use
            (only needed for ``obs_source="ee_asset"`` or the default
            Earth Engine ingestion — not needed at all if you only use
            ``obs_source="csv"``). If omitted, uses whatever
            is already configured for the environment (see
            ``savana.ee_init.initialize``) — set this explicitly if
            you have more than one Google Cloud project and the wrong
            one keeps getting picked up.

    Returns:
        A fully populated :class:`RainfallAssessment` — inspect
        ``.validation_by_zone_df`` / ``.validation_overall_df``
        directly, or call ``.summarize()``, ``.answer("...")``,
        ``.show()``, ``.export_workbook(...)`` on it, same as building
        one by hand.

    Example::

        from savana.rainfall import validate_against_gpcc

        # Reproduce the paper exactly:
        result = validate_against_gpcc()

        # One station, two products, a shorter recent period:
        result = validate_against_gpcc(
            stations=(-1.5, 12.4),
            products=["CHIRPS", "GPM_IMERG"],
            start_year=2018, end_year=2023,
        )
        print(result.validation_overall_df)
    """
    from pathlib import Path

    from . import stations as _stations

    stations_df = _stations.load_stations_any(stations)
    is_default_wa_network = stations is None

    if products is None:
        selected_products = config.DEFAULT_PRODUCTS
    else:
        unknown = [p for p in products if p not in config.DEFAULT_PRODUCTS]
        if unknown:
            raise ValueError(
                f"Unknown product(s) {unknown}. Known products: "
                f"{sorted(config.DEFAULT_PRODUCTS)}."
            )
        selected_products = {p: config.DEFAULT_PRODUCTS[p] for p in products}

    if obs_source is None:
        obs_source = "ee_asset" if is_default_wa_network else "download"

    cache_path = Path(cache_dir) if cache_dir else None

    ra = RainfallAssessment(
        products=selected_products,
        stations=stations_df,
        zones_gdf=zones_gdf,
        zones_fc=zones_fc,
        rain_threshold=rain_threshold,
        cache_dir=cache_path,
        ee_project=ee_project,
    )

    obs_kwargs = {}
    if obs_source == "csv":
        if obs_csv is None:
            raise ValueError('obs_source="csv" requires obs_csv=<path>.')
        obs_kwargs = {"obs_csv": obs_csv}
    elif obs_source == "download":
        obs_kwargs = {"start_year": start_year, "end_year": end_year}
        if cache_path:
            obs_kwargs["data_dir"] = cache_path
    # "ee_asset" takes no year kwargs — filtered to the requested range below instead

    ra.get_observations(source=obs_source, **obs_kwargs)
    ra.obs_df = ra.obs_df[
        (ra.obs_df["year"] >= start_year) & (ra.obs_df["year"] <= end_year)
    ].reset_index(drop=True)
    if ra.obs_df.empty:
        raise ValueError(
            f"No GPCC observations found for {start_year}-{end_year} with "
            f"obs_source={obs_source!r}. Check the year range against what "
            f"that source actually covers."
        )
    ra.ingest(start=f"{start_year}-01-01", end=f"{end_year}-12-31")
    ra.extract(cache_dir=cache_path)
    if zones_gdf is not None or zones_fc is not None:
        ra.assign_zones()
    ra.merge()
    ra.validate()
    ra.analyze_thresholds()
    ra.score()

    if cache_path:
        cache_path.mkdir(parents=True, exist_ok=True)
        ra.merged_df.to_csv(
            cache_path
            / (
                "merged_obs_grid_zoned.csv"
                if "zone" in ra.merged_df.columns
                else "merged_obs_grid.csv"
            ),
            index=False,
        )
        if ra.validation_by_zone_df is not None:
            ra.validation_by_zone_df.to_csv(
                cache_path / "validation_by_zone.csv", index=False
            )
        ra.validation_overall_df.to_csv(
            cache_path / "validation_overall.csv", index=False
        )
        ra.ranking_df.to_csv(cache_path / "product_ranking.csv", index=False)
        ra.threshold_df.to_csv(cache_path / "threshold_sensitivity.csv", index=False)
        print(f"  Result CSVs written to: {cache_path}")

    return ra

savana.rainfall.config

Default configuration for global precipitation product assessment.

Everything here is a default, every public function in savana.rainfall accepts overrides, so a user assessing a different region, a different subset of products, their own gauge network, or their own application weights is not locked into the West Africa study configuration. The WA study (16 GPCC FDD v2022 stations, 5 ecological zones, 6 global precipitation products, 7 conservation/water-management applications) ships as the default so savana.rainfall is useful out of the box and reproduces the original manuscript, but nothing here is required to use the package on a different AOI.

Nothing in this module touches Earth Engine or hits the network, it is pure data, safe to import eagerly. EE objects (ee.FeatureCollection, ee.Geometry) are built lazily, inside functions in :mod:.stations and :mod:.ingestion, exactly where the JS/EE equivalents built them.

default_stations_wa()

The 16 West Africa GPCC gauge stations, as a pandas.DataFrame.

This is the manuscript's real validation network, not a synthetic placeholder, the default stations_df used throughout savana.rainfall when no stations_df is supplied. Any function accepting stations_df= accepts a DataFrame with this same shape (station_id, station_name, lon, lat, elevation_m, source) for a different gauge network.

Source code in savana/rainfall/config.py

def default_stations_wa():
    """The 16 West Africa GPCC gauge stations, as a pandas.DataFrame.

    This is the manuscript's real validation network, not a synthetic
    placeholder, the default ``stations_df`` used throughout
    ``savana.rainfall`` when no ``stations_df`` is supplied. Any function
    accepting ``stations_df=`` accepts a DataFrame with this same shape
    (station_id, station_name, lon, lat, elevation_m, source) for a
    different gauge network.
    """
    import pandas as pd

    return pd.DataFrame(DEFAULT_STATIONS_WA_RAW, columns=DEFAULT_STATION_COLUMNS)

savana.rainfall.stations

Gauge station metadata and precipitation observation loading.

Every function here works on an arbitrary stations_df — a pandas.DataFrame with at minimum station_id, lon, lat columns (station_name, elevation_m, source are recommended but not required). savana.rainfall.config.default_stations_wa() supplies the 16-station West Africa GPCC network as a convenient default so the package works out of the box, but nothing here assumes those specific stations. Point this module at your own gauge network by building a DataFrame in that shape and passing it as stations_df= throughout.

Four ways to get observations, in increasing order of "how much can this handle a station set that isn't the WA 16":

1. :func:load_stations_from_csv — you already have your own station metadata + observation CSVs. Fully general.
2. :func:download_gpcc — downloads the public GPCC Full Data Daily v2022 archive and extracts at whatever station coordinates you give it. Fully general, works for any station anywhere GPCC has coverage, but downloads ~440 MB and is slow the first time.
3. :func:load_gpcc_obs_from_asset — fast, but only returns rows for station_ids that exist in the given EE table asset. The packaged default asset (:data:config.DEFAULT_GPCC_ASSET_WA) covers only the 16 WA stations; point asset_id at your own pre-extracted table for a different network, or use option 1/2 instead.

All three return real gauge observations. There is deliberately no synthetic/simulated observation option: GPCC gauge data is the ground truth this package validates against, so fabricating it would make every resulting metric meaningless.

download_gpcc(stations_df=None, start_year=2001, end_year=2020, data_dir=None, keep_raw=False)

Download the public GPCC archive and extract at any station set.

Works for any stations_df (defaults to the WA 16), anywhere the GPCC 1.0-degree grid has coverage — this is the fully general path, unlike :func:load_gpcc_obs_from_asset which only covers whatever stations happen to already be in an EE asset.

Downloads are cached: files already present in data_dir are skipped, so re-running after a partial failure only fetches what's missing. Requires requests, xarray, netCDF4 (installed with pip install "savana[rainfall]").

Parameters:

Name	Type	Description	Default
stations_df		defaults to :func:config.default_stations_wa.	None
start_year, end_year		inclusive year range.	required
data_dir	str | Path | None	local cache directory. Defaults to ./savana_rainfall_data in the current working directory.	None
keep_raw	bool	if False (default), deletes the raw yearly NetCDF files after extraction to save disk space (~20 MB/year).	False

Returns:

Type	Description
	pandas.DataFrame with columns station_id, year, month, obs_mm_day.

Source code in savana/rainfall/stations.py

def download_gpcc(
    stations_df=None,
    start_year: int = 2001,
    end_year: int = 2020,
    data_dir: str | Path | None = None,
    keep_raw: bool = False,
):
    """Download the public GPCC archive and extract at any station set.

    Works for any ``stations_df`` (defaults to the WA 16), anywhere the
    GPCC 1.0-degree grid has coverage — this is the fully general path,
    unlike :func:`load_gpcc_obs_from_asset` which only covers whatever
    stations happen to already be in an EE asset.

    Downloads are cached: files already present in ``data_dir`` are
    skipped, so re-running after a partial failure only fetches what's
    missing. Requires ``requests``, ``xarray``, ``netCDF4`` (installed
    with ``pip install "savana[rainfall]"``).

    Args:
        stations_df: defaults to :func:`config.default_stations_wa`.
        start_year, end_year: inclusive year range.
        data_dir: local cache directory. Defaults to
            ``./savana_rainfall_data`` in the current working directory.
        keep_raw: if False (default), deletes the raw yearly NetCDF
            files after extraction to save disk space (~20 MB/year).

    Returns:
        pandas.DataFrame with columns station_id, year, month, obs_mm_day.
    """
    import pandas as pd

    stations_df = (
        stations_df if stations_df is not None else config.default_stations_wa()
    )
    _validate_stations_df(stations_df)

    data_dir = Path(data_dir) if data_dir else Path.cwd() / "savana_rainfall_data"
    raw_dir = data_dir / "gpcc_raw"
    raw_dir.mkdir(parents=True, exist_ok=True)

    print(
        f"  GPCC download: {len(stations_df)} station(s), "
        f"{start_year}-{end_year}, cache: {raw_dir}"
    )

    all_rows = []
    for year in range(start_year, end_year + 1):
        try:
            nc_path = _download_gpcc_year(year, raw_dir)
        except Exception as exc:  # noqa: BLE001
            print(f"  \u26a0  {year}: download failed ({type(exc).__name__}: {exc})")
            continue

        try:
            df_yr = _extract_monthly_means(nc_path, stations_df)
            df_yr = df_yr[df_yr["year"] == year]
            all_rows.append(df_yr)
            expected = len(stations_df) * 12
            flag = "" if len(df_yr) >= expected * 0.9 else " \u26a0"
            print(f"  {year}: {len(df_yr)} rows (expected {expected}){flag}")
        except Exception as exc:  # noqa: BLE001
            print(f"  \u26a0  {year}: extraction failed ({type(exc).__name__}: {exc})")
            continue
        finally:
            if not keep_raw:
                nc_path.unlink(missing_ok=True)

    if not all_rows:
        raise RuntimeError("No GPCC data could be extracted for any year.")

    combined = pd.concat(all_rows, ignore_index=True)
    combined = combined.sort_values(["station_id", "year", "month"]).reset_index(
        drop=True
    )

    out_path = data_dir / f"gpcc_obs_{start_year}_{end_year}.csv"
    combined.to_csv(out_path, index=False)
    print(f"  Done: {len(combined):,} rows -> {out_path}")
    return combined

get_observations(stations_df=None, source='download', **kwargs)

Single entry point for getting gauge observations, any station set.

Parameters:

Name	Type	Description	Default
stations_df		defaults to :func:config.default_stations_wa.	None
source	str	one of - "download" (default): :func:download_gpcc — fully general, works for any station, slow on first run. - "ee_asset": :func:load_gpcc_obs_from_asset — fast, limited to whatever stations are already in the asset. - "csv": :func:load_stations_from_csv's obs half — requires obs_csv= in kwargs.	'download'
**kwargs		forwarded to the selected loader.	{}

Returns:

Type	Description
	pandas.DataFrame with columns station_id, year, month, obs_mm_day.

Source code in savana/rainfall/stations.py

def get_observations(stations_df=None, source: str = "download", **kwargs):
    """Single entry point for getting gauge observations, any station set.

    Args:
        stations_df: defaults to :func:`config.default_stations_wa`.
        source: one of
            - ``"download"`` (default): :func:`download_gpcc` — fully
              general, works for any station, slow on first run.
            - ``"ee_asset"``: :func:`load_gpcc_obs_from_asset` — fast,
              limited to whatever stations are already in the asset.
            - ``"csv"``: :func:`load_stations_from_csv`'s obs half —
              requires ``obs_csv=`` in ``kwargs``.
        **kwargs: forwarded to the selected loader.

    Returns:
        pandas.DataFrame with columns station_id, year, month, obs_mm_day.
    """
    stations_df = (
        stations_df if stations_df is not None else config.default_stations_wa()
    )

    if source == "download":
        return download_gpcc(stations_df, **kwargs)
    if source == "ee_asset":
        return load_gpcc_obs_from_asset(stations_df, **kwargs)
    if source == "csv":
        obs_csv = kwargs.pop("obs_csv", None)
        if obs_csv is None:
            raise ValueError('source="csv" requires obs_csv=<path> in kwargs.')
        import pandas as pd

        obs_df = pd.read_csv(obs_csv)
        missing = REQUIRED_OBS_COLUMNS - set(obs_df.columns)
        if missing:
            raise ValueError(f"obs_csv is missing column(s): {sorted(missing)}")
        return obs_df

    raise ValueError(
        f"Unknown source {source!r}. Expected one of: "
        f'"download", "ee_asset", "csv".'
    )

load_gpcc_obs_from_asset(stations_df=None, asset_id=None)

Load pre-extracted GPCC observations from an Earth Engine table asset.

Fast (no download, no NetCDF processing) but only returns rows for station_id values that already exist in the asset. Any station in stations_df not found in the asset is reported via a printed warning, not silently dropped without explanation — use :func:download_gpcc for those instead, or build your own asset with station_id, year, month, obs_mm_day columns and pass its ID here.

Parameters:

Name	Type	Description	Default
stations_df		defaults to :func:config.default_stations_wa.	None
asset_id	str | None	EE table asset ID. Defaults to :data:config.DEFAULT_GPCC_ASSET_WA, which only covers the 16 default WA stations — pass your own asset_id for any other station set.	None

Returns:

Type	Description
	pandas.DataFrame with columns station_id, year, month, obs_mm_day.

Source code in savana/rainfall/stations.py

def load_gpcc_obs_from_asset(stations_df=None, asset_id: str | None = None):
    """Load pre-extracted GPCC observations from an Earth Engine table asset.

    Fast (no download, no NetCDF processing) but only returns rows for
    ``station_id`` values that already exist in the asset. Any station in
    ``stations_df`` not found in the asset is reported via a printed
    warning, not silently dropped without explanation — use
    :func:`download_gpcc` for those instead, or build your own asset with
    ``station_id, year, month, obs_mm_day`` columns and pass its ID here.

    Args:
        stations_df: defaults to :func:`config.default_stations_wa`.
        asset_id: EE table asset ID. Defaults to
            :data:`config.DEFAULT_GPCC_ASSET_WA`, which only covers the
            16 default WA stations — pass your own asset_id for any
            other station set.

    Returns:
        pandas.DataFrame with columns station_id, year, month, obs_mm_day.
    """
    import ee
    import pandas as pd

    stations_df = (
        stations_df if stations_df is not None else config.default_stations_wa()
    )
    asset_id = asset_id or config.DEFAULT_GPCC_ASSET_WA
    _validate_stations_df(stations_df)

    fc = ee.FeatureCollection(asset_id)
    info = fc.getInfo()
    rows = [f["properties"] for f in info["features"]]
    obs_df = pd.DataFrame(rows)
    if obs_df.empty:
        raise ValueError(f"EE asset {asset_id!r} returned no features.")

    missing_cols = REQUIRED_OBS_COLUMNS - set(obs_df.columns)
    if missing_cols:
        raise ValueError(
            f"EE asset {asset_id!r} is missing expected column(s) "
            f"{sorted(missing_cols)}. Expected {sorted(REQUIRED_OBS_COLUMNS)}."
        )

    obs_df["year"] = obs_df["year"].astype(int)
    obs_df["month"] = obs_df["month"].astype(int)
    obs_df["obs_mm_day"] = obs_df["obs_mm_day"].astype(float)

    wanted_ids = set(stations_df["station_id"])
    found_ids = set(obs_df["station_id"].unique())
    obs_df = obs_df[obs_df["station_id"].isin(wanted_ids)].reset_index(drop=True)

    not_found = wanted_ids - found_ids
    if not_found:
        print(
            f"  \u26a0  {len(not_found)} station(s) in stations_df have no data "
            f"in asset {asset_id!r}: {sorted(not_found)}\n"
            f"     Use download_gpcc() for these, or point asset_id at an "
            f"asset that includes them."
        )

    print(
        f"  GPCC observations loaded from asset: {len(obs_df):,} rows "
        f"({obs_df['station_id'].nunique()} station(s))"
    )
    return obs_df

load_stations_any(stations=None)

Turn almost anything describing station locations into a proper stations_df — the single entry point every high-level function (:func:savana.rainfall.pipeline.validate_against_gpcc) uses so a user never has to hand-build a DataFrame just to try one station.

Accepts

• None -> :func:savana.rainfall.config.default_stations_wa (the 16 WA GPCC stations).
• an existing stations_df (DataFrame with station_id, lon, lat) -> validated and returned as-is.
• a path to a .geojson/.json file of Point features -> one station per feature; station_id/station_name are read from feature properties if present, else auto-generated.
• a path to a .csv file -> loaded via :func:load_stations_from_csv's station-table shape.
• a list of (lon, lat) or (station_id, lon, lat) tuples, or a list of dicts with at least lon/lat keys.
• a single station as (lon, lat) or [lon, lat] — both a tuple and a plain 2-element list work.

Returns:

Type	Description
	A validated stations_df.

Source code in savana/rainfall/stations.py

def load_stations_any(stations=None):
    """Turn almost anything describing station locations into a proper
    ``stations_df`` — the single entry point every high-level function
    (:func:`savana.rainfall.pipeline.validate_against_gpcc`) uses so a
    user never has to hand-build a DataFrame just to try one station.

    Accepts:
        - ``None`` -> :func:`savana.rainfall.config.default_stations_wa`
          (the 16 WA GPCC stations).
        - an existing ``stations_df`` (DataFrame with
          ``station_id, lon, lat``) -> validated and returned as-is.
        - a path to a ``.geojson``/``.json`` file of Point features ->
          one station per feature; ``station_id``/``station_name`` are
          read from feature properties if present, else auto-generated.
        - a path to a ``.csv`` file -> loaded via
          :func:`load_stations_from_csv`'s station-table shape.
        - a list of ``(lon, lat)`` or ``(station_id, lon, lat)`` tuples,
          or a list of dicts with at least ``lon``/``lat`` keys.
        - a single station as ``(lon, lat)`` or ``[lon, lat]`` — both a
          tuple and a plain 2-element list work.

    Returns:
        A validated ``stations_df``.
    """
    import pandas as pd

    if stations is None:
        return config.default_stations_wa()

    if hasattr(stations, "columns"):  # already a DataFrame
        _validate_stations_df(stations)
        return stations

    if isinstance(stations, (str, Path)):
        # A plain string can be either a file path or typed coordinates
        # ("-0.17, 5.56", or several stations as "lon,lat; lon,lat").
        # Typed coordinates are common from GUIs and command lines, and
        # silently treating "-0.17, 5.56" as a filename produces a very
        # confusing "unknown '.56' file type" error, so try parsing as
        # coordinates first whenever the text can't be a path.
        parsed = _stations_from_coord_string(stations)
        if parsed is not None:
            stations = parsed
        else:
            path = Path(stations)
            if path.suffix.lower() in (".geojson", ".json"):
                return _stations_from_geojson(path)
            if path.suffix.lower() == ".csv":
                df, _ = load_stations_from_csv(path)
                return df
            raise ValueError(
                f"Don't know how to load stations from {str(stations)!r}. "
                f"Expected a .geojson/.json/.csv path, or coordinates "
                f'like "-0.17, 5.56" (lon, lat) or '
                f'"-0.17,5.56; 2.1,6.5" for several stations.'
            )

    # Single station shorthand: (lon, lat) or [lon, lat] -- two plain
    # numbers, not a list of multiple stations. Deliberately checked
    # before the general list/tuple-of-stations branch below, and
    # deliberately accepts both tuple and list (parentheses around a
    # single list literal in Python don't make it a tuple -- ([-1.5,
    # 12.4]) is just [-1.5, 12.4] -- so both forms need to work).
    if (
        isinstance(stations, (tuple, list))
        and len(stations) == 2
        and all(isinstance(x, (int, float)) for x in stations)
    ):
        stations = [tuple(stations)]

    if isinstance(stations, (list, tuple)):
        rows = []
        for i, item in enumerate(stations):
            if isinstance(item, dict):
                row = dict(item)
                row.setdefault("station_id", f"S{i + 1:03d}")
            elif len(item) == 2:
                row = {"station_id": f"S{i + 1:03d}", "lon": item[0], "lat": item[1]}
            elif len(item) == 3:
                row = {"station_id": item[0], "lon": item[1], "lat": item[2]}
            else:
                raise ValueError(f"Can't parse station entry: {item!r}")
            rows.append(row)
        df = pd.DataFrame(rows)
        _validate_stations_df(df)
        return df

    raise ValueError(
        f"Don't know how to interpret stations={stations!r} (type "
        f"{type(stations).__name__}). See load_stations_any() docstring "
        f"for accepted formats."
    )

load_stations_from_csv(stations_csv, obs_csv=None)

Load your own station metadata (and optionally observations) from CSV.

stations_csv must have columns station_id, station_name, lon, lat, elevation_m, source (matching :data:config.DEFAULT_STATION_COLUMNS). obs_csv, if given, must have columns station_id, year, month, obs_mm_day.

This is the fully general entry point for a station network that isn't West Africa's 16 GPCC stations at all — bring your own gauge metadata and (optionally) your own already-extracted observations.

Source code in savana/rainfall/stations.py

def load_stations_from_csv(stations_csv: str | Path, obs_csv: str | Path | None = None):
    """Load your own station metadata (and optionally observations) from CSV.

    ``stations_csv`` must have columns
    ``station_id, station_name, lon, lat, elevation_m, source``
    (matching :data:`config.DEFAULT_STATION_COLUMNS`).
    ``obs_csv``, if given, must have columns
    ``station_id, year, month, obs_mm_day``.

    This is the fully general entry point for a station network that
    isn't West Africa's 16 GPCC stations at all — bring your own gauge
    metadata and (optionally) your own already-extracted observations.
    """
    import pandas as pd

    stations_df = pd.read_csv(stations_csv)
    missing_stn = set(config.DEFAULT_STATION_COLUMNS) - set(stations_df.columns)
    if missing_stn:
        raise ValueError(f"stations_csv is missing column(s): {sorted(missing_stn)}")
    _validate_stations_df(stations_df)
    print(f"  Stations loaded: {len(stations_df)} from {stations_csv}")

    if obs_csv is None:
        return stations_df, None

    obs_df = pd.read_csv(obs_csv)
    missing_obs = REQUIRED_OBS_COLUMNS - set(obs_df.columns)
    if missing_obs:
        raise ValueError(f"obs_csv is missing column(s): {sorted(missing_obs)}")
    obs_df["year"] = obs_df["year"].astype(int)
    obs_df["month"] = obs_df["month"].astype(int)
    obs_df["obs_mm_day"] = obs_df["obs_mm_day"].astype(float)
    print(f"  Observations loaded: {len(obs_df):,} rows from {obs_csv}")

    return stations_df, obs_df

preview_map(stations_df=None, m=None, zoom=5)

A quick interactive map of station locations — the first thing to check before extracting or validating anything: "are these actually where I think they are?"

Parameters:

Name	Type	Description	Default
stations_df		defaults to :func:config.default_stations_wa.	None
m		an existing geemap.Map to add to, or a new one is created.	None
zoom	int	zoom level when centering on the stations.	5

Returns:

Type	Description
	A geemap.Map with one styled point layer for the stations.
	Click a point on the map to see its properties
	(station_id, station_name, lon, lat) in
	geemap's built-in inspector panel.

Source code in savana/rainfall/stations.py

def preview_map(stations_df=None, m=None, zoom: int = 5):
    """A quick interactive map of station locations — the first thing to
    check before extracting or validating anything: "are these actually
    where I think they are?"

    Args:
        stations_df: defaults to :func:`config.default_stations_wa`.
        m: an existing ``geemap.Map`` to add to, or a new one is created.
        zoom: zoom level when centering on the stations.

    Returns:
        A ``geemap.Map`` with one styled point layer for the stations.
        Click a point on the map to see its properties
        (``station_id``, ``station_name``, ``lon``, ``lat``) in
        geemap's built-in inspector panel.
    """
    import geemap

    stations_df = (
        stations_df if stations_df is not None else config.default_stations_wa()
    )
    _validate_stations_df(stations_df)

    if m is None:
        m = geemap.Map()

    fc = stations_to_ee_fc(stations_df)
    center_lon = float(stations_df.lon.mean())
    center_lat = float(stations_df.lat.mean())
    m.set_center(center_lon, center_lat, zoom)
    m.add_layer(fc.style(**{"color": "FFEB3B", "pointSize": 6}), {}, "Gauge Stations")
    return m

stations_to_ee_fc(stations_df)

Convert any stations DataFrame to an ee.FeatureCollection of points.

Works for any stations_df meeting :data:REQUIRED_STATION_COLUMNS — not specific to the WA network. Extra columns are copied through as feature properties.

Source code in savana/rainfall/stations.py

def stations_to_ee_fc(stations_df):
    """Convert any stations DataFrame to an ``ee.FeatureCollection`` of points.

    Works for any ``stations_df`` meeting :data:`REQUIRED_STATION_COLUMNS`
    — not specific to the WA network. Extra columns are copied through
    as feature properties.
    """
    import ee

    _validate_stations_df(stations_df)

    features = []
    extra_cols = [c for c in stations_df.columns if c not in ("lon", "lat")]
    for _, row in stations_df.iterrows():
        props = {c: row[c] for c in extra_cols}
        features.append(
            ee.Feature(ee.Geometry.Point([float(row.lon), float(row.lat)]), props)
        )
    return ee.FeatureCollection(features)

savana.rainfall.zones

Ecological/climatic zone construction and station assignment.

No packaged shapefile, zones are built, from whatever base regions and split logic you give :func:build_zones_from_bands. This is a direct, generalised port of the author's GEE zone-delineation script: 3 named base climatic regions, each optionally split by a latitude band into one or more final ecological zones. The West Africa 5-zone scheme (:data:config.DEFAULT_ZONE_DEFS_WA + :data:config.DEFAULT_ZONE_BASE_ASSETS_WA) is just one configuration of this generic builder, not special-cased logic, a different region, a different number of base regions, or no latitude splitting at all, all use the same function.

Three ways to get zone geometry, in order of generality:

1. :func:build_zones_from_bands, build zones yourself from any named base regions (EE assets, or ee.FeatureCollection/ee.Geometry objects you already have) plus your own split rules. Fully general, works for any region, any number of zones, any split logic (or none).
2. :func:load_zones_from_file, you already have a zone boundary file (shapefile, GeoJSON, GeoPackage, ...) from QGIS or elsewhere. Fully general, no Earth Engine involved at all.
3. :func:single_region_zone, you don't want zone stratification at all, just one study-area boundary. Wraps any boundary (a file path, a geojson dict, an ee.Geometry, or a bounding box) into a one-zone table so the rest of the package (which only ever asks "is this station's zone-name X") doesn't need a special no-zone code path.

:func:assign_zones then attaches a zone label to a stations DataFrame from any of the above, or falls back to a documented latitude-band heuristic if no zone geometry is available at all.

assign_zones(stations_df, zones_gdf=None, zones_fc=None, name_field=None, use_default_if_none=False)

Add a zone column to stations_df.

Parameters:

Name	Type	Description	Default
stations_df		any DataFrame with station_id, lon, lat.	required
zones_gdf		a local geopandas.GeoDataFrame (from :func:load_zones_from_file, :func:single_region_zone, or :func:zones_fc_to_gdf), joined locally via geopandas.	None
zones_fc		a live ee.FeatureCollection (from :func:build_zones_from_bands or :func:default_wa_zones) , joined via Earth Engine (filterBounds per station), no geopandas required.	None
name_field	str | None	property/column holding the zone name. Defaults to :data:config.DEFAULT_ZONE_NAME_FIELD ("zone_name").	None
use_default_if_none	bool	if True and neither zones_gdf nor zones_fc is given, attempts :func:default_wa_zones, only useful if you're the author (or have access to the same EE assets). False by default, since that default is not portable to other users/regions, pass your own zones_gdf/zones_fc instead, or accept the coarse latitude fallback.	False

Returns:

Type	Description
	Copy of stations_df with a new zone column. Any station
	that can't be matched to a real zone polygon falls back to the
	latitude-band heuristic, with a printed warning, this always
	returns a usable zone column, never leaves it null.

Source code in savana/rainfall/zones.py

def assign_zones(
    stations_df,
    zones_gdf=None,
    zones_fc=None,
    name_field: str | None = None,
    use_default_if_none: bool = False,
):
    """Add a ``zone`` column to ``stations_df``.

    Args:
        stations_df: any DataFrame with ``station_id, lon, lat``.
        zones_gdf: a local ``geopandas.GeoDataFrame`` (from
            :func:`load_zones_from_file`, :func:`single_region_zone`, or
            :func:`zones_fc_to_gdf`), joined locally via geopandas.
        zones_fc: a live ``ee.FeatureCollection`` (from
            :func:`build_zones_from_bands` or :func:`default_wa_zones`)
           , joined via Earth Engine (``filterBounds`` per station), no
            geopandas required.
        name_field: property/column holding the zone name. Defaults to
            :data:`config.DEFAULT_ZONE_NAME_FIELD` (``"zone_name"``).
        use_default_if_none: if True and neither ``zones_gdf`` nor
            ``zones_fc`` is given, attempts :func:`default_wa_zones`,
            only useful if you're the author (or have access to the
            same EE assets). False by default, since that default is
            not portable to other users/regions, pass your own
            ``zones_gdf``/``zones_fc`` instead, or accept the coarse
            latitude fallback.

    Returns:
        Copy of ``stations_df`` with a new ``zone`` column. Any station
        that can't be matched to a real zone polygon falls back to the
        latitude-band heuristic, with a printed warning, this always
        returns a usable ``zone`` column, never leaves it null.
    """
    stations_df = stations_df.copy()
    name_field = name_field or config.DEFAULT_ZONE_NAME_FIELD

    if zones_gdf is None and zones_fc is None and use_default_if_none:
        try:
            zones_fc = default_wa_zones()
        except Exception as exc:  # noqa: BLE001
            print(
                f"  \u26a0  Could not build default WA zones "
                f"({type(exc).__name__}: {exc}). Falling back to latitude "
                f"bands for all stations."
            )

    if zones_fc is not None:
        try:
            import ee

            from .stations import stations_to_ee_fc

            station_fc = stations_to_ee_fc(stations_df)

            def _tag(feature):
                pt = feature.geometry()
                match = zones_fc.filterBounds(pt).first()
                return feature.set(
                    "zone",
                    ee.Algorithms.If(match, ee.Feature(match).get(name_field), None),
                )

            tagged = station_fc.map(_tag).getInfo()
            zone_by_id = {
                f["properties"]["station_id"]: f["properties"].get("zone")
                for f in tagged["features"]
            }
            stations_df["zone"] = stations_df["station_id"].map(zone_by_id)
        except Exception as exc:  # noqa: BLE001
            print(
                f"  \u26a0  EE zone join failed ({type(exc).__name__}: {exc}), "
                f"falling back to latitude bands."
            )
            stations_df["zone"] = None

    elif zones_gdf is not None:
        try:
            import geopandas as gpd
            from shapely.geometry import Point

            pts = gpd.GeoDataFrame(
                stations_df,
                geometry=[Point(xy) for xy in zip(stations_df.lon, stations_df.lat)],
                crs=zones_gdf.crs or "EPSG:4326",
            )
            joined = gpd.sjoin(pts, zones_gdf[[name_field, "geometry"]], how="left")
            stations_df["zone"] = joined[name_field].values
        except ImportError:
            print(
                "  \u26a0  geopandas/shapely not installed, using latitude-band "
                'fallback. Install with `pip install "savana[rainfall]"`.'
            )
            stations_df["zone"] = None
    else:
        stations_df["zone"] = None

    unmatched = stations_df["zone"].isna()
    if unmatched.any():
        ids = stations_df.loc[unmatched, "station_id"].tolist()
        if zones_gdf is not None or zones_fc is not None:
            print(
                f"  \u26a0  {len(ids)} station(s) had no real zone match, using "
                f"latitude fallback: {ids}"
            )
        stations_df.loc[unmatched, "zone"] = stations_df.loc[unmatched, "lat"].apply(
            _assign_zone_by_latitude
        )

    return stations_df

build_zones_from_bands(base_zones, zone_defs, bounds=None)

Build a final zones ee.FeatureCollection from named base regions, each optionally split by a latitude band.

This is the generic version of the GEE script's whole zone-building pipeline (its Sections 1, 3, 4), nothing here is specific to West Africa or to exactly 3 base regions / 5 output zones.

Parameters:

Name	Type	Description	Default
base_zones	dict	{name: source} where each source is an EE asset ID (str), an already-loaded ee.FeatureCollection, or an ee.Geometry. These are your raw regions before any splitting, e.g. {"Sahelian": "projects/x/assets/y", ...} for the WA case, or e.g. {"my_watershed": my_geometry} for a single custom region.	required
zone_defs	list[dict]	list of dicts, each describing one output zone: - zone_name (required): the label written to the zone_name property (or whatever config.DEFAULT_ZONE_NAME_FIELD is). - source_zone (required): key into base_zones this output zone is derived from. - lat_min, lat_max (optional): latitude band to clip to. Omit both (or span the full region) to pass the source region through unmodified, the pattern for "one base region = one output zone, no further splitting". - any other keys (zone_id, color_hex, rainfall_mm_yr, notes, ...) are copied through as feature properties, same as the GEE script's ZONE_DEFS.	required
bounds	tuple[float, float, float, float] | None	(min_lon, min_lat, max_lon, max_lat) used only to clip latitude bands to a sensible extent. Defaults to a generous global-ish box if not given, set this to your own study area's bounds for anything other than West Africa.	None

Returns:

Type	Description
	ee.FeatureCollection of the final zone polygons, one
	feature per zone_defs entry, with all of that entry's keys
	as properties.

Source code in savana/rainfall/zones.py

def build_zones_from_bands(
    base_zones: dict,
    zone_defs: list[dict],
    bounds: tuple[float, float, float, float] | None = None,
):
    """Build a final zones ``ee.FeatureCollection`` from named base
    regions, each optionally split by a latitude band.

    This is the generic version of the GEE script's whole zone-building
    pipeline (its Sections 1, 3, 4), nothing here is specific to West
    Africa or to exactly 3 base regions / 5 output zones.

    Args:
        base_zones: ``{name: source}`` where each ``source`` is an EE
            asset ID (str), an already-loaded ``ee.FeatureCollection``,
            or an ``ee.Geometry``. These are your raw regions before any
            splitting, e.g. ``{"Sahelian": "projects/x/assets/y", ...}``
            for the WA case, or e.g. ``{"my_watershed": my_geometry}``
            for a single custom region.
        zone_defs: list of dicts, each describing one output zone:
            - ``zone_name`` (required): the label written to the
              ``zone_name`` property (or whatever ``config.DEFAULT_ZONE_NAME_FIELD``
              is).
            - ``source_zone`` (required): key into ``base_zones`` this
              output zone is derived from.
            - ``lat_min``, ``lat_max`` (optional): latitude band to
              clip to. Omit both (or span the full region) to pass the
              source region through unmodified, the pattern for "one
              base region = one output zone, no further splitting".
            - any other keys (``zone_id``, ``color_hex``,
              ``rainfall_mm_yr``, notes, ...) are copied through as
              feature properties, same as the GEE script's ZONE_DEFS.
        bounds: ``(min_lon, min_lat, max_lon, max_lat)`` used only to
            clip latitude bands to a sensible extent. Defaults to a
            generous global-ish box if not given, set this to your own
            study area's bounds for anything other than West Africa.

    Returns:
        ``ee.FeatureCollection`` of the final zone polygons, one
        feature per ``zone_defs`` entry, with all of that entry's keys
        as properties.
    """
    import ee

    bounds = bounds or (-180.0, -85.0, 180.0, 85.0)

    loaded_bases = {}
    for name, source in base_zones.items():
        if isinstance(source, str):
            loaded_bases[name] = (
                ee.FeatureCollection(source).geometry().dissolve(ee.ErrorMargin(100))
            )
        elif isinstance(source, ee.FeatureCollection):
            loaded_bases[name] = source.geometry().dissolve(ee.ErrorMargin(100))
        else:  # already an ee.Geometry
            loaded_bases[name] = source

    features = []
    for zdef in zone_defs:
        if "zone_name" not in zdef or "source_zone" not in zdef:
            raise ValueError(
                f"zone_def is missing required key(s): {zdef}. "
                f"Every zone_def needs at least 'zone_name' and 'source_zone'."
            )
        if zdef["source_zone"] not in loaded_bases:
            raise ValueError(
                f"zone_def {zdef['zone_name']!r} references source_zone "
                f"{zdef['source_zone']!r}, not found in base_zones "
                f"({sorted(base_zones)})."
            )

        source_geom = loaded_bases[zdef["source_zone"]]
        lat_min, lat_max = zdef.get("lat_min"), zdef.get("lat_max")
        if lat_min is not None and lat_max is not None:
            band = _lat_band(lat_min, lat_max, bounds)
            geometry = source_geom.intersection(band, ee.ErrorMargin(100))
        else:
            geometry = source_geom

        props = {k: v for k, v in zdef.items() if k not in ("lat_min", "lat_max")}
        features.append(ee.Feature(geometry, props))

    return ee.FeatureCollection(features)

default_wa_zones(bounds=None)

Build the West Africa 5-zone scheme from the author's own base EE assets (:data:config.DEFAULT_ZONE_BASE_ASSETS_WA) using :data:config.DEFAULT_ZONE_DEFS_WA.

This is just the WA study's configuration of :func:build_zones_from_bands, call that function directly with your own base_zones/zone_defs for a different region.

Requires the 3 base assets to actually exist and be readable by the caller's EE account, they're the author's own uploaded shapefiles, not a public dataset. If you're not the author, either ask for read access, upload your own copies and pass your own base_zones dict, or use a completely different region's data.

Source code in savana/rainfall/zones.py

def default_wa_zones(bounds=None):
    """Build the West Africa 5-zone scheme from the author's own base
    EE assets (:data:`config.DEFAULT_ZONE_BASE_ASSETS_WA`) using
    :data:`config.DEFAULT_ZONE_DEFS_WA`.

    This is just the WA study's *configuration* of
    :func:`build_zones_from_bands`, call that function directly with
    your own ``base_zones``/``zone_defs`` for a different region.

    Requires the 3 base assets to actually exist and be readable by the
    caller's EE account, they're the author's own uploaded shapefiles,
    not a public dataset. If you're not the author, either ask for read
    access, upload your own copies and pass your own
    ``base_zones`` dict, or use a completely different region's data.
    """
    return build_zones_from_bands(
        config.DEFAULT_ZONE_BASE_ASSETS_WA,
        config.DEFAULT_ZONE_DEFS_WA,
        bounds=bounds or config.DEFAULT_ZONE_BOUNDS_WA,
    )

export_zones(zones_fc, asset_id=None, drive_folder=None, drive_description='ecological_zones', file_format='GeoJSON')

Export a built zones FeatureCollection, port of the GEE script's three export buttons (asset / GeoJSON / Shapefile), as background ee.batch tasks rather than a UI panel.

Parameters:

Name	Type	Description	Default
zones_fc		from :func:build_zones_from_bands (or any ee.FeatureCollection).	required
asset_id	str | None	if given, submits an Export.table.toAsset task so the zones can be reloaded quickly later via ee.FeatureCollection(asset_id) instead of rebuilding from base regions every time.	None
drive_folder, drive_description, file_format		if drive_folder is given, submits an Export.table.toDrive task (file_format one of "GeoJSON", "SHP", "CSV", ...).	required

Returns:

Type	Description
	list of submitted ee.batch.Task objects (already started,
	check task.status() for progress, same as any other GEE
	batch export).

Source code in savana/rainfall/zones.py

def export_zones(
    zones_fc,
    asset_id: str | None = None,
    drive_folder: str | None = None,
    drive_description: str = "ecological_zones",
    file_format: str = "GeoJSON",
):
    """Export a built zones FeatureCollection, port of the GEE script's
    three export buttons (asset / GeoJSON / Shapefile), as background
    ``ee.batch`` tasks rather than a UI panel.

    Args:
        zones_fc: from :func:`build_zones_from_bands` (or any
            ``ee.FeatureCollection``).
        asset_id: if given, submits an ``Export.table.toAsset`` task so
            the zones can be reloaded quickly later via
            ``ee.FeatureCollection(asset_id)`` instead of rebuilding
            from base regions every time.
        drive_folder, drive_description, file_format: if
            ``drive_folder`` is given, submits an
            ``Export.table.toDrive`` task (``file_format`` one of
            ``"GeoJSON"``, ``"SHP"``, ``"CSV"``, ...).

    Returns:
        list of submitted ``ee.batch.Task`` objects (already started,
        check ``task.status()`` for progress, same as any other GEE
        batch export).
    """
    import ee

    tasks = []
    if asset_id:
        task = ee.batch.Export.table.toAsset(
            collection=zones_fc, description=drive_description, assetId=asset_id
        )
        task.start()
        tasks.append(task)
        print(f"  Submitted asset export: {asset_id}")

    if drive_folder:
        task = ee.batch.Export.table.toDrive(
            collection=zones_fc,
            description=drive_description,
            folder=drive_folder,
            fileFormat=file_format,
        )
        task.start()
        tasks.append(task)
        print(f"  Submitted Drive export ({file_format}) to folder {drive_folder!r}")

    if not tasks:
        print("  Nothing submitted, pass asset_id and/or drive_folder.")
    return tasks

load_zones_from_file(path)

Load your own zone boundaries from any vector file geopandas can read (shapefile, GeoJSON, GeoPackage, ...). Fully general, for a user who already has zone geometry from QGIS or elsewhere and doesn't need :func:build_zones_from_bands at all.

Source code in savana/rainfall/zones.py

def load_zones_from_file(path):
    """Load your own zone boundaries from any vector file geopandas can
    read (shapefile, GeoJSON, GeoPackage, ...). Fully general, for a
    user who already has zone geometry from QGIS or elsewhere and
    doesn't need :func:`build_zones_from_bands` at all.
    """
    import geopandas as gpd

    return gpd.read_file(path)

single_region_zone(boundary, zone_name='Study Area')

Wrap one boundary as a one-row zones table, for a user who wants a specific area of interest but no zone stratification.

Parameters:

Name	Type	Description	Default
boundary		any of, a local vector file path (shapefile, GeoJSON, ...), a GeoJSON-like dict, an ee.Geometry, or a (min_lon, min_lat, max_lon, max_lat) bounding box.	required
zone_name	str	the single zone label everything in boundary will be assigned.	'Study Area'

Returns:

Type	Description
	A geopandas.GeoDataFrame with one row (usable with
	func:assign_zones's zones_gdf=), unless boundary is
	an ee.Geometry, in which case an ee.FeatureCollection is
	returned instead (usable with zones_fc=).

Source code in savana/rainfall/zones.py

def single_region_zone(boundary, zone_name: str = "Study Area"):
    """Wrap one boundary as a one-row zones table, for a user who
    wants a specific area of interest but no zone stratification.

    Args:
        boundary: any of, a local vector file path (shapefile,
            GeoJSON, ...), a GeoJSON-like dict, an ``ee.Geometry``, or
            a ``(min_lon, min_lat, max_lon, max_lat)`` bounding box.
        zone_name: the single zone label everything in ``boundary``
            will be assigned.

    Returns:
        A ``geopandas.GeoDataFrame`` with one row (usable with
        :func:`assign_zones`'s ``zones_gdf=``), unless ``boundary`` is
        an ``ee.Geometry``, in which case an ``ee.FeatureCollection`` is
        returned instead (usable with ``zones_fc=``).
    """
    try:
        import ee

        if isinstance(boundary, ee.Geometry):
            props = {config.DEFAULT_ZONE_NAME_FIELD: zone_name}
            return ee.FeatureCollection([ee.Feature(boundary, props)])
    except ImportError:
        pass  # ee not installed -> boundary can't be an ee.Geometry, fall through

    import geopandas as gpd

    if isinstance(boundary, (str,)):
        gdf = gpd.read_file(boundary)
        geom = gdf.geometry.unary_union
    elif isinstance(boundary, dict):
        from shapely.geometry import shape

        geom = shape(boundary)
    elif isinstance(boundary, (tuple, list)) and len(boundary) == 4:
        from shapely.geometry import box

        geom = box(*boundary)
    else:
        raise ValueError(
            f"Unrecognised boundary type: {type(boundary)}. Expected a file "
            f"path, GeoJSON dict, ee.Geometry, or (min_lon, min_lat, max_lon, "
            f"max_lat) bounding box."
        )

    return gpd.GeoDataFrame(
        {config.DEFAULT_ZONE_NAME_FIELD: [zone_name]}, geometry=[geom], crs="EPSG:4326"
    )

zone_areas_km2(zones_fc, name_field=None)

Add an area_km2 property to every feature, port of the GEE script's area-reporting section. Returns the FeatureCollection with the extra property; call .getInfo() or use :func:zones_fc_to_gdf to inspect it locally.

Source code in savana/rainfall/zones.py

def zone_areas_km2(zones_fc, name_field: str | None = None):
    """Add an ``area_km2`` property to every feature, port of the GEE
    script's area-reporting section. Returns the FeatureCollection with
    the extra property; call ``.getInfo()`` or use
    :func:`zones_fc_to_gdf` to inspect it locally.
    """
    import ee

    name_field = name_field or config.DEFAULT_ZONE_NAME_FIELD

    def _add_area(f):
        area_km2 = f.geometry().area(ee.ErrorMargin(100)).divide(1e6).round()
        return f.set("area_km2", area_km2)

    zones_with_area = zones_fc.map(_add_area)
    info = zones_with_area.select([name_field, "area_km2"]).getInfo()
    for f in info["features"]:
        p = f["properties"]
        print(f"  {p.get(name_field)}: {p.get('area_km2'):,.0f} km2")
    return zones_with_area

zones_fc_to_gdf(zones_fc, name_field=None)

Pull a (typically small, a handful of zone polygons) EE FeatureCollection down to a local geopandas.GeoDataFrame, for use with :func:assign_zones's local-join path, or for saving to a file yourself.

Source code in savana/rainfall/zones.py

def zones_fc_to_gdf(zones_fc, name_field: str | None = None):
    """Pull a (typically small, a handful of zone polygons) EE
    FeatureCollection down to a local ``geopandas.GeoDataFrame``, for
    use with :func:`assign_zones`'s local-join path, or for saving to a
    file yourself.
    """
    import geopandas as gpd
    from shapely.geometry import shape

    info = zones_fc.getInfo()
    rows = [f["properties"] for f in info["features"]]
    geoms = [shape(f["geometry"]) for f in info["features"]]
    return gpd.GeoDataFrame(rows, geometry=geoms, crs="EPSG:4326")

savana.rainfall.ingestion

Precipitation product ingestion, harmonise any product catalogue to a common monthly mean mm/day ImageCollection.

Every function accepts a products dict (see :data:savana.rainfall.config.DEFAULT_PRODUCTS for the required shape) and a roi, so this works for a different product catalogue or a different region, not just the WA six-product/study-area default.

Nothing here calls ee.Initialize(), that's the caller's responsibility (see :mod:savana.ee_init, reused as-is), consistent with the rest of savana never initialising EE as a side effect of import.

build_roi(stations_df=None, bounds=None, buffer_deg=2.0)

Build an ee.Geometry region of interest.

Parameters:

Name	Type	Description	Default
stations_df		if given (and bounds is None), the ROI is the bounding box of the stations plus buffer_deg on each side, works for any station set, anywhere.	None
bounds		explicit (min_lon, min_lat, max_lon, max_lat), takes priority over stations_df if given.	None
buffer_deg	float	degrees of padding added around the station bbox.	2.0

If neither is given, falls back to the West Africa study bounds used in the manuscript (5-25N, 20W-15E).

Source code in savana/rainfall/ingestion.py

def build_roi(stations_df=None, bounds=None, buffer_deg: float = 2.0):
    """Build an ``ee.Geometry`` region of interest.

    Args:
        stations_df: if given (and ``bounds`` is None), the ROI is the
            bounding box of the stations plus ``buffer_deg`` on each
            side, works for any station set, anywhere.
        bounds: explicit ``(min_lon, min_lat, max_lon, max_lat)``, takes
            priority over ``stations_df`` if given.
        buffer_deg: degrees of padding added around the station bbox.

    If neither is given, falls back to the West Africa study bounds
    used in the manuscript (5-25N, 20W-15E).
    """
    import ee

    if bounds is not None:
        min_lon, min_lat, max_lon, max_lat = bounds
    elif stations_df is not None and len(stations_df) > 0:
        min_lon = float(stations_df.lon.min()) - buffer_deg
        max_lon = float(stations_df.lon.max()) + buffer_deg
        min_lat = float(stations_df.lat.min()) - buffer_deg
        max_lat = float(stations_df.lat.max()) + buffer_deg
    else:
        min_lon, min_lat, max_lon, max_lat = (-20.0, 5.0, 15.0, 25.0)

    return ee.Geometry.Rectangle([min_lon, min_lat, max_lon, max_lat])

export_merra2_yearly_assets(years, roi, asset_folder, products=None)

Export one daily-aggregated MERRA-2 asset per year to asset_folder, as a background EE batch task per year.

Splitting the export by year keeps each task well under GEE's per-request compute/timeout limits. Returns the list of submitted ee.batch.Task objects, check task.status() for progress; this can take hours for a full 20-year run and is meant to be fire-and-forget, not awaited synchronously.

Source code in savana/rainfall/ingestion.py

def export_merra2_yearly_assets(
    years: list[int], roi, asset_folder: str, products: dict | None = None
):
    """Export one daily-aggregated MERRA-2 asset per year to
    ``asset_folder``, as a background EE batch task per year.

    Splitting the export by year keeps each task well under GEE's
    per-request compute/timeout limits. Returns the list of submitted
    ``ee.batch.Task`` objects, check ``task.status()`` for progress;
    this can take hours for a full 20-year run and is meant to be
    fire-and-forget, not awaited synchronously.
    """
    import ee

    products = products if products is not None else config.DEFAULT_PRODUCTS
    spec = products["MERRA2"]
    tasks = []
    for year in years:
        start, end = f"{year}-01-01", f"{year + 1}-01-01"
        ic = (
            ee.ImageCollection(spec["collection"])
            .filterDate(start, end)
            .filterBounds(roi)
        )
        days = ee.List.sequence(
            0, ee.Date(end).difference(ee.Date(start), "day").subtract(1)
        )

        def _day_mean(d, _ic=ic, _spec=spec, _start=start):
            d0 = ee.Date(_start).advance(d, "day")
            d1 = d0.advance(1, "day")
            img = (
                _ic.filterDate(d0, d1)
                .select(_spec["band"])
                .mean()
                .multiply(_spec["scale_factor"])
                .rename("precip_mm_day")
            )
            return img.set("system:time_start", d0.millis())

        daily_stack = ee.ImageCollection(days.map(_day_mean)).toBands()
        asset_id = f"{asset_folder}/merra2_daily_{year}"
        task = ee.batch.Export.image.toAsset(
            image=daily_stack,
            description=f"merra2_daily_{year}",
            assetId=asset_id,
            region=roi,
            scale=config.DEFAULT_TARGET_RESOLUTION_M,
            maxPixels=1e10,
        )
        task.start()
        tasks.append(task)
        print(f"  Submitted export task: {asset_id}")
    return tasks

load_all_products(start, end, roi=None, products=None, stations_df=None)

Load every product in products (default: all 6) as monthly mm/day ImageCollections.

Returns {product_name: ee.ImageCollection}. roi defaults to :func:build_roi from stations_df if not given.

Source code in savana/rainfall/ingestion.py

def load_all_products(
    start: str, end: str, roi=None, products: dict | None = None, stations_df=None
):
    """Load every product in ``products`` (default: all 6) as monthly
    mm/day ImageCollections.

    Returns ``{product_name: ee.ImageCollection}``. ``roi`` defaults to
    :func:`build_roi` from ``stations_df`` if not given.
    """
    products = products if products is not None else config.DEFAULT_PRODUCTS
    if roi is None:
        roi = build_roi(stations_df)

    out = {}
    for name in products:
        print(f"  Loading {name} ...")
        out[name] = load_product(name, start, end, roi, products)
    return out

load_merra2_from_assets(asset_folder, start, end, roi)

Rebuild the monthly mm/day MERRA-2 ImageCollection from pre-exported yearly daily-aggregate assets (see :func:export_merra2_yearly_assets).

Source code in savana/rainfall/ingestion.py

def load_merra2_from_assets(asset_folder: str, start: str, end: str, roi):
    """Rebuild the monthly mm/day MERRA-2 ImageCollection from
    pre-exported yearly daily-aggregate assets (see
    :func:`export_merra2_yearly_assets`).
    """
    import ee

    start_year = int(start[:4])
    end_year = int(end[:4])
    daily_images = []
    for year in range(start_year, end_year + 1):
        asset_id = f"{asset_folder}/merra2_daily_{year}"
        try:
            stack = ee.Image(asset_id)
        except Exception as exc:  # noqa: BLE001
            print(f"  \u26a0  Missing MERRA-2 asset for {year}: {exc}")
            continue
        band_names = stack.bandNames()
        n_bands = band_names.size().getInfo()
        for i in range(n_bands):
            band = ee.String(band_names.get(i))
            day_offset = ee.Number.parse(band.slice(0, band.index("_")))
            date = ee.Date(f"{year}-01-01").advance(day_offset, "day")
            img = (
                stack.select([band])
                .rename("precip_mm_day")
                .set("system:time_start", date.millis())
            )
            daily_images.append(img)

    daily_ic = ee.ImageCollection(daily_images)
    return _monthly_from_daily(daily_ic, "precip_mm_day", roi, start, end)

load_product(name, start, end, roi, products=None)

Load one precipitation product as a monthly mean mm/day ee.ImageCollection, dispatched by its conversion type.

Parameters:

Name	Type	Description	Default
name	str	key into products (default :data:config.DEFAULT_PRODUCTS).	required
start, end		'YYYY-MM-DD', clipped to the product's real availability window (see :data:config.DEFAULT_PRODUCT_DATE_RANGES).	required
roi		ee.Geometry (see :func:build_roi).	required
products	dict | None	catalogue dict. Bring your own for a different product set, must have the shape of :data:config.DEFAULT_PRODUCTS.	None

Source code in savana/rainfall/ingestion.py

def load_product(name: str, start: str, end: str, roi, products: dict | None = None):
    """Load one precipitation product as a monthly mean mm/day
    ``ee.ImageCollection``, dispatched by its ``conversion`` type.

    Args:
        name: key into ``products`` (default
            :data:`config.DEFAULT_PRODUCTS`).
        start, end: 'YYYY-MM-DD', clipped to the product's real
            availability window (see
            :data:`config.DEFAULT_PRODUCT_DATE_RANGES`).
        roi: ``ee.Geometry`` (see :func:`build_roi`).
        products: catalogue dict. Bring your own for a different
            product set, must have the shape of
            :data:`config.DEFAULT_PRODUCTS`.
    """
    products = products if products is not None else config.DEFAULT_PRODUCTS
    if name not in products:
        raise ValueError(f"Unknown product {name!r}. Known: {sorted(products)}")

    spec = products[name]
    clipped_start, clipped_end = _clip_dates(name, start, end, products)

    if spec["native_temporal"] == "hourly":
        return _load_hourly_to_monthly(name, spec, roi, clipped_start, clipped_end)

    loader = _LOADERS.get(spec["conversion"])
    if loader is None:
        raise ValueError(f"{name}: unknown conversion type {spec['conversion']!r}")
    return loader(name, spec, roi, clipped_start, clipped_end)

savana.rainfall.extraction

Point-sample precipitation products at gauge stations, and merge with observations into the long-format table :mod:.validation consumes.

Works for any stations_df, not just the WA 16, extraction is purely a function of whatever station coordinates you give it.

extract_all_products(products_ic, stations_df, cache_dir=None)

Extract every product in products_ic (from :func:savana.rainfall.ingestion.load_all_products) at every station, and stack into one long-format DataFrame.

Parameters:

Name	Type	Description	Default
products_ic	dict	{name: ee.ImageCollection}.	required
stations_df		any stations DataFrame.	required
cache_dir		if given, each product's extraction is cached to cache_dir/precip_extraction_<NAME>.csv (matching the original per-product CSV workflow), a re-run with the same cache_dir reuses whatever's already there instead of re-extracting from Earth Engine, and any product missing from the cache is extracted and added to it. Delete the relevant CSV (or the whole folder) to force a fresh pull.	None

Returns:

Type	Description
	Long-format DataFrame: ``station_id, year, month, product,
	sim_mm_day``.

Source code in savana/rainfall/extraction.py

def extract_all_products(products_ic: dict, stations_df, cache_dir=None):
    """Extract every product in ``products_ic`` (from
    :func:`savana.rainfall.ingestion.load_all_products`) at every
    station, and stack into one long-format DataFrame.

    Args:
        products_ic: ``{name: ee.ImageCollection}``.
        stations_df: any stations DataFrame.
        cache_dir: if given, each product's extraction is cached to
            ``cache_dir/precip_extraction_<NAME>.csv`` (matching the
            original per-product CSV workflow), a re-run with the same
            ``cache_dir`` reuses whatever's already there instead of
            re-extracting from Earth Engine, and any product missing
            from the cache is extracted and added to it. Delete the
            relevant CSV (or the whole folder) to force a fresh pull.

    Returns:
        Long-format DataFrame: ``station_id, year, month, product,
        sim_mm_day``.
    """
    import pandas as pd

    frames = []
    for name, ic in products_ic.items():
        cache_path = (
            Path(cache_dir) / f"precip_extraction_{name}.csv" if cache_dir else None
        )
        if cache_path is not None and cache_path.exists():
            print(f"  Using cached extraction: {cache_path}")
            frames.append(pd.read_csv(cache_path))
            continue

        df = extract_product_at_stations(ic, stations_df, name)
        if cache_path is not None:
            cache_path.parent.mkdir(parents=True, exist_ok=True)
            df.to_csv(cache_path, index=False)
            print(f"  Cached: {cache_path}")
        frames.append(df)

    return pd.concat(frames, ignore_index=True)

extract_product_at_stations(product_ic, stations_df, product_name)

Point-sample one monthly mm/day ee.ImageCollection at every station in stations_df.

Returns a long-format pandas.DataFrame: station_id, year, month, product, sim_mm_day.

Source code in savana/rainfall/extraction.py

def extract_product_at_stations(product_ic, stations_df, product_name: str):
    """Point-sample one monthly mm/day ``ee.ImageCollection`` at every
    station in ``stations_df``.

    Returns a long-format ``pandas.DataFrame``:
    ``station_id, year, month, product, sim_mm_day``.
    """
    import ee
    import pandas as pd

    from .stations import stations_to_ee_fc

    fc = stations_to_ee_fc(stations_df)

    def _sample_image(img):
        date = ee.Date(img.get("system:time_start"))
        sampled = img.reduceRegions(
            collection=fc,
            reducer=ee.Reducer.first(),
            scale=config.DEFAULT_TARGET_RESOLUTION_M,
        )
        return sampled.map(
            lambda f: f.set(
                {
                    "year": date.get("year"),
                    "month": date.get("month"),
                }
            )
        )

    sampled_fc = product_ic.map(_sample_image).flatten()
    info = sampled_fc.getInfo()

    rows = []
    for f in info["features"]:
        p = f["properties"]
        rows.append(
            {
                "station_id": p.get("station_id"),
                "year": int(p["year"]),
                "month": int(p["month"]),
                "product": product_name,
                "sim_mm_day": p.get("first"),
            }
        )
    df = pd.DataFrame(rows)
    print(f"  Extracted {product_name}: {len(df):,} station-months")
    return df

merge_with_observations(sim_long_df, obs_df, stations_df=None)

Merge long-format simulated values with observations, optionally attaching station metadata (including zone if already assigned).

Parameters:

Name	Type	Description	Default
sim_long_df		from :func:extract_all_products, columns station_id, year, month, product, sim_mm_day.	required
obs_df		from :mod:.stations, columns station_id, year, month, obs_mm_day.	required
stations_df		optional, to bring along zone (from :func:savana.rainfall.zones.assign_zones) or any other station attribute, joined on station_id.	None

Returns:

Type	Description
	Long-format DataFrame ready for :mod:.validation:
	station_id, year, month, product, sim_mm_day, obs_mm_day
	(+ any joined station columns).

Source code in savana/rainfall/extraction.py

def merge_with_observations(sim_long_df, obs_df, stations_df=None):
    """Merge long-format simulated values with observations, optionally
    attaching station metadata (including ``zone`` if already assigned).

    Args:
        sim_long_df: from :func:`extract_all_products`, columns
            ``station_id, year, month, product, sim_mm_day``.
        obs_df: from :mod:`.stations`, columns
            ``station_id, year, month, obs_mm_day``.
        stations_df: optional, to bring along ``zone`` (from
            :func:`savana.rainfall.zones.assign_zones`) or any other
            station attribute, joined on ``station_id``.

    Returns:
        Long-format DataFrame ready for :mod:`.validation`:
        ``station_id, year, month, product, sim_mm_day, obs_mm_day``
        (+ any joined station columns).
    """
    merged = sim_long_df.merge(obs_df, on=["station_id", "year", "month"], how="inner")

    if stations_df is not None:
        extra_cols = [c for c in stations_df.columns if c != "station_id"]
        merged = merged.merge(
            stations_df[["station_id"] + extra_cols], on="station_id", how="left"
        )

    n_before = len(sim_long_df)
    n_after = len(merged)
    if n_after < n_before:
        print(
            f"  Merge: {n_after:,}/{n_before:,} simulated station-months matched "
            f"an observation ({n_before - n_after:,} unmatched, check obs "
            f"coverage for those station-months)."
        )
    return merged

savana.rainfall.validation

Continuous and categorical validation metrics.

Two metric classes, matching the manuscript's dual-class framework:

• Continuous (:func:compute_continuous): bias, pbias, mae, rmse, r, r2, nse, kge, how well magnitude and pattern agree.
• Categorical (:func:compute_categorical): pod, far, csi, ets, freq_bias, how well wet/dry events are detected above a threshold.

Both take plain obs/sim array-likes, so they work regardless of which stations, products, or zones produced them. The four levels of spatial/temporal aggregation used in the manuscript (per-station, per-zone, per-season, pooled) are all just different group_cols to the single :func:validate_grouped function, there's no separate per-station/per-zone/per-season implementation to keep in sync.

add_season_column(merged_df, month_col='month')

Add a season column (DJF/MAM/JJA/SON) derived from month_col.

Source code in savana/rainfall/validation.py

def add_season_column(merged_df, month_col: str = "month"):
    """Add a ``season`` column (DJF/MAM/JJA/SON) derived from ``month_col``."""
    merged_df = merged_df.copy()
    merged_df["season"] = merged_df[month_col].map(_SEASON_MONTHS)
    return merged_df

compute_all_metrics(obs, sim, threshold=None)

Continuous + categorical metrics for one obs/sim pair, merged.

Source code in savana/rainfall/validation.py

def compute_all_metrics(obs, sim, threshold: float | None = None) -> dict:
    """Continuous + categorical metrics for one obs/sim pair, merged."""
    out = compute_continuous(obs, sim)
    out.update(compute_categorical(obs, sim, threshold))
    return out

compute_categorical(obs, sim, threshold=None)

Categorical wet/dry detection metrics from a 2x2 contingency table.

A record is "wet" if the value >= threshold (default :data:config.DEFAULT_RAIN_THRESHOLD_MM_DAY, the WMO standard of 1.0 mm/day), else "dry".

Source code in savana/rainfall/validation.py

def compute_categorical(obs, sim, threshold: float | None = None) -> dict:
    """Categorical wet/dry detection metrics from a 2x2 contingency table.

    A record is "wet" if the value >= ``threshold`` (default
    :data:`config.DEFAULT_RAIN_THRESHOLD_MM_DAY`, the WMO standard of
    1.0 mm/day), else "dry".
    """
    import numpy as np
    import pandas as pd

    threshold = (
        threshold if threshold is not None else config.DEFAULT_RAIN_THRESHOLD_MM_DAY
    )

    o = pd.Series(obs).astype(float)
    s = pd.Series(sim).astype(float)
    mask = o.notna() & s.notna()
    o, s = o[mask].to_numpy(), s[mask].to_numpy()
    n = len(o)

    keys = ["pod", "far", "csi", "ets", "freq_bias"]
    count_keys = ["hits", "misses", "false_alarms", "correct_negatives"]
    if n < 1:
        return {"n": 0, **{k: float("nan") for k in keys + count_keys}}

    obs_wet, sim_wet = o >= threshold, s >= threshold
    hits = int(np.sum(obs_wet & sim_wet))
    misses = int(np.sum(obs_wet & ~sim_wet))
    false_alarms = int(np.sum(~obs_wet & sim_wet))
    correct_negatives = int(np.sum(~obs_wet & ~sim_wet))

    pod = hits / (hits + misses) if (hits + misses) > 0 else float("nan")
    far = (
        false_alarms / (hits + false_alarms)
        if (hits + false_alarms) > 0
        else float("nan")
    )
    csi_denom = hits + misses + false_alarms
    csi = hits / csi_denom if csi_denom > 0 else float("nan")

    total = hits + misses + false_alarms + correct_negatives
    hits_random = (
        (hits + misses) * (hits + false_alarms) / total if total > 0 else float("nan")
    )
    ets_denom = hits + misses + false_alarms - hits_random
    ets = (
        (hits - hits_random) / ets_denom
        if ets_denom not in (0, float("nan")) and not np.isnan(ets_denom)
        else float("nan")
    )

    freq_bias = (
        (hits + false_alarms) / (hits + misses) if (hits + misses) > 0 else float("nan")
    )

    return {
        "n": n,
        "threshold": threshold,
        "pod": pod,
        "far": far,
        "csi": csi,
        "ets": ets,
        "freq_bias": freq_bias,
        "hits": hits,
        "misses": misses,
        "false_alarms": false_alarms,
        "correct_negatives": correct_negatives,
    }

compute_continuous(obs, sim)

Continuous performance metrics for one obs/sim pair.

Rows where either obs or sim is missing are dropped before computation, matching the manuscript's paired-record approach. Returns {"n": 0, ...all-NaN} if fewer than 2 valid pairs remain (metrics like r/NSE/KGE are undefined below that).

Source code in savana/rainfall/validation.py

def compute_continuous(obs, sim) -> dict:
    """Continuous performance metrics for one obs/sim pair.

    Rows where either ``obs`` or ``sim`` is missing are dropped before
    computation, matching the manuscript's paired-record approach.
    Returns ``{"n": 0, ...all-NaN}`` if fewer than 2 valid pairs remain
    (metrics like r/NSE/KGE are undefined below that).
    """
    import numpy as np
    import pandas as pd

    o = pd.Series(obs).astype(float)
    s = pd.Series(sim).astype(float)
    mask = o.notna() & s.notna()
    o, s = o[mask].to_numpy(), s[mask].to_numpy()
    n = len(o)

    keys = ["bias", "pbias", "mae", "rmse", "r", "r2", "nse", "kge"]
    if n < 2:
        return {"n": n, **{k: float("nan") for k in keys}}

    obar, sbar = o.mean(), s.mean()
    diff = s - o

    bias = float(diff.mean())
    pbias = float(100.0 * diff.sum() / o.sum()) if o.sum() != 0 else float("nan")
    mae = float(np.abs(diff).mean())
    rmse = float(np.sqrt((diff**2).mean()))

    o_std, s_std = o.std(), s.std()
    if o_std == 0 or s_std == 0:
        r = float("nan")
    else:
        r = float(np.corrcoef(o, s)[0, 1])
    r2 = float(r**2) if not np.isnan(r) else float("nan")

    denom_nse = ((o - obar) ** 2).sum()
    nse = (
        float(1 - ((s - o) ** 2).sum() / denom_nse) if denom_nse != 0 else float("nan")
    )

    if np.isnan(r) or obar == 0 or o_std == 0:
        kge = float("nan")
    else:
        alpha = s_std / o_std
        beta = sbar / obar
        kge = float(1 - np.sqrt((r - 1) ** 2 + (alpha - 1) ** 2 + (beta - 1) ** 2))

    return {
        "n": n,
        "bias": bias,
        "pbias": pbias,
        "mae": mae,
        "rmse": rmse,
        "r": r,
        "r2": r2,
        "nse": nse,
        "kge": kge,
    }

rank_products(validation_df, metric='kge', group_cols=None)

Rank products within each group by a single metric (default KGE, the manuscript's primary ranking metric, see methods 2.4.1).

group_cols defaults to every column in validation_df except "product" and the metric columns, i.e. whatever grouping level the input DataFrame already represents (zone, station, season...). An empty result (e.g. pooled/overall validation with no zone column) ranks across the whole table as a single group, rather than failing.

Source code in savana/rainfall/validation.py

def rank_products(
    validation_df, metric: str = "kge", group_cols: list[str] | None = None
):
    """Rank products within each group by a single metric (default KGE,
    the manuscript's primary ranking metric, see methods 2.4.1).

    ``group_cols`` defaults to every column in ``validation_df`` except
    ``"product"`` and the metric columns, i.e. whatever grouping level
    the input DataFrame already represents (zone, station, season...).
    An empty result (e.g. pooled/overall validation with no zone column)
    ranks across the whole table as a single group, rather than failing.
    """
    df = validation_df.copy()
    if group_cols is None:
        non_metric = set(config.DEFAULT_METRICS_FLAT) | {
            "n",
            "threshold",
            "hits",
            "misses",
            "false_alarms",
            "correct_negatives",
            "product",
        }
        group_cols = [c for c in df.columns if c not in non_metric]

    ascending = metric in ("far", "rmse", "mae", "bias")
    if group_cols:
        df["rank"] = (
            df.groupby(group_cols)[metric]
            .rank(ascending=ascending, method="min")
            .astype(int)
        )
    else:
        df["rank"] = df[metric].rank(ascending=ascending, method="min").astype(int)
    sort_cols = group_cols + ["rank"]
    return df.sort_values(sort_cols).reset_index(drop=True)

validate_by_season(merged_df, threshold=None)

Metrics per (zone, season, product).

Source code in savana/rainfall/validation.py

def validate_by_season(merged_df, threshold=None):
    """Metrics per (zone, season, product)."""
    if "season" not in merged_df.columns:
        merged_df = add_season_column(merged_df)
    return validate_grouped(
        merged_df, ["zone", "season", "product"], threshold=threshold
    )

validate_by_station(merged_df, threshold=None)

Metrics per (station_id, product), manuscript's finest level.

Source code in savana/rainfall/validation.py

def validate_by_station(merged_df, threshold=None):
    """Metrics per (station_id, product), manuscript's finest level."""
    return validate_grouped(merged_df, ["station_id", "product"], threshold=threshold)

validate_by_zone(merged_df, threshold=None)

Metrics per (zone, product), the manuscript's primary analytical lens.

merged_df must already have a zone column (see :func:savana.rainfall.zones.assign_zones).

Source code in savana/rainfall/validation.py

def validate_by_zone(merged_df, threshold=None):
    """Metrics per (zone, product), the manuscript's primary analytical lens.

    ``merged_df`` must already have a ``zone`` column
    (see :func:`savana.rainfall.zones.assign_zones`).
    """
    if "zone" not in merged_df.columns:
        raise ValueError(
            "merged_df has no 'zone' column, run zones.assign_zones() on your "
            "stations_df and merge it in before calling validate_by_zone()."
        )
    return validate_grouped(merged_df, ["zone", "product"], threshold=threshold)

validate_grouped(merged_df, group_cols, obs_col='obs_mm_day', sim_col='sim_mm_day', threshold=None)

Compute continuous + categorical metrics for each group in merged_df.

merged_df must be long-format with one row per (station, year, month, product) and columns obs_col/sim_col (see :func:savana.rainfall.extraction.merge_with_observations).

This single function implements all four aggregation levels used in the manuscript, pass the group_cols that define the level:

• per-station: ["station_id", "product"]
• per-zone: ["zone", "product"]
• per-season: ["zone", "season", "product"]
• pooled/overall: ["product"]

Returns a DataFrame with one row per group, group_cols + all metrics.

Source code in savana/rainfall/validation.py

def validate_grouped(
    merged_df,
    group_cols: list[str],
    obs_col: str = "obs_mm_day",
    sim_col: str = "sim_mm_day",
    threshold: float | None = None,
):
    """Compute continuous + categorical metrics for each group in ``merged_df``.

    ``merged_df`` must be long-format with one row per
    (station, year, month, product) and columns ``obs_col``/``sim_col``
    (see :func:`savana.rainfall.extraction.merge_with_observations`).

    This single function implements all four aggregation levels used in
    the manuscript, pass the ``group_cols`` that define the level:

    - per-station:  ``["station_id", "product"]``
    - per-zone:     ``["zone", "product"]``
    - per-season:   ``["zone", "season", "product"]``
    - pooled/overall: ``["product"]``

    Returns a DataFrame with one row per group, group_cols + all metrics.
    """
    import pandas as pd

    rows = []
    for key, sub in merged_df.groupby(group_cols, dropna=False):
        key = key if isinstance(key, tuple) else (key,)
        metrics = compute_all_metrics(sub[obs_col], sub[sim_col], threshold)
        rows.append(dict(zip(group_cols, key), **metrics))
    return pd.DataFrame(rows)

validate_overall(merged_df, threshold=None)

Metrics per product, pooled across all stations/zones.

Source code in savana/rainfall/validation.py

def validate_overall(merged_df, threshold=None):
    """Metrics per product, pooled across all stations/zones."""
    return validate_grouped(merged_df, ["product"], threshold=threshold)

savana.rainfall.thresholds

Rain-detection threshold sensitivity analysis.

Categorical metrics (POD, FAR, CSI, ETS) depend on the wet/dry threshold used to classify a record, this matters most in dryland zones where near-zero rainfall makes categorical detection structurally unstable (see manuscript sections 2.5 and the Saharian zone note in :data:savana.rainfall.config.DEFAULT_ZONE_NOTES). This module sweeps :func:savana.rainfall.validation.compute_categorical across multiple thresholds rather than duplicating its logic.

stability_summary(threshold_df, metric='csi', group_cols=None)

Rank product/group robustness across the threshold sweep.

Returns one row per group with the metric's mean, std, and coefficient of variation across all swept thresholds, a low CV means the product's performance on that metric is stable regardless of exactly where the wet/dry line is drawn; a high CV flags threshold-sensitive rankings (per manuscript 2.5's structural instability finding in near-zero-rainfall zones).

Source code in savana/rainfall/thresholds.py

def stability_summary(
    threshold_df, metric: str = "csi", group_cols: list[str] | None = None
):
    """Rank product/group robustness across the threshold sweep.

    Returns one row per group with the metric's mean, std, and
    coefficient of variation across all swept thresholds, a low CV
    means the product's performance on that metric is stable regardless
    of exactly where the wet/dry line is drawn; a high CV flags
    threshold-sensitive rankings (per manuscript 2.5's structural
    instability finding in near-zero-rainfall zones).
    """
    if group_cols is None:
        group_cols = [c for c in ("zone", "product") if c in threshold_df.columns]

    agg = threshold_df.groupby(group_cols)[metric].agg(["mean", "std"]).reset_index()
    agg["cv"] = agg["std"] / agg["mean"].replace(0, float("nan"))
    return agg.sort_values(group_cols[:-1] + ["cv"] if len(group_cols) > 1 else "cv")

threshold_sensitivity(merged_df, thresholds=None, group_cols=None, obs_col='obs_mm_day', sim_col='sim_mm_day')

Categorical metrics at each of several rain-detection thresholds.

Parameters:

Name	Type	Description	Default
merged_df		long-format obs/sim data (see :func:savana.rainfall.extraction.merge_with_observations).	required
thresholds	list[float] | None	mm/day values to sweep. Defaults to :data:config.DEFAULT_THRESHOLD_SWEEP_MM_DAY (0.1, 0.5, 1.0, 2.0, 5.0 mm/day).	None
group_cols	list[str] | None	grouping level, e.g. ["zone", "product"] (default) or ["product"] for pooled.	None

Returns:

Type	Description
	DataFrame with one row per (group, threshold) combination.

Source code in savana/rainfall/thresholds.py

def threshold_sensitivity(
    merged_df,
    thresholds: list[float] | None = None,
    group_cols: list[str] | None = None,
    obs_col: str = "obs_mm_day",
    sim_col: str = "sim_mm_day",
):
    """Categorical metrics at each of several rain-detection thresholds.

    Args:
        merged_df: long-format obs/sim data (see
            :func:`savana.rainfall.extraction.merge_with_observations`).
        thresholds: mm/day values to sweep. Defaults to
            :data:`config.DEFAULT_THRESHOLD_SWEEP_MM_DAY`
            (0.1, 0.5, 1.0, 2.0, 5.0 mm/day).
        group_cols: grouping level, e.g. ``["zone", "product"]``
            (default) or ``["product"]`` for pooled.

    Returns:
        DataFrame with one row per (group, threshold) combination.
    """
    import pandas as pd

    thresholds = (
        thresholds if thresholds is not None else config.DEFAULT_THRESHOLD_SWEEP_MM_DAY
    )
    if group_cols is None:
        group_cols = ["zone", "product"] if "zone" in merged_df.columns else ["product"]

    rows = []
    for threshold in thresholds:
        for key, sub in merged_df.groupby(group_cols, dropna=False):
            key = key if isinstance(key, tuple) else (key,)
            metrics = validation.compute_categorical(
                sub[obs_col], sub[sim_col], threshold
            )
            rows.append(dict(zip(group_cols, key), **metrics))

    df = pd.DataFrame(rows)
    return df.sort_values(group_cols + ["threshold"]).reset_index(drop=True)

savana.rainfall.spatial

Pixel-wise spatial diagnostics, the scripted counterpart to the interactive GEE app's on-demand map layers (bias/correlation/trend/ agreement), so those ~20 exploratory analyses produce real exportable outputs instead of only living as live map interaction.

Every function takes an explicit reference ImageCollection rather than assuming GPCC gridded data is available, the manuscript's GPCC reference is point-station only (see :mod:.stations), so by default these functions compare products against each other (inter-product agreement) or use whichever gridded product you designate as the reference for a given call.

agreement_map(products_ic)

Inter-product agreement: pixel-wise standard deviation across all products' mean-period image, normalised by the ensemble mean (coefficient of variation). Low CV = products agree spatially.

Source code in savana/rainfall/spatial.py

def agreement_map(products_ic: dict):
    """Inter-product agreement: pixel-wise standard deviation across all
    products' mean-period image, normalised by the ensemble mean
    (coefficient of variation). Low CV = products agree spatially.
    """
    import ee

    means = [ic.select("precip_mm_day").mean() for ic in products_ic.values()]
    stack = ee.ImageCollection(means)
    mean_img = stack.mean().rename("ensemble_mean")
    std_img = stack.reduce(ee.Reducer.stdDev()).rename("ensemble_std")
    cv = std_img.divide(mean_img.max(0.01)).rename("agreement_cv")
    return mean_img, std_img, cv

bias_map(product_ic, reference_ic)

Mean pixel-wise bias (product - reference) over the full period, in mm/day.

Source code in savana/rainfall/spatial.py

def bias_map(product_ic, reference_ic):
    """Mean pixel-wise bias (product - reference) over the full period,
    in mm/day.
    """
    p_mean = product_ic.select("precip_mm_day").mean().rename("product_mean")
    r_mean = reference_ic.select("precip_mm_day").mean().rename("reference_mean")
    return p_mean.subtract(r_mean).rename("bias_mm_day")

correlation_map(product_ic, reference_ic)

Pixel-wise Pearson correlation between two monthly ImageCollections over the full period, via ee.Reducer.pearsonsCorrelation on a time-matched image pair stack.

Source code in savana/rainfall/spatial.py

def correlation_map(product_ic, reference_ic):
    """Pixel-wise Pearson correlation between two monthly ImageCollections
    over the full period, via ``ee.Reducer.pearsonsCorrelation`` on a
    time-matched image pair stack.
    """
    import ee

    def _pair(img):
        date = img.get("system:time_start")
        match = reference_ic.filter(ee.Filter.eq("system:time_start", date)).first()
        return (
            img.select("precip_mm_day")
            .rename("product")
            .addBands(ee.Image(match).select("precip_mm_day").rename("reference"))
            .set("system:time_start", date)
        )

    paired = product_ic.map(_pair)
    corr = paired.select(["product", "reference"]).reduce(
        ee.Reducer.pearsonsCorrelation()
    )
    return corr.select("correlation").rename("r")

preview_bias_map(product_ic, reference_ic, product_name='', reference_name='', region=None, m=None)

A quick INTER-PRODUCT bias map (product minus another gridded product), before running formal validation, "roughly where do these two products disagree spatially?" Direct port of the GEE app's Bias Map button.

IMPORTANT: this is product vs. product, never product vs. GPCC. GPCC exists in this package only as point gauge observations (see :mod:.stations), there's no gridded GPCC raster to difference a product against pixel-by-pixel. For the actual "does this product agree with real GPCC ground truth" spatial check, use :func:preview_station_bias_map instead, which plots true bias at each gauge location. This function is for a different, valid question, "how much do CHIRPS and GPM-IMERG disagree with each other spatially", not a validation check.

Parameters:

Name	Type	Description	Default
product_ic, reference_ic		monthly mm/day ImageCollections (both gridded products, neither is GPCC).	required
product_name, reference_name		labels for the map layers.	required
region		clip to this ee.Geometry if given.	None
m		an existing geemap.Map to add to, or a new one is created.	None

Returns:

Type	Description
	The geemap.Map with bias and percent-bias layers added.

Source code in savana/rainfall/spatial.py

def preview_bias_map(
    product_ic, reference_ic, product_name="", reference_name="", region=None, m=None
):
    """A quick INTER-PRODUCT bias map (product minus another gridded
    product), before running formal validation, "roughly where do
    these two products disagree spatially?" Direct port of the GEE
    app's Bias Map button.

    IMPORTANT: this is product vs. product, never product vs. GPCC.
    GPCC exists in this package only as point gauge observations (see
    :mod:`.stations`), there's no gridded GPCC raster to difference a
    product against pixel-by-pixel. For the actual "does this product
    agree with real GPCC ground truth" spatial check, use
    :func:`preview_station_bias_map` instead, which plots true bias at
    each gauge location. This function is for a different, valid
    question, "how much do CHIRPS and GPM-IMERG disagree with each
    other spatially", not a validation check.

    Args:
        product_ic, reference_ic: monthly mm/day ImageCollections (both
            gridded products, neither is GPCC).
        product_name, reference_name: labels for the map layers.
        region: clip to this ``ee.Geometry`` if given.
        m: an existing ``geemap.Map`` to add to, or a new one is created.

    Returns:
        The ``geemap.Map`` with bias and percent-bias layers added.
    """
    import geemap

    prod_mean = product_ic.select("precip_mm_day").mean()
    ref_mean = reference_ic.select("precip_mm_day").mean()
    bias = prod_mean.subtract(ref_mean).rename("bias_mm_day")
    pbias = (
        prod_mean.subtract(ref_mean)
        .divide(ref_mean.add(1e-6))
        .multiply(100)
        .clamp(-80, 80)
        .rename("pbias_pct")
    )
    if region is not None:
        bias, pbias = bias.clip(region), pbias.clip(region)

    if m is None:
        m = geemap.Map()
        if region is not None:
            m.centerObject(region, 6)

    suffix = f"{product_name} vs {reference_name}" if product_name else "Bias"
    m.add_layer(bias, config.DEFAULT_VIS_PARAMS["bias"], f"Bias (mm/d), {suffix}")
    m.add_layer(pbias, config.DEFAULT_VIS_PARAMS["pbias"], f"% Bias, {suffix}")
    return m

preview_mean_map(product_ic, product_name='', region=None, kind='daily', m=None, obs_df=None, stations_df=None)

A quick look at one product's long-term mean rainfall on an interactive map, before running any validation, just "does this product's spatial pattern look sane over my area?"

Direct port of the GEE app's Annual Total / Mean Daily Rate map buttons.

Parameters:

Name	Type	Description	Default
product_ic		a monthly mm/day ee.ImageCollection (from :func:savana.rainfall.ingestion.load_product).	required
product_name	str	label for the map layer.	''
region		clip to this ee.Geometry if given.	None
kind	str	"daily" (mm/day, default) or "annual" (mm/yr, mean daily rate x 365.25), selects which :data:config.DEFAULT_VIS_PARAMS entry is used.	'daily'
m		an existing geemap.Map to add to, or a new one is created.	None
obs_df, stations_df		if BOTH given, overlays real GPCC station values as colored markers on top of the raster, NOT a rasterized/interpolated GPCC surface (GPCC stays point data throughout this package), just each gauge's true mean observed value, plotted at its real location, colored on the same scale as the raster underneath it, so you can eyeball whether the raster's color at a station roughly matches that station's actual marker color. Click a marker (with geemap's Inspector tool active) to see both the exact GPCC value and the raster's pixel value at that same point side by side.	required

Returns:

Type	Description
	The geemap.Map with the mean layer added (and the GPCC
	overlay, if requested).

Source code in savana/rainfall/spatial.py

def preview_mean_map(
    product_ic,
    product_name: str = "",
    region=None,
    kind: str = "daily",
    m=None,
    obs_df=None,
    stations_df=None,
):
    """A quick look at one product's long-term mean rainfall on an
    interactive map, before running any validation, just "does this
    product's spatial pattern look sane over my area?"

    Direct port of the GEE app's Annual Total / Mean Daily Rate map
    buttons.

    Args:
        product_ic: a monthly mm/day ``ee.ImageCollection`` (from
            :func:`savana.rainfall.ingestion.load_product`).
        product_name: label for the map layer.
        region: clip to this ``ee.Geometry`` if given.
        kind: ``"daily"`` (mm/day, default) or ``"annual"`` (mm/yr,
            mean daily rate x 365.25), selects which
            :data:`config.DEFAULT_VIS_PARAMS` entry is used.
        m: an existing ``geemap.Map`` to add to, or a new one is created.
        obs_df, stations_df: if BOTH given, overlays real GPCC station
            values as colored markers on top of the raster, NOT a
            rasterized/interpolated GPCC surface (GPCC stays point data
            throughout this package), just each gauge's true mean
            observed value, plotted at its real location, colored on the
            same scale as the raster underneath it, so you can eyeball
            whether the raster's color at a station roughly matches that
            station's actual marker color. Click a marker (with geemap's
            Inspector tool active) to see both the exact GPCC value and
            the raster's pixel value at that same point side by side.

    Returns:
        The ``geemap.Map`` with the mean layer added (and the GPCC
        overlay, if requested).
    """
    import geemap

    if kind not in ("daily", "annual"):
        raise ValueError(f"kind must be 'daily' or 'annual', got {kind!r}")

    mean_img = product_ic.select("precip_mm_day").mean()
    if kind == "annual":
        mean_img = mean_img.multiply(365.25)
    if region is not None:
        mean_img = mean_img.clip(region)

    if m is None:
        m = geemap.Map()
        if region is not None:
            m.centerObject(region, 6)

    vis = config.DEFAULT_VIS_PARAMS[kind]
    label = f"{'Annual Total' if kind == 'annual' else 'Mean Daily'}, {product_name}"
    m.add_layer(mean_img, vis, label)

    if obs_df is not None and stations_df is not None:
        _add_gpcc_overlay(m, obs_df, stations_df, vis, kind)

    return m

preview_station_bias_map(merged_df, product, m=None, zoom=5)

Per-station mean bias against REAL GPCC observations, plotted as colored markers, the spatial check that's actually anchored to ground truth, unlike :func:preview_bias_map (which can only ever compare two gridded products against each other, since GPCC has no gridded form in this package).

Parameters:

Name	Type	Description	Default
merged_df		long-format obs/sim table with station coordinates already joined in, i.e. from :func:savana.rainfall.extraction.merge_with_observations called with stations_df= (which :meth:savana.rainfall.pipeline.RainfallAssessment.merge always does), so lon/lat columns are present.	required
product	str	which product's bias to show.	required
m		an existing geemap.Map to add to, or a new one is created.	None

Returns:

Type	Description
	A geemap.Map with one marker per station, blue if that
	product overestimates GPCC there, red if it underestimates.
	Click a marker to see the exact bias value.

Source code in savana/rainfall/spatial.py

def preview_station_bias_map(merged_df, product: str, m=None, zoom: int = 5):
    """Per-station mean bias against REAL GPCC observations, plotted as
    colored markers, the spatial check that's actually anchored to
    ground truth, unlike :func:`preview_bias_map` (which can only ever
    compare two gridded products against each other, since GPCC has no
    gridded form in this package).

    Args:
        merged_df: long-format obs/sim table with station coordinates
            already joined in, i.e. from
            :func:`savana.rainfall.extraction.merge_with_observations`
            called with ``stations_df=`` (which
            :meth:`savana.rainfall.pipeline.RainfallAssessment.merge`
            always does), so ``lon``/``lat`` columns are present.
        product: which product's bias to show.
        m: an existing ``geemap.Map`` to add to, or a new one is created.

    Returns:
        A ``geemap.Map`` with one marker per station, blue if that
        product overestimates GPCC there, red if it underestimates.
        Click a marker to see the exact bias value.
    """
    import ee
    import geemap

    sub = merged_df[merged_df["product"] == product]
    if sub.empty:
        raise ValueError(f"No rows for product={product!r} in merged_df.")
    missing = {"lon", "lat"} - set(sub.columns)
    if missing:
        raise ValueError(
            f"merged_df is missing {sorted(missing)}, station coordinates "
            f"weren't joined in. Call merge_with_observations(..., "
            f"stations_df=your_stations_df), or just use "
            f"RainfallAssessment.merge(), which does this automatically."
        )

    per_station = (
        sub.assign(_bias=sub["sim_mm_day"] - sub["obs_mm_day"])
        .groupby(["station_id", "lon", "lat"], as_index=False)["_bias"]
        .mean()
        .rename(columns={"_bias": "bias_mm_day"})
    )

    if m is None:
        m = geemap.Map()
        m.set_center(float(per_station.lon.mean()), float(per_station.lat.mean()), zoom)

    features = [
        ee.Feature(
            ee.Geometry.Point([float(r.lon), float(r.lat)]),
            {"station_id": r.station_id, "bias_mm_day": round(float(r.bias_mm_day), 3)},
        )
        for r in per_station.itertuples()
    ]
    fc = ee.FeatureCollection(features)
    over = fc.filter(ee.Filter.gte("bias_mm_day", 0))
    under = fc.filter(ee.Filter.lt("bias_mm_day", 0))
    m.add_layer(
        over.style(**{"color": "0D47A1", "pointSize": 8}),
        {},
        f"{product} overestimates GPCC",
    )
    m.add_layer(
        under.style(**{"color": "B71C1C", "pointSize": 8}),
        {},
        f"{product} underestimates GPCC",
    )
    return m

threshold_sensitivity_map(product_ic, reference_ic, thresholds=None)

Pixel-wise CSI at each of several thresholds, spatial counterpart to :func:savana.rainfall.thresholds.threshold_sensitivity.

Returns {threshold: ee.Image} of pixel-wise CSI.

Source code in savana/rainfall/spatial.py

def threshold_sensitivity_map(
    product_ic, reference_ic, thresholds: list[float] | None = None
):
    """Pixel-wise CSI at each of several thresholds, spatial counterpart
    to :func:`savana.rainfall.thresholds.threshold_sensitivity`.

    Returns ``{threshold: ee.Image}`` of pixel-wise CSI.
    """
    import ee

    thresholds = (
        thresholds if thresholds is not None else config.DEFAULT_THRESHOLD_SWEEP_MM_DAY
    )
    out = {}
    for tau in thresholds:

        def _pair(img, _tau=tau):
            date = img.get("system:time_start")
            match = ee.Image(
                reference_ic.filter(ee.Filter.eq("system:time_start", date)).first()
            )
            obs_wet = match.select("precip_mm_day").gte(_tau)
            sim_wet = img.select("precip_mm_day").gte(_tau)
            hit = obs_wet.And(sim_wet).rename("hit")
            miss = obs_wet.And(sim_wet.Not()).rename("miss")
            fa = obs_wet.Not().And(sim_wet).rename("fa")
            return hit.addBands([miss, fa]).set("system:time_start", date)

        counts = product_ic.map(_pair).select(["hit", "miss", "fa"]).sum()
        csi = (
            counts.select("hit")
            .divide(
                counts.select("hit")
                .add(counts.select("miss"))
                .add(counts.select("fa"))
                .max(1)
            )
            .rename("csi")
        )
        out[tau] = csi
    return out

trend_map(product_ic, unit='mm/day/year')

Pixel-wise linear trend over time via ee.Reducer.linearFit (time in years since first image).

Source code in savana/rainfall/spatial.py

def trend_map(product_ic, unit: str = "mm/day/year"):
    """Pixel-wise linear trend over time via
    ``ee.Reducer.linearFit`` (time in years since first image).
    """
    import ee

    first_date = ee.Date(product_ic.first().get("system:time_start"))

    def _add_time_band(img):
        t = ee.Date(img.get("system:time_start")).difference(first_date, "year")
        return img.addBands(ee.Image.constant(t).rename("t").float())

    stacked = product_ic.map(_add_time_band).select(["t", "precip_mm_day"])
    fit = stacked.reduce(ee.Reducer.linearFit())
    slope = fit.select("scale").rename(f"trend_{unit.replace('/', '_')}")
    return slope

zonal_rank_table(products_ic, zones_gdf, reference_ic=None, name_field=None)

Zone-mean bias/trend per product, as a flat table, the scripted counterpart to the GEE app's zonal-ranking map layer.

Requires zones_gdf (see :func:savana.rainfall.zones.default_zones_wa). If reference_ic is given, also includes zone-mean bias against it.

Source code in savana/rainfall/spatial.py

def zonal_rank_table(
    products_ic: dict, zones_gdf, reference_ic=None, name_field: str | None = None
):
    """Zone-mean bias/trend per product, as a flat table, the scripted
    counterpart to the GEE app's zonal-ranking map layer.

    Requires ``zones_gdf`` (see
    :func:`savana.rainfall.zones.default_zones_wa`). If ``reference_ic``
    is given, also includes zone-mean bias against it.
    """
    import ee
    import pandas as pd

    name_field = name_field or config.DEFAULT_ZONE_NAME_FIELD
    zones_fc = ee.FeatureCollection(
        [
            ee.Feature(ee.Geometry(g.__geo_interface__), {name_field: n})
            for g, n in zip(zones_gdf.geometry, zones_gdf[name_field])
        ]
    )

    rows = []
    for name, ic in products_ic.items():
        mean_img = ic.select("precip_mm_day").mean().rename("mean_mm_day")
        trend_img = trend_map(ic)
        bands = mean_img.addBands(trend_img)
        if reference_ic is not None:
            bands = bands.addBands(bias_map(ic, reference_ic))

        stats = bands.reduceRegions(
            collection=zones_fc,
            reducer=ee.Reducer.mean(),
            scale=config.DEFAULT_TARGET_RESOLUTION_M,
        ).getInfo()

        for f in stats["features"]:
            p = f["properties"]
            row = {"product": name, "zone": p.get(name_field)}
            row.update({k: v for k, v in p.items() if k not in (name_field,)})
            rows.append(row)

    return pd.DataFrame(rows)

savana.rainfall.decision

Application-weighted product scoring and the interactive decision support workbook.

:func:score_products implements two normalisation modes:

• "fixed" (default): each metric normalised against a fixed plausible range (KGE against -1..1, NSE against -5..1, |PBIAS| against 0..60, etc — see :data:savana.rainfall.config.DEFAULT_NORMALIZATION_BOUNDS). This is what the shipped decision workbook and reported figures actually use — verified by reproducing WA_Precipitation_Decision_Tool_v2.xlsx's SELECTOR-sheet scores exactly (Fire-risk x Saharian x CHIRPS = 0.7551).
• "zone_relative": per-zone min-max across whatever products are being compared, matching the manuscript's written formula (section 2.6). Scale-invariant — recommended if you add/remove products from the default six, since fixed bounds were tuned for the original product set's plausible range.

:func:build_workbook writes a live spreadsheet tool (data-validation dropdowns + INDEX/MATCH formulas that recompute instantly), not a static report — the same design as the current WA_Precipitation_Decision_Tool_v2.xlsx, generalised to any products/zones/apps rather than hardcoded to the WA six.

best_product(scores_df, app, zone=None)

The top-scoring product for a given application (and optionally zone). Returns (product, score) or (None, None) if no match.

Source code in savana/rainfall/decision.py

def best_product(scores_df, app: str, zone: str | None = None):
    """The top-scoring product for a given application (and optionally
    zone). Returns ``(product, score)`` or ``(None, None)`` if no match.
    """
    sub = scores_df[scores_df["app"] == app]
    if zone is not None and "zone" in sub.columns:
        sub = sub[sub["zone"] == zone]
    if sub.empty:
        return None, None
    top = sub.sort_values("score", ascending=False).iloc[0]
    return top["product"], float(top["score"])

build_workbook(out_path, validation_by_zone_df, validation_overall_df=None, ranking_df=None, threshold_df=None, scores_df=None, app_weights=None, zone_notes=None)

Write the interactive decision-support workbook.

Sheet layout matches the current (v2) design: flat DATA_* sheets holding the real numbers, APP_WEIGHTS as a visible reference table, SCORES (flat app/zone/product/score — restores compatibility with fig_application_rankings_v4.py, which reads this exact sheet name/shape), and two live sheets driven by data-validation dropdowns + INDEX/MATCH formulas: SELECTOR (pick an application + zone, see every product ranked) and SCORECARD (pick a zone + product, see its raw metrics).

Parameters:

Name	Type	Description	Default
out_path		destination .xlsx path.	required
validation_by_zone_df		from :func:savana.rainfall.validation.validate_by_zone.	required
validation_overall_df, ranking_df, threshold_df		optional companion tables (validate_overall, rank_products, threshold_sensitivity outputs) — written as-is if given.	required
scores_df		from :func:score_products. Computed automatically from validation_by_zone_df + app_weights if not given.	None
app_weights, zone_notes		default to :data:config.DEFAULT_APP_WEIGHTS / :data:config.DEFAULT_ZONE_NOTES.	required

Source code in savana/rainfall/decision.py

def build_workbook(
    out_path,
    validation_by_zone_df,
    validation_overall_df=None,
    ranking_df=None,
    threshold_df=None,
    scores_df=None,
    app_weights: dict | None = None,
    zone_notes: dict | None = None,
):
    """Write the interactive decision-support workbook.

    Sheet layout matches the current (v2) design: flat ``DATA_*`` sheets
    holding the real numbers, ``APP_WEIGHTS`` as a visible reference
    table, ``SCORES`` (flat app/zone/product/score — restores
    compatibility with ``fig_application_rankings_v4.py``, which reads
    this exact sheet name/shape), and two live sheets driven by
    data-validation dropdowns + ``INDEX``/``MATCH`` formulas:
    ``SELECTOR`` (pick an application + zone, see every product ranked)
    and ``SCORECARD`` (pick a zone + product, see its raw metrics).

    Args:
        out_path: destination .xlsx path.
        validation_by_zone_df: from
            :func:`savana.rainfall.validation.validate_by_zone`.
        validation_overall_df, ranking_df, threshold_df: optional
            companion tables (validate_overall, rank_products,
            threshold_sensitivity outputs) — written as-is if given.
        scores_df: from :func:`score_products`. Computed automatically
            from ``validation_by_zone_df`` + ``app_weights`` if not
            given.
        app_weights, zone_notes: default to
            :data:`config.DEFAULT_APP_WEIGHTS` /
            :data:`config.DEFAULT_ZONE_NOTES`.
    """
    import openpyxl
    from openpyxl.styles import Font, PatternFill
    from openpyxl.utils import get_column_letter
    from openpyxl.worksheet.datavalidation import DataValidation

    app_weights = app_weights if app_weights is not None else config.DEFAULT_APP_WEIGHTS
    zone_notes = zone_notes if zone_notes is not None else config.DEFAULT_ZONE_NOTES

    # Pooled-only runs (no zones assigned) have validation_by_zone_df =
    # None. The workbook still works fine in that case -- we synthesize
    # a single "pooled" zone from the overall table so every sheet
    # (DATA_by_zone, SCORES, and the SELECTOR/SCORECARD dropdowns) has
    # one consistent zone to key on, instead of crashing on None.
    if validation_by_zone_df is None:
        if validation_overall_df is None:
            raise ValueError(
                "build_workbook needs validation_by_zone_df or "
                "validation_overall_df (run validation first)."
            )
        validation_by_zone_df = validation_overall_df.copy()
        if "zone" not in validation_by_zone_df.columns:
            validation_by_zone_df.insert(0, "zone", "pooled")

    if scores_df is None:
        scores_df = score_products(validation_by_zone_df, weights=app_weights)
    elif "zone" not in scores_df.columns:
        scores_df = scores_df.copy()
        scores_df.insert(scores_df.columns.get_loc("product"), "zone", "pooled")

    wb = openpyxl.Workbook()
    wb.remove(wb.active)

    header_fill = PatternFill("solid", fgColor="1A6B1A")
    header_font = Font(color="FFFFFF", bold=True)

    def _write_df(ws, df):
        ws.append(list(df.columns))
        for cell in ws[1]:
            cell.fill = header_fill
            cell.font = header_font
        for row in df.itertuples(index=False):
            ws.append(list(row))
        for i, col in enumerate(df.columns, start=1):
            width = max(10, min(28, int(df[col].astype(str).str.len().max() or 10) + 2))
            ws.column_dimensions[get_column_letter(i)].width = width

    ws_zone = wb.create_sheet("DATA_by_zone")
    _write_df(ws_zone, validation_by_zone_df)

    if validation_overall_df is not None:
        _write_df(wb.create_sheet("DATA_overall"), validation_overall_df)
    if ranking_df is not None:
        _write_df(wb.create_sheet("DATA_ranking"), ranking_df)
    if threshold_df is not None:
        _write_df(wb.create_sheet("DATA_threshold"), threshold_df)

    ws_weights = wb.create_sheet("APP_WEIGHTS")
    metric_cols = sorted({m for w in app_weights.values() for m in w})
    ws_weights.append(["Application"] + metric_cols + ["Primary focus"])
    for cell in ws_weights[1]:
        cell.fill = header_fill
        cell.font = header_font
    focus = config.DEFAULT_APP_FOCUS
    for app, w in app_weights.items():
        ws_weights.append(
            [app] + [w.get(m, 0.0) for m in metric_cols] + [focus.get(app, "")]
        )

    ws_scores = wb.create_sheet("SCORES")
    _write_df(
        ws_scores,
        (
            scores_df[["app", "zone", "product", "score"]]
            if "zone" in scores_df.columns
            else scores_df[["app", "product", "score"]]
        ),
    )

    # ── SELECTOR: dropdown App + Zone -> ranked product table ──
    ws_sel = wb.create_sheet("SELECTOR")
    ws_sel["B2"] = "Select Application:"
    ws_sel["B2"].font = Font(bold=True)
    ws_sel["C2"] = list(app_weights.keys())[0]
    ws_sel["E2"] = "Select Zone:"
    ws_sel["E2"].font = Font(bold=True)
    zones_available = (
        sorted(scores_df["zone"].unique()) if "zone" in scores_df.columns else []
    )
    ws_sel["F2"] = zones_available[0] if zones_available else ""

    dv_app = DataValidation(type="list", formula1=f'"{",".join(app_weights.keys())}"')
    ws_sel.add_data_validation(dv_app)
    dv_app.add(ws_sel["C2"])
    if zones_available:
        dv_zone = DataValidation(type="list", formula1=f'"{",".join(zones_available)}"')
        ws_sel.add_data_validation(dv_zone)
        dv_zone.add(ws_sel["F2"])

    ws_sel["B4"] = "Zone note:"
    ws_sel["B4"].font = Font(bold=True)
    ws_sel["C4"] = (
        "=IFERROR(VLOOKUP(F2, {"
        + ",".join(f'"{z}","{n}"' for z, n in zone_notes.items())
        + '}, 2, FALSE), "")'
    )

    header_row = 6
    ws_sel.cell(header_row, 2, "Rank").font = header_font
    ws_sel.cell(header_row, 3, "Product").font = header_font
    ws_sel.cell(header_row, 4, "Score").font = header_font
    for c in (2, 3, 4):
        ws_sel.cell(header_row, c).fill = header_fill

    # Helper lookup table (hidden columns J:M): app, zone, product, score,
    # plus a concatenated key — same shape as the SCORES sheet, written
    # again here so SELECTOR's formulas don't depend on sheet order.
    key_col, app_col, zone_col, prod_col, score_col = "J", "K", "L", "M", "N"
    ws_sel[f"{app_col}1"], ws_sel[f"{zone_col}1"] = "app", "zone"
    ws_sel[f"{prod_col}1"], ws_sel[f"{score_col}1"] = "product", "score"
    ws_sel[f"{key_col}1"] = "key"
    for i, r in enumerate(scores_df.itertuples(index=False), start=2):
        zone_val = getattr(r, "zone", "")
        ws_sel[f"{app_col}{i}"] = r.app
        ws_sel[f"{zone_col}{i}"] = zone_val
        ws_sel[f"{prod_col}{i}"] = r.product
        ws_sel[f"{score_col}{i}"] = r.score
        ws_sel[f"{key_col}{i}"] = f'={app_col}{i}&"|"&{zone_col}{i}&"|"&{prod_col}{i}'
    last_row = len(scores_df) + 1

    products = sorted(scores_df["product"].unique())
    for i, prod in enumerate(products, start=1):
        r = header_row + i
        ws_sel.cell(r, 2, i)
        ws_sel.cell(r, 3, prod)
        formula = (
            f"=IFERROR(INDEX(${score_col}$2:${score_col}${last_row},"
            f'MATCH($C$2&"|"&$F$2&"|"&"{prod}",'
            f'${key_col}$2:${key_col}${last_row},0)),"")'
        )
        ws_sel.cell(r, 4, formula)

    for col, width in zip("BCDEF", (8, 22, 10, 16, 20)):
        ws_sel.column_dimensions[col].width = width
    for col in (key_col, app_col, zone_col, prod_col, score_col):
        ws_sel.column_dimensions[col].width = 14

    # ── SCORECARD: dropdown Zone + Product -> raw metrics ──
    ws_card = wb.create_sheet("SCORECARD")
    ws_card["B2"] = "Zone:"
    ws_card["B2"].font = Font(bold=True)
    ws_card["C2"] = zones_available[0] if zones_available else ""
    ws_card["D2"] = "Product:"
    ws_card["D2"].font = Font(bold=True)
    ws_card["E2"] = products[0] if products else ""

    if zones_available:
        dv_zone2 = DataValidation(
            type="list", formula1=f'"{",".join(zones_available)}"'
        )
        ws_card.add_data_validation(dv_zone2)
        dv_zone2.add(ws_card["C2"])
    dv_prod = DataValidation(type="list", formula1=f'"{",".join(products)}"')
    ws_card.add_data_validation(dv_prod)
    dv_prod.add(ws_card["E2"])

    metric_report_cols = [
        c for c in validation_by_zone_df.columns if c in config.DEFAULT_METRICS_FLAT
    ]
    n_data_rows = len(validation_by_zone_df) + 1
    for i, m in enumerate(metric_report_cols, start=4):
        ws_card.cell(i, 2, m)
        col_letter = get_column_letter(validation_by_zone_df.columns.get_loc(m) + 1)
        formula = (
            f"=IFERROR(INDEX(DATA_by_zone!${col_letter}$2:${col_letter}${n_data_rows},"
            f"MATCH(1,(DATA_by_zone!$A$2:$A${n_data_rows}=$C$2)*"
            f"(DATA_by_zone!$B$2:$B${n_data_rows}=$E$2),0)))"
        )
        ws_card.cell(i, 3, formula)
    ws_card.column_dimensions["B"].width = 14
    ws_card.column_dimensions["C"].width = 14

    wb.save(out_path)
    print(f"  Decision workbook written: {out_path}")
    return out_path

score_products(validation_df, weights=None, normalization='fixed', bounds=None, group_cols=None)

Application-weighted composite score for every (zone, product) combination, for every application in weights.

Parameters:

Name	Type	Description	Default
validation_df		per-zone (or per-station/pooled) metrics table, e.g. from :func:savana.rainfall.validation.validate_by_zone. Must have a product column and the metric columns referenced by weights (kge, r, nse, pod, far, csi, pbias by default).	required
weights	dict | None	{application_name: {metric: weight, ...}}, weights summing to 1.0 per application. Defaults to :data:config.DEFAULT_APP_WEIGHTS (the 7 conservation/ water-management applications) — pass your own for different applications or priorities.	None
normalization	str	"fixed" (default, matches shipped results) or "zone_relative" (matches the manuscript's written formula — see module docstring).	'fixed'
bounds	dict | None	only used when normalization="fixed". Defaults to :data:config.DEFAULT_NORMALIZATION_BOUNDS.	None
group_cols	list[str] | None	columns identifying each row's context (default: ["zone"] if present, else none — i.e. pooled).	None

Returns:

Type	Description
	Long-format DataFrame: group_cols + ["app", "product", "score"].

Source code in savana/rainfall/decision.py

def score_products(
    validation_df,
    weights: dict | None = None,
    normalization: str = "fixed",
    bounds: dict | None = None,
    group_cols: list[str] | None = None,
):
    """Application-weighted composite score for every (zone, product)
    combination, for every application in ``weights``.

    Args:
        validation_df: per-zone (or per-station/pooled) metrics table,
            e.g. from :func:`savana.rainfall.validation.validate_by_zone`.
            Must have a ``product`` column and the metric columns
            referenced by ``weights`` (kge, r, nse, pod, far, csi, pbias
            by default).
        weights: ``{application_name: {metric: weight, ...}}``, weights
            summing to 1.0 per application. Defaults to
            :data:`config.DEFAULT_APP_WEIGHTS` (the 7 conservation/
            water-management applications) — pass your own for
            different applications or priorities.
        normalization: ``"fixed"`` (default, matches shipped results) or
            ``"zone_relative"`` (matches the manuscript's written
            formula — see module docstring).
        bounds: only used when ``normalization="fixed"``. Defaults to
            :data:`config.DEFAULT_NORMALIZATION_BOUNDS`.
        group_cols: columns identifying each row's context (default:
            ``["zone"]`` if present, else none — i.e. pooled).

    Returns:
        Long-format DataFrame: ``group_cols + ["app", "product", "score"]``.
    """
    import pandas as pd

    weights = weights if weights is not None else config.DEFAULT_APP_WEIGHTS
    bounds = bounds if bounds is not None else config.DEFAULT_NORMALIZATION_BOUNDS
    if group_cols is None:
        group_cols = ["zone"] if "zone" in validation_df.columns else []

    metric_cols = sorted({m for w in weights.values() for m in w})
    missing = [m for m in metric_cols if m not in validation_df.columns]
    if missing:
        raise ValueError(
            f"validation_df is missing metric column(s) required by weights: "
            f"{missing}"
        )

    df = validation_df.copy()

    if normalization == "fixed":
        for m in metric_cols:
            if m not in bounds:
                raise ValueError(f"No normalization bounds given for metric {m!r}")
            df[f"_norm_{m}"] = df[m].apply(
                lambda v, mm=m: _normalise_fixed(v, bounds[mm])
            )
    elif normalization == "zone_relative":
        for m in metric_cols:
            invert = bounds.get(m, (None, None, False))[2]
            grp = df.groupby(group_cols)[m] if group_cols else df[m]

            def _rel(s, _invert=invert):
                lo, hi = s.min(), s.max()
                if hi == lo:
                    return s.apply(lambda _: float("nan"))
                frac = (s - lo) / (hi - lo)
                return 1 - frac if _invert else frac

            if group_cols:
                df[f"_norm_{m}"] = grp.transform(_rel)
            else:
                df[f"_norm_{m}"] = _rel(df[m])
    else:
        raise ValueError(
            f"Unknown normalization {normalization!r}: expected "
            f'"fixed" or "zone_relative".'
        )

    rows = []
    for app, app_weights in weights.items():
        score = sum(df[f"_norm_{m}"] * w for m, w in app_weights.items())
        for idx, s in score.items():
            row = {c: df.at[idx, c] for c in group_cols}
            row["app"] = app
            row["product"] = df.at[idx, "product"]
            row["score"] = round(float(s), 4) if s == s else float("nan")  # NaN-safe
            rows.append(row)

    scores_df = pd.DataFrame(rows)
    sort_cols = group_cols + ["app", "score"]
    return scores_df.sort_values(
        sort_cols, ascending=[True] * len(group_cols) + [True, False]
    ).reset_index(drop=True)

savana.rainfall.insights

Grounded facts, summary, and Q&A for a rainfall assessment.

Mirrors :mod:savana.insights's pattern exactly: :func:compute_facts computes everything once, with every section independently wrapped so one failure (e.g. no threshold data was run) doesn't take down the rest; :func:summarize turns facts into readable prose; :func:answer does grounded keyword-based retrieval against the facts dict. Every number in the output traces back to something actually computed, never fabricated, the same rule as savana.insights.

answer(facts, question)

Grounded keyword-based answer to a question about the assessment.

Matches application names and zone names appearing (case-insensitive, substring) in question against facts, and reports only what was actually computed. Falls back to :func:summarize if nothing specific matches.

Source code in savana/rainfall/insights.py

def answer(facts: dict, question: str) -> str:
    """Grounded keyword-based answer to a question about the assessment.

    Matches application names and zone names appearing (case-insensitive,
    substring) in ``question`` against ``facts``, and reports only what
    was actually computed. Falls back to :func:`summarize` if nothing
    specific matches.
    """
    q = question.lower()

    matched_apps = [a for a in facts.get("apps", []) if a.lower() in q]
    matched_zones = [z for z in facts.get("zones", []) if z.lower() in q]

    if matched_apps and matched_zones:
        lines = []
        for app in matched_apps:
            for zone in matched_zones:
                pair = facts.get("best_by_app_zone", {}).get(app, {}).get(zone)
                if pair:
                    prod, score = pair
                    note = facts.get("zone_notes", {}).get(zone, "")
                    lines.append(
                        f"For {app} in the {zone} zone, {prod} scores highest "
                        f"({score:.3f})." + (f" Note: {note}" if note else "")
                    )
        if lines:
            return "\n".join(lines)

    if matched_apps:
        lines = []
        for app in matched_apps:
            pair = facts.get("best_by_app_pooled", {}).get(app)
            if pair:
                prod, score = pair
                lines.append(
                    f"For {app} (pooled across all zones), {prod} scores highest "
                    f"({score:.3f}). Zone-specific selection is usually preferred "
                    f"— ask about a specific zone for a more precise answer."
                )
        if lines:
            return "\n".join(lines)

    if matched_zones:
        lines = []
        for zone in matched_zones:
            pair = facts.get("top_kge_by_zone", {}).get(zone)
            note = facts.get("zone_notes", {}).get(zone, "")
            if pair:
                prod, kge = pair
                lines.append(
                    f"In the {zone} zone, {prod} has the best overall KGE "
                    f"({kge:.3f})." + (f" Note: {note}" if note else "")
                )
        if lines:
            return "\n".join(lines)

    if "bias" in q or "overestimat" in q or "underestimat" in q:
        worst = facts.get("worst_pbias_by_zone", {})
        if worst:
            lines = [f"  - {z}: {p} ({b:+.1f}%)" for z, (p, b) in worst.items()]
            return "Largest zone-level biases (PBIAS):\n" + "\n".join(lines)

    return (
        "I couldn't match that to a specific application or zone in this "
        "assessment. Here's the full summary instead:\n\n" + summarize(facts)
    )

compute_facts(scores_df, validation_df, ranking_df=None, threshold_df=None, zone_notes=None)

Compute every grounded fact available from a rainfall assessment.

Parameters:

Name	Type	Description	Default
scores_df		from :func:savana.rainfall.decision.score_products.	required
validation_df		from :func:savana.rainfall.validation.validate_by_zone (or any grouped validation output with a product column).	required
ranking_df		optional, from :func:savana.rainfall.validation.rank_products.	None
threshold_df		optional, from :func:savana.rainfall.thresholds.threshold_sensitivity.	None
zone_notes	dict | None	defaults to :data:config.DEFAULT_ZONE_NOTES.	None

Returns:

Type	Description
dict	dict with keys: ``n_products, n_zones, apps, zones, products,
dict	best_by_app_zone, best_by_app_pooled, top_kge_by_zone,
dict	threshold_stability, zone_notes, warnings``.

Source code in savana/rainfall/insights.py

def compute_facts(
    scores_df,
    validation_df,
    ranking_df=None,
    threshold_df=None,
    zone_notes: dict | None = None,
) -> dict:
    """Compute every grounded fact available from a rainfall assessment.

    Args:
        scores_df: from :func:`savana.rainfall.decision.score_products`.
        validation_df: from
            :func:`savana.rainfall.validation.validate_by_zone` (or any
            grouped validation output with a ``product`` column).
        ranking_df: optional, from
            :func:`savana.rainfall.validation.rank_products`.
        threshold_df: optional, from
            :func:`savana.rainfall.thresholds.threshold_sensitivity`.
        zone_notes: defaults to :data:`config.DEFAULT_ZONE_NOTES`.

    Returns:
        dict with keys: ``n_products, n_zones, apps, zones, products,
        best_by_app_zone, best_by_app_pooled, top_kge_by_zone,
        threshold_stability, zone_notes, warnings``.
    """
    zone_notes = zone_notes if zone_notes is not None else config.DEFAULT_ZONE_NOTES
    facts: dict = {"warnings": []}

    try:
        facts["products"] = sorted(scores_df["product"].unique())
        facts["apps"] = sorted(scores_df["app"].unique())
        facts["zones"] = (
            sorted(scores_df["zone"].unique()) if "zone" in scores_df.columns else []
        )
        facts["n_products"] = len(facts["products"])
        facts["n_zones"] = len(facts["zones"])
    except Exception as exc:  # noqa: BLE001
        facts["warnings"].append(f"Could not read scores_df structure: {exc}")
        facts["products"], facts["apps"], facts["zones"] = [], [], []
        facts["n_products"], facts["n_zones"] = 0, 0

    try:
        best_by_app_zone = {}
        for app in facts["apps"]:
            best_by_app_zone[app] = {}
            for zone in facts["zones"]:
                prod, score = decision.best_product(scores_df, app, zone)
                if prod is not None:
                    best_by_app_zone[app][zone] = (prod, score)
        facts["best_by_app_zone"] = best_by_app_zone
    except Exception as exc:  # noqa: BLE001
        facts["warnings"].append(f"Could not compute best_by_app_zone: {exc}")
        facts["best_by_app_zone"] = {}

    try:
        best_by_app_pooled = {}
        for app in facts["apps"]:
            prod, score = decision.best_product(scores_df, app, zone=None)
            if prod is not None:
                best_by_app_pooled[app] = (prod, score)
        facts["best_by_app_pooled"] = best_by_app_pooled
    except Exception as exc:  # noqa: BLE001
        facts["warnings"].append(f"Could not compute best_by_app_pooled: {exc}")
        facts["best_by_app_pooled"] = {}

    try:
        top_kge_by_zone = {}
        if "kge" in validation_df.columns and "zone" in validation_df.columns:
            for zone, sub in validation_df.groupby("zone"):
                top = sub.sort_values("kge", ascending=False).iloc[0]
                top_kge_by_zone[zone] = (top["product"], round(float(top["kge"]), 3))
        facts["top_kge_by_zone"] = top_kge_by_zone
    except Exception as exc:  # noqa: BLE001
        facts["warnings"].append(f"Could not compute top_kge_by_zone: {exc}")
        facts["top_kge_by_zone"] = {}

    try:
        worst_bias = {}
        if "pbias" in validation_df.columns and "zone" in validation_df.columns:
            for zone, sub in validation_df.groupby("zone"):
                worst = sub.iloc[sub["pbias"].abs().idxmax() - sub.index[0]]
                worst_bias[zone] = (worst["product"], round(float(worst["pbias"]), 2))
        facts["worst_pbias_by_zone"] = worst_bias
    except Exception as exc:  # noqa: BLE001
        facts["warnings"].append(f"Could not compute worst_pbias_by_zone: {exc}")
        facts["worst_pbias_by_zone"] = {}

    facts["threshold_stability"] = {}
    if threshold_df is not None:
        try:
            from . import thresholds as _thresholds

            stab = _thresholds.stability_summary(threshold_df, metric="csi")
            facts["threshold_stability"] = {
                (row.get("zone", "pooled"), row["product"]): round(float(row["cv"]), 3)
                for row in stab.to_dict("records")
                if row["cv"] == row["cv"]  # drop NaN
            }
        except Exception as exc:  # noqa: BLE001
            facts["warnings"].append(f"Could not compute threshold_stability: {exc}")

    facts["ranking"] = ranking_df
    facts["zone_notes"] = {z: n for z, n in zone_notes.items() if z in facts["zones"]}

    return facts

summarize(facts)

Turn compute_facts() output into a readable text summary.

Source code in savana/rainfall/insights.py

def summarize(facts: dict) -> str:
    """Turn ``compute_facts()`` output into a readable text summary."""
    lines = []
    lines.append(
        f"Assessed {facts.get('n_products', 0)} precipitation product(s) "
        f"across {facts.get('n_zones', 0)} zone(s): "
        f"{', '.join(facts.get('zones', [])) or 'pooled only'}."
    )

    best_pooled = facts.get("best_by_app_pooled", {})
    if best_pooled:
        lines.append("\nBest product per application (pooled across zones):")
        for app, (prod, score) in best_pooled.items():
            lines.append(f"  - {app}: {prod} (score {score:.3f})")

    top_kge = facts.get("top_kge_by_zone", {})
    if top_kge:
        lines.append("\nBest KGE (primary ranking metric) per zone:")
        for zone, (prod, kge) in top_kge.items():
            lines.append(f"  - {zone}: {prod} (KGE {kge:.3f})")

    worst_bias = facts.get("worst_pbias_by_zone", {})
    if worst_bias:
        lines.append("\nLargest |PBIAS| per zone (worth a closer look before use):")
        for zone, (prod, pbias) in worst_bias.items():
            lines.append(f"  - {zone}: {prod} ({pbias:+.1f}%)")

    zone_notes = facts.get("zone_notes", {})
    if zone_notes:
        lines.append("\nZone notes:")
        for zone, note in zone_notes.items():
            lines.append(f"  - {zone}: {note}")

    if facts.get("warnings"):
        lines.append("\nWarnings (some facts could not be computed):")
        for w in facts["warnings"]:
            lines.append(f"  - {w}")

    return "\n".join(lines)

savana.rainfall.viz

Static matplotlib figures for a rainfall assessment.

Reads directly from the DataFrames produced by :mod:.validation and :mod:.decision, no Excel round-trip required, though :func:recommendation_heatmap also happily reads a SCORES sheet exported by :func:savana.rainfall.decision.build_workbook if that's more convenient (same shape either way: app, zone, product, score).

application_ranking_bars(scores_df, app)

Bar chart of every product's score for one application, one bar group per zone.

Source code in savana/rainfall/viz.py

def application_ranking_bars(scores_df, app: str):
    """Bar chart of every product's score for one application, one bar
    group per zone.
    """
    import numpy as np

    sub = scores_df[scores_df["app"] == app]
    zones = sorted(sub["zone"].unique()) if "zone" in sub.columns else ["pooled"]
    products = sorted(sub["product"].unique())

    fig, ax = _get_fig_ax(figsize=(1.5 * len(zones) + 2, 5))
    width = 0.8 / len(products)
    x = np.arange(len(zones))
    for i, prod in enumerate(products):
        vals = [
            sub[(sub.get("zone", "pooled") == z) & (sub["product"] == prod)][
                "score"
            ].mean()
            for z in zones
        ]
        ax.bar(x + i * width, vals, width, label=prod)
    ax.set_xticks(x + width * (len(products) - 1) / 2)
    ax.set_xticklabels(zones, rotation=30, ha="right")
    ax.set_ylabel("Application-weighted score")
    ax.set_title(f"Product ranking, {app}")
    ax.legend(fontsize=8, ncol=2)
    fig.tight_layout()
    return fig

metric_heatmap(validation_df, metric='kge', group_col='zone')

Zone x product heatmap of one metric.

Source code in savana/rainfall/viz.py

def metric_heatmap(validation_df, metric: str = "kge", group_col: str = "zone"):
    """Zone x product heatmap of one metric."""
    pivot = validation_df.pivot_table(index=group_col, columns="product", values=metric)
    fig, ax = _get_fig_ax(
        figsize=(1.2 * len(pivot.columns) + 2, 0.6 * len(pivot.index) + 2)
    )
    im = ax.imshow(pivot.values, cmap="RdYlGn", aspect="auto")
    ax.set_xticks(range(len(pivot.columns)))
    ax.set_xticklabels(pivot.columns, rotation=45, ha="right")
    ax.set_yticks(range(len(pivot.index)))
    ax.set_yticklabels(pivot.index)
    for i in range(len(pivot.index)):
        for j in range(len(pivot.columns)):
            v = pivot.values[i, j]
            if v == v:  # not NaN
                ax.text(j, i, f"{v:.2f}", ha="center", va="center", fontsize=8)
    ax.set_title(f"{metric.upper()} by {group_col} x product")
    fig.colorbar(im, ax=ax, shrink=0.8, label=metric.upper())
    fig.tight_layout()
    return fig

preview_comparison(merged_df, station_id=None, product=None)

A quick "does this look right?" comparison of observed vs simulated values, before computing formal validation metrics. Scatter with a 1:1 reference line, one color per product (or filtered to one product/station if given). Mirrors the GEE app's per-station validation scatter chart.

Source code in savana/rainfall/viz.py

def preview_comparison(
    merged_df, station_id: str | None = None, product: str | None = None
):
    """A quick "does this look right?" comparison of observed vs
    simulated values, before computing formal validation metrics.
    Scatter with a 1:1 reference line, one color per product (or
    filtered to one product/station if given). Mirrors the GEE app's
    per-station validation scatter chart.
    """
    df = merged_df
    if station_id is not None:
        df = df[df["station_id"] == station_id]
    if product is not None:
        df = df[df["product"] == product]
    if df.empty:
        raise ValueError("No rows match the given station_id/product filter.")

    fig, ax = _get_fig_ax(figsize=(6, 6))
    for prod, sub in df.groupby("product"):
        ax.scatter(sub["obs_mm_day"], sub["sim_mm_day"], s=10, alpha=0.5, label=prod)

    lim = max(df["obs_mm_day"].max(), df["sim_mm_day"].max()) * 1.05
    ax.plot([0, lim], [0, lim], "k--", linewidth=1, label="1:1")
    ax.set_xlim(0, lim)
    ax.set_ylim(0, lim)
    ax.set_xlabel("Observed (mm/day)")
    ax.set_ylabel("Simulated (mm/day)")
    title = "Obs vs Sim"
    if station_id:
        title += f", {station_id}"
    if product:
        title += f", {product}"
    ax.set_title(title)
    ax.legend(fontsize=8)
    fig.tight_layout()
    return fig

preview_observations(obs_df, station_id=None)

A quick time-series look at raw GPCC observations, before extracting or validating any product, "does the reference data itself look sane?" One line per station, or a single station if station_id is given.

Source code in savana/rainfall/viz.py

def preview_observations(obs_df, station_id: str | None = None):
    """A quick time-series look at raw GPCC observations, before
    extracting or validating any product, "does the reference data
    itself look sane?" One line per station, or a single station if
    ``station_id`` is given.
    """
    import pandas as pd

    df = obs_df if station_id is None else obs_df[obs_df["station_id"] == station_id]
    if df.empty:
        raise ValueError(f"No observations found for station_id={station_id!r}")

    fig, ax = _get_fig_ax(figsize=(10, 4))
    for sid, sub in df.groupby("station_id"):
        sub = sub.sort_values(["year", "month"])
        t = pd.to_datetime(
            sub["year"].astype(str) + "-" + sub["month"].astype(str) + "-01"
        )
        ax.plot(t, sub["obs_mm_day"], label=sid, linewidth=1)
    ax.set_ylabel("Observed (mm/day)")
    ax.set_title(
        "GPCC observations" + (f", {station_id}" if station_id else ", all stations")
    )
    if df["station_id"].nunique() > 1:
        ax.legend(fontsize=7, ncol=4)
    fig.tight_layout()
    return fig

recommendation_heatmap(scores_df)

App x zone grid, each cell showing the single best product + its score, the "decision matrix" view, folding in fig_application_rankings_v4.py's recommendation-heatmap figure.

Source code in savana/rainfall/viz.py

def recommendation_heatmap(scores_df):
    """App x zone grid, each cell showing the single best product +
    its score, the "decision matrix" view, folding in
    fig_application_rankings_v4.py's recommendation-heatmap figure.
    """
    from . import decision

    apps = sorted(scores_df["app"].unique())
    zones = (
        sorted(scores_df["zone"].unique())
        if "zone" in scores_df.columns
        else ["pooled"]
    )

    best_scores = [[float("nan")] * len(zones) for _ in apps]
    labels = [[""] * len(zones) for _ in apps]
    for i, app in enumerate(apps):
        for j, zone in enumerate(zones):
            prod, score = decision.best_product(
                scores_df, app, zone if zone != "pooled" else None
            )
            if prod is not None:
                best_scores[i][j] = score
                labels[i][j] = f"{prod}\n{score:.2f}"

    fig, ax = _get_fig_ax(figsize=(1.5 * len(zones) + 3, 0.6 * len(apps) + 2))
    im = ax.imshow(best_scores, cmap="RdYlGn", vmin=0, vmax=1, aspect="auto")
    ax.set_xticks(range(len(zones)))
    ax.set_xticklabels(zones, rotation=30, ha="right")
    ax.set_yticks(range(len(apps)))
    ax.set_yticklabels(apps)
    for i in range(len(apps)):
        for j in range(len(zones)):
            if labels[i][j]:
                ax.text(j, i, labels[i][j], ha="center", va="center", fontsize=7)
    ax.set_title("Recommended product by application x zone")
    fig.colorbar(im, ax=ax, shrink=0.8, label="Application-weighted score")
    fig.tight_layout()
    return fig

taylor_diagram(validation_df, zone=None, ref_std=1.0)

Simplified Taylor diagram (correlation vs normalised std dev) for every product, optionally filtered to one zone.

Source code in savana/rainfall/viz.py

def taylor_diagram(validation_df, zone: str | None = None, ref_std: float = 1.0):
    """Simplified Taylor diagram (correlation vs normalised std dev)
    for every product, optionally filtered to one zone.
    """
    import matplotlib.pyplot as plt
    import numpy as np

    df = validation_df if zone is None else validation_df[validation_df["zone"] == zone]

    fig = plt.figure(figsize=(7, 7))
    ax = fig.add_subplot(111, polar=True)
    ax.set_thetalim(0, np.pi / 2)
    ax.set_xticks(np.arccos([1, 0.9, 0.7, 0.5, 0.3, 0])[::-1])
    ax.set_xticklabels(["1.0", "0.9", "0.7", "0.5", "0.3", "0"][::-1])

    for _, row in df.iterrows():
        r = row.get("r", float("nan"))
        if r != r:
            continue
        theta = np.arccos(max(-1, min(1, r)))
        # r2 as a stand-in radial "spread" proxy when true sim/obs std ratio
        # isn't in the table; callers with std ratios can pass their own.
        radius = row.get("std_ratio", 1.0)
        ax.plot(theta, radius, "o", label=row["product"], markersize=8)

    ax.set_title(f"Taylor diagram{f', {zone}' if zone else ''}")
    ax.legend(loc="upper left", bbox_to_anchor=(1.05, 1.0), fontsize=8)
    fig.tight_layout()
    return fig

zonal_boxplot(validation_df, metric='kge')

Distribution of one metric across zones, one box per product.

Source code in savana/rainfall/viz.py

def zonal_boxplot(validation_df, metric: str = "kge"):
    """Distribution of one metric across zones, one box per product."""
    products = sorted(validation_df["product"].unique())
    fig, ax = _get_fig_ax(figsize=(1.3 * len(products) + 2, 5))
    data = [
        validation_df.loc[validation_df["product"] == p, metric].dropna()
        for p in products
    ]
    ax.boxplot(data, labels=products)
    ax.set_ylabel(metric.upper())
    ax.set_title(f"{metric.upper()} distribution across zones, by product")
    ax.tick_params(axis="x", rotation=45)
    fig.tight_layout()
    return fig