Classification API Reference

This page documents the land-system classification API, the top-level savana.* modules. For precipitation product assessment, see the separate Precipitation Assessment API Reference, which documents the savana.rainfall.* modules.

Same names, different modules

A few module names appear in both packages: config, thresholds, and viz. On this page they always mean the classification versions (savana.config, savana.thresholds, savana.viz). Their rainfall namesakes (savana.rainfall.config, etc.) are entirely separate and documented on the rainfall API page.

High-level entry point

savana.pipeline

High-level orchestration: the "few lines of code" entry point.

Mirrors the workflow in mainrun.js step for step, but wrapped as a single Python object so a Jupyter user can go from an AOI to a validated, multi-epoch classified savanna land-system map in a handful of calls instead of hand-assembling every module.

Example

import savana clf = savana.SavanaClassifier( ... aoi="projects/ee-desmond/assets/NewParkMerged", ... name_filter="Kogyae", ... park_name="Kogyae", ... epochs=[2017, 2019, 2021, 2024], ... ) clf.run() clf.maps[2024] # ee.Image, classified 2024 land systems clf.accuracy_summary() # pandas.DataFrame, one row per model clf.show(2024) # interactive geemap.Map in the notebook

SavanaClassifier

End-to-end adaptive savanna land-system classifier for one AOI.

All parameters have sane defaults matching the original manuscript methodology; override any of them for a different landscape, class scheme, or sensor configuration.

Source code in savana/pipeline.py

class SavanaClassifier:
    """End-to-end adaptive savanna land-system classifier for one AOI.

    All parameters have sane defaults matching the original manuscript
    methodology; override any of them for a different landscape,
    class scheme, or sensor configuration.
    """

    def __init__(
        self,
        aoi,
        name_filter: str | None = None,
        park_name: str = "AOI",
        epochs: list[int] | None = None,
        reference_year: int | None = None,
        class_info: dict | None = None,
        class_property: str = config.CLASS_PROPERTY,
        n_clusters: int = config.DEFAULT_N_CLUSTERS,
        points_per_class: int = config.DEFAULT_POINTS_PER_CLASS,
        candidates_per_cluster: int = config.DEFAULT_CANDIDATES_PER_CLUSTER,
        confidence_margin: float = config.DEFAULT_CONFIDENCE_MARGIN,
        random_seed: int = config.DEFAULT_RANDOM_SEED,
        phenology_min_year: int = config.DEFAULT_PHENOLOGY_MIN_YEAR,
        rf_trees: int = 150,
        scale: int = config.DEFAULT_EXPORT_SCALE,
        crs: str = config.DEFAULT_CRS,
        ee_project: str | None = None,
    ):
        ee_init.initialize(project=ee_project)

        self.region = ee_init.load_aoi(aoi, name_filter=name_filter)
        self.park_name = park_name
        self.epochs = sorted(epochs or [2024])
        self.reference_year = reference_year or self.epochs[-1]
        self.class_info = class_info or config.DEFAULT_CLASS_INFO
        self.class_property = class_property
        self.n_clusters = n_clusters
        self.points_per_class = points_per_class
        self.candidates_per_cluster = candidates_per_cluster
        self.confidence_margin = confidence_margin
        self.random_seed = random_seed
        self.phenology_min_year = phenology_min_year
        self.rf_trees = rf_trees
        self.scale = scale
        self.crs = crs

        # Populated by .run()
        self.idx: dict | None = None
        self.T: dict | None = None
        self.masks: dict | None = None
        self.embedding = None
        self.gcps = None
        self.models: dict | None = None
        self.maps: dict = {}
        self.change: dict | None = None

    # -- pipeline stages, callable individually or via .run() ----------

    def build_features(self):
        """Build composites, indices, RUE, thresholds, and masks for the reference year."""
        s2_annual = composites.sentinel2_annual(self.reference_year, self.region)
        s2_dry = composites.seasonal_composite(
            f"{self.reference_year}-03-01", f"{self.reference_year}-05-15", self.region
        )
        s2_wet = composites.seasonal_composite(
            f"{self.reference_year}-05-01", f"{self.reference_year}-07-15", self.region
        )
        pcts = composites.percentile_composites(self.reference_year, self.region)
        self.embedding = composites.embedding_image(self.reference_year, self.region)

        self.idx = indices.compute(s2_annual, s2_dry, s2_wet, pcts["p10"], pcts["p90"])
        self.rue_annual = rue_mod.compute_annual(self.reference_year, self.region)
        self.T = thresholds_mod.compute(self.idx, self.region)
        self.masks = masks_mod.compute(self.idx, self.T)
        return self

    def sample_training_points(self):
        """Unsupervised clustering + rule-based labelling + class balancing."""
        cluster_result = sampling.cluster_embedding(
            self.embedding,
            self.region,
            n_clusters=self.n_clusters,
            seed=self.random_seed,
        )
        self.gcps = sampling.build_gcps(
            self.embedding,
            self.idx,
            cluster_result["clusters"],
            self.T,
            self.region,
            n_classes=len(self.class_info),
            points_per_class=self.points_per_class,
            scale=self.scale,
            class_property=self.class_property,
            n_clusters=self.n_clusters,
            candidates_per_cluster=self.candidates_per_cluster,
            confidence_margin=self.confidence_margin,
        )
        return self

    def train(self):
        """Train the 4-model ablation + master classifiers."""
        self.models = classifiers.train_all_models(
            self.gcps,
            self.embedding,
            self.idx,
            self.rue_annual["rue"],
            self.region,
            class_property=self.class_property,
            n_trees=self.rf_trees,
            seed=self.random_seed,
            class_order=sorted(self.class_info.keys()),
        )
        return self

    def classify(self):
        """Classify every requested epoch year."""
        self.maps = classifiers.classify_all_epochs(
            self.epochs,
            self.models,
            self.region,
            park_name=self.park_name,
            embedding_current_year=self.reference_year,
            embedding_current_image=self.embedding,
            phenology_min_year=self.phenology_min_year,
        )
        return self

    def analyse_change(self):
        """Run conservative + RUE-validated change detection across epochs."""
        if len(self.epochs) >= 2:
            self.change = change_mod.analyse(
                self.maps, self.epochs, self.region, park_name=self.park_name
            )
        return self

    def run(self):
        """Run the full pipeline: features -> sampling -> training -> classification -> change."""
        return (
            self.build_features()
            .sample_training_points()
            .train()
            .classify()
            .analyse_change()
        )

    # -- results & reporting --------------------------------------------

    def accuracy_summary(self):
        """One row per model (A/B/C/D) with overall accuracy, kappa, PA/UA."""
        matrices = {
            "a": self.models["matrix_a"],
            "b": self.models["matrix_b"],
            "c": self.models["matrix_c"],
            "d": self.models["matrix_d"],
        }
        return accuracy_mod.summary_dataframe(
            matrices, park_name=self.park_name, class_info=self.class_info
        )

    def confusion_matrices(self):
        """Full per-class confusion matrix table across all 4 models."""
        matrices = {
            "a": self.models["matrix_a"],
            "b": self.models["matrix_b"],
            "c": self.models["matrix_c"],
            "d": self.models["matrix_d"],
        }
        return accuracy_mod.confusion_matrix_dataframe(
            matrices, park_name=self.park_name, class_info=self.class_info
        )

    def class_areas(self):
        """Per-epoch class area statistics (km2) as a pandas DataFrame."""
        stats_scale = self.change["stats_scale"] if self.change else self.scale
        return exports_mod.class_areas_dataframe(
            self.maps, self.epochs, self.region, stats_scale
        )

    def export(self, drive_folder: str | None = None, asset_folder: str | None = None):
        """Export classified maps (+ change products, if computed) to Drive/Assets."""
        tasks = exports_mod.export_classified_maps(
            self.maps,
            self.epochs,
            self.region,
            park_name=self.park_name,
            drive_folder=drive_folder,
            asset_folder=asset_folder,
            scale=self.scale,
            crs=self.crs,
        )
        if self.change is not None and (drive_folder or asset_folder):
            tasks += exports_mod.export_change_products(
                self.change,
                self.region,
                park_name=self.park_name,
                drive_folder=drive_folder,
                asset_folder=asset_folder,
                scale=self.scale,
                crs=self.crs,
            )
        return tasks

    def show(self, year: int | None = None, m=None):
        """Display a classified epoch (default: reference year) on an interactive map."""
        year = year or self.reference_year
        return viz.show_classified_map(
            self.maps[year], region=self.region, class_info=self.class_info, m=m
        )

    def show_gcps(self, with_background: bool = True, m=None):
        """Display the ground control points on the map, colored by class.

        A sanity check on the sampling/labelling step — where the
        training points actually landed and whether their classes look
        spatially sensible — before trusting the classifier they train.
        Requires .sample_training_points() (or .run()) to have completed.

        Args:
            with_background: If True (default), shows the reference
                year's classified map underneath the points, dimmed, so
                you can visually compare point placement against the
                result. If False, points are shown alone.
        """
        if self.gcps is None:
            raise RuntimeError("Call .sample_training_points() (or .run()) first.")
        background = None
        if with_background and self.reference_year in self.maps:
            background = self.maps[self.reference_year]
        return viz.show_gcps(
            self.gcps,
            region=self.region,
            class_info=self.class_info,
            class_property=self.class_property,
            background=background,
            m=m,
        )

    def show_years(self, years: list[int] | None = None, m=None):
        """Display several classified epochs as toggleable layers on one map.

        Uses geemap's layer panel — check/uncheck each year's checkbox to
        flip between them. Defaults to all epochs the classifier ran.

        >>> clf.show_years()             # all epochs
        >>> clf.show_years([2019, 2024]) # just these two
        """
        return viz.show_multi_year_map(
            self.maps, years=years, region=self.region, class_info=self.class_info, m=m
        )

    def show_geolibre(self, year: int | None = None, m=None):
        """Display a classified epoch inside the GeoLibre Jupyter widget.

        Alternative to .show() — same idea, different map backend.
        Requires: pip install "savana[geolibre]" (Python >= 3.11).
        """
        from . import viz_geolibre

        year = year or self.reference_year
        return viz_geolibre.show_classified_map(
            self.maps[year], region=self.region, class_info=self.class_info, m=m
        )

    def show_years_geolibre(self, years: list[int] | None = None, m=None):
        """Display several classified epochs as toggleable layers in GeoLibre.

        Alternative to .show_years() — same idea, different map backend.
        Requires: pip install "savana[geolibre]" (Python >= 3.11).
        """
        from . import viz_geolibre

        return viz_geolibre.show_multi_year_map(
            self.maps, years=years, region=self.region, class_info=self.class_info, m=m
        )

    def compare(self, left=2024, right="SATELLITE", m=None):
        """Side-by-side swipe comparison between two years, or a year vs. a basemap.

        ``left``/``right`` each accept either an epoch year (int, must be
        in ``self.maps``) or a basemap name string (e.g. ``"SATELLITE"``,
        ``"HYBRID"``, ``"ROADMAP"``, ``"Esri.WorldImagery"``). Drag the
        handle in the middle of the resulting map to swipe.

        >>> clf.compare(2019, 2024)              # two classified years
        >>> clf.compare(2024, "SATELLITE")        # classified year vs. basemap
        """
        left_img = self.maps[left] if isinstance(left, int) else left
        right_img = self.maps[right] if isinstance(right, int) else right
        left_label = str(left) if isinstance(left, int) else left
        right_label = str(right) if isinstance(right, int) else right
        return viz.compare_split_map(
            left_img,
            right_img,
            left_label=left_label,
            right_label=right_label,
            region=self.region,
            class_info=self.class_info,
            m=m,
        )

    def show_change(self, m=None):
        """Display change-detection layers on an interactive map."""
        if self.change is None:
            raise RuntimeError(
                "Call .analyse_change() (or .run()) with >= 2 epochs first."
            )
        return viz.show_change_map(self.change, region=self.region, m=m)

    def facts(self) -> dict:
        """Compute the grounded facts dict — real numbers from your actual results.

        This is the single source of truth for .summarize() and .answer();
        call it directly if you want the raw structured data instead of text.
        """
        from . import insights

        return insights.compute_facts(self)

    def summarize(self) -> str:
        """Plain-English report generated entirely from real computed results.

        No AI, no invented numbers — every figure here traces back to
        .class_areas() / .accuracy_summary() / the change-detection stats.
        """
        from . import insights

        return insights.summarize(self.facts())

    def answer(self, question: str) -> str:
        """Answer a question about your results using only computed facts.

        Simple keyword matching, not an LLM — it can only ever report
        numbers the pipeline actually produced, so it can't hallucinate.
        Try asking about a class's area, the dominant class, accuracy,
        or change between years.

        >>> clf.answer("how much core woodland is there in 2024?")
        >>> clf.answer("what changed between the years?")
        >>> clf.answer("how accurate is the model?")
        """
        from . import insights

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

accuracy_summary()

One row per model (A/B/C/D) with overall accuracy, kappa, PA/UA.

Source code in savana/pipeline.py

def accuracy_summary(self):
    """One row per model (A/B/C/D) with overall accuracy, kappa, PA/UA."""
    matrices = {
        "a": self.models["matrix_a"],
        "b": self.models["matrix_b"],
        "c": self.models["matrix_c"],
        "d": self.models["matrix_d"],
    }
    return accuracy_mod.summary_dataframe(
        matrices, park_name=self.park_name, class_info=self.class_info
    )

analyse_change()

Run conservative + RUE-validated change detection across epochs.

Source code in savana/pipeline.py

def analyse_change(self):
    """Run conservative + RUE-validated change detection across epochs."""
    if len(self.epochs) >= 2:
        self.change = change_mod.analyse(
            self.maps, self.epochs, self.region, park_name=self.park_name
        )
    return self

answer(question)

Answer a question about your results using only computed facts.

Simple keyword matching, not an LLM — it can only ever report numbers the pipeline actually produced, so it can't hallucinate. Try asking about a class's area, the dominant class, accuracy, or change between years.

clf.answer("how much core woodland is there in 2024?") clf.answer("what changed between the years?") clf.answer("how accurate is the model?")

Source code in savana/pipeline.py

def answer(self, question: str) -> str:
    """Answer a question about your results using only computed facts.

    Simple keyword matching, not an LLM — it can only ever report
    numbers the pipeline actually produced, so it can't hallucinate.
    Try asking about a class's area, the dominant class, accuracy,
    or change between years.

    >>> clf.answer("how much core woodland is there in 2024?")
    >>> clf.answer("what changed between the years?")
    >>> clf.answer("how accurate is the model?")
    """
    from . import insights

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

build_features()

Build composites, indices, RUE, thresholds, and masks for the reference year.

Source code in savana/pipeline.py

def build_features(self):
    """Build composites, indices, RUE, thresholds, and masks for the reference year."""
    s2_annual = composites.sentinel2_annual(self.reference_year, self.region)
    s2_dry = composites.seasonal_composite(
        f"{self.reference_year}-03-01", f"{self.reference_year}-05-15", self.region
    )
    s2_wet = composites.seasonal_composite(
        f"{self.reference_year}-05-01", f"{self.reference_year}-07-15", self.region
    )
    pcts = composites.percentile_composites(self.reference_year, self.region)
    self.embedding = composites.embedding_image(self.reference_year, self.region)

    self.idx = indices.compute(s2_annual, s2_dry, s2_wet, pcts["p10"], pcts["p90"])
    self.rue_annual = rue_mod.compute_annual(self.reference_year, self.region)
    self.T = thresholds_mod.compute(self.idx, self.region)
    self.masks = masks_mod.compute(self.idx, self.T)
    return self

class_areas()

Per-epoch class area statistics (km2) as a pandas DataFrame.

Source code in savana/pipeline.py

def class_areas(self):
    """Per-epoch class area statistics (km2) as a pandas DataFrame."""
    stats_scale = self.change["stats_scale"] if self.change else self.scale
    return exports_mod.class_areas_dataframe(
        self.maps, self.epochs, self.region, stats_scale
    )

classify()

Classify every requested epoch year.

Source code in savana/pipeline.py

def classify(self):
    """Classify every requested epoch year."""
    self.maps = classifiers.classify_all_epochs(
        self.epochs,
        self.models,
        self.region,
        park_name=self.park_name,
        embedding_current_year=self.reference_year,
        embedding_current_image=self.embedding,
        phenology_min_year=self.phenology_min_year,
    )
    return self

compare(left=2024, right='SATELLITE', m=None)

Side-by-side swipe comparison between two years, or a year vs. a basemap.

left/right each accept either an epoch year (int, must be in self.maps) or a basemap name string (e.g. "SATELLITE", "HYBRID", "ROADMAP", "Esri.WorldImagery"). Drag the handle in the middle of the resulting map to swipe.

clf.compare(2019, 2024) # two classified years clf.compare(2024, "SATELLITE") # classified year vs. basemap

Source code in savana/pipeline.py

def compare(self, left=2024, right="SATELLITE", m=None):
    """Side-by-side swipe comparison between two years, or a year vs. a basemap.

    ``left``/``right`` each accept either an epoch year (int, must be
    in ``self.maps``) or a basemap name string (e.g. ``"SATELLITE"``,
    ``"HYBRID"``, ``"ROADMAP"``, ``"Esri.WorldImagery"``). Drag the
    handle in the middle of the resulting map to swipe.

    >>> clf.compare(2019, 2024)              # two classified years
    >>> clf.compare(2024, "SATELLITE")        # classified year vs. basemap
    """
    left_img = self.maps[left] if isinstance(left, int) else left
    right_img = self.maps[right] if isinstance(right, int) else right
    left_label = str(left) if isinstance(left, int) else left
    right_label = str(right) if isinstance(right, int) else right
    return viz.compare_split_map(
        left_img,
        right_img,
        left_label=left_label,
        right_label=right_label,
        region=self.region,
        class_info=self.class_info,
        m=m,
    )

confusion_matrices()

Full per-class confusion matrix table across all 4 models.

Source code in savana/pipeline.py

def confusion_matrices(self):
    """Full per-class confusion matrix table across all 4 models."""
    matrices = {
        "a": self.models["matrix_a"],
        "b": self.models["matrix_b"],
        "c": self.models["matrix_c"],
        "d": self.models["matrix_d"],
    }
    return accuracy_mod.confusion_matrix_dataframe(
        matrices, park_name=self.park_name, class_info=self.class_info
    )

export(drive_folder=None, asset_folder=None)

Export classified maps (+ change products, if computed) to Drive/Assets.

Source code in savana/pipeline.py

def export(self, drive_folder: str | None = None, asset_folder: str | None = None):
    """Export classified maps (+ change products, if computed) to Drive/Assets."""
    tasks = exports_mod.export_classified_maps(
        self.maps,
        self.epochs,
        self.region,
        park_name=self.park_name,
        drive_folder=drive_folder,
        asset_folder=asset_folder,
        scale=self.scale,
        crs=self.crs,
    )
    if self.change is not None and (drive_folder or asset_folder):
        tasks += exports_mod.export_change_products(
            self.change,
            self.region,
            park_name=self.park_name,
            drive_folder=drive_folder,
            asset_folder=asset_folder,
            scale=self.scale,
            crs=self.crs,
        )
    return tasks

facts()

Compute the grounded facts dict — real numbers from your actual results.

This is the single source of truth for .summarize() and .answer(); call it directly if you want the raw structured data instead of text.

Source code in savana/pipeline.py

def facts(self) -> dict:
    """Compute the grounded facts dict — real numbers from your actual results.

    This is the single source of truth for .summarize() and .answer();
    call it directly if you want the raw structured data instead of text.
    """
    from . import insights

    return insights.compute_facts(self)

run()

Run the full pipeline: features -> sampling -> training -> classification -> change.

Source code in savana/pipeline.py

def run(self):
    """Run the full pipeline: features -> sampling -> training -> classification -> change."""
    return (
        self.build_features()
        .sample_training_points()
        .train()
        .classify()
        .analyse_change()
    )

sample_training_points()

Unsupervised clustering + rule-based labelling + class balancing.

Source code in savana/pipeline.py

def sample_training_points(self):
    """Unsupervised clustering + rule-based labelling + class balancing."""
    cluster_result = sampling.cluster_embedding(
        self.embedding,
        self.region,
        n_clusters=self.n_clusters,
        seed=self.random_seed,
    )
    self.gcps = sampling.build_gcps(
        self.embedding,
        self.idx,
        cluster_result["clusters"],
        self.T,
        self.region,
        n_classes=len(self.class_info),
        points_per_class=self.points_per_class,
        scale=self.scale,
        class_property=self.class_property,
        n_clusters=self.n_clusters,
        candidates_per_cluster=self.candidates_per_cluster,
        confidence_margin=self.confidence_margin,
    )
    return self

show(year=None, m=None)

Display a classified epoch (default: reference year) on an interactive map.

Source code in savana/pipeline.py

def show(self, year: int | None = None, m=None):
    """Display a classified epoch (default: reference year) on an interactive map."""
    year = year or self.reference_year
    return viz.show_classified_map(
        self.maps[year], region=self.region, class_info=self.class_info, m=m
    )

show_change(m=None)

Display change-detection layers on an interactive map.

Source code in savana/pipeline.py

def show_change(self, m=None):
    """Display change-detection layers on an interactive map."""
    if self.change is None:
        raise RuntimeError(
            "Call .analyse_change() (or .run()) with >= 2 epochs first."
        )
    return viz.show_change_map(self.change, region=self.region, m=m)

show_gcps(with_background=True, m=None)

Display the ground control points on the map, colored by class.

A sanity check on the sampling/labelling step — where the training points actually landed and whether their classes look spatially sensible — before trusting the classifier they train. Requires .sample_training_points() (or .run()) to have completed.

Parameters:

Name	Type	Description	Default
with_background	bool	If True (default), shows the reference year's classified map underneath the points, dimmed, so you can visually compare point placement against the result. If False, points are shown alone.	True

Source code in savana/pipeline.py

def show_gcps(self, with_background: bool = True, m=None):
    """Display the ground control points on the map, colored by class.

    A sanity check on the sampling/labelling step — where the
    training points actually landed and whether their classes look
    spatially sensible — before trusting the classifier they train.
    Requires .sample_training_points() (or .run()) to have completed.

    Args:
        with_background: If True (default), shows the reference
            year's classified map underneath the points, dimmed, so
            you can visually compare point placement against the
            result. If False, points are shown alone.
    """
    if self.gcps is None:
        raise RuntimeError("Call .sample_training_points() (or .run()) first.")
    background = None
    if with_background and self.reference_year in self.maps:
        background = self.maps[self.reference_year]
    return viz.show_gcps(
        self.gcps,
        region=self.region,
        class_info=self.class_info,
        class_property=self.class_property,
        background=background,
        m=m,
    )

show_geolibre(year=None, m=None)

Display a classified epoch inside the GeoLibre Jupyter widget.

Alternative to .show() — same idea, different map backend. Requires: pip install "savana[geolibre]" (Python >= 3.11).

Source code in savana/pipeline.py

def show_geolibre(self, year: int | None = None, m=None):
    """Display a classified epoch inside the GeoLibre Jupyter widget.

    Alternative to .show() — same idea, different map backend.
    Requires: pip install "savana[geolibre]" (Python >= 3.11).
    """
    from . import viz_geolibre

    year = year or self.reference_year
    return viz_geolibre.show_classified_map(
        self.maps[year], region=self.region, class_info=self.class_info, m=m
    )

show_years(years=None, m=None)

Display several classified epochs as toggleable layers on one map.

Uses geemap's layer panel — check/uncheck each year's checkbox to flip between them. Defaults to all epochs the classifier ran.

clf.show_years() # all epochs clf.show_years([2019, 2024]) # just these two

Source code in savana/pipeline.py

def show_years(self, years: list[int] | None = None, m=None):
    """Display several classified epochs as toggleable layers on one map.

    Uses geemap's layer panel — check/uncheck each year's checkbox to
    flip between them. Defaults to all epochs the classifier ran.

    >>> clf.show_years()             # all epochs
    >>> clf.show_years([2019, 2024]) # just these two
    """
    return viz.show_multi_year_map(
        self.maps, years=years, region=self.region, class_info=self.class_info, m=m
    )

show_years_geolibre(years=None, m=None)

Display several classified epochs as toggleable layers in GeoLibre.

Alternative to .show_years() — same idea, different map backend. Requires: pip install "savana[geolibre]" (Python >= 3.11).

Source code in savana/pipeline.py

def show_years_geolibre(self, years: list[int] | None = None, m=None):
    """Display several classified epochs as toggleable layers in GeoLibre.

    Alternative to .show_years() — same idea, different map backend.
    Requires: pip install "savana[geolibre]" (Python >= 3.11).
    """
    from . import viz_geolibre

    return viz_geolibre.show_multi_year_map(
        self.maps, years=years, region=self.region, class_info=self.class_info, m=m
    )

summarize()

Plain-English report generated entirely from real computed results.

No AI, no invented numbers — every figure here traces back to .class_areas() / .accuracy_summary() / the change-detection stats.

Source code in savana/pipeline.py

def summarize(self) -> str:
    """Plain-English report generated entirely from real computed results.

    No AI, no invented numbers — every figure here traces back to
    .class_areas() / .accuracy_summary() / the change-detection stats.
    """
    from . import insights

    return insights.summarize(self.facts())

train()

Train the 4-model ablation + master classifiers.

Source code in savana/pipeline.py

def train(self):
    """Train the 4-model ablation + master classifiers."""
    self.models = classifiers.train_all_models(
        self.gcps,
        self.embedding,
        self.idx,
        self.rue_annual["rue"],
        self.region,
        class_property=self.class_property,
        n_trees=self.rf_trees,
        seed=self.random_seed,
        class_order=sorted(self.class_info.keys()),
    )
    return self

classify_landscape(aoi, epochs=None, park_name='AOI', **kwargs)

One-call convenience wrapper: build, run, and return a fitted classifier.

clf = savana.classify_landscape( ... "path/to/my_park.geojson", epochs=[2020, 2024], park_name="My Park" ... ) clf.show()

Source code in savana/pipeline.py

def classify_landscape(
    aoi, epochs: list[int] | None = None, park_name: str = "AOI", **kwargs
) -> SavanaClassifier:
    """One-call convenience wrapper: build, run, and return a fitted classifier.

    >>> clf = savana.classify_landscape(
    ...     "path/to/my_park.geojson", epochs=[2020, 2024], park_name="My Park"
    ... )
    >>> clf.show()
    """
    clf = SavanaClassifier(aoi, epochs=epochs, park_name=park_name, **kwargs)
    clf.run()
    return clf

Feature engineering

savana.composites

Sentinel-2, seasonal, percentile, and AlphaEarth embedding composites.

Ported from kogyae.js Phase 2 (Sentinel-2 reference imagery) and Phase 3 (AlphaEarth embeddings), and from the year-loop composite logic in Phase 8, the standalone composites.js module you use in mainrun.js was empty, so this reconstructs it from the working monolithic script with the exact same signatures mainrun.js expects (getSentinel2Annual, getSeasonalComposite, getPercentileComposites, getEmbeddingImage), plus a graceful fallback for seasons/years with sparse Cloud Score+ coverage.

embedding_image(year, region)

Annual AlphaEarth satellite embedding image (64 bands: A01..A64).

Equivalent to getEmbeddingImage / C.getEmbeddingImage.

Source code in savana/composites.py

def embedding_image(year: int, region):
    """Annual AlphaEarth satellite embedding image (64 bands: A01..A64).

    Equivalent to ``getEmbeddingImage`` / ``C.getEmbeddingImage``.
    """
    import ee

    start = ee.Date.fromYMD(year, 1, 1)
    end = start.advance(1, "year")
    return (
        ee.ImageCollection(config.EMBEDDING_COLLECTION)
        .filter(ee.Filter.date(start, end))
        .filter(ee.Filter.bounds(region))
        .mosaic()
        .clip(region)
    )

percentile_composites(year, region, percentiles=(10, 90), cloud_score_threshold=0.65)

Per-band percentile composites (default p10/p90) for a calendar year.

Used to build the seasonal-amplitude / stability indices. Equivalent to C.getPercentileComposites, returns a dict keyed "p{percentile}" (e.g. {"p10": image, "p90": image}).

Source code in savana/composites.py

def percentile_composites(
    year: int,
    region,
    percentiles: tuple[int, int] = (10, 90),
    cloud_score_threshold: float = 0.65,
):
    """Per-band percentile composites (default p10/p90) for a calendar year.

    Used to build the seasonal-amplitude / stability indices.
    Equivalent to ``C.getPercentileComposites``, returns a dict keyed
    ``"p{percentile}"`` (e.g. ``{"p10": image, "p90": image}``).
    """
    import ee

    start = ee.Date.fromYMD(year, 1, 1)
    end = start.advance(1, "year")
    cs_plus = (
        ee.ImageCollection(config.CLOUD_SCORE_COLLECTION)
        .filter(ee.Filter.date(start, end))
        .filter(ee.Filter.bounds(region))
    )
    s2_masked = (
        ee.ImageCollection(config.S2_COLLECTION)
        .filter(ee.Filter.date(start, end))
        .filter(ee.Filter.bounds(region))
        .linkCollection(cs_plus, cs_plus.first().bandNames())
        .map(lambda img: img.updateMask(img.select("cs").gte(cloud_score_threshold)))
        .select("B.*")
    )
    band_names = s2_masked.first().bandNames()
    out = {}
    for p in percentiles:
        img = (
            s2_masked.reduce(ee.Reducer.percentile([p])).rename(band_names).clip(region)
        )
        out[f"p{p}"] = img
    return out

seasonal_composite(start_date, end_date, region, cloud_score_threshold=0.55, cloud_pct_fallback=50)

Cloud-masked median Sentinel-2 composite for an arbitrary date range.

Falls back to a simple CLOUDY_PIXEL_PERCENTAGE filter if Cloud Score+ has no coverage for the window, and to a zero-filled image if there are no scenes at all (keeps downstream band math from failing on sparse early-record years). Equivalent to getSeasonalComposite / C.getSeasonalComposite.

Source code in savana/composites.py

def seasonal_composite(
    start_date: str,
    end_date: str,
    region,
    cloud_score_threshold: float = 0.55,
    cloud_pct_fallback: float = 50,
):
    """Cloud-masked median Sentinel-2 composite for an arbitrary date range.

    Falls back to a simple ``CLOUDY_PIXEL_PERCENTAGE`` filter if Cloud
    Score+ has no coverage for the window, and to a zero-filled image
    if there are no scenes at all (keeps downstream band math from
    failing on sparse early-record years). Equivalent to
    ``getSeasonalComposite`` / ``C.getSeasonalComposite``.
    """
    import ee

    cs_plus = (
        ee.ImageCollection(config.CLOUD_SCORE_COLLECTION)
        .filter(ee.Filter.date(start_date, end_date))
        .filter(ee.Filter.bounds(region))
    )
    s2 = (
        ee.ImageCollection(config.S2_COLLECTION)
        .filter(ee.Filter.date(start_date, end_date))
        .filter(ee.Filter.bounds(region))
    )
    has_scenes = s2.size().gt(0)
    has_cs = cs_plus.size().gt(0)

    with_cs = (
        s2.linkCollection(cs_plus, cs_plus.first().bandNames())
        .map(lambda img: img.updateMask(img.select("cs").gte(cloud_score_threshold)))
        .select("B.*")
        .median()
        .clip(region)
    )
    with_cloud_pct = (
        s2.filter(ee.Filter.lt("CLOUDY_PIXEL_PERCENTAGE", cloud_pct_fallback))
        .select("B.*")
        .median()
        .clip(region)
    )
    empty = (
        ee.Image.constant(0)
        .rename("B8")
        .addBands(ee.Image.constant(0).rename("B4"))
        .addBands(ee.Image.constant(0).rename("B11"))
        .addBands(ee.Image.constant(0).rename("B3"))
        .clip(region)
    )
    return ee.Image(
        ee.Algorithms.If(
            has_scenes,
            ee.Image(ee.Algorithms.If(has_cs, with_cs, with_cloud_pct)),
            empty,
        )
    )

sentinel2_annual(year, region, cloud_score_threshold=0.65)

Cloud-masked median Sentinel-2 SR composite for a calendar year.

Equivalent to getSentinel2Composite / C.getSentinel2Annual.

Source code in savana/composites.py

def sentinel2_annual(year: int, region, cloud_score_threshold: float = 0.65):
    """Cloud-masked median Sentinel-2 SR composite for a calendar year.

    Equivalent to ``getSentinel2Composite`` / ``C.getSentinel2Annual``.
    """
    import ee

    start = ee.Date.fromYMD(year, 1, 1)
    end = start.advance(1, "year")
    s2 = (
        ee.ImageCollection(config.S2_COLLECTION)
        .filter(ee.Filter.date(start, end))
        .filter(ee.Filter.bounds(region))
    )
    cs_plus = ee.ImageCollection(config.CLOUD_SCORE_COLLECTION)
    return (
        s2.linkCollection(cs_plus, cs_plus.first().bandNames())
        .map(lambda img: img.updateMask(img.select("cs").gte(cloud_score_threshold)))
        .select("B.*")
        .median()
        .clip(region)
    )

savana.indices

Spectral index computation and the 14-band phenological stack.

Direct port of indices.js.

build_pheno_stack(idx, rue_img)

Build the 14-band phenological stack used in Model D training.

rue_img must be a single band named 'RUE'.

Source code in savana/indices.py

def build_pheno_stack(idx: dict, rue_img):
    """Build the 14-band phenological stack used in Model D training.

    ``rue_img`` must be a single band named ``'RUE'``.
    """
    import ee

    return ee.Image.cat(
        [
            idx["ndvi_dry"],  # 1
            idx["ndmi_dry"],  # 2
            idx["ndvi_amp"],  # 3
            idx["ndmi"],  # 4
            idx["ndbi"],  # 5
            idx["ndvi"],  # 6
            idx["ndvi_wet"],  # 7
            idx["ndmi_wet"],  # 8
            idx["ndvi_p10"],  # 9
            idx["ndmi_p10"],  # 10
            idx["ndvi_p90"],  # 11
            idx["ndmi_p90"],  # 12
            idx["ndvi_p_amp"],  # 13
            rue_img,  # 14, RUE
        ]
    )

compute(s2_annual, s2_dry, s2_wet, p10, p90)

Compute all spectral indices from seasonal/percentile composites.

Returns a dict, access as idx["ndvi"], idx["ndvi_dry"], etc.

Source code in savana/indices.py

def compute(s2_annual, s2_dry, s2_wet, p10, p90) -> dict:
    """Compute all spectral indices from seasonal/percentile composites.

    Returns a dict, access as ``idx["ndvi"]``, ``idx["ndvi_dry"]``, etc.
    """
    ndvi = s2_annual.normalizedDifference(["B8", "B4"]).rename("NDVI")
    ndmi = s2_annual.normalizedDifference(["B8", "B11"]).rename("NDMI")
    mndwi = s2_annual.normalizedDifference(["B3", "B11"]).rename("MNDWI")
    ndbi = s2_annual.normalizedDifference(["B11", "B8"]).rename("NDBI")
    ndvi_dry = s2_dry.normalizedDifference(["B8", "B4"]).rename("NDVI_dry")
    ndmi_dry = s2_dry.normalizedDifference(["B8", "B11"]).rename("NDMI_dry")
    ndbi_dry = s2_dry.normalizedDifference(["B11", "B8"]).rename("NDBI_dry")
    ndvi_wet = s2_wet.normalizedDifference(["B8", "B4"]).rename("NDVI_wet")
    ndmi_wet = s2_wet.normalizedDifference(["B8", "B11"]).rename("NDMI_wet")
    ndvi_amp = ndvi_wet.subtract(ndvi_dry).rename("NDVI_amp")
    ndvi_p10 = p10.normalizedDifference(["B8", "B4"]).rename("NDVI_p10")
    ndmi_p10 = p10.normalizedDifference(["B8", "B11"]).rename("NDMI_p10")
    ndvi_p90 = p90.normalizedDifference(["B8", "B4"]).rename("NDVI_p90")
    ndmi_p90 = p90.normalizedDifference(["B8", "B11"]).rename("NDMI_p90")
    ndvi_p_amp = ndvi_p90.subtract(ndvi_p10).rename("NDVI_p_amp")

    return {
        "ndvi": ndvi,
        "ndmi": ndmi,
        "mndwi": mndwi,
        "ndbi": ndbi,
        "ndvi_dry": ndvi_dry,
        "ndmi_dry": ndmi_dry,
        "ndbi_dry": ndbi_dry,
        "ndvi_wet": ndvi_wet,
        "ndmi_wet": ndmi_wet,
        "ndvi_amp": ndvi_amp,
        "ndvi_p10": ndvi_p10,
        "ndmi_p10": ndmi_p10,
        "ndvi_p90": ndvi_p90,
        "ndmi_p90": ndmi_p90,
        "ndvi_p_amp": ndvi_p_amp,
    }

savana.rue

Rain Use Efficiency (RUE): integrated NDVI normalised by rainfall.

Tile-boundary bias is removed by normalising integrated NDVI by the count of valid (cloud-free) months before dividing by annual rainfall. Direct port of rue.js.

compute_annual(year, region, cloud_score_threshold=0.65)

Full integrated-NDVI / CHIRPS RUE for a training year.

Returns {"chirps": image, "indvi": image, "rue": image} where rue is a single band named 'RUE' (used as training feature #14).

Source code in savana/rue.py

def compute_annual(year: int, region, cloud_score_threshold: float = 0.65) -> dict:
    """Full integrated-NDVI / CHIRPS RUE for a training year.

    Returns ``{"chirps": image, "indvi": image, "rue": image}`` where
    ``rue`` is a single band named ``'RUE'`` (used as training feature #14).
    """
    import ee

    yr = str(year)
    chirps = (
        ee.ImageCollection(config.CHIRPS_COLLECTION)
        .filter(ee.Filter.date(f"{yr}-01-01", f"{yr}-12-31"))
        .filter(ee.Filter.bounds(region))
        .sum()
        .clip(region)
        .rename("annual_rainfall_mm")
    )

    months = ee.List.sequence(1, 12)

    def _monthly(m):
        start = ee.Date.fromYMD(year, m, 1)
        end = start.advance(1, "month")
        cs_plus_m = (
            ee.ImageCollection(config.CLOUD_SCORE_COLLECTION)
            .filter(ee.Filter.date(start, end))
            .filter(ee.Filter.bounds(region))
        )
        s2_m = (
            ee.ImageCollection(config.S2_COLLECTION)
            .filter(ee.Filter.date(start, end))
            .filter(ee.Filter.bounds(region))
            .linkCollection(cs_plus_m, cs_plus_m.first().bandNames())
            .map(
                lambda img: img.updateMask(img.select("cs").gte(cloud_score_threshold))
            )
        )
        monthly_img = ee.Image(
            ee.Algorithms.If(
                s2_m.size().gt(0),
                s2_m.select(["B8", "B4"])
                .median()
                .normalizedDifference()
                .rename("NDVI")
                .multiply(30),
                ee.Image.constant(0).rename("NDVI").selfMask(),
            )
        )
        return monthly_img.clip(region)

    monthly_ndvi = ee.ImageCollection(months.map(_monthly))
    valid_month_count = monthly_ndvi.count().rename("valid_months")
    indvi = monthly_ndvi.sum().divide(valid_month_count).multiply(12).rename("iNDVI")
    rue = indvi.divide(chirps.add(ee.Image.constant(1))).rename("RUE").clip(region)

    return {"chirps": chirps, "indvi": indvi, "rue": rue}

epoch_rue(year, region, cloud_score_threshold=0.65)

RUE for a specific epoch year. Returns a single band named RUE_{year}.

Source code in savana/rue.py

def epoch_rue(year: int, region, cloud_score_threshold: float = 0.65):
    """RUE for a specific epoch year. Returns a single band named ``RUE_{year}``."""
    import ee

    yr = str(year)
    rainfall = (
        ee.ImageCollection(config.CHIRPS_COLLECTION)
        .filter(ee.Filter.date(f"{yr}-01-01", f"{yr}-12-31"))
        .filter(ee.Filter.bounds(region))
        .sum()
        .rename("rainfall")
        .clip(region)
    )
    cs_plus = (
        ee.ImageCollection(config.CLOUD_SCORE_COLLECTION)
        .filter(ee.Filter.date(f"{yr}-01-01", f"{yr}-12-31"))
        .filter(ee.Filter.bounds(region))
    )
    s2 = (
        ee.ImageCollection(config.S2_COLLECTION)
        .filter(ee.Filter.date(f"{yr}-01-01", f"{yr}-12-31"))
        .filter(ee.Filter.bounds(region))
    )
    annual_ndvi = ee.Image(
        ee.Algorithms.If(
            s2.size().gt(0),
            s2.linkCollection(cs_plus, cs_plus.first().bandNames())
            .map(
                lambda img: img.normalizedDifference(["B8", "B4"])
                .rename("NDVI")
                .updateMask(img.select("cs").gte(cloud_score_threshold))
            )
            .mean()
            .rename("NDVI")
            .clip(region),
            ee.Image.constant(0).rename("NDVI").clip(region),
        )
    )
    return (
        annual_ndvi.multiply(1000)
        .divide(rainfall.add(ee.Image.constant(1)))
        .rename(f"RUE_{yr}")
        .clip(region)
    )

savana.thresholds

Adaptive percentile-based threshold derivation.

All thresholds are derived from this AOI's own index distribution, nothing is hardcoded, so the same code transfers to any savanna landscape automatically. Direct port of thresholds.js.

compute(idx, region)

Compute all classification thresholds from index percentiles.

Returns a dict of ee.Number, access as T["CORE_NDVI_DRY"] etc.

Source code in savana/thresholds.py

def compute(idx: dict, region) -> dict:
    """Compute all classification thresholds from index percentiles.

    Returns a dict of ``ee.Number``, access as ``T["CORE_NDVI_DRY"]`` etc.
    """
    import ee

    pct_image = ee.Image.cat(
        [
            idx["ndvi"],
            idx["ndmi"],
            idx["ndbi"],
            idx["ndvi_dry"],
            idx["ndmi_dry"],
            idx["ndvi_wet"],
            idx["ndmi_wet"],
            idx["ndvi_amp"],
        ]
    )
    pct_dict = pct_image.reduceRegion(
        reducer=ee.Reducer.percentile([5, 10, 25, 50, 75, 90, 95]),
        geometry=region,
        scale=100,
        maxPixels=1e9,
        tileScale=8,
    )

    def p(band, pct):
        return ee.Number(pct_dict.get(f"{band}_p{pct}"))

    return {
        # Anthropogenic, high NDBI or very sparse vegetation
        "ANTHRO_NDBI": p("NDBI", 90),
        "ANTHRO_NDVI_MAX": ee.Number(0.20),
        # Riparian, evergreen gallery forest, high NDMI both seasons
        "RIPARIAN_NDMI": p("NDMI", 95),
        "RIPARIAN_NDMI_DRY": p("NDMI_dry", 75).add(
            p("NDMI_dry", 90).subtract(p("NDMI_dry", 75)).multiply(0.5)
        ),
        "RIPARIAN_NDVI_DRY": p("NDVI_dry", 50).add(
            p("NDVI_dry", 75).subtract(p("NDVI_dry", 50)).multiply(0.5)
        ),
        # Core Woodland, dense closed canopy, top 25% dry-season NDVI
        "CORE_NDVI_DRY": p("NDVI_dry", 75),
        "CORE_NDMI": p("NDMI", 50).add(
            p("NDMI", 75).subtract(p("NDMI", 50)).multiply(0.5)
        ),
        # Grassland, loses canopy in dry season, high amplitude
        "GRASS_NDVI_DRY_MAX": p("NDVI_dry", 25),
        "GRASS_AMP_MIN": p("NDVI_amp", 10).add(
            p("NDVI_amp", 25).subtract(p("NDVI_amp", 10)).multiply(0.75)
        ),
        "GRASS_NDMI_DRY_MAX": p("NDMI_dry", 25).add(
            p("NDMI_dry", 50).subtract(p("NDMI_dry", 25)).multiply(0.75)
        ),
        # Shrub-Transition, intermediate dry NDVI, low dry moisture
        "SHRUB_NDVI_DRY_MIN": p("NDVI_dry", 10),
        "SHRUB_NDVI_DRY_MAX": p("NDVI_dry", 50),
        "SHRUB_NDMI_DRY_MIN": p("NDMI_dry", 25),
        "SHRUB_NDMI_DRY_MAX": p("NDMI_dry", 50).add(
            p("NDMI_dry", 75).subtract(p("NDMI_dry", 50)).multiply(0.20)
        ),
        "SHRUB_AMP_MIN": p("NDVI_amp", 10).add(
            p("NDVI_amp", 25).subtract(p("NDVI_amp", 10)).multiply(0.30)
        ),
        "SHRUB_AMP_MAX": p("NDVI_amp", 90),
        # Open Woodland, intermediate canopy, sub-populations split
        "OPEN_NDVI_DRY_MIN": p("NDVI_dry", 25),
        "OPEN_NDVI_DRY_MAX": p("NDVI_dry", 75),
        "OPEN_NDMI_MIN": p("NDMI", 25),
        "OPEN_NDMI_MAX": p("NDMI", 75),
        "OPEN_NDVI_DRY_MID": p("NDVI_dry", 50),
        "OPEN_NDMI_SPLIT_LOW": p("NDMI", 50),
        "OPEN_NDMI_SPLIT_HIGH": p("NDMI", 50).add(
            p("NDMI", 75).subtract(p("NDMI", 50)).multiply(0.5)
        ),
    }

savana.masks

Mutually exclusive land system masks.

Priority order: Anthropogenic -> Riparian -> Core -> Grassland -> Shrub-Transition -> Open Woodland. Each mask explicitly excludes all higher-priority classes. Direct port of masks.js.

add_layers(masks, m=None, class_info=None)

Add each mask as a layer to a geemap/folium Map (or ee.Map in Code Editor context).

m should be a geemap.Map instance in a notebook. If omitted, a new one is created and returned.

Source code in savana/masks.py

def add_layers(masks: dict, m=None, class_info: dict | None = None):
    """Add each mask as a layer to a geemap/folium Map (or ``ee.Map`` in Code Editor context).

    ``m`` should be a ``geemap.Map`` instance in a notebook. If omitted,
    a new one is created and returned.
    """
    from . import config

    if m is None:
        import geemap

        m = geemap.Map()

    info = class_info or config.DEFAULT_CLASS_INFO
    color_by_name = {
        "core": info[1]["color"],
        "open": info[2]["color"],
        "shrub": info[3]["color"],
        "grass": info[4]["color"],
        "riparian": info[5]["color"],
        "anthro": info[6]["color"],
    }
    labels = {
        "core": "MASK: Core Woodland",
        "open": "MASK: Open Woodland",
        "shrub": "MASK: Shrub-Transition",
        "grass": "MASK: Grassland",
        "riparian": "MASK: Riparian",
        "anthro": "MASK: Anthropogenic",
    }
    for key in ["core", "open", "shrub", "grass", "riparian", "anthro"]:
        m.addLayer(
            masks[key].selfMask(),
            {"palette": [color_by_name[key]]},
            labels[key],
            False,
        )
    return m

compute(idx, T)

Compute all 6 land system masks from indices and thresholds.

Returns {"anthro", "riparian", "core", "grass", "shrub", "open"}.

Source code in savana/masks.py

def compute(idx: dict, T: dict) -> dict:
    """Compute all 6 land system masks from indices and thresholds.

    Returns ``{"anthro", "riparian", "core", "grass", "shrub", "open"}``.
    """
    mask_anthro = (
        idx["ndbi"].gt(T["ANTHRO_NDBI"]).Or(idx["ndvi"].lt(T["ANTHRO_NDVI_MAX"]))
    )

    mask_riparian = (
        idx["ndmi"]
        .gt(T["RIPARIAN_NDMI"])
        .And(idx["ndmi_dry"].gt(T["RIPARIAN_NDMI_DRY"]))
        .And(idx["ndvi_dry"].gt(T["RIPARIAN_NDVI_DRY"]))
        .And(mask_anthro.Not())
    )

    mask_core = (
        idx["ndvi_dry"]
        .gt(T["CORE_NDVI_DRY"])
        .And(idx["ndmi"].gt(T["CORE_NDMI"]))
        .And(mask_anthro.Not())
        .And(mask_riparian.Not())
    )

    mask_grass = (
        idx["ndvi_dry"]
        .lt(T["GRASS_NDVI_DRY_MAX"])
        .And(
            idx["ndvi_amp"]
            .gt(T["GRASS_AMP_MIN"])
            .Or(idx["ndmi_dry"].lt(T["GRASS_NDMI_DRY_MAX"]))
        )
        .And(mask_anthro.Not())
        .And(mask_riparian.Not())
        .And(mask_core.Not())
    )

    mask_shrub = (
        idx["ndvi_dry"]
        .gte(T["SHRUB_NDVI_DRY_MIN"])
        .And(idx["ndvi_dry"].lt(T["SHRUB_NDVI_DRY_MAX"]))
        .And(idx["ndmi_dry"].gte(T["SHRUB_NDMI_DRY_MIN"]))
        .And(idx["ndmi_dry"].lt(T["SHRUB_NDMI_DRY_MAX"]))
        .And(idx["ndvi_amp"].gte(T["SHRUB_AMP_MIN"]))
        .And(idx["ndvi_amp"].lt(T["SHRUB_AMP_MAX"]))
        .And(mask_anthro.Not())
        .And(mask_riparian.Not())
        .And(mask_core.Not())
        .And(mask_grass.Not())
    )

    mask_open_dry = (
        idx["ndvi_dry"]
        .gte(T["OPEN_NDVI_DRY_MIN"])
        .And(idx["ndvi_dry"].lt(T["OPEN_NDVI_DRY_MID"]))
        .And(idx["ndmi"].gte(T["OPEN_NDMI_MIN"]))
        .And(idx["ndmi"].lte(T["OPEN_NDMI_SPLIT_HIGH"]))
        .And(mask_anthro.Not())
        .And(mask_riparian.Not())
        .And(mask_core.Not())
        .And(mask_grass.Not())
        .And(mask_shrub.Not())
    )

    mask_open_moist = (
        idx["ndvi_dry"]
        .gte(T["OPEN_NDVI_DRY_MID"])
        .And(idx["ndvi_dry"].lte(T["OPEN_NDVI_DRY_MAX"]))
        .And(idx["ndmi"].gt(T["OPEN_NDMI_SPLIT_LOW"]))
        .And(idx["ndmi"].lte(T["OPEN_NDMI_MAX"]))
        .And(mask_anthro.Not())
        .And(mask_riparian.Not())
        .And(mask_core.Not())
        .And(mask_grass.Not())
        .And(mask_shrub.Not())
    )

    return {
        "anthro": mask_anthro,
        "riparian": mask_riparian,
        "core": mask_core,
        "grass": mask_grass,
        "shrub": mask_shrub,
        "open": mask_open_dry.Or(mask_open_moist),
    }

coverage(masks, region)

Fraction of AOI covered by each mask (target range ~0.05-0.30 each).

Source code in savana/masks.py

def coverage(masks: dict, region) -> dict:
    """Fraction of AOI covered by each mask (target range ~0.05-0.30 each)."""
    import ee

    out = {}
    for name, mask in masks.items():
        stats = mask.unmask(0).reduceRegion(
            reducer=ee.Reducer.mean(),
            geometry=region,
            scale=100,
            maxPixels=1e9,
            tileScale=8,
        )
        out[name] = stats
    return out

Training and classification

savana.sampling

GCP (ground control point) generation.

Unsupervised clustering -> cluster-stratified candidate sampling -> rule-based provisional labelling with a confidence-margin filter -> class balancing -> embedding extraction. Direct port of sampling.js (and the equivalent Phase 4 logic in kogyae.js).

assign_labels(candidates, T, confidence_margin=config.DEFAULT_CONFIDENCE_MARGIN, class_property=config.CLASS_PROPERTY)

Assign rule-based provisional labels; filter out low-confidence and invalid points.

Returns the filtered, valid-labelled ee.FeatureCollection.

Source code in savana/sampling.py

def assign_labels(
    candidates,
    T: dict,
    confidence_margin: float = config.DEFAULT_CONFIDENCE_MARGIN,
    class_property: str = config.CLASS_PROPERTY,
):
    """Assign rule-based provisional labels; filter out low-confidence and invalid points.

    Returns the filtered, valid-labelled ``ee.FeatureCollection``.
    """
    import ee

    def _label(f):
        ndvi_v = ee.Number(f.get("NDVI"))
        ndmi_v = ee.Number(f.get("NDMI"))
        ndbi_v = ee.Number(f.get("NDBI"))
        ndvi_dry = ee.Number(f.get("NDVI_dry"))
        ndmi_dry = ee.Number(f.get("NDMI_dry"))
        ndvi_amp = ee.Number(ee.Algorithms.If(f.get("NDVI_amp"), f.get("NDVI_amp"), 0))

        is_ambiguous = (
            ndvi_dry.gt(T["CORE_NDVI_DRY"].subtract(confidence_margin))
            .And(ndvi_dry.lt(T["CORE_NDVI_DRY"].add(confidence_margin)))
            .Or(
                ndvi_dry.gt(T["GRASS_NDVI_DRY_MAX"].subtract(confidence_margin)).And(
                    ndvi_dry.lt(T["GRASS_NDVI_DRY_MAX"].add(confidence_margin))
                )
            )
            .Or(
                ndvi_dry.gt(T["OPEN_NDVI_DRY_MIN"].subtract(confidence_margin)).And(
                    ndvi_dry.lt(T["OPEN_NDVI_DRY_MIN"].add(confidence_margin))
                )
            )
        )

        label = ee.Number(
            ee.Algorithms.If(
                ndbi_v.gt(T["ANTHRO_NDBI"]).Or(ndvi_v.lt(T["ANTHRO_NDVI_MAX"])),
                6,
                ee.Algorithms.If(
                    ndmi_v.gt(T["RIPARIAN_NDMI"])
                    .And(ndmi_dry.gt(T["RIPARIAN_NDMI_DRY"]))
                    .And(ndvi_dry.gt(T["RIPARIAN_NDVI_DRY"])),
                    5,
                    ee.Algorithms.If(
                        ndvi_dry.gt(T["CORE_NDVI_DRY"]).And(ndmi_v.gt(T["CORE_NDMI"])),
                        1,
                        ee.Algorithms.If(
                            ndvi_dry.lt(T["GRASS_NDVI_DRY_MAX"]).And(
                                ndvi_amp.gt(T["GRASS_AMP_MIN"]).Or(
                                    ndmi_dry.lt(T["GRASS_NDMI_DRY_MAX"])
                                )
                            ),
                            4,
                            ee.Algorithms.If(
                                ndvi_dry.gte(T["SHRUB_NDVI_DRY_MIN"])
                                .And(ndvi_dry.lt(T["SHRUB_NDVI_DRY_MAX"]))
                                .And(ndmi_dry.gte(T["SHRUB_NDMI_DRY_MIN"]))
                                .And(ndmi_dry.lt(T["SHRUB_NDMI_DRY_MAX"]))
                                .And(ndvi_amp.gte(T["SHRUB_AMP_MIN"]))
                                .And(ndvi_amp.lt(T["SHRUB_AMP_MAX"])),
                                3,
                                ee.Algorithms.If(
                                    ndvi_dry.gte(T["OPEN_NDVI_DRY_MIN"])
                                    .And(ndvi_dry.lte(T["OPEN_NDVI_DRY_MAX"]))
                                    .And(ndmi_v.gte(T["OPEN_NDMI_MIN"]))
                                    .And(ndmi_v.lte(T["OPEN_NDMI_MAX"])),
                                    2,
                                    -1,
                                ),
                            ),
                        ),
                    ),
                ),
            )
        )

        is_rare = label.eq(3).Or(label.eq(5))
        final_label = ee.Number(
            ee.Algorithms.If(is_rare, label, ee.Algorithms.If(is_ambiguous, 99, label))
        )
        return f.set(class_property, final_label)

    labelled = candidates.map(_label)
    valid = labelled.filter(ee.Filter.gt(class_property, 0)).filter(
        ee.Filter.neq(class_property, 99)
    )
    return valid

build_gcps(embedding, idx, clusters, T, region, n_classes=6, points_per_class=config.DEFAULT_POINTS_PER_CLASS, scale=config.DEFAULT_EXPORT_SCALE, class_property=config.CLASS_PROPERTY, n_clusters=config.DEFAULT_N_CLUSTERS, candidates_per_cluster=config.DEFAULT_CANDIDATES_PER_CLUSTER, confidence_margin=config.DEFAULT_CONFIDENCE_MARGIN)

End-to-end: sample candidates -> label -> balance -> extract embeddings.

Returns the final ee.FeatureCollection of ground control points with embedding bands attached, ready for classifier training.

Source code in savana/sampling.py

def build_gcps(
    embedding,
    idx: dict,
    clusters,
    T: dict,
    region,
    n_classes: int = 6,
    points_per_class: int = config.DEFAULT_POINTS_PER_CLASS,
    scale: int = config.DEFAULT_EXPORT_SCALE,
    class_property: str = config.CLASS_PROPERTY,
    n_clusters: int = config.DEFAULT_N_CLUSTERS,
    candidates_per_cluster: int = config.DEFAULT_CANDIDATES_PER_CLUSTER,
    confidence_margin: float = config.DEFAULT_CONFIDENCE_MARGIN,
):
    """End-to-end: sample candidates -> label -> balance -> extract embeddings.

    Returns the final ``ee.FeatureCollection`` of ground control points
    with embedding bands attached, ready for classifier training.
    """
    import ee

    candidates = sample_candidates(
        idx,
        clusters,
        region,
        n_clusters=n_clusters,
        candidates_per_cluster=candidates_per_cluster,
        scale=scale,
    )
    labelled = assign_labels(
        candidates,
        T,
        confidence_margin=confidence_margin,
        class_property=class_property,
    )

    def take(fc, class_val, seed):
        sub = fc.filter(ee.Filter.eq(class_property, class_val))
        return ee.FeatureCollection(
            ee.Algorithms.If(
                sub.size().gte(points_per_class),
                sub.randomColumn("pick", seed).sort("pick").limit(points_per_class),
                sub,
            )
        )

    balanced = None
    for i, class_val in enumerate(range(1, n_classes + 1)):
        chunk = take(labelled, class_val, 101 + i)
        balanced = chunk if balanced is None else balanced.merge(chunk)

    gcps = embedding.sampleRegions(
        collection=balanced.filterBounds(region),
        properties=[
            class_property,
            "cluster6",
            "NDVI",
            "NDMI",
            "NDBI",
            "NDVI_dry",
            "NDMI_dry",
            "NDVI_wet",
            "NDMI_wet",
            "NDVI_amp",
        ],
        scale=scale,
        geometries=True,
        tileScale=8,
    ).filter(ee.Filter.notNull(embedding.bandNames()))

    return gcps

cluster_embedding(embedding, region, n_clusters=config.DEFAULT_N_CLUSTERS, seed=config.DEFAULT_RANDOM_SEED)

Unsupervised k-means clustering in AlphaEarth embedding space.

Returns {"clusters": image, "samples": feature_collection}.

Source code in savana/sampling.py

def cluster_embedding(
    embedding,
    region,
    n_clusters: int = config.DEFAULT_N_CLUSTERS,
    seed: int = config.DEFAULT_RANDOM_SEED,
):
    """Unsupervised k-means clustering in AlphaEarth embedding space.

    Returns ``{"clusters": image, "samples": feature_collection}``.
    """
    import ee

    samples = embedding.sample(
        region=region, scale=30, numPixels=3000, seed=seed, tileScale=8
    )
    kmeans = ee.Clusterer.wekaKMeans(nClusters=n_clusters, seed=seed).train(samples)
    clusters = embedding.cluster(kmeans).toInt().rename("cluster6")
    return {"clusters": clusters, "samples": samples}

sample_candidates(idx, clusters, region, n_clusters=config.DEFAULT_N_CLUSTERS, candidates_per_cluster=config.DEFAULT_CANDIDATES_PER_CLUSTER, scale=config.DEFAULT_EXPORT_SCALE)

Cluster-stratified candidate point sampling.

Source code in savana/sampling.py

def sample_candidates(
    idx: dict,
    clusters,
    region,
    n_clusters: int = config.DEFAULT_N_CLUSTERS,
    candidates_per_cluster: int = config.DEFAULT_CANDIDATES_PER_CLUSTER,
    scale: int = config.DEFAULT_EXPORT_SCALE,
):
    """Cluster-stratified candidate point sampling."""
    import ee

    candidate_image = ee.Image.cat(
        [
            idx["ndvi"],
            idx["ndmi"],
            idx["mndwi"],
            idx["ndbi"],
            idx["ndvi_dry"],
            idx["ndmi_dry"],
            idx["ndbi_dry"],
            idx["ndvi_wet"],
            idx["ndmi_wet"],
            idx["ndvi_amp"],
            clusters,
        ]
    ).clip(region)

    keep_props = [
        "cluster6",
        "NDVI",
        "NDMI",
        "NDBI",
        "NDVI_dry",
        "NDMI_dry",
        "NDBI_dry",
        "NDVI_wet",
        "NDMI_wet",
        "NDVI_amp",
    ]

    def sample_cluster(cluster_id, seed):
        return (
            candidate_image.updateMask(clusters.eq(cluster_id))
            .sample(
                region=region,
                scale=scale,
                numPixels=candidates_per_cluster * 10,
                seed=seed,
                geometries=True,
                tileScale=8,
            )
            .filter(
                ee.Filter.notNull(
                    [
                        "cluster6",
                        "NDVI",
                        "NDMI",
                        "NDBI",
                        "NDVI_dry",
                        "NDMI_dry",
                        "NDVI_wet",
                        "NDMI_wet",
                    ]
                )
            )
            .randomColumn("pick", seed)
            .sort("pick")
            .limit(candidates_per_cluster)
            .map(lambda f: ee.Feature(f.geometry(), {p: f.get(p) for p in keep_props}))
        )

    seeds = [100 + 10 * i for i in range(n_clusters)]
    candidates = ee.FeatureCollection(
        [sample_cluster(i, seeds[i]) for i in range(n_clusters)]
    ).flatten()
    return candidates

savana.classifiers

Four-model ablation training and multi-epoch classification.

Ported from kogyae.js Phase 5 (training), Phase 7 (accuracy split), and Phase 8 (multi-epoch classification), the standalone classifers.js module referenced by mainrun.js was empty, so this reconstructs it with the exact signatures mainrun.js expects (trainAllModels, classifyAllEpochs).

Models

A: KNN (k=3) | AlphaEarth embeddings only [baseline] B: Random Forest | AlphaEarth embeddings only C: Random Forest | Phenology indices only [CIRCULAR, ablation only, not used operationally, since labels were derived from the same indices] D: Random Forest | Embeddings + Phenology [PRIMARY, used for mapping]

classify_all_epochs(epochs, models, region, park_name='AOI', embedding_current_year=None, embedding_current_image=None, phenology_min_year=config.DEFAULT_PHENOLOGY_MIN_YEAR, smooth_radius_px=1)

Classify every epoch year with the appropriate master classifier.

Years >= phenology_min_year use Model D (embeddings + phenology, recomputed for that specific year). Earlier years, where seasonal Sentinel-2 coverage is typically too sparse for reliable phenology , fall back to Model B (embeddings only). This generalises the hardcoded "if year === 2017" special case in the original script.

embedding_current_year/embedding_current_image: if one of the epochs is the same year the embedding used for training was already computed for, pass it in to avoid recomputing it.

Returns {year: classified_image}.

Source code in savana/classifiers.py

def classify_all_epochs(
    epochs: list[int],
    models: dict,
    region,
    park_name: str = "AOI",
    embedding_current_year: int | None = None,
    embedding_current_image=None,
    phenology_min_year: int = config.DEFAULT_PHENOLOGY_MIN_YEAR,
    smooth_radius_px: int = 1,
) -> dict:
    """Classify every epoch year with the appropriate master classifier.

     Years >= ``phenology_min_year`` use Model D (embeddings + phenology,
     recomputed for that specific year). Earlier years, where seasonal
     Sentinel-2 coverage is typically too sparse for reliable phenology
    , fall back to Model B (embeddings only). This generalises the
     hardcoded "if year === 2017" special case in the original script.

     ``embedding_current_year``/``embedding_current_image``: if one of
     the epochs is the same year the embedding used for training was
     already computed for, pass it in to avoid recomputing it.

     Returns ``{year: classified_image}``.
    """
    import ee

    class_property = models["class_property"]
    classified_maps = {}

    for year in epochs:
        if embedding_current_year is not None and year == embedding_current_year:
            emb_img = embedding_current_image
        else:
            emb_img = composites.embedding_image(year, region)

        if year < phenology_min_year:
            classified = (
                emb_img.classify(models["master_b"])
                .rename(class_property)
                .clip(region)
                .toByte()
                .focal_mode(
                    radius=smooth_radius_px, units="pixels", kernelType="square"
                )
                .rename(class_property)
                .clip(region)
                .toByte()
                .set("year", year)
                .set("park", park_name)
                .set("classifier", "ModelB_RF_Embeddings")
                .set("system:time_start", ee.Date.fromYMD(year, 1, 1).millis())
            )
        else:
            yr = str(year)
            s2_annual_yr = composites.sentinel2_annual(year, region)
            s2_dry_yr = composites.seasonal_composite(
                f"{yr}-03-01", f"{yr}-05-15", region
            )
            s2_wet_yr = composites.seasonal_composite(
                f"{yr}-05-01", f"{yr}-07-15", region
            )
            pcts_yr = composites.percentile_composites(year, region)

            idx_yr = indices.compute(
                s2_annual_yr, s2_dry_yr, s2_wet_yr, pcts_yr["p10"], pcts_yr["p90"]
            )
            idx_yr = {k: v.unmask(0) for k, v in idx_yr.items()}
            rue_yr = rue_mod.epoch_rue(year, region).select([f"RUE_{yr}"]).rename("RUE")

            pheno_stack_yr = indices.build_pheno_stack(idx_yr, rue_yr)
            combined_img_yr = ee.Image.cat([emb_img, pheno_stack_yr])

            classified = (
                combined_img_yr.classify(models["master_d"])
                .rename(class_property)
                .clip(region)
                .toByte()
                .focal_mode(
                    radius=smooth_radius_px, units="pixels", kernelType="square"
                )
                .rename(class_property)
                .clip(region)
                .toByte()
                .set("year", year)
                .set("park", park_name)
                .set("classifier", "ModelD_RF_Embeddings_Phenology")
                .set("system:time_start", ee.Date.fromYMD(year, 1, 1).millis())
            )

        classified_maps[year] = classified

    return classified_maps

train_all_models(gcps, embedding, idx, rue_img, region, class_property=config.CLASS_PROPERTY, n_trees=150, split_fraction=0.7, seed=config.DEFAULT_RANDOM_SEED, class_order=None)

Train the 4-model ablation and the two "master" classifiers used for mapping (Model B for years without reliable phenology, Model D for all other years).

Returns a dict with the trained classifiers, error matrices, and the exact band lists each classifier expects (so downstream classification always matches training feature order).

Source code in savana/classifiers.py

def train_all_models(
    gcps,
    embedding,
    idx: dict,
    rue_img,
    region,
    class_property: str = config.CLASS_PROPERTY,
    n_trees: int = 150,
    split_fraction: float = 0.7,
    seed: int = config.DEFAULT_RANDOM_SEED,
    class_order: list[int] | None = None,
) -> dict:
    """Train the 4-model ablation and the two "master" classifiers used
    for mapping (Model B for years without reliable phenology, Model D
    for all other years).

    Returns a dict with the trained classifiers, error matrices, and
    the exact band lists each classifier expects (so downstream
    classification always matches training feature order).
    """
    import ee

    class_order = class_order or [1, 2, 3, 4, 5, 6]

    training_data = embedding.sampleRegions(
        collection=gcps,
        properties=[class_property],
        scale=config.DEFAULT_EXPORT_SCALE,
        tileScale=8,
    ).filter(ee.Filter.notNull(embedding.bandNames()))

    with_random = gcps.randomColumn("split", seed)
    train_set = with_random.filter(ee.Filter.lt("split", split_fraction))
    valid_set = with_random.filter(ee.Filter.gte("split", split_fraction))

    pheno_stack = indices.build_pheno_stack(idx, rue_img)
    pheno_bands = pheno_stack.bandNames()
    combined_stack = ee.Image.cat([embedding, pheno_stack])
    combined_bands = combined_stack.bandNames()

    train_full = combined_stack.sampleRegions(
        collection=train_set,
        properties=[class_property],
        scale=config.DEFAULT_EXPORT_SCALE,
        tileScale=8,
    ).filter(ee.Filter.notNull(combined_bands))
    valid_full = combined_stack.sampleRegions(
        collection=valid_set,
        properties=[class_property],
        scale=config.DEFAULT_EXPORT_SCALE,
        tileScale=8,
    ).filter(ee.Filter.notNull(combined_bands))

    emb_bands = embedding.bandNames()
    class_prop_list = ee.List([class_property])
    train_emb = train_full.select(emb_bands.cat(class_prop_list))
    valid_emb = valid_full.select(emb_bands.cat(class_prop_list))
    train_pheno = train_full.select(pheno_bands.cat(class_prop_list))
    valid_pheno = valid_full.select(pheno_bands.cat(class_prop_list))

    # Model A, KNN (k=3) | Embeddings only [baseline]
    model_a = ee.Classifier.smileKNN(3).train(
        features=train_emb, classProperty=class_property, inputProperties=emb_bands
    )
    matrix_a = valid_emb.classify(model_a).errorMatrix(
        actual=class_property, predicted="classification", order=class_order
    )

    # Model B, RF (150 trees) | Embeddings only
    model_b = ee.Classifier.smileRandomForest(
        numberOfTrees=n_trees,
        variablesPerSplit=8,
        minLeafPopulation=1,
        bagFraction=0.632,
        seed=seed,
    ).train(features=train_emb, classProperty=class_property, inputProperties=emb_bands)
    matrix_b = valid_emb.classify(model_b).errorMatrix(
        actual=class_property, predicted="classification", order=class_order
    )

    # Model C, RF (150 trees) | Phenology only [CIRCULAR, ablation diagnostic only]
    model_c = ee.Classifier.smileRandomForest(
        numberOfTrees=n_trees,
        variablesPerSplit=4,
        minLeafPopulation=1,
        bagFraction=0.632,
        seed=seed,
    ).train(
        features=train_pheno, classProperty=class_property, inputProperties=pheno_bands
    )
    matrix_c = valid_pheno.classify(model_c).errorMatrix(
        actual=class_property, predicted="classification", order=class_order
    )

    # Model D, RF (150 trees) | Embeddings + Phenology [PRIMARY]
    model_d = ee.Classifier.smileRandomForest(
        numberOfTrees=n_trees,
        variablesPerSplit=9,
        minLeafPopulation=1,
        bagFraction=0.632,
        seed=seed,
    ).train(
        features=train_full,
        classProperty=class_property,
        inputProperties=combined_bands,
    )
    matrix_d = valid_full.classify(model_d).errorMatrix(
        actual=class_property, predicted="classification", order=class_order
    )

    # "Master" classifiers used for actual epoch mapping, trained on the
    # FULL gcps/trainFull sets (not just the 70% split) for best final quality.
    master_b = ee.Classifier.smileRandomForest(
        numberOfTrees=n_trees,
        variablesPerSplit=8,
        minLeafPopulation=1,
        bagFraction=0.632,
        seed=seed,
    ).train(
        features=training_data,
        classProperty=class_property,
        inputProperties=embedding.bandNames(),
    )

    model_d_band_names = embedding.bandNames().cat(ee.List(indices.PHENO_BAND_ORDER))
    master_d = ee.Classifier.smileRandomForest(
        numberOfTrees=n_trees,
        variablesPerSplit=9,
        minLeafPopulation=1,
        bagFraction=0.632,
        seed=seed,
    ).train(
        features=train_full,
        classProperty=class_property,
        inputProperties=model_d_band_names,
    )

    return {
        "model_a": model_a,
        "matrix_a": matrix_a,
        "model_b": model_b,
        "matrix_b": matrix_b,
        "model_c": model_c,
        "matrix_c": matrix_c,
        "model_d": model_d,
        "matrix_d": matrix_d,
        "master_b": master_b,
        "master_d": master_d,
        "embedding_bands": embedding.bandNames(),
        "model_d_band_names": model_d_band_names,
        "class_property": class_property,
    }

savana.change

Conservative change detection + Rain-Use-Efficiency inter-annual variability.

Direct port of change.js. Distinguishes genuine structural change from rainfall-driven apparent change (a common false-positive source in savanna change detection) using RUE coefficient-of-variation as a filter.

analyse(classified_maps, epochs, region, park_name='AOI', rue_cv_threshold=0.15)

Run conservative change detection across all epochs.

Requires exactly the epochs present as keys in classified_maps; the "conservative" and RUE checks specifically use the first and last epoch, plus stability through any provided middle epochs.

Returns a dict with the change stack, conservative/genuine/variable change masks, transition codes, RUE-CV image, and per-epoch RUE images (keyed rue_{year}).

Source code in savana/change.py

def analyse(
    classified_maps: dict,
    epochs: list[int],
    region,
    park_name: str = "AOI",
    rue_cv_threshold: float = 0.15,
) -> dict:
    """Run conservative change detection across all epochs.

    Requires exactly the epochs present as keys in ``classified_maps``;
    the "conservative" and RUE checks specifically use the first and
    last epoch, plus stability through any provided middle epochs.

    Returns a dict with the change stack, conservative/genuine/variable
    change masks, transition codes, RUE-CV image, and per-epoch RUE
    images (keyed ``rue_{year}``).
    """
    import ee

    epochs_sorted = sorted(epochs)
    first_year, last_year = epochs_sorted[0], epochs_sorted[-1]
    stats_scale = _auto_stats_scale(region)

    def band_name(y):
        return f"ls_{y}"

    change_stack = ee.Image.cat(
        [classified_maps[y].rename(band_name(y)) for y in epochs_sorted]
    )

    # Conservative change: stable at both ends of the sequence, but
    # different overall (requires >= 3 epochs to be meaningful; with
    # exactly 2 epochs this reduces to a simple pairwise change mask).
    if len(epochs_sorted) >= 4:
        second_year, second_last_year = epochs_sorted[1], epochs_sorted[-2]
        stable_early = change_stack.select(band_name(first_year)).eq(
            change_stack.select(band_name(second_year))
        )
        stable_late = change_stack.select(band_name(second_last_year)).eq(
            change_stack.select(band_name(last_year))
        )
        conservative_change = (
            stable_early.And(stable_late)
            .And(
                change_stack.select(band_name(first_year)).neq(
                    change_stack.select(band_name(last_year))
                )
            )
            .rename("conservative_change")
        )
        stable_throughout = ee.Image(1).clip(region)
        for a, b in zip(epochs_sorted[:-1], epochs_sorted[1:]):
            stable_throughout = stable_throughout.And(
                change_stack.select(band_name(a)).eq(change_stack.select(band_name(b)))
            )
        stable_throughout = stable_throughout.rename("stable_all_epochs")
    else:
        conservative_change = (
            change_stack.select(band_name(first_year))
            .neq(change_stack.select(band_name(last_year)))
            .rename("conservative_change")
        )
        stable_throughout = conservative_change.Not().rename("stable_all_epochs")

    conservative_transition = (
        change_stack.select(band_name(first_year))
        .multiply(10)
        .add(change_stack.select(band_name(last_year)))
        .updateMask(conservative_change)
        .rename("conservative_transition")
    )

    # RUE inter-annual variability across the same epochs.
    rue_images = {y: rue_mod.epoch_rue(y, region) for y in epochs_sorted}
    rue_stack = ee.Image.cat(list(rue_images.values()))
    rue_cv = (
        rue_stack.reduce(ee.Reducer.stdDev())
        .divide(rue_stack.reduce(ee.Reducer.mean()))
        .rename("RUE_CV")
    )

    genuine_change = conservative_change.And(rue_cv.lt(rue_cv_threshold)).rename(
        f"genuine_change_{first_year}_{last_year}"
    )
    variable_change = conservative_change.And(rue_cv.gte(rue_cv_threshold)).rename(
        f"variable_change_{first_year}_{last_year}"
    )

    result = {
        "change_stack": change_stack,
        "conservative_change": conservative_change,
        "conservative_transition": conservative_transition,
        "stable_throughout": stable_throughout,
        "genuine_change": genuine_change,
        "variable_change": variable_change,
        "rue_cv": rue_cv,
        "stats_scale": stats_scale,
        "first_year": first_year,
        "last_year": last_year,
    }
    for y, img in rue_images.items():
        result[f"rue_{y}"] = img
    return result

class_area_stats(classified_maps, epochs, region, scale)

Per-epoch class area statistics (km2), grouped by class code.

Source code in savana/change.py

def class_area_stats(classified_maps: dict, epochs: list[int], region, scale) -> dict:
    """Per-epoch class area statistics (km2), grouped by class code."""
    import ee

    out = {}
    for year in epochs:
        groups = (
            ee.Image.pixelArea()
            .divide(1e6)
            .addBands(classified_maps[year])
            .reduceRegion(
                reducer=ee.Reducer.sum().group(groupField=1, groupName="landSystem"),
                geometry=region,
                scale=scale,
                maxPixels=1e10,
                tileScale=8,
            )
        )
        out[year] = groups
    return out

Accuracy and export

savana.accuracy

Confusion matrix and accuracy summary reporting.

Ports the CSV structure from accuracy.js (24-row full confusion matrix across 4 models x 6 classes, plus a 4-row per-model summary), but returns pandas.DataFrame directly for notebook use, Drive/CSV export is available separately via :mod:savana.exports.

confusion_matrix_dataframe(matrices, park_name='AOI', class_info=None)

Build the full 24-row (4 models x N classes) confusion matrix table.

matrices maps model key ("a","b","c","d") to an ee.ConfusionMatrix.

Source code in savana/accuracy.py

def confusion_matrix_dataframe(
    matrices: dict, park_name: str = "AOI", class_info: dict | None = None
):
    """Build the full 24-row (4 models x N classes) confusion matrix table.

    ``matrices`` maps model key ("a","b","c","d") to an ``ee.ConfusionMatrix``.
    """
    import pandas as pd

    info = class_info or config.DEFAULT_CLASS_INFO
    classes = _class_labels(info)

    rows = []
    for key, cm in matrices.items():
        model_name, _, _ = MODEL_LABELS.get(key, (key, key, key))
        oa = cm.accuracy().getInfo()
        kappa = cm.kappa().getInfo()
        pa = cm.producersAccuracy().getInfo()
        ua = cm.consumersAccuracy().getInfo()
        arr = cm.array().getInfo()
        for i, (code, label) in enumerate(classes):
            row = {
                "park": park_name,
                "model": model_name,
                "actual_class_code": code,
                "actual_class_label": label,
                "overall_accuracy": oa,
                "kappa": kappa,
                "producer_accuracy": pa[i][0] if i < len(pa) else None,
                "user_accuracy": ua[0][i] if i < len(ua[0]) else None,
            }
            for j, (_, pred_label) in enumerate(classes):
                row[f"pred_{pred_label}"] = (
                    arr[i][j] if i < len(arr) and j < len(arr[i]) else None
                )
            rows.append(row)
    return pd.DataFrame(rows)

print_summary(matrices, park_name='AOI')

Console summary mirroring the ablation-comparison prints in kogyai.js.

Source code in savana/accuracy.py

def print_summary(matrices: dict, park_name: str = "AOI") -> None:
    """Console summary mirroring the ablation-comparison prints in kogyai.js."""
    print(f"--- ACCURACY: {park_name} ---")
    for key in ["a", "b", "c", "d"]:
        if key not in matrices:
            continue
        cm = matrices[key]
        note = " [CIRCULAR, inflated, diagnostic only]" if key == "c" else ""
        oa = cm.accuracy().getInfo()
        kappa = cm.kappa().getInfo()
        print(f"Model {key.upper()} | OA: {oa:.4f} | Kappa: {kappa:.4f}{note}")
    print("")
    print("Interpretation guide:")
    print("  B > A  -> RF outperforms KNN on the same embeddings")
    print("  D > B  -> phenology adds value beyond embeddings alone")
    print("  D > C  -> embeddings add value beyond indices alone")
    print("  C is circular (labels derived from the same indices), diagnostic only")
    print("  Model D (primary) is used for all epoch mapping.")

summary_dataframe(matrices, park_name='AOI', class_info=None)

One row per model with overall accuracy, kappa, and per-class PA/UA.

Source code in savana/accuracy.py

def summary_dataframe(
    matrices: dict, park_name: str = "AOI", class_info: dict | None = None
):
    """One row per model with overall accuracy, kappa, and per-class PA/UA."""
    import pandas as pd

    info = class_info or config.DEFAULT_CLASS_INFO
    classes = _class_labels(info)

    rows = []
    for key, cm in matrices.items():
        model_code, model_desc, feature_space = MODEL_LABELS.get(key, (key, key, key))
        oa = cm.accuracy().getInfo()
        kappa = cm.kappa().getInfo()
        pa = cm.producersAccuracy().getInfo()
        ua = cm.consumersAccuracy().getInfo()
        row = {
            "park": park_name,
            "model_code": key.upper(),
            "model_description": model_desc,
            "feature_space": feature_space,
            "overall_accuracy": oa,
            "kappa": kappa,
        }
        for i, (_, label) in enumerate(classes):
            row[f"PA_{label}"] = pa[i][0] if i < len(pa) else None
            row[f"UA_{label}"] = ua[0][i] if i < len(ua[0]) else None
        rows.append(row)
    return pd.DataFrame(rows)

savana.exports

Export helpers: Google Drive, Earth Engine Assets, and local CSV/GeoTIFF.

Ports the export logic in exports.js, generalised so the folder/asset path, CRS, and scale are all parameters instead of hardcoded to one researcher's Drive folder and asset project.

class_areas_dataframe(classified_maps, epochs, region, scale)

Client-side (getInfo) per-epoch class area table as a pandas DataFrame.

Source code in savana/exports.py

def class_areas_dataframe(classified_maps: dict, epochs: list[int], region, scale: int):
    """Client-side (getInfo) per-epoch class area table as a pandas DataFrame."""
    import pandas as pd

    from . import change as change_mod

    stats = change_mod.class_area_stats(classified_maps, epochs, region, scale)
    rows = []
    for year, groups_ee in stats.items():
        groups = groups_ee.getInfo().get("groups", [])
        for g in groups:
            rows.append(
                {"year": year, "landSystem": g["landSystem"], "area_km2": g["sum"]}
            )
    return pd.DataFrame(rows)

export_change_products(chg, region, park_name='AOI', drive_folder=None, asset_folder=None, scale=config.DEFAULT_EXPORT_SCALE, crs=config.DEFAULT_CRS, start=True)

Export change-detection and RUE products to Drive and/or Assets.

Source code in savana/exports.py

def export_change_products(
    chg: dict,
    region,
    park_name: str = "AOI",
    drive_folder: str | None = None,
    asset_folder: str | None = None,
    scale: int = config.DEFAULT_EXPORT_SCALE,
    crs: str = config.DEFAULT_CRS,
    start: bool = True,
) -> list:
    """Export change-detection and RUE products to Drive and/or Assets."""
    import ee

    tasks = []

    def _drive(image, name_suffix, sub_scale=None):
        if not drive_folder:
            return
        task = ee.batch.Export.image.toDrive(
            image=image,
            description=f"{park_name}_{name_suffix}",
            folder=drive_folder,
            fileNamePrefix=f"{park_name}_{name_suffix.lower()}",
            region=region,
            scale=sub_scale or scale,
            crs=crs,
            maxPixels=1e10,
        )
        if start:
            task.start()
        tasks.append(task)

    def _asset(image, name_suffix, sub_scale=None):
        if not asset_folder:
            return
        task = ee.batch.Export.image.toAsset(
            image=image,
            description=f"{park_name}_{name_suffix}_Asset",
            assetId=f"{asset_folder}/{park_name}/{name_suffix.lower()}",
            region=region,
            scale=sub_scale or scale,
            crs=crs,
            maxPixels=1e10,
            pyramidingPolicy={".default": "MODE"},
        )
        if start:
            task.start()
        tasks.append(task)

    _drive(chg["conservative_change"].toByte(), "ConservativeChange")
    _drive(chg["conservative_transition"].toByte(), "ConservativeTransition")
    _asset(chg["conservative_change"].toByte(), "ConservativeChange")

    _drive(chg["genuine_change"].toByte(), "GenuineChange")
    _asset(chg["genuine_change"].toByte(), "GenuineChange")

    _drive(chg["rue_cv"].toFloat(), "RUE_CV", sub_scale=100)
    _asset(chg["rue_cv"].toFloat(), "RUE_CV", sub_scale=100)

    return tasks

export_classified_maps(classified_maps, epochs, region, park_name='AOI', drive_folder=None, asset_folder=None, scale=config.DEFAULT_EXPORT_SCALE, crs=config.DEFAULT_CRS, start=True)

Export each epoch's classified map to Drive and/or an EE Asset folder.

At least one of drive_folder / asset_folder should be given, or nothing will be exported. Returns the list of started EE tasks.

Source code in savana/exports.py

def export_classified_maps(
    classified_maps: dict,
    epochs: list[int],
    region,
    park_name: str = "AOI",
    drive_folder: str | None = None,
    asset_folder: str | None = None,
    scale: int = config.DEFAULT_EXPORT_SCALE,
    crs: str = config.DEFAULT_CRS,
    start: bool = True,
) -> list:
    """Export each epoch's classified map to Drive and/or an EE Asset folder.

    At least one of ``drive_folder`` / ``asset_folder`` should be given,
    or nothing will be exported. Returns the list of started EE tasks.
    """
    import ee

    tasks = []
    for year in epochs:
        img = classified_maps[year]
        if drive_folder:
            task = ee.batch.Export.image.toDrive(
                image=img,
                description=f"{park_name}_LandSystem_{year}",
                folder=drive_folder,
                fileNamePrefix=f"{park_name}_land_system_{year}",
                region=region,
                scale=scale,
                crs=crs,
                maxPixels=1e10,
            )
            if start:
                task.start()
            tasks.append(task)
        if asset_folder:
            task = ee.batch.Export.image.toAsset(
                image=img,
                description=f"{park_name}_LandSystem_Asset_{year}",
                assetId=f"{asset_folder}/{park_name}/land_system_{year}",
                region=region,
                scale=scale,
                crs=crs,
                maxPixels=1e10,
                pyramidingPolicy={".default": "MODE"},
            )
            if start:
                task.start()
            tasks.append(task)
    return tasks

Insights and visualisation

savana.insights

Grounded analytical insights: knowledge derived only from real results.

The design principle here is deliberate: every number in summarize() and answer() traces back to something actually computed by the pipeline (class_areas(), accuracy_summary(), the change-detection stats), never invented, interpolated, or guessed. compute_facts() is the single source of truth; both text-producing functions only ever read from it. This keeps savana's reporting honest even as it grows, if a future version adds LLM-phrased summaries, that layer should sit on top of these same facts, never replace them.

answer(facts, question)

Answer a natural-language question using only precomputed facts.

This is deliberately simple keyword matching, not an LLM, it can only ever report numbers that are actually in facts, so it cannot hallucinate a result the pipeline didn't produce. Questions it doesn't recognise get an honest "don't know" rather than a guess.

Source code in savana/insights.py

def answer(facts: dict, question: str) -> str:
    """Answer a natural-language question using only precomputed facts.

    This is deliberately simple keyword matching, not an LLM, it can
    only ever report numbers that are actually in ``facts``, so it
    cannot hallucinate a result the pipeline didn't produce. Questions
    it doesn't recognise get an honest "don't know" rather than a guess.
    """
    q = question.lower()
    epochs = facts["epochs"]
    class_names = list(facts["class_names"].values())

    # Which epoch is being asked about? Default to the most recent.
    year = next((y for y in epochs if str(y) in q), epochs[-1])

    # Does the question mention a specific class?
    matched_class = next((name for name in class_names if name.lower() in q), None)

    if any(w in q for w in ["change", "lost", "gained", "convert"]) and facts.get(
        "change"
    ):
        chg = facts["change"]
        return (
            f"Between {chg['first_year']} and {chg['last_year']}, "
            f"{chg['genuine_change_pct_of_area']}% of {facts['park_name']}'s area shows genuine "
            f"structural land-system change; a further {chg['variable_change_pct_of_area']}% shows "
            f"apparent change that's more likely rainfall-driven variability than real conversion."
        )

    if any(w in q for w in ["accura", "model", "reliab", "confidence"]) and facts.get(
        "accuracy"
    ):
        acc = facts["accuracy"]
        return (
            f"The primary classification model (AlphaEarth embeddings + phenology) achieved "
            f"{acc['primary_model_accuracy']:.1%} overall accuracy "
            f"(kappa {acc['primary_model_kappa']:.3f}) for {facts['park_name']}."
        )

    if matched_class:
        area = facts["area_by_epoch"].get(year, {}).get(matched_class)
        pct = facts["pct_by_epoch"].get(year, {}).get(matched_class)
        if area is not None:
            return (
                f"In {year}, {matched_class} covered {area:.1f} km2 "
                f"({pct}% of {facts['park_name']})."
            )
        return f"No {matched_class} area was found for {year} in the computed results."

    if any(w in q for w in ["dominant", "most", "largest", "majority"]):
        dominant = facts["dominant_class"].get(year)
        pct = facts["pct_by_epoch"].get(year, {}).get(dominant)
        return (
            f"The dominant land-system class in {year} was {dominant} "
            f"({pct}% of {facts['park_name']})."
        )

    if any(w in q for w in ["total", "area", "size", "how big"]):
        total = facts["total_area_km2"].get(year)
        return f"The total classified area of {facts['park_name']} in {year} was {total:.1f} km2."

    return (
        "I can only answer from what the pipeline actually computed, try asking about "
        "a specific class's area, the dominant class, overall accuracy, or change between years."
    )

compute_facts(clf)

Extract a structured dict of real, computed facts from a fitted classifier.

Requires clf.run() (or at least .classify()) to have completed. This is the "knowledge base", everything else in this module reads from its output, never from the raw ee.Image objects directly.

Each section (area, accuracy, change) is computed independently and guarded against Earth Engine timeouts, a slow/large AOI causing one section to time out will not prevent the others from returning. Any section that fails is set to None and noted in facts["warnings"] rather than raising, since a partial, honest answer is better than a hard crash on results that mostly did compute successfully.

Source code in savana/insights.py

def compute_facts(clf) -> dict:
    """Extract a structured dict of real, computed facts from a fitted classifier.

    Requires ``clf.run()`` (or at least ``.classify()``) to have completed.
    This is the "knowledge base", everything else in this module reads
    from its output, never from the raw ee.Image objects directly.

    Each section (area, accuracy, change) is computed independently and
    guarded against Earth Engine timeouts, a slow/large AOI causing one
    section to time out will not prevent the others from returning. Any
    section that fails is set to ``None`` and noted in ``facts["warnings"]``
    rather than raising, since a partial, honest answer is better than a
    hard crash on results that mostly did compute successfully.
    """
    info = clf.class_info
    class_names = {code: v["name"] for code, v in info.items()}
    epochs = sorted(clf.epochs)

    facts: dict = {
        "park_name": clf.park_name,
        "epochs": epochs,
        "class_names": class_names,
        "area_by_epoch": {},  # {year: {class_name: km2}}
        "pct_by_epoch": {},  # {year: {class_name: pct_of_total}}
        "total_area_km2": {},  # {year: total_km2}
        "dominant_class": {},  # {year: class_name}
        "accuracy": None,
        "change": None,
        "warnings": [],
    }

    # --- Area stats (per-epoch class areas) ---
    try:
        areas_df = clf.class_areas()
        for year in epochs:
            year_rows = areas_df[areas_df["year"] == year]
            by_class = {}
            for _, row in year_rows.iterrows():
                code = int(row["landSystem"])
                name = class_names.get(code, f"class_{code}")
                by_class[name] = float(row["area_km2"])
            total = sum(by_class.values())
            facts["area_by_epoch"][year] = by_class
            facts["total_area_km2"][year] = total
            pct = {}
            if total > 0:
                for name, area in by_class.items():
                    pct[name] = round(100 * area / total, 1)
            facts["pct_by_epoch"][year] = pct
            facts["dominant_class"][year] = (
                max(by_class, key=by_class.get) if by_class else None
            )
    except (
        Exception
    ) as exc:  # noqa: BLE001 - deliberately broad: any EE failure here is non-fatal
        facts["warnings"].append(
            f"Area statistics unavailable ({type(exc).__name__}: {exc})."
        )

    # --- Accuracy, best model by overall accuracy, plus the primary model (D) specifically ---
    try:
        acc_df = clf.accuracy_summary()
        if acc_df is not None and len(acc_df) > 0:
            best_row = acc_df.loc[acc_df["overall_accuracy"].idxmax()]
            primary_row = acc_df[acc_df["model_code"] == "D"]
            primary_oa = None
            primary_kappa = None
            if len(primary_row):
                primary_oa = round(float(primary_row.iloc[0]["overall_accuracy"]), 4)
                primary_kappa = round(float(primary_row.iloc[0]["kappa"]), 4)
            facts["accuracy"] = {
                "best_model_code": best_row["model_code"],
                "best_model_description": best_row["model_description"],
                "best_overall_accuracy": round(float(best_row["overall_accuracy"]), 4),
                "best_kappa": round(float(best_row["kappa"]), 4),
                "primary_model_accuracy": primary_oa,
                "primary_model_kappa": primary_kappa,
            }
    except Exception as exc:  # noqa: BLE001
        facts["warnings"].append(
            f"Accuracy statistics unavailable ({type(exc).__name__}: {exc})."
        )

    # --- Change detection, only attempted if >= 2 epochs were run ---
    if clf.change is not None:
        try:
            chg = clf.change
            first_year, last_year = chg["first_year"], chg["last_year"]
            # Single combined image + single reduceRegion call instead of two
            # separate ones, halves the round trips to Earth Engine for this section.
            combined = (
                chg["conservative_change"]
                .unmask(0)
                .rename("conservative")
                .addBands(chg["genuine_change"].unmask(0).rename("genuine"))
            )
            stats = combined.reduceRegion(
                reducer=_mean_reducer(),
                geometry=clf.region,
                scale=chg["stats_scale"],
                maxPixels=1e10,
                tileScale=8,
                bestEffort=True,
            ).getInfo()
            conservative_frac = stats.get("conservative", 0) or 0
            genuine_frac = stats.get("genuine", 0) or 0
            facts["change"] = {
                "first_year": first_year,
                "last_year": last_year,
                "conservative_change_pct_of_area": round(100 * conservative_frac, 2),
                "genuine_change_pct_of_area": round(100 * genuine_frac, 2),
                "variable_change_pct_of_area": round(
                    100 * (conservative_frac - genuine_frac), 2
                ),
            }
        except Exception as exc:  # noqa: BLE001
            facts["warnings"].append(
                f"Change statistics unavailable ({type(exc).__name__}: {exc})."
            )

    return facts

summarize(facts)

Turn a facts dict into a plain-English narrative report.

Source code in savana/insights.py

def summarize(facts: dict) -> str:
    """Turn a facts dict into a plain-English narrative report."""
    lines = []
    park = facts["park_name"]
    epochs = facts["epochs"]

    lines.append(
        f"Land-system classification summary for {park} ({', '.join(map(str, epochs))}):"
    )
    lines.append("")

    for year in epochs:
        total = facts["total_area_km2"].get(year, 0)
        dominant = facts["dominant_class"].get(year)
        pct = facts["pct_by_epoch"].get(year, {})
        lines.append(f"{year}: total classified area {total:.1f} km2.")
        if dominant:
            lines.append(
                f"  Dominant class: {dominant} ({pct.get(dominant, 0)}% of the area)."
            )
        for name, p in sorted(pct.items(), key=lambda kv: -kv[1]):
            area = facts["area_by_epoch"][year].get(name, 0)
            lines.append(f"  - {name}: {area:.1f} km2 ({p}%)")
        lines.append("")

    if facts.get("accuracy"):
        acc = facts["accuracy"]
        if acc["primary_model_accuracy"] is not None:
            lines.append(
                f"Model accuracy: the primary model (embeddings + phenology) reached "
                f"{acc['primary_model_accuracy']:.1%} overall accuracy "
                f"(kappa {acc['primary_model_kappa']:.3f})."
            )
        else:
            lines.append("Model accuracy: primary model results unavailable.")
        lines.append(
            f"Best-performing model overall: {acc['best_model_description']} "
            f"({acc['best_overall_accuracy']:.1%} accuracy, kappa {acc['best_kappa']:.3f})."
        )
        lines.append("")

    if facts.get("change"):
        chg = facts["change"]
        lines.append(
            f"Change {chg['first_year']} to {chg['last_year']}: "
            f"{chg['conservative_change_pct_of_area']}% of the area shows conservative "
            f"(stable-to-stable) land-system change."
        )
        lines.append(
            f"  Of that, {chg['genuine_change_pct_of_area']}% of the total area is genuine "
            f"structural change (validated against rainfall variability), while "
            f"{chg['variable_change_pct_of_area']}% appears to be rainfall-driven apparent "
            f"change rather than true land-system conversion."
        )

    if facts.get("warnings"):
        lines.append("")
        lines.append("Note: some sections could not be computed and are omitted above:")
        for w in facts["warnings"]:
            lines.append(f"  - {w}")

    return "\n".join(lines)

savana.viz

Interactive visualization helpers built on geemap/leafmap for Jupyter.

add_legend(m, class_info=None, title='Land System Classes')

Add a class legend to a geemap.Map.

Source code in savana/viz.py

def add_legend(m, class_info: dict | None = None, title: str = "Land System Classes"):
    """Add a class legend to a geemap.Map."""
    info = class_info or config.DEFAULT_CLASS_INFO
    legend_dict = {v["name"]: f"#{v['color']}" for v in info.values()}
    m.add_legend(title=title, legend_dict=legend_dict)
    return m

compare_split_map(left, right, left_label='Left', right_label='Right', region=None, class_info=None, m=None, zoom=12)

Side-by-side swipe comparison between two layers.

Each of left/right can be either: - an ee.Image (e.g. a classified year, or clf.maps[2019]), rendered with the land-system palette/legend - a basemap name string (e.g. "SATELLITE", "HYBRID", "ROADMAP", "Esri.WorldImagery"), passed straight to geemap

Drag the handle in the middle of the map to swipe between them, this also works for "classified year vs. underlying satellite imagery" by passing a basemap name string as one side.

Source code in savana/viz.py

def compare_split_map(
    left,
    right,
    left_label: str = "Left",
    right_label: str = "Right",
    region=None,
    class_info: dict | None = None,
    m=None,
    zoom: int = 12,
):
    """Side-by-side swipe comparison between two layers.

    Each of ``left``/``right`` can be either:
        - an ``ee.Image`` (e.g. a classified year, or ``clf.maps[2019]``),
          rendered with the land-system palette/legend
        - a basemap name string (e.g. ``"SATELLITE"``, ``"HYBRID"``,
          ``"ROADMAP"``, ``"Esri.WorldImagery"``), passed straight to geemap

    Drag the handle in the middle of the map to swipe between them,
    this also works for "classified year vs. underlying satellite
    imagery" by passing a basemap name string as one side.
    """
    import ee
    import geemap

    info = class_info or config.DEFAULT_CLASS_INFO

    def _to_layer(side, label):
        if isinstance(side, str):
            return side  # basemap name, geemap resolves this itself
        if isinstance(side, ee.Image):
            return geemap.ee_tile_layer(side, config.class_vis_params(info), label)
        raise TypeError(
            "compare_split_map sides must be an ee.Image or a basemap "
            f"name string, got {type(side)!r}"
        )

    left_layer = _to_layer(left, left_label)
    right_layer = _to_layer(right, right_label)

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

    m.split_map(left_layer=left_layer, right_layer=right_layer)
    return m

show_change_map(chg, region=None, m=None, zoom=12)

Display conservative/genuine/variable change layers on a geemap.Map.

Source code in savana/viz.py

def show_change_map(chg: dict, region=None, m=None, zoom: int = 12):
    """Display conservative/genuine/variable change layers on a geemap.Map."""
    import geemap

    if m is None:
        m = geemap.Map()
    if region is not None:
        m.centerObject(region, zoom)
    m.add_layer(
        chg["conservative_change"].selfMask(),
        {"palette": ["8b0000"]},
        "Conservative change",
        True,
    )
    m.add_layer(
        chg["genuine_change"].selfMask(),
        {"palette": ["d73027"]},
        "Genuine structural change",
        False,
    )
    m.add_layer(
        chg["variable_change"].selfMask(),
        {"palette": ["fc8d59"]},
        "Rainfall-driven apparent change",
        False,
    )
    m.add_layer(
        chg["rue_cv"],
        {"min": 0, "max": 0.3, "palette": ["1a9641", "ffffbf", "d73027"]},
        "RUE coefficient of variation",
        False,
    )
    return m

show_classified_map(classified_image, region=None, class_info=None, m=None, zoom=12)

Display a classified land-system image on an interactive geemap.Map.

Source code in savana/viz.py

def show_classified_map(
    classified_image,
    region=None,
    class_info: dict | None = None,
    m=None,
    zoom: int = 12,
):
    """Display a classified land-system image on an interactive geemap.Map."""
    import geemap

    info = class_info or config.DEFAULT_CLASS_INFO
    if m is None:
        m = geemap.Map()
    if region is not None:
        m.centerObject(region, zoom)
    m.add_layer(
        classified_image, config.class_vis_params(info), "Land System Classification"
    )
    add_legend(m, info)
    return m

show_gcps(gcps, region=None, class_info=None, class_property=config.CLASS_PROPERTY, background=None, m=None, zoom=12, point_size=5)

Display ground control points on a map, colored by assigned class.

Lets you visually sanity-check the sampling/labelling step, where the training points actually landed, and whether their classes look spatially sensible, before trusting the classifier trained on them.

Parameters:

Name	Type	Description	Default
gcps		The ee.FeatureCollection of ground control points (e.g. clf.gcps), with class_property set on each feature.	required
region		AOI to center the map on (e.g. clf.region).	None
class_info	dict | None	Class scheme (defaults to the standard 6-class one).	None
class_property	str	Property name holding the class code on each point (defaults to savana's standard "landSystem").	CLASS_PROPERTY
background		Optional ee.Image to show underneath the points (e.g. a classified year, or a Sentinel-2 composite), makes it easier to judge whether points look correctly placed.	None
m		Existing geemap.Map to add to, or a new one is created.	None
point_size	int	Marker size in pixels.	5

Source code in savana/viz.py

def show_gcps(
    gcps,
    region=None,
    class_info: dict | None = None,
    class_property: str = config.CLASS_PROPERTY,
    background=None,
    m=None,
    zoom: int = 12,
    point_size: int = 5,
):
    """Display ground control points on a map, colored by assigned class.

    Lets you visually sanity-check the sampling/labelling step, where
    the training points actually landed, and whether their classes look
    spatially sensible, before trusting the classifier trained on them.

    Args:
        gcps: The ``ee.FeatureCollection`` of ground control points
            (e.g. ``clf.gcps``), with ``class_property`` set on each
            feature.
        region: AOI to center the map on (e.g. ``clf.region``).
        class_info: Class scheme (defaults to the standard 6-class one).
        class_property: Property name holding the class code on each
            point (defaults to savana's standard ``"landSystem"``).
        background: Optional ee.Image to show underneath the points
            (e.g. a classified year, or a Sentinel-2 composite), makes
            it easier to judge whether points look correctly placed.
        m: Existing geemap.Map to add to, or a new one is created.
        point_size: Marker size in pixels.
    """
    import ee
    import geemap

    info = class_info or config.DEFAULT_CLASS_INFO
    if m is None:
        m = geemap.Map()
    if region is not None:
        m.centerObject(region, zoom)

    if background is not None:
        m.add_layer(background, config.class_vis_params(info), "Background", True, 0.6)

    for code, entry in sorted(info.items()):
        class_points = gcps.filter(ee.Filter.eq(class_property, code))
        styled = class_points.style(color=entry["color"], pointSize=point_size)
        m.add_layer(styled, {}, f"GCPs: {entry['name']}")

    add_legend(m, info, title="Ground Control Points")
    return m

show_multi_year_map(classified_maps, years=None, region=None, class_info=None, m=None, zoom=12)

Add every requested epoch as its own toggleable layer on one map.

Uses geemap's built-in layer panel, each year gets its own checkbox, so you can flip between them (or view several at once with opacity sliders) without re-running anything.

Source code in savana/viz.py

def show_multi_year_map(
    classified_maps: dict,
    years: list[int] | None = None,
    region=None,
    class_info: dict | None = None,
    m=None,
    zoom: int = 12,
):
    """Add every requested epoch as its own toggleable layer on one map.

    Uses geemap's built-in layer panel, each year gets its own checkbox,
    so you can flip between them (or view several at once with opacity
    sliders) without re-running anything.
    """
    import geemap

    info = class_info or config.DEFAULT_CLASS_INFO
    years = sorted(years or classified_maps.keys())
    if m is None:
        m = geemap.Map()
    if region is not None:
        m.centerObject(region, zoom)
    for year in years:
        m.add_layer(
            classified_maps[year], config.class_vis_params(info), f"Land System {year}"
        )
    add_legend(m, info)
    return m

savana.viz_geolibre

Display savana classification results inside the GeoLibre Jupyter widget.

GeoLibre (https://geolibre.app, MIT-licensed, by the same author as geoai) ships a Jupyter-native Python package with a leafmap-style API. This module is a thin adapter: it converts savana's ee.Image outputs into XYZ tile URLs (via Earth Engine's own tile server) that GeoLibre's add_tile_layer() can display, the same way you'd add any other raster tile source.

This is intentionally a lightweight integration, it depends only on the public geolibre PyPI package, versioned and released the same way as every other savana dependency. It is not a GeoLibre plugin (that would be TypeScript code living inside GeoLibre's own repo/build system); see the project roadmap for that as a possible future, separate effort.

Requires: pip install "savana[geolibre]" (needs Python >= 3.11, since that is GeoLibre's own minimum, this is stricter than savana's core Python >= 3.10 requirement).

Known limitation: GeoLibre's swipe/compare tool is currently a UI-only plugin (Plugins menu > Swipe) with no scriptable Python entry point yet. Once GeoLibre exposes one, a compare_geolibre() will be added here to match savana.viz.compare_split_map().

show_classified_map(classified_image, region=None, class_info=None, m=None, zoom=10, layout='embed')

Display one classified land-system image inside the GeoLibre widget.

Returns a geolibre.Map, display it in a notebook cell by putting it as the last expression, same as any other Jupyter widget.

Source code in savana/viz_geolibre.py

def show_classified_map(
    classified_image,
    region=None,
    class_info: dict | None = None,
    m=None,
    zoom: int = 10,
    layout: str = "embed",
):
    """Display one classified land-system image inside the GeoLibre widget.

    Returns a ``geolibre.Map``, display it in a notebook cell by putting
    it as the last expression, same as any other Jupyter widget.
    """
    from geolibre import Map

    info = class_info or config.DEFAULT_CLASS_INFO
    tile_url = _ee_image_to_tile_url(classified_image, config.class_vis_params(info))

    if m is None:
        if region is not None:
            lng, lat = _region_center(region)
            m = Map(center=(lng, lat), zoom=zoom, layout=layout)
        else:
            m = Map(layout=layout)
    elif region is not None:
        lng, lat = _region_center(region)
        m.set_center(lng, lat, zoom=zoom)

    m.add_tile_layer(tile_url, name="Land System Classification")
    return m

show_multi_year_map(classified_maps, years=None, region=None, class_info=None, m=None, zoom=10, layout='embed')

Add every requested epoch as its own toggleable tile layer in GeoLibre.

Each year appears as its own entry in GeoLibre's Layers panel, with its own visibility checkbox and opacity slider (the same panel you already saw in the app), no extra code needed on your end to toggle between them once this cell has run.

Source code in savana/viz_geolibre.py

def show_multi_year_map(
    classified_maps: dict,
    years: list[int] | None = None,
    region=None,
    class_info: dict | None = None,
    m=None,
    zoom: int = 10,
    layout: str = "embed",
):
    """Add every requested epoch as its own toggleable tile layer in GeoLibre.

    Each year appears as its own entry in GeoLibre's Layers panel, with its
    own visibility checkbox and opacity slider (the same panel you already
    saw in the app), no extra code needed on your end to toggle between
    them once this cell has run.
    """
    from geolibre import Map

    info = class_info or config.DEFAULT_CLASS_INFO
    years = sorted(years or classified_maps.keys())

    if m is None:
        if region is not None:
            lng, lat = _region_center(region)
            m = Map(center=(lng, lat), zoom=zoom, layout=layout)
        else:
            m = Map(layout=layout)
    elif region is not None:
        lng, lat = _region_center(region)
        m.set_center(lng, lat, zoom=zoom)

    vis_params = config.class_vis_params(info)
    for year in years:
        tile_url = _ee_image_to_tile_url(classified_maps[year], vis_params)
        m.add_tile_layer(tile_url, name=f"Land System {year}")
    return m

Natural-language agent

savana.agents

SavanaGeoAgent: the one agent class for savana - grounded Q&A, map control, and a chat UI, all in one place.

This is built directly ON TOP of geoai's own agent infrastructure (geoai.agents) rather than a parallel reimplementation: the map is a real geoai.Map (leafmap/MapLibre-based, the same class geoai's own demos use), map control comes from geoai's real, full-featured MapTools (fly_to, add_basemap, add_vector, add_raster, add_cog_layer, remove_layer, and more), and model creation reuses geoai's own create_anthropic_model/create_openai_model/create_gemini_model. savana adds its own grounded Q&A tools (summarize, class_area, accuracy, change, etc.) alongside geoai's map tools on one combined Strands agent.

from savana.agents import SavanaGeoAgent

agent = SavanaGeoAgent(clf, model="anthropic")
agent.ask("How much core woodland is there in 2024?")   # savana Q&A
agent.ask("Fly to the study area and add a satellite basemap")  # geoai map tools
agent.show_ui()                                          # chat UI + live map, inline

Requires: pip install "savana[agents]" (installs geoai-py[agents], which brings in strands-agents, leafmap, and the LLM provider SDKs).

SavanaGeoAgent

The one agent class for savana: grounded Q&A + full map control + chat UI.

Built on geoai's real Map/MapTools/model-factory infrastructure (see module docstring), savana adds its own grounded query tools alongside geoai's map-control tools on one combined agent.

Deliberately one class, not two, pass clf, rainfall, or both. Whichever you pass determines which grounded tool set(s) get loaded, so someone working on both a land-system classification and a rainfall assessment for the same study area gets one agent and one ask(), not two agents to keep track of.

Parameters:

Name	Type	Description	Default
clf		A SavanaClassifier that has already run. Optional if rainfall is given.	None
rainfall		A savana.rainfall.pipeline.RainfallAssessment that has already been scored (.score() called). Optional if clf is given.	None
model	str	Either a provider name ("anthropic", "openai", "gemini", "ollama", uses that provider's env-var API key, or a local Ollama server, and a sensible default model id) or an already-built Strands model instance.	'anthropic'
model_id	Optional[str]	Optional explicit model id, used only when model is a provider name string.	None
map_instance		Optional existing geoai.Map (leafmap/MapLibre) to control. If omitted, geoai creates a default one.	None
max_tokens	int	Explicit max output tokens for the Anthropic provider specifically, always set explicitly here (not left to provider defaults), since omitting it is what causes a bare KeyError: 'max_tokens' in some Strands/Anthropic version combinations.	4096
**model_kwargs	Any	Passed through to geoai's model factory.	{}

Source code in savana/agents.py

class SavanaGeoAgent:
    """The one agent class for savana: grounded Q&A + full map control + chat UI.

    Built on geoai's real ``Map``/``MapTools``/model-factory infrastructure
    (see module docstring), savana adds its own grounded query tools
    alongside geoai's map-control tools on one combined agent.

    Deliberately one class, not two, pass ``clf``, ``rainfall``, or
    both. Whichever you pass determines which grounded tool set(s) get
    loaded, so someone working on both a land-system classification and
    a rainfall assessment for the same study area gets one agent and
    one ``ask()``, not two agents to keep track of.

    Args:
        clf: A ``SavanaClassifier`` that has already run. Optional if
            ``rainfall`` is given.
        rainfall: A ``savana.rainfall.pipeline.RainfallAssessment`` that
            has already been scored (``.score()`` called). Optional if
            ``clf`` is given.
        model: Either a provider name (``"anthropic"``, ``"openai"``,
            ``"gemini"``, ``"ollama"``, uses that provider's env-var API
            key, or a local Ollama server, and a sensible default model
            id) or an already-built Strands model instance.
        model_id: Optional explicit model id, used only when ``model`` is
            a provider name string.
        map_instance: Optional existing ``geoai.Map`` (leafmap/MapLibre)
            to control. If omitted, geoai creates a default one.
        max_tokens: Explicit max output tokens for the Anthropic provider
            specifically, always set explicitly here (not left to
            provider defaults), since omitting it is what causes a bare
            ``KeyError: 'max_tokens'`` in some Strands/Anthropic version
            combinations.
        **model_kwargs: Passed through to geoai's model factory.
    """

    def __init__(
        self,
        clf=None,
        rainfall=None,
        model: str = "anthropic",
        model_id: Optional[str] = None,
        map_instance=None,
        max_tokens: int = 4096,
        **model_kwargs: Any,
    ):
        if clf is None and rainfall is None:
            raise ValueError(
                "SavanaGeoAgent needs at least one of clf= (a fitted "
                "SavanaClassifier) or rainfall= (a scored RainfallAssessment)."
            )
        _require_geoai_agents()
        from geoai.agents import MapTools
        from geoai.agents.map_tools import MapSession
        from strands import Agent

        # Shared chat state, a plain agent.ask("...") call in any cell
        # and typing into show_ui()'s own text box both write here, so
        # whichever is currently displayed stays in sync with the other.
        self._history: list[str] = []
        self._chat_output = None  # set by show_ui() once displayed

        # Import each provider's model factory individually, not every
        # installed geoai version has every provider (e.g. some older
        # versions lack create_gemini_model), so a missing one shouldn't
        # block using a provider that IS available.
        factories: dict = {}
        try:
            from geoai.agents import create_anthropic_model

            factories["anthropic"] = lambda **kw: create_anthropic_model(
                max_tokens=max_tokens, **kw
            )
        except ImportError:
            pass
        try:
            from geoai.agents import create_openai_model

            factories["openai"] = create_openai_model
        except ImportError:
            pass
        try:
            from geoai.agents import create_gemini_model

            factories["gemini"] = create_gemini_model
        except ImportError:
            pass
        try:
            from geoai.agents import create_ollama_model

            factories["ollama"] = create_ollama_model
        except ImportError:
            pass

        # Real geoai map + map tools, not a savana-specific reimplementation.
        self._session = MapSession(map_instance)
        self._map_tools = MapTools(self._session)

        # Add the layer-toggle panel ONCE, up front, it's a live,
        # reactive MapLibre control that automatically tracks every layer
        # added afterward by any tool. Calling it again per-layer (an
        # earlier version of this code did that) risks stacking duplicate
        # panels instead of just staying in sync.
        try:
            self._session.m.add_layer_control()
        except Exception as exc:  # noqa: BLE001
            import warnings

            warnings.warn(
                f"Could not add the layer toggle panel: {type(exc).__name__}: {exc}",
                stacklevel=2,
            )

        self._qa_tools = _SavanaQATools(clf) if clf is not None else None
        self._savana_map_tools = (
            _SavanaMapTools(clf, self._session) if clf is not None else None
        )
        self._rainfall_qa_tools = (
            _RainfallQATools(rainfall) if rainfall is not None else None
        )

        if isinstance(model, str) and model.lower() in factories:
            kwargs = dict(model_kwargs)
            if model_id:
                kwargs["model_id"] = model_id
            model_instance = factories[model.lower()](**kwargs)
        elif isinstance(model, str):
            raise ValueError(
                f"Provider {model!r} is not available (either unknown, or your "
                f"installed geoai version doesn't export its model factory). "
                f"Available in this environment: {list(factories)}, "
                "or pass an already-built Strands model instance."
            )
        else:
            model_instance = model  # assume caller passed a real Strands model

        map_tool_names = [
            "fly_to",
            "create_map",
            "zoom_to",
            "jump_to",
            "add_basemap",
            "add_vector",
            "add_raster",
            "add_cog_layer",
            "remove_layer",
            "get_layer_names",
            "set_terrain",
            "remove_terrain",
            "add_overture_3d_buildings",
            "set_paint_property",
            "set_layout_property",
            "set_color",
            "set_opacity",
            "set_visibility",
            "add_marker",
            "set_pitch",
        ]
        map_tools = [getattr(self._map_tools, name) for name in map_tool_names]

        system_prompt = _GEOAI_MAP_SYSTEM_PROMPT
        combined_tools = list(map_tools)
        if self._qa_tools is not None:
            system_prompt += _SAVANA_PROMPT_ADDENDUM
            combined_tools += (
                self._qa_tools.build_tools() + self._savana_map_tools.build_tools()
            )
        if self._rainfall_qa_tools is not None:
            system_prompt += _RAINFALL_PROMPT_ADDENDUM
            combined_tools += self._rainfall_qa_tools.build_tools()

        agent_name = "Savana Agent"
        if clf is not None and rainfall is not None:
            agent_name = "Savana Land-System + Rainfall Agent"
        elif rainfall is not None:
            agent_name = "Savana Rainfall Agent"
        else:
            agent_name = "Savana Land-System Agent"

        self._agent = Agent(
            name=agent_name,
            model=model_instance,
            system_prompt=system_prompt,
            tools=combined_tools,
            callback_handler=None,
        )

    @property
    def map(self):
        """The live geoai.Map (leafmap/MapLibre) this agent controls."""
        return self._session.m

    def ask(self, prompt: str) -> str:
        """Send a single-turn question, get a plain-text answer back.

        If show_ui() is currently displayed, this also updates it,
        asking from a plain cell and typing into the UI box both write
        to the same visible chat log.
        """
        self._history.append(f"You: {prompt}")
        self._history.append("Agent is thinking...")
        self._render_chat()
        try:
            result = self._agent(prompt)
            answer = getattr(result, "final_text", str(result))
        except Exception as exc:  # noqa: BLE001
            self._history.pop()
            self._history.append(f"Agent error: {type(exc).__name__}: {exc}")
            self._history.append("")
            self._render_chat()
            raise
        self._history.pop()
        self._history.append(f"Agent: {answer}")
        self._history.append("")
        self._render_chat()
        return answer

    def _render_chat(self):
        """Redraw the show_ui() chat panel, if one is currently displayed."""
        if self._chat_output is None:
            return
        self._chat_output.clear_output(wait=True)
        with self._chat_output:
            for line in self._history:
                print(line)

    def __call__(self, prompt: str):
        """Full Strands result object (same as calling the agent directly)."""
        return self._agent(prompt)

    def show_ui(self, height: int = 600):
        """Display the live geoai map + a chat box side by side, inline.

        Calling ``agent.ask(...)`` in a separate cell also updates this
        panel, if it's currently displayed, they share the same chat log.

        Requires: ``ipywidgets`` (installed with the ``agents`` extra).
        """
        try:
            import ipywidgets as widgets
            from IPython.display import display
        except ImportError as exc:  # pragma: no cover
            raise ImportError(
                "show_ui() requires ipywidgets. Install with: pip install ipywidgets"
            ) from exc

        map_panel = widgets.VBox(
            [widgets.HTML("<b>Map</b>"), self.map],
            layout=widgets.Layout(
                flex="1 1 0%", min_width="480px", height=f"{height}px"
            ),
        )

        self._chat_output = widgets.Output(
            layout=widgets.Layout(
                border="1px solid #ccc",
                padding="8px",
                height=f"{height - 60}px",
                overflow_y="auto",
            )
        )
        self._render_chat()  # show anything already asked before show_ui() was called

        text_box = widgets.Text(
            placeholder="Ask about your results, or ask to fly/add basemap/compare...",
            layout=widgets.Layout(width="80%"),
        )
        send_button = widgets.Button(description="Send", button_style="primary")

        def _send(_=None):
            question = text_box.value.strip()
            if not question:
                return
            text_box.value = ""
            try:
                self.ask(question)
            except Exception:  # noqa: BLE001
                pass  # already recorded to chat by ask() itself

        send_button.on_click(_send)
        text_box.on_submit(_send)

        chat_panel = widgets.VBox(
            [
                widgets.HTML("<b>Chat</b>"),
                self._chat_output,
                widgets.HBox([text_box, send_button]),
            ],
            layout=widgets.Layout(flex="1 1 0%", min_width="360px"),
        )

        display(widgets.HBox([map_panel, chat_panel]))

map property

The live geoai.Map (leafmap/MapLibre) this agent controls.

__call__(prompt)

Full Strands result object (same as calling the agent directly).

Source code in savana/agents.py

def __call__(self, prompt: str):
    """Full Strands result object (same as calling the agent directly)."""
    return self._agent(prompt)

ask(prompt)

Send a single-turn question, get a plain-text answer back.

If show_ui() is currently displayed, this also updates it, asking from a plain cell and typing into the UI box both write to the same visible chat log.

Source code in savana/agents.py

def ask(self, prompt: str) -> str:
    """Send a single-turn question, get a plain-text answer back.

    If show_ui() is currently displayed, this also updates it,
    asking from a plain cell and typing into the UI box both write
    to the same visible chat log.
    """
    self._history.append(f"You: {prompt}")
    self._history.append("Agent is thinking...")
    self._render_chat()
    try:
        result = self._agent(prompt)
        answer = getattr(result, "final_text", str(result))
    except Exception as exc:  # noqa: BLE001
        self._history.pop()
        self._history.append(f"Agent error: {type(exc).__name__}: {exc}")
        self._history.append("")
        self._render_chat()
        raise
    self._history.pop()
    self._history.append(f"Agent: {answer}")
    self._history.append("")
    self._render_chat()
    return answer

show_ui(height=600)

Display the live geoai map + a chat box side by side, inline.

Calling agent.ask(...) in a separate cell also updates this panel, if it's currently displayed, they share the same chat log.

Requires: ipywidgets (installed with the agents extra).

Source code in savana/agents.py

def show_ui(self, height: int = 600):
    """Display the live geoai map + a chat box side by side, inline.

    Calling ``agent.ask(...)`` in a separate cell also updates this
    panel, if it's currently displayed, they share the same chat log.

    Requires: ``ipywidgets`` (installed with the ``agents`` extra).
    """
    try:
        import ipywidgets as widgets
        from IPython.display import display
    except ImportError as exc:  # pragma: no cover
        raise ImportError(
            "show_ui() requires ipywidgets. Install with: pip install ipywidgets"
        ) from exc

    map_panel = widgets.VBox(
        [widgets.HTML("<b>Map</b>"), self.map],
        layout=widgets.Layout(
            flex="1 1 0%", min_width="480px", height=f"{height}px"
        ),
    )

    self._chat_output = widgets.Output(
        layout=widgets.Layout(
            border="1px solid #ccc",
            padding="8px",
            height=f"{height - 60}px",
            overflow_y="auto",
        )
    )
    self._render_chat()  # show anything already asked before show_ui() was called

    text_box = widgets.Text(
        placeholder="Ask about your results, or ask to fly/add basemap/compare...",
        layout=widgets.Layout(width="80%"),
    )
    send_button = widgets.Button(description="Send", button_style="primary")

    def _send(_=None):
        question = text_box.value.strip()
        if not question:
            return
        text_box.value = ""
        try:
            self.ask(question)
        except Exception:  # noqa: BLE001
            pass  # already recorded to chat by ask() itself

    send_button.on_click(_send)
    text_box.on_submit(_send)

    chat_panel = widgets.VBox(
        [
            widgets.HTML("<b>Chat</b>"),
            self._chat_output,
            widgets.HBox([text_box, send_button]),
        ],
        layout=widgets.Layout(flex="1 1 0%", min_width="360px"),
    )

    display(widgets.HBox([map_panel, chat_panel]))

Configuration and Earth Engine

savana.config

Default configuration for savanna landscape classification.

Everything here is a default, every public function in the package accepts overrides, so a user classifying a different savanna system with a different class scheme is not locked into these values.

class_palette(class_info=None)

Ordered hex palette (by ascending class code) for map visualisation.

Source code in savana/config.py

def class_palette(class_info: dict[int, dict[str, str]] | None = None) -> list[str]:
    """Ordered hex palette (by ascending class code) for map visualisation."""
    info = class_info or DEFAULT_CLASS_INFO
    return [info[k]["color"] for k in sorted(info)]

class_vis_params(class_info=None)

Earth Engine visualization params for a classified land-system image.

Source code in savana/config.py

def class_vis_params(class_info: dict[int, dict[str, str]] | None = None) -> dict:
    """Earth Engine visualization params for a classified land-system image."""
    info = class_info or DEFAULT_CLASS_INFO
    keys = sorted(info)
    return {"min": min(keys), "max": max(keys), "palette": class_palette(info)}

savana.ee_init

Earth Engine session setup and flexible AOI (area-of-interest) loading.

The original GEE script hardcoded a single park asset (projects/ee-desmond/assets/NewParkMerged) and filtered it by name. This module generalises that so any user can classify their study area, supplied as an Earth Engine asset ID, a GeoJSON/Shapefile path, a geopandas.GeoDataFrame, or an ee.Geometry/ee.FeatureCollection directly.

initialize(project=None, force=False)

Initialize the Earth Engine Python API (auth if needed).

Parameters:

Name	Type	Description	Default
project	str | None	Google Cloud project registered for Earth Engine use. If omitted, uses whatever is already configured for the environment (EARTHENGINE_PROJECT env var or prior ee.Authenticate() state).	None
force	bool	Re-initialize even if already initialized this session.	False

Source code in savana/ee_init.py

def initialize(project: str | None = None, force: bool = False) -> None:
    """Initialize the Earth Engine Python API (auth if needed).

    Args:
        project: Google Cloud project registered for Earth Engine use.
            If omitted, uses whatever is already configured for the
            environment (``EARTHENGINE_PROJECT`` env var or prior
            ``ee.Authenticate()`` state).
        force: Re-initialize even if already initialized this session.
    """
    global _EE_INITIALIZED
    import ee

    if _EE_INITIALIZED and not force:
        return

    project = project or os.environ.get("EARTHENGINE_PROJECT")
    try:
        if project:
            ee.Initialize(project=project)
        else:
            ee.Initialize()
    except Exception:
        ee.Authenticate()
        if project:
            ee.Initialize(project=project)
        else:
            ee.Initialize()
    _EE_INITIALIZED = True

load_aoi(source, name_filter=None)

Resolve any of several AOI input types into a single ee.Geometry.

Parameters:

Name	Type	Description	Default
source	Any	One of: - ee.Geometry or ee.FeatureCollection (used directly) - str Earth Engine asset ID, e.g. "projects/x/assets/parks" - str path to a local GeoJSON / Shapefile / GeoPackage - geopandas.GeoDataFrame / geopandas.GeoSeries - dict GeoJSON geometry or feature	required
name_filter	str | None	If source resolves to a FeatureCollection with multiple features (e.g. a parks database), filter to the feature(s) whose NAME property equals this value before dissolving to a single geometry. If the asset uses a different property name, filter it yourself beforehand and pass an ee.Geometry instead.	None

Returns:

Type	Description
Geometry	ee.Geometry: a single (possibly multi-part) dissolved geometry.

Source code in savana/ee_init.py

def load_aoi(source: Any, name_filter: str | None = None) -> ee.Geometry:
    """Resolve any of several AOI input types into a single ``ee.Geometry``.

    Args:
        source: One of:
            - ``ee.Geometry`` or ``ee.FeatureCollection`` (used directly)
            - str Earth Engine asset ID, e.g. ``"projects/x/assets/parks"``
            - str path to a local GeoJSON / Shapefile / GeoPackage
            - ``geopandas.GeoDataFrame`` / ``geopandas.GeoSeries``
            - dict GeoJSON geometry or feature
        name_filter: If ``source`` resolves to a FeatureCollection with
            multiple features (e.g. a parks database), filter to the
            feature(s) whose ``NAME`` property equals this value before
            dissolving to a single geometry. If the asset uses a
            different property name, filter it yourself beforehand and
            pass an ``ee.Geometry`` instead.

    Returns:
        ee.Geometry: a single (possibly multi-part) dissolved geometry.
    """
    import ee

    initialize()

    if isinstance(source, ee.Geometry):
        return source

    if isinstance(source, ee.FeatureCollection):
        fc = source
        if name_filter is not None:
            fc = fc.filter(ee.Filter.eq("NAME", name_filter))
        return fc.geometry().dissolve(maxError=1)

    if isinstance(source, dict):
        geom = source.get("geometry", source)
        return ee.Geometry(geom)

    if isinstance(source, str):
        # Earth Engine asset IDs don't have file extensions and aren't
        # local paths; anything with a recognised geo file extension
        # (or that exists on disk) is treated as a local vector file.
        _, ext = os.path.splitext(source)
        if os.path.exists(source) or ext.lower() in (
            ".geojson",
            ".json",
            ".shp",
            ".gpkg",
        ):
            return _load_local_vector(source, name_filter)
        # Otherwise assume it's an EE asset ID.
        fc = ee.FeatureCollection(source)
        if name_filter is not None:
            fc = fc.filter(ee.Filter.eq("NAME", name_filter))
        return fc.geometry().dissolve(maxError=1)

    # geopandas GeoDataFrame / GeoSeries duck-typed via __geo_interface__
    if hasattr(source, "__geo_interface__"):
        return _geodataframe_to_ee_geometry(source)

    raise TypeError(
        f"Unsupported AOI source type: {type(source)!r}. Pass an "
        "ee.Geometry, ee.FeatureCollection, EE asset ID string, local "
        "vector file path, or a geopandas GeoDataFrame."
    )