API reference

An upstream river basin and everything open that can be clipped to it.

Create one from an outlet coordinate and every layer method afterwards is masked to the polygon, not to its bounding box::

import basinkit as bk

basin = bk.Basin.from_point(26.87, 87.15)    # Sapta Koshi at Chatara
basin.area_km2
dem = basin.dem()                            # xarray, clipped + masked
lc  = basin.landcover()
rain = basin.precipitation(2000, 2024)       # basin-mean monthly series
basin.download_all("koshi/")                 # the whole default stack

Attributes:

Name	Type	Description
geometry	shapely geometry	Basin polygon in EPSG:4326.
provenance	dict	Which backend and which dataset version produced the polygon. This travels with the basin so that a result is always attributable, and it is written into every export.

Source code in basinkit/basin.py

class Basin:
    """An upstream river basin and everything open that can be clipped to it.

    Create one from an outlet coordinate and every layer method afterwards is
    masked to the polygon, not to its bounding box::

        import basinkit as bk

        basin = bk.Basin.from_point(26.87, 87.15)    # Sapta Koshi at Chatara
        basin.area_km2
        dem = basin.dem()                            # xarray, clipped + masked
        lc  = basin.landcover()
        rain = basin.precipitation(2000, 2024)       # basin-mean monthly series
        basin.download_all("koshi/")                 # the whole default stack

    Attributes
    ----------
    geometry : shapely geometry
        Basin polygon in EPSG:4326.
    provenance : dict
        Which backend and which dataset version produced the polygon. This
        travels with the basin so that a result is always attributable, and it
        is written into every export.
    """

    def __init__(self, geometry, provenance: dict | None = None) -> None:
        self.geometry = geometry
        self.provenance = provenance or {}
        #: How far into the cache ledger this basin has already attributed.
        self._cache_mark = 0
        self._cache: dict[str, Any] = {}

    # -- constructors ------------------------------------------------------
    @classmethod
    def from_point(
        cls, lat: float, lon: float, *, backend: str = "auto", **kwargs
    ) -> Basin:
        """Delineate the basin upstream of an outlet coordinate.

        Parameters
        ----------
        backend
            ``'auto'`` (default), ``'hydrobasins'``, ``'dem'`` or ``'api'``.
            See :mod:`basinkit.delineate` for what each one is good at.
        """
        from .delineate import delineate

        # Coerced before it is compared, because a string reached the
        # comparison and raised "'<=' not supported between str and int" from
        # inside this check -- a message about Python's type system in place
        # of the one written here about coordinates.
        try:
            lat = float(lat)
            lon = float(lon)
        except (TypeError, ValueError):
            raise ValueError(
                f"({lat!r}, {lon!r}) is not a coordinate pair. Pass numbers: "
                "from_point(26.87, 87.15)."
            ) from None
        if not (lat == lat and lon == lon) or lat in (float("inf"), float("-inf")) \
                or lon in (float("inf"), float("-inf")):
            raise ValueError(
                f"({lat}, {lon}) is not a coordinate pair: a latitude or "
                "longitude cannot be NaN or infinite."
            )
        if not -90 <= lat <= 90 or not -180 <= lon <= 180:
            raise ValueError(
                f"({lat}, {lon}) is not a valid lat/lon. Note the order is "
                "(lat, lon) -- swapping them is the usual cause."
            )
        geom, prov = delineate(lat, lon, backend=backend, **kwargs)
        return cls(geom, prov)

    @classmethod
    def from_geometry(cls, geometry, provenance: dict | None = None,
                      crs: str | None = None) -> Basin:
        """Wrap a polygon you already have (a gauge basin, an official boundary).

        ``crs`` is the CRS the geometry arrives in. Anything other than
        EPSG:4326 is reprojected, because every step downstream of here works
        in lon/lat. Leaving it out keeps the old behaviour and takes the
        geometry as EPSG:4326 already, so existing callers are unaffected.
        """
        if crs is not None:
            import geopandas as gpd

            gs = gpd.GeoSeries([geometry], crs=crs)
            if gs.crs is not None and gs.crs.to_epsg() != 4326:
                geometry = gs.to_crs("EPSG:4326").iloc[0]
        return cls(geometry, provenance or {"backend": "user-supplied"})

    @classmethod
    def from_file(cls, path: str | Path) -> Basin:
        """Load a basin from any vector file geopandas can read."""
        import geopandas as gpd

        gdf = gpd.read_file(path)
        if gdf.crs and gdf.crs.to_epsg() != 4326:
            gdf = gdf.to_crs("EPSG:4326")
        return cls(gdf.union_all(), {"backend": "file", "path": str(path)})

    # -- properties --------------------------------------------------------
    @property
    def area_km2(self) -> float:
        """Basin area via an equal-area projection centred on the basin itself."""
        if "area" not in self._cache:
            from .clip import basin_area_km2

            self._cache["area"] = basin_area_km2(self.geometry)
        return self._cache["area"]

    @property
    def bounds(self) -> tuple[float, float, float, float]:
        return self.geometry.bounds

    @property
    def centroid(self) -> tuple[float, float]:
        c = self.geometry.centroid
        return (c.y, c.x)

    @property
    def bbox_efficiency(self) -> float:
        """Basin area as a fraction of its bounding-box area.

        This is the number that justifies polygon clipping. A compact basin
        scores near 0.7; a long dendritic one can drop below 0.25, meaning a
        bbox-based download wastes three quarters of everything it transfers
        and biases every basin average with a neighbour's pixels.
        """
        from shapely.geometry import box

        from .clip import basin_area_km2

        return self.area_km2 / basin_area_km2(box(*self.geometry.bounds))

    def __repr__(self) -> str:
        backend = self.provenance.get("backend", "?")
        lat, lon = self.centroid
        return (
            f"<Basin area={self.area_km2:,.0f} km2 "
            f"centroid=({lat:.3f}, {lon:.3f}) backend={backend!r}>"
        )

    # -- what was fetched, and when ---------------------------------------
    #: Keys worth keeping out of a layer's attrs. The sources already tag
    #: every array they return; what was missing was a place to keep the tags
    #: once the array had been handed on and forgotten.
    _ACQUISITION_KEYS = ("basinkit_product", "license", "citation",
                         "basinkit_native_res_m", "basinkit_output_res_m",
                         "basinkit_resampling", "basinkit_coarsen_factor",
                         "basinkit_tiles_merged", "basinkit_sources")

    def _acquired(self, out, call: str, **extra):
        """Record what a layer call fetched, then hand the layer back unchanged.

        A product name carries an epoch -- "ESA WorldCover 2021" -- but not a
        retrieval date, and the two answer different questions. The epoch says
        which edition was asked for; the date says which bytes arrived, which
        is what tells you whether a rerun months later differs because you
        changed something or because the provider reissued a tile.

        Appended rather than replaced: fetching the same layer twice in a run
        is a fact about the run, not something to tidy away.
        """
        import datetime

        attrs = getattr(out, "attrs", None) or {}
        record = {"call": call,
                  "retrieved_utc": datetime.datetime.now(
                      datetime.timezone.utc).isoformat(timespec="seconds")}
        for key in self._ACQUISITION_KEYS:
            if key in attrs:
                record[key.replace("basinkit_", "")] = attrs[key]
        # The attrs say what arrived; the extras say what was asked for. Where
        # both name the same thing the arrival wins, or a report prints the key
        # a caller typed -- "cop30" -- in place of the product that answered
        # it, "Copernicus DEM GLO-30".
        for key, value in extra.items():
            if value is not None:
                record.setdefault(key, value)

        # Which files this call actually read, and what was in them. The
        # product name and the date are what the run says about itself; the
        # digest is what a later run can check. Without it "the inputs had not
        # moved" is an assertion, and the record exists so that it need not be.
        from . import cache

        seen = cache.ledger()
        if len(seen) > self._cache_mark:
            # One call that opened the same file twice read the same bytes
            # twice, which says nothing a reader needs, so each file appears
            # once per layer. Across layers they repeat, because a second
            # layer reading the same tile is a fact about the run.
            fresh, sources = seen[self._cache_mark:], []
            for s in fresh:
                entry = {k: s[k] for k in ("url", "bytes", "sha256", "from_cache")}
                if entry not in sources:
                    sources.append(entry)
            record["sources"] = sources
            self._cache_mark = len(seen)

        self.provenance.setdefault("layers", []).append(record)
        return out

    def acquisitions(self):
        """Every layer fetched on this basin, in the order it was fetched.

        This is the row set a batch run should keep beside each basin: rerun
        the same outlet in six months and a diff over these records says
        whether the inputs moved.
        """
        return list(self.provenance.get("layers", []))

    # -- layers ------------------------------------------------------------
    def dem(self, product: str = "cop30", **kwargs):
        """Elevation, clipped and masked to the basin."""
        from .sources.dem import dem

        out = dem(self.geometry, product=product, **kwargs)
        # Some GDAL/PROJ combinations -- QGIS's bundled stack among them --
        # open a GeoTIFF and hand back no CRS at all. Everything downstream
        # then fails with rioxarray's "CRS not found", naming the symptom
        # rather than the cause. The source is lon/lat, so say so.
        if getattr(getattr(out, "rio", None), "crs", None) is None:
            out = out.rio.write_crs("EPSG:4326")
        return self._acquired(out, "dem", requested=product)

    def landcover(self, year: int | None = None, source: str = "worldcover",
                  **kwargs):
        """Land cover, as a 2-D array with its class legend in ``.attrs``.

        ``worldcover`` is ESA WorldCover at 10 m (2020 or 2021); ``esri`` is the
        ESRI / Impact Observatory annual series. Both return the same shape of
        object, and each carries its own legend -- the two number their classes
        differently, and code 10 is tree cover in one and cloud in the other.

        ``year=None`` means 2021 for WorldCover and the latest year published
        for this location for ESRI.
        """
        from .sources.landcover import esri_lulc, worldcover

        if source == "worldcover":
            return self._acquired(
                worldcover(self.geometry, year=year or 2021, **kwargs), "landcover")
        if source == "esri":
            return self._acquired(
                esri_lulc(self.geometry, year=year, **kwargs), "landcover")
        raise ValueError(f"Unknown land cover source {source!r}: use 'worldcover' or 'esri'.")

    def soil(self, prop: str = "clay", depth: str = "0-5cm", **kwargs):
        """A SoilGrids property. See :data:`basinkit.sources.soil.PROPERTIES`."""
        from .sources.soil import soilgrids

        return self._acquired(
            soilgrids(self.geometry, prop=prop, depth=depth, **kwargs),
            "soil", property=prop, depth=depth)

    def available_water_capacity(self, depth: str = "0-5cm"):
        """Plant-available water capacity (field capacity minus wilting point)."""
        from .sources.soil import available_water_capacity

        return available_water_capacity(self.geometry, depth=depth)

    def precipitation(self, start=2000, end=None, source: str = "chirps", **kwargs):
        """Basin-mean rainfall time series. ``chirps``, ``persiann`` or ``terraclimate``.

        Every source returns the same kind of object: a named
        ``xarray.DataArray`` of precipitation indexed by time. It did not used
        to -- ``terraclimate`` handed back a whole ``Dataset``, so code that
        worked on one source raised an attribute error on another, and the
        function's own name promised one quantity while returning thirteen.
        For the rest of the TerraClimate variables call
        :meth:`water_balance`, or ``sources.climate.terraclimate`` directly.
        """
        from .sources.climate import chirps, persiann, terraclimate

        if source == "chirps":
            return self._acquired(
                chirps(self.geometry, start, end, **kwargs),
                "precipitation", source="chirps", years=f"{start}-{end or 'latest'}")
        if source == "persiann":
            return self._acquired(
                persiann(self.geometry, str(start), end, **kwargs),
                "precipitation", source="persiann", years=f"{start}-{end or 'latest'}")
        if source == "terraclimate":
            ds = terraclimate(self.geometry, ("ppt",), int(start), end, **kwargs)
            out = ds["ppt"]
            out.name = "precipitation"
            out.attrs.setdefault("units", "mm/month")
            # The dataset-level provenance -- which years are missing, whether
            # the clip held -- belongs to this series too, and is the only
            # record that it is incomplete.
            for key, value in ds.attrs.items():
                out.attrs.setdefault(key, value)
            return self._acquired(out, "precipitation",
                                  source="terraclimate",
                                  years=f"{start}-{end or 'latest'}")
        raise ValueError(
            f"Unknown precipitation source {source!r}: use 'chirps', 'persiann' "
            "or 'terraclimate'."
        )

    def water_balance(self, start: int = 2000, end: int | None = None):
        """Monthly P / AET / PET / Q / soil-moisture balance from TerraClimate."""
        from .sources.climate import water_balance

        return water_balance(self.geometry, start, end)

    def surface_water(self, layer: str = "occurrence", **kwargs):
        """JRC Global Surface Water: a pre-reduced 37-year Landsat water stack."""
        from .sources.water import global_surface_water

        return self._acquired(
            global_surface_water(self.geometry, layer=layer, **kwargs),
            "surface_water", layer=layer)

    def attributes(self, prefixes: tuple[str, ...] | None = None, **kwargs):
        """281 pre-computed BasinATLAS attributes for this basin.

        The row returned belongs to the outlet's HydroBASINS unit, and its
        ``_u`` columns are already aggregated over everything upstream -- so
        this characterises the whole catchment without touching a raster.

        Costs one 2.7 GB download the first time, then nothing.
        """
        from .sources.attributes import describe, hydroatlas

        hybas_id = self.provenance.get("outlet_hybas_id")
        looked_up = None
        if hybas_id is None:
            # Refusing here was the wrong answer three times over. The
            # HydroBASINS unit containing an outlet can be looked up for any
            # basin, whichever backend drew the boundary; this call already
            # costs a 2.7 GB download, so the regional file it needs is not
            # the expensive part; and the QGIS algorithm for the same thing
            # delineates with backend='hydrobasins' itself, so the plugin
            # answered a question the library declined. The id is looked up
            # and recorded as looked up, which is the part that matters: the
            # row describes the unit at the outlet, not the polygon in hand.
            from .delineate.hydrobasins import reported_upland_km2

            outlet = self.provenance.get("outlet")
            if outlet is not None:
                lat, lon = float(outlet[0]), float(outlet[1])
            else:
                lat, lon = self.centroid
            found = reported_upland_km2(lat, lon, allow_download=True,
                                        progress=kwargs.get("progress", True))
            if not found or not found.get("hybas_id"):
                raise ValueError(
                    "BasinATLAS is keyed by HydroBASINS id, and no "
                    f"HydroBASINS unit was found at ({lat:.5f}, {lon:.5f}) to "
                    "key this basin by. Pass a geometry to "
                    "basinkit.sources.attributes.hydroatlas() directly."
                )
            hybas_id = found["hybas_id"]
            looked_up = found

        gdf = hydroatlas(hybas_id=hybas_id, prefixes=prefixes, **kwargs)
        out = describe(gdf.iloc[0])
        if looked_up is not None and isinstance(out, dict):
            out["_basinkit_note"] = (
                f"This basin came from the "
                f"{self.provenance.get('backend', 'unknown')!r} backend, which "
                "does not carry a HydroBASINS id, so unit "
                f"{hybas_id} was looked up as the one containing the outlet. "
                "Its upstream columns describe what HydroBASINS says drains to "
                f"that unit's outlet -- {looked_up['up_area_km2']:,.0f} km2 -- "
                f"against the {self.area_km2:,.1f} km2 of the polygon in hand. "
                "Where those two differ, the attributes describe the former."
            )
        return out

    def rivers(self, min_order: int = 0, **kwargs):
        """HydroRIVERS reaches inside the basin, with discharge and stream order."""
        from .sources.vectors import hydrorivers

        return self._acquired(
            hydrorivers(self.geometry, min_order=min_order, **kwargs), "rivers")

    def lakes(self, min_area_km2: float = 0.0, **kwargs):
        """HydroLAKES water bodies inside the basin."""
        from .sources.vectors import hydrolakes

        return self._acquired(
            hydrolakes(self.geometry, min_area_km2=min_area_km2, **kwargs), "lakes")

    def sentinel2(self, start: str, end: str, *, cloud_cover: float = 20,
                  bands: list[str] | None = None, composite: str | None = "median",
                  **kwargs):
        """Sentinel-2 L2A over the basin, cloud-filtered and optionally composited."""
        from .sources.stac import composite as reduce_time
        from .sources.stac import stac_search, stac_stack

        stack_kw = _stack_kwargs(kwargs)
        items = stac_search(
            "sentinel2", geometry=self.geometry, start=start, end=end,
            cloud_cover=cloud_cover, **kwargs
        )
        ds = stac_stack(items, self.geometry, bands=bands or ["blue", "green", "red", "nir"],
                        collection="sentinel2", **stack_kw)
        return reduce_time(ds, composite) if composite else ds

    def landsat(self, start: str, end: str, *, cloud_cover: float = 20,
                bands: list[str] | None = None, composite: str | None = "median",
                **kwargs):
        """Landsat Collection 2 Level-2 over the basin (1982 to present)."""
        from .sources.stac import composite as reduce_time
        from .sources.stac import stac_search, stac_stack

        stack_kw = _stack_kwargs(kwargs)
        items = stac_search(
            "landsat", geometry=self.geometry, start=start, end=end,
            cloud_cover=cloud_cover, **kwargs
        )
        ds = stac_stack(items, self.geometry, bands=bands or ["blue", "green", "red", "nir08"],
                        collection="landsat", **stack_kw)
        return reduce_time(ds, composite) if composite else ds

    def sentinel1(self, start: str, end: str, *, bands: list[str] | None = None,
                  composite: str | None = "median", **kwargs):
        """Sentinel-1 RTC: terrain-corrected radar, so it works through cloud."""
        from .sources.stac import composite as reduce_time
        from .sources.stac import stac_search, stac_stack

        stack_kw = _stack_kwargs(kwargs)
        items = stac_search(
            "sentinel1_rtc", geometry=self.geometry, start=start, end=end, **kwargs
        )
        ds = stac_stack(items, self.geometry, bands=bands or ["vv", "vh"],
                        collection="sentinel1_rtc", **stack_kw)
        return reduce_time(ds, composite) if composite else ds

    # -- summaries ---------------------------------------------------------
    # -- terrain surfaces -------------------------------------------------
    #
    # Each of these is computed from the elevation array, so passing ``dem=``
    # lets one download serve all of them.

    def slope(self, *, dem=None, degrees: bool = True, **kwargs):
        """Steepest descent at every cell, in degrees by default."""
        from .terrain import slope

        return slope(self.dem(**kwargs) if dem is None else dem, degrees=degrees)

    def aspect(self, *, dem=None, **kwargs):
        """Which way each slope faces, clockwise from north. Flat is NaN."""
        from .terrain import aspect

        return aspect(self.dem(**kwargs) if dem is None else dem)

    def hillshade(self, *, dem=None, azimuth: float = 315.0,
                  altitude: float = 45.0, z_factor: float = 1.0, **kwargs):
        """Shaded relief, 0 to 1, for a figure that reads as terrain."""
        from .terrain import hillshade

        return hillshade(self.dem(**kwargs) if dem is None else dem,
                         azimuth=azimuth, altitude=altitude, z_factor=z_factor)

    def curvature(self, *, dem=None, **kwargs):
        """Profile curvature: positive where the ground sheds, negative where
        it collects."""
        from .terrain import curvature

        return curvature(self.dem(**kwargs) if dem is None else dem)

    def flow_accumulation(self, *, dem=None, **kwargs):
        """Cells draining through each cell, routed over this basin's own DEM."""
        from .terrain import flow_accumulation

        return flow_accumulation(self.dem(**kwargs) if dem is None else dem)

    def twi(self, *, dem=None, **kwargs):
        """Topographic wetness index: where water gathers and ground saturates."""
        from .terrain import twi

        return twi(self.dem(**kwargs) if dem is None else dem)

    def hand(self, *, dem=None, min_area_km2: float = 1.0, **kwargs):
        """Height above the nearest drainage: the terrain layer flood work wants."""
        from .terrain import hand

        return hand(self.dem(**kwargs) if dem is None else dem,
                    min_area_km2=min_area_km2)

    def tpi(self, *, dem=None, window: int = 11, **kwargs):
        """Topographic position index: height above the surrounding ground."""
        from .terrain import tpi

        return tpi(self.dem(**kwargs) if dem is None else dem, window=window)

    def tri(self, *, dem=None, **kwargs):
        """Terrain ruggedness index: how far a cell sits from its neighbours."""
        from .terrain import tri

        return tri(self.dem(**kwargs) if dem is None else dem)

    def roughness(self, *, dem=None, **kwargs):
        """Local relief within each cell's 3x3 neighbourhood."""
        from .terrain import roughness

        return roughness(self.dem(**kwargs) if dem is None else dem)

    def landform(self, *, dem=None, window: int = 11, **kwargs):
        """Six slope-position classes: valley, slopes, flat, ridge."""
        from .terrain import landform

        return landform(self.dem(**kwargs) if dem is None else dem, window=window)

    def drainage_density(self, *, dem=None, thresholds=None,
                         chosen_km2=None, **kwargs):
        """Drainage density across the range of defensible thresholds.

        Not one number: the density moves by a large factor across channel
        initiation thresholds that are all defensible, so the curve and the
        threshold used are what make the figure comparable with anyone else's.
        """
        from .terrain import drainage_density

        return drainage_density(self.dem(**kwargs) if dem is None else dem,
                                thresholds=thresholds, chosen_km2=chosen_km2)

    def report(self, path, *, title=None, dem=None, rivers=None, **kwargs):
        """Write the eight-page PDF report for this basin.

        Cover with the elevation data suitability grade, elevation and
        hypsometry, slope and aspect, the shape of the ground, the channel
        network and Horton's laws, the morphometric table with symbols and
        references, what the answer rests on, and a methods page.

        ``max_pixels=`` sets the pixel budget for the elevation grid the
        report is built from; the cover states the resolution it worked out
        to. The default is :data:`basinkit.report.REPORT_MAX_PIXELS` and not
        the package-wide one, because every page derives arrays of its own
        from this grid.
        """
        from .report import REPORT_MAX_PIXELS, report

        # The budget has to be applied here, because this method fetched the
        # grid itself and handed report() a finished one -- so report()'s own
        # max_pixels was never reached, and every caller that goes through
        # this method, which is all of them including the QGIS and ArcGIS
        # algorithms, still read at the package-wide 100 Mpx default. That is
        # what reached 6.1 GB on the Koshi and was killed by the kernel. The
        # earlier fix was real and was in the wrong function.
        if dem is None:
            kwargs.setdefault("max_pixels", REPORT_MAX_PIXELS)
            dem = self.dem(**kwargs)
        return report(self, path, title=title, dem=dem, rivers=rivers)

    def data_quality(self, **kwargs):
        """A quality indicator for every layer, each against something independent.

        The elevation grade split by terrain class; two independently produced
        land cover maps compared; SoilGrids' own 90 per cent interval; CHIRPS
        against TerraClimate; permanent against seasonal water; and the checks
        that exist for the boundary. Each entry carries the measurement, the
        threshold it was judged against, and what was compared with what.
        """
        from .quality import data_quality

        return data_quality(self, **kwargs)

    def dem_suitability(self, *, dem=None, support_map: bool = False, **kwargs):
        """Whether the elevation model supports terrain analysis in this basin.

        Five measured tests -- depression filling, the slope noise floor,
        relief against the model's vertical error, level water surfaces and
        missing coverage -- combined into ``HIGH``, ``MODERATE`` or
        ``LIMITED``. It grades the terrain products, not the basin boundary.
        """
        from .suitability import suitability

        return suitability(self, dem=self.dem(**kwargs) if dem is None else dem,
                           support_map=support_map)

    def validation_envelope(self, *, suitability=None):
        """Which part of the validation set this basin belongs to.

        The accuracy figures published with basinkit are medians over 2,740
        delineations at 2,550 gauges, and they differ by two orders of
        magnitude between the largest and the smallest band. This reports the
        error that was measured at *this* basin's size, and names the
        conditions the validation never covered, so that a figure which does
        not apply here is not quoted as though it did.

        Reads only the area, so it downloads nothing. Pass a
        :meth:`dem_suitability` result to compare the terrain descriptors too.
        """
        from .envelope import validation_envelope

        return validation_envelope(self, suitability=suitability)

    def streams(self, *, dem=None, min_area_km2: float = 1.0, **kwargs):
        """The channel network routed from this basin's own elevation."""
        from .terrain import streams

        return streams(self.dem(**kwargs) if dem is None else dem,
                       min_area_km2=min_area_km2)

    def zonal(self, values, zones=None, *, bins=None, labels=None, **kwargs):
        """Summarise one layer inside the classes of another.

        ``basin.zonal(basin.precipitation_grid, basin.landcover())`` answers
        how much rain falls on each land-cover class; passing ``bins`` cuts a
        continuous layer such as elevation into bands instead. Area is summed
        from the true size of every cell, so a basin spanning several degrees
        of latitude is not weighted towards its southern edge.
        """
        from .zonal import zonal

        if zones is None:
            zones = self.landcover(**kwargs)
        return zonal(values, zones, labels=labels, bins=bins)

    def landcover_change(self, start: int, end: int, *, source: str = "esri",
                         **kwargs):
        """What became what, between two years, in square kilometres.

        Only the ESRI annual series carries a year for every year, so that is
        the default here even though ``landcover()`` defaults to WorldCover.
        The diagonal is the ground that did not change; everything off it is a
        transition, and the table is the honest form of a deforestation or
        urban-growth figure because it shows what the loss became.
        """
        import numpy as np
        import pandas as pd

        from .sources.landcover import _declared_classes
        from .terrain import cell_area_km2

        first = self.landcover(year=start, source=source, **kwargs)
        second = self.landcover(year=end, source=source, **kwargs)
        if first.shape != second.shape:
            second = second.rio.reproject_match(first)

        legend = _declared_classes(first) or {}
        a = np.asarray(first.values, dtype="float64")
        b = np.asarray(second.values, dtype="float64")
        areas = cell_area_km2(first)
        usable = np.isfinite(a) & np.isfinite(b) & (a > 0) & (b > 0)

        rows = []
        for code_from in np.unique(a[usable]):
            for code_to in np.unique(b[usable]):
                picked = usable & (a == code_from) & (b == code_to)
                if not picked.any():
                    continue
                rows.append({
                    "from": legend.get(int(code_from), int(code_from)),
                    "to": legend.get(int(code_to), int(code_to)),
                    "area_km2": round(float(areas[picked].sum()), 3),
                    "changed": bool(code_from != code_to),
                })
        frame = pd.DataFrame(rows).sort_values("area_km2", ascending=False)
        frame.attrs["years"] = (start, end)
        return frame.reset_index(drop=True)

    def precipitation_trend(self, *, start: int = 1985, end: int | None = None,
                            alpha: float = 0.05, **kwargs):
        """Is rainfall over this basin trending, and by how much a year."""
        from .climate import trend

        series = self.precipitation(start=start, end=end, **kwargs)
        s = series.to_series() if hasattr(series, "to_series") else series
        yearly = s.groupby(s.index.year).sum()
        yearly.index = __import__("pandas").to_datetime(
            [f"{int(y)}-01-01" for y in yearly.index]
        )
        return trend(yearly, alpha=alpha)

    def spi(self, *, scale: int = 3, start: int = 1985, end: int | None = None,
            **kwargs):
        """Standardized Precipitation Index over this basin, month by month."""
        from .climate import spi

        series = self.precipitation(start=start, end=end, **kwargs)
        s = series.to_series() if hasattr(series, "to_series") else series
        return spi(s, scale=scale)

    def subbasins(self, **kwargs):
        """The sub-catchments this basin is assembled from, with their routing.

        ``delineate`` dissolves the HydroBASINS units it walked into a single
        polygon. This hands back the pieces instead, each with ``NEXT_DOWN``,
        which is the routing graph itself: every distributed model wants
        sub-catchments and the links between them, and nothing here has to be
        inferred from geometry afterwards.

        Available for the ``hydrobasins`` backend, which is the default.
        """
        from .delineate.hydrobasins import upstream_units

        prov = self.provenance
        if prov.get("backend") != "hydrobasins":
            raise ValueError(
                "subbasins() reads the units the HydroBASINS traversal walked, "
                f"and this basin came from the {prov.get('backend', 'unknown')!r} "
                "backend. Delineate with backend='hydrobasins' to get them."
            )
        return upstream_units(
            prov["region"], prov["outlet_hybas_id"], prov.get("level", 12), **kwargs
        )

    def curve_number(self, *, cn_table=None, **kwargs) -> dict:
        """The SCS curve number, composed from this basin's cover and soil.

            cn = basin.curve_number()
            cn["composite_cn"], cn["cn_amc_iii_wet"]

        Nothing new is downloaded: the land cover and the soil texture are the
        ones this package already fetches. Two judgements go into it -- the
        hydrologic soil group from texture rather than a soil survey, and ESA
        WorldCover classes matched to TR-55 cover types -- and both are
        returned in ``assumptions`` and ``cover_mapping`` so the number can be
        reproduced or disagreed with.
        """
        from .runoff import curve_number

        return curve_number(self, cn_table=cn_table, **kwargs)

    def concentration_time(self, *, morphometry=None, **kwargs) -> dict:
        """How long the basin takes to respond to rain, by several formulas.

            t = basin.concentration_time()
            t["estimates"]["kirpich_min"]
            t["spread_factor"]

        Time of concentration is not measured; every value is an empirical
        formula fitted to a particular landscape, and on one basin they
        routinely disagree by a factor of two. All of them are returned, with
        what each was fitted to, so the choice stays with the person making
        it. Pass ``morphometry=`` the dict you already have to avoid
        recomputing it.
        """
        from .response import concentration_time

        m = morphometry if morphometry is not None else self.morphometry(**kwargs)
        return concentration_time(m)

    def terrain_stats(self, *, dem=None, **kwargs) -> dict:
        """Elevation, relief and mean slope: the standard morphometry.

        Pass ``dem=`` the elevation you already hold. Fetching a second one
        here is not only wasteful: if the first was coarsened to a pixel
        budget and this one is not, the two disagree about the highest point
        of the same basin, and a report that prints both looks wrong whichever
        number the reader believes.
        """
        import numpy as np

        elev = dem if dem is not None else self.dem(**kwargs)
        vals = np.asarray(elev.values, dtype="float64")
        vals = vals[np.isfinite(vals)]
        if vals.size == 0:
            return {}

        # Slope comes from the same routine the slope map uses, which keeps
        # nodata as nodata. Filling the space around the polygon with zeros
        # before differencing puts a plain around the basin and a cliff at its
        # edge, and the mean then answers to those instead: on the Helmand
        # above the Kajaki Dam, which fills 43% of its bounding box, it read 7
        # degrees where the basin's own slopes average 16. That is the
        # bounding-box mistake this package exists to avoid, made against
        # itself.
        from .terrain import slope as _slope_of

        slope = np.asarray(_slope_of(elev).values, dtype="float64")
        finite_slope = slope[np.isfinite(slope)]

        return {
            "area_km2": round(self.area_km2, 2),
            "elev_min_m": round(float(vals.min()), 1),
            "elev_max_m": round(float(vals.max()), 1),
            "elev_mean_m": round(float(vals.mean()), 1),
            "relief_m": round(float(vals.max() - vals.min()), 1),
            "slope_mean_deg": (round(float(finite_slope.mean()), 2)
                               if finite_slope.size else None),
            "bbox_efficiency": round(self.bbox_efficiency, 3),
        }

    def morphometry(self, **kwargs) -> dict:
        """The classical Horton-Strahler-Schumm morphometric parameters.

            m = basin.morphometry()
            m["areal"]["drainage_density_km_per_km2"]
            m["network"]        # one row per Strahler order

        Counts Strahler *streams*, not the reaches a river dataset splits them
        into, and measures area, perimeter and every length in one equal-area
        projection. Both matter: on the Koshi, counting reaches turns the
        bifurcation ratios into values that are not physically possible.

        See :mod:`basinkit.morphometry` for what each parameter is and for what
        the numbers can and cannot be compared against.
        """
        from .morphometry import morphometry

        return morphometry(self, **kwargs)

    def summary(self, *, terrain: bool = True, landcover: bool = True) -> dict:
        """A one-call characterisation of the basin."""
        out: dict[str, Any] = {
            "area_km2": round(self.area_km2, 2),
            "centroid_lat_lon": [round(v, 5) for v in self.centroid],
            "bounds": [round(v, 5) for v in self.bounds],
            "bbox_efficiency": round(self.bbox_efficiency, 3),
            "provenance": self.provenance,
        }
        if terrain:
            try:
                out["terrain"] = self.terrain_stats()
            except Exception as exc:
                out["terrain"] = {"error": str(exc)}
        if landcover:
            try:
                from .sources.landcover import class_fractions

                out["landcover_fractions"] = class_fractions(self.landcover())
            except Exception as exc:
                out["landcover_fractions"] = {"error": str(exc)}
        return out

    # -- licensing ---------------------------------------------------------
    def license_report(self, layers: tuple[str, ...] | None = None) -> str:
        """Attribution and licence text for the layers you used.

        Print this into your methods section. Every layer basinkit fetches by
        default is CC BY 4.0 or more permissive, which means it can be
        redistributed and used commercially -- but only if it is attributed.
        """
        keys = layers or catalog.DEFAULT_STACK
        lines = ["Data sources and licences", "=" * 26, ""]
        for key in keys:
            try:
                ds = catalog.get(key)
            except KeyError:
                continue
            lines.append(f"{ds.name}")
            lines.append(f"    Licence : {ds.license}")
            lines.append(f"    Access  : {ds.route}")
            if ds.citation:
                lines.append(f"    Cite    : {ds.citation}")
            if not ds.commercial_ok:
                lines.append("    WARNING : commercial use not permitted")
            if not ds.redistributable:
                lines.append("    WARNING : redistribution not permitted")
            lines.append("")
        return "\n".join(lines)

    @staticmethod
    def check_license(key: str, *, commercial: bool = False,
                      redistribute: bool = False) -> None:
        """Raise if a dataset's licence forbids the intended use."""
        ds = catalog.get(key)
        if commercial and not ds.commercial_ok:
            raise LicenseError(
                f"{ds.name} is licensed {ds.license}, which forbids commercial use. "
                f"{ds.notes}"
            )
        if redistribute and not ds.redistributable:
            raise LicenseError(
                f"{ds.name} may not be redistributed under {ds.license}. {ds.notes}"
            )

    # -- export ------------------------------------------------------------
    def to_geojson(self, path: str | Path | None = None) -> str:
        import geopandas as gpd

        gdf = gpd.GeoDataFrame(
            {"area_km2": [self.area_km2],
             "backend": [self.provenance.get("backend", "")],
             "source": [self.provenance.get("source_dataset", "")]},
            geometry=[self.geometry], crs="EPSG:4326",
        )
        if path:
            gdf.to_file(path, driver="GeoJSON")
            return str(path)
        return gdf.to_json()

    def download_all(
        self,
        outdir: str | Path,
        layers: tuple[str, ...] = ("dem", "landcover", "soil", "surface_water",
                                   "precipitation", "rivers"),
        *,
        start: int = 2000,
        end: int | None = None,
        progress: bool = True,
    ) -> dict:
        """Fetch the default stack and write it to ``outdir``.

        This is the "give me everything" button. Each layer is attempted
        independently, so one failing source (a polar basin with no CHIRPS, say)
        does not abort the rest -- failures are recorded in the manifest
        alongside the successes.
        """
        outdir = Path(outdir)
        outdir.mkdir(parents=True, exist_ok=True)

        manifest: dict[str, Any] = {
            "basin": {
                "area_km2": round(self.area_km2, 2),
                "bounds": list(self.bounds),
                "centroid_lat_lon": list(self.centroid),
            },
            "provenance": self.provenance,
            "layers": {},
            "failed": {},
        }

        self.to_geojson(outdir / "basin.geojson")
        manifest["layers"]["basin"] = "basin.geojson"

        def _write_raster(da, name: str) -> str:
            fp = outdir / f"{name}.tif"
            # Declare nodata in the file header. The array is already masked
            # outside the basin, but without a declared nodata value QGIS and
            # ArcGIS paint that area solid black instead of transparent, and
            # rasterio's masked read returns no mask at all -- so a correctly
            # clipped raster looks and behaves like an unclipped one the moment
            # it leaves Python.
            import numpy as np

            if da.rio.nodata is None:
                if np.issubdtype(da.dtype, np.floating):
                    da = da.rio.write_nodata(np.nan, encoded=False)
                else:
                    da = da.rio.write_nodata(0, encoded=False)
            da.rio.to_raster(fp, compress="deflate", tiled=True)
            return fp.name

        jobs = {
            "dem": lambda: _write_raster(self.dem(progress=progress), "dem"),
            "landcover": lambda: _write_raster(
                self.landcover(progress=progress), "landcover"
            ),
            "soil": lambda: _write_raster(self.soil("clay"), "soil_clay_0-5cm"),
            "surface_water": lambda: _write_raster(
                self.surface_water(progress=progress), "surface_water_occurrence"
            ),
            "precipitation": lambda: self._write_series(
                self.precipitation(start, end), outdir, "precipitation_chirps"
            ),
            "rivers": lambda: self._write_vector(
                self.rivers(progress=progress), outdir, "rivers"
            ),
            "lakes": lambda: self._write_vector(
                self.lakes(progress=progress), outdir, "lakes"
            ),
        }

        for name in layers:
            if name not in jobs:
                manifest["failed"][name] = f"unknown layer {name!r}"
                continue
            try:
                manifest["layers"][name] = jobs[name]()
            except Exception as exc:
                manifest["failed"][name] = f"{type(exc).__name__}: {exc}"

        # Explicit UTF-8: the licence text and the manifest both carry
        # non-ASCII, and Windows would otherwise write them as cp1252.
        (outdir / "LICENSES.txt").write_text(self.license_report(),
                                             encoding="utf-8")
        manifest["layers"]["licenses"] = "LICENSES.txt"
        (outdir / "manifest.json").write_text(
            json.dumps(manifest, indent=2, default=str), encoding="utf-8")
        return manifest

    @staticmethod
    def _write_series(da, outdir: Path, name: str) -> str:
        fp = outdir / f"{name}.csv"
        da.to_dataframe().to_csv(fp)
        return fp.name

    @staticmethod
    def _write_vector(gdf, outdir: Path, name: str) -> str:
        fp = outdir / f"{name}.gpkg"
        if len(gdf) == 0:
            return f"{name}: none within basin"
        gdf.to_file(fp, driver="GPKG")
        return fp.name

    # -- viz ---------------------------------------------------------------
    def export_3d(self, path: str | Path, **kwargs):
        """Write an interactive 3D page for this basin: terrain, imagery, rivers.

        One self-contained HTML file with everything embedded, so it opens with
        no network. See :func:`basinkit.viz3d.export_3d` for the options.

            basin.export_3d("koshi.html")
            basin.export_3d("koshi.html", texture=None)   # elevation only, small

        This is a way of looking at the layers this package fetches. It makes no
        claim the other methods do not already make.
        """
        from .viz3d import export_3d

        return export_3d(self, path, **kwargs)

    def explore(self, **kwargs):
        """Interactive map of the basin. Needs ``pip install 'basinkit[viz]'``."""
        from .viz import explore

        return explore(self, **kwargs)

    def plot(self, **kwargs):
        """Static matplotlib figure: hypsometry, boundary and river network."""
        from .viz import plot

        return plot(self, **kwargs)

area_km2 property

Basin area via an equal-area projection centred on the basin itself.

bbox_efficiency property

Basin area as a fraction of its bounding-box area.

This is the number that justifies polygon clipping. A compact basin scores near 0.7; a long dendritic one can drop below 0.25, meaning a bbox-based download wastes three quarters of everything it transfers and biases every basin average with a neighbour's pixels.

acquisitions()

Every layer fetched on this basin, in the order it was fetched.

This is the row set a batch run should keep beside each basin: rerun the same outlet in six months and a diff over these records says whether the inputs moved.

Source code in basinkit/basin.py

def acquisitions(self):
    """Every layer fetched on this basin, in the order it was fetched.

    This is the row set a batch run should keep beside each basin: rerun
    the same outlet in six months and a diff over these records says
    whether the inputs moved.
    """
    return list(self.provenance.get("layers", []))

aspect(*, dem=None, **kwargs)

Which way each slope faces, clockwise from north. Flat is NaN.

Source code in basinkit/basin.py

def aspect(self, *, dem=None, **kwargs):
    """Which way each slope faces, clockwise from north. Flat is NaN."""
    from .terrain import aspect

    return aspect(self.dem(**kwargs) if dem is None else dem)

attributes(prefixes=None, **kwargs)

281 pre-computed BasinATLAS attributes for this basin.

The row returned belongs to the outlet's HydroBASINS unit, and its _u columns are already aggregated over everything upstream -- so this characterises the whole catchment without touching a raster.

Costs one 2.7 GB download the first time, then nothing.

Source code in basinkit/basin.py

def attributes(self, prefixes: tuple[str, ...] | None = None, **kwargs):
    """281 pre-computed BasinATLAS attributes for this basin.

    The row returned belongs to the outlet's HydroBASINS unit, and its
    ``_u`` columns are already aggregated over everything upstream -- so
    this characterises the whole catchment without touching a raster.

    Costs one 2.7 GB download the first time, then nothing.
    """
    from .sources.attributes import describe, hydroatlas

    hybas_id = self.provenance.get("outlet_hybas_id")
    looked_up = None
    if hybas_id is None:
        # Refusing here was the wrong answer three times over. The
        # HydroBASINS unit containing an outlet can be looked up for any
        # basin, whichever backend drew the boundary; this call already
        # costs a 2.7 GB download, so the regional file it needs is not
        # the expensive part; and the QGIS algorithm for the same thing
        # delineates with backend='hydrobasins' itself, so the plugin
        # answered a question the library declined. The id is looked up
        # and recorded as looked up, which is the part that matters: the
        # row describes the unit at the outlet, not the polygon in hand.
        from .delineate.hydrobasins import reported_upland_km2

        outlet = self.provenance.get("outlet")
        if outlet is not None:
            lat, lon = float(outlet[0]), float(outlet[1])
        else:
            lat, lon = self.centroid
        found = reported_upland_km2(lat, lon, allow_download=True,
                                    progress=kwargs.get("progress", True))
        if not found or not found.get("hybas_id"):
            raise ValueError(
                "BasinATLAS is keyed by HydroBASINS id, and no "
                f"HydroBASINS unit was found at ({lat:.5f}, {lon:.5f}) to "
                "key this basin by. Pass a geometry to "
                "basinkit.sources.attributes.hydroatlas() directly."
            )
        hybas_id = found["hybas_id"]
        looked_up = found

    gdf = hydroatlas(hybas_id=hybas_id, prefixes=prefixes, **kwargs)
    out = describe(gdf.iloc[0])
    if looked_up is not None and isinstance(out, dict):
        out["_basinkit_note"] = (
            f"This basin came from the "
            f"{self.provenance.get('backend', 'unknown')!r} backend, which "
            "does not carry a HydroBASINS id, so unit "
            f"{hybas_id} was looked up as the one containing the outlet. "
            "Its upstream columns describe what HydroBASINS says drains to "
            f"that unit's outlet -- {looked_up['up_area_km2']:,.0f} km2 -- "
            f"against the {self.area_km2:,.1f} km2 of the polygon in hand. "
            "Where those two differ, the attributes describe the former."
        )
    return out

available_water_capacity(depth='0-5cm')

Plant-available water capacity (field capacity minus wilting point).

Source code in basinkit/basin.py

def available_water_capacity(self, depth: str = "0-5cm"):
    """Plant-available water capacity (field capacity minus wilting point)."""
    from .sources.soil import available_water_capacity

    return available_water_capacity(self.geometry, depth=depth)

check_license(key, *, commercial=False, redistribute=False) staticmethod

Raise if a dataset's licence forbids the intended use.

Source code in basinkit/basin.py

@staticmethod
def check_license(key: str, *, commercial: bool = False,
                  redistribute: bool = False) -> None:
    """Raise if a dataset's licence forbids the intended use."""
    ds = catalog.get(key)
    if commercial and not ds.commercial_ok:
        raise LicenseError(
            f"{ds.name} is licensed {ds.license}, which forbids commercial use. "
            f"{ds.notes}"
        )
    if redistribute and not ds.redistributable:
        raise LicenseError(
            f"{ds.name} may not be redistributed under {ds.license}. {ds.notes}"
        )

concentration_time(*, morphometry=None, **kwargs)

How long the basin takes to respond to rain, by several formulas.

t = basin.concentration_time()
t["estimates"]["kirpich_min"]
t["spread_factor"]

Time of concentration is not measured; every value is an empirical formula fitted to a particular landscape, and on one basin they routinely disagree by a factor of two. All of them are returned, with what each was fitted to, so the choice stays with the person making it. Pass morphometry= the dict you already have to avoid recomputing it.

Source code in basinkit/basin.py

def concentration_time(self, *, morphometry=None, **kwargs) -> dict:
    """How long the basin takes to respond to rain, by several formulas.

        t = basin.concentration_time()
        t["estimates"]["kirpich_min"]
        t["spread_factor"]

    Time of concentration is not measured; every value is an empirical
    formula fitted to a particular landscape, and on one basin they
    routinely disagree by a factor of two. All of them are returned, with
    what each was fitted to, so the choice stays with the person making
    it. Pass ``morphometry=`` the dict you already have to avoid
    recomputing it.
    """
    from .response import concentration_time

    m = morphometry if morphometry is not None else self.morphometry(**kwargs)
    return concentration_time(m)

curvature(*, dem=None, **kwargs)

Profile curvature: positive where the ground sheds, negative where it collects.

Source code in basinkit/basin.py

def curvature(self, *, dem=None, **kwargs):
    """Profile curvature: positive where the ground sheds, negative where
    it collects."""
    from .terrain import curvature

    return curvature(self.dem(**kwargs) if dem is None else dem)

curve_number(*, cn_table=None, **kwargs)

The SCS curve number, composed from this basin's cover and soil.

cn = basin.curve_number()
cn["composite_cn"], cn["cn_amc_iii_wet"]

Nothing new is downloaded: the land cover and the soil texture are the ones this package already fetches. Two judgements go into it -- the hydrologic soil group from texture rather than a soil survey, and ESA WorldCover classes matched to TR-55 cover types -- and both are returned in assumptions and cover_mapping so the number can be reproduced or disagreed with.

Source code in basinkit/basin.py

def curve_number(self, *, cn_table=None, **kwargs) -> dict:
    """The SCS curve number, composed from this basin's cover and soil.

        cn = basin.curve_number()
        cn["composite_cn"], cn["cn_amc_iii_wet"]

    Nothing new is downloaded: the land cover and the soil texture are the
    ones this package already fetches. Two judgements go into it -- the
    hydrologic soil group from texture rather than a soil survey, and ESA
    WorldCover classes matched to TR-55 cover types -- and both are
    returned in ``assumptions`` and ``cover_mapping`` so the number can be
    reproduced or disagreed with.
    """
    from .runoff import curve_number

    return curve_number(self, cn_table=cn_table, **kwargs)

data_quality(**kwargs)

A quality indicator for every layer, each against something independent.

The elevation grade split by terrain class; two independently produced land cover maps compared; SoilGrids' own 90 per cent interval; CHIRPS against TerraClimate; permanent against seasonal water; and the checks that exist for the boundary. Each entry carries the measurement, the threshold it was judged against, and what was compared with what.

Source code in basinkit/basin.py

def data_quality(self, **kwargs):
    """A quality indicator for every layer, each against something independent.

    The elevation grade split by terrain class; two independently produced
    land cover maps compared; SoilGrids' own 90 per cent interval; CHIRPS
    against TerraClimate; permanent against seasonal water; and the checks
    that exist for the boundary. Each entry carries the measurement, the
    threshold it was judged against, and what was compared with what.
    """
    from .quality import data_quality

    return data_quality(self, **kwargs)

dem(product='cop30', **kwargs)

Elevation, clipped and masked to the basin.

Source code in basinkit/basin.py

def dem(self, product: str = "cop30", **kwargs):
    """Elevation, clipped and masked to the basin."""
    from .sources.dem import dem

    out = dem(self.geometry, product=product, **kwargs)
    # Some GDAL/PROJ combinations -- QGIS's bundled stack among them --
    # open a GeoTIFF and hand back no CRS at all. Everything downstream
    # then fails with rioxarray's "CRS not found", naming the symptom
    # rather than the cause. The source is lon/lat, so say so.
    if getattr(getattr(out, "rio", None), "crs", None) is None:
        out = out.rio.write_crs("EPSG:4326")
    return self._acquired(out, "dem", requested=product)

dem_suitability(*, dem=None, support_map=False, **kwargs)

Whether the elevation model supports terrain analysis in this basin.

Five measured tests -- depression filling, the slope noise floor, relief against the model's vertical error, level water surfaces and missing coverage -- combined into HIGH, MODERATE or LIMITED. It grades the terrain products, not the basin boundary.

Source code in basinkit/basin.py

def dem_suitability(self, *, dem=None, support_map: bool = False, **kwargs):
    """Whether the elevation model supports terrain analysis in this basin.

    Five measured tests -- depression filling, the slope noise floor,
    relief against the model's vertical error, level water surfaces and
    missing coverage -- combined into ``HIGH``, ``MODERATE`` or
    ``LIMITED``. It grades the terrain products, not the basin boundary.
    """
    from .suitability import suitability

    return suitability(self, dem=self.dem(**kwargs) if dem is None else dem,
                       support_map=support_map)

download_all(outdir, layers=('dem', 'landcover', 'soil', 'surface_water', 'precipitation', 'rivers'), *, start=2000, end=None, progress=True)

Fetch the default stack and write it to outdir.

This is the "give me everything" button. Each layer is attempted independently, so one failing source (a polar basin with no CHIRPS, say) does not abort the rest -- failures are recorded in the manifest alongside the successes.

Source code in basinkit/basin.py

def download_all(
    self,
    outdir: str | Path,
    layers: tuple[str, ...] = ("dem", "landcover", "soil", "surface_water",
                               "precipitation", "rivers"),
    *,
    start: int = 2000,
    end: int | None = None,
    progress: bool = True,
) -> dict:
    """Fetch the default stack and write it to ``outdir``.

    This is the "give me everything" button. Each layer is attempted
    independently, so one failing source (a polar basin with no CHIRPS, say)
    does not abort the rest -- failures are recorded in the manifest
    alongside the successes.
    """
    outdir = Path(outdir)
    outdir.mkdir(parents=True, exist_ok=True)

    manifest: dict[str, Any] = {
        "basin": {
            "area_km2": round(self.area_km2, 2),
            "bounds": list(self.bounds),
            "centroid_lat_lon": list(self.centroid),
        },
        "provenance": self.provenance,
        "layers": {},
        "failed": {},
    }

    self.to_geojson(outdir / "basin.geojson")
    manifest["layers"]["basin"] = "basin.geojson"

    def _write_raster(da, name: str) -> str:
        fp = outdir / f"{name}.tif"
        # Declare nodata in the file header. The array is already masked
        # outside the basin, but without a declared nodata value QGIS and
        # ArcGIS paint that area solid black instead of transparent, and
        # rasterio's masked read returns no mask at all -- so a correctly
        # clipped raster looks and behaves like an unclipped one the moment
        # it leaves Python.
        import numpy as np

        if da.rio.nodata is None:
            if np.issubdtype(da.dtype, np.floating):
                da = da.rio.write_nodata(np.nan, encoded=False)
            else:
                da = da.rio.write_nodata(0, encoded=False)
        da.rio.to_raster(fp, compress="deflate", tiled=True)
        return fp.name

    jobs = {
        "dem": lambda: _write_raster(self.dem(progress=progress), "dem"),
        "landcover": lambda: _write_raster(
            self.landcover(progress=progress), "landcover"
        ),
        "soil": lambda: _write_raster(self.soil("clay"), "soil_clay_0-5cm"),
        "surface_water": lambda: _write_raster(
            self.surface_water(progress=progress), "surface_water_occurrence"
        ),
        "precipitation": lambda: self._write_series(
            self.precipitation(start, end), outdir, "precipitation_chirps"
        ),
        "rivers": lambda: self._write_vector(
            self.rivers(progress=progress), outdir, "rivers"
        ),
        "lakes": lambda: self._write_vector(
            self.lakes(progress=progress), outdir, "lakes"
        ),
    }

    for name in layers:
        if name not in jobs:
            manifest["failed"][name] = f"unknown layer {name!r}"
            continue
        try:
            manifest["layers"][name] = jobs[name]()
        except Exception as exc:
            manifest["failed"][name] = f"{type(exc).__name__}: {exc}"

    # Explicit UTF-8: the licence text and the manifest both carry
    # non-ASCII, and Windows would otherwise write them as cp1252.
    (outdir / "LICENSES.txt").write_text(self.license_report(),
                                         encoding="utf-8")
    manifest["layers"]["licenses"] = "LICENSES.txt"
    (outdir / "manifest.json").write_text(
        json.dumps(manifest, indent=2, default=str), encoding="utf-8")
    return manifest

drainage_density(*, dem=None, thresholds=None, chosen_km2=None, **kwargs)

Drainage density across the range of defensible thresholds.

Not one number: the density moves by a large factor across channel initiation thresholds that are all defensible, so the curve and the threshold used are what make the figure comparable with anyone else's.

Source code in basinkit/basin.py

def drainage_density(self, *, dem=None, thresholds=None,
                     chosen_km2=None, **kwargs):
    """Drainage density across the range of defensible thresholds.

    Not one number: the density moves by a large factor across channel
    initiation thresholds that are all defensible, so the curve and the
    threshold used are what make the figure comparable with anyone else's.
    """
    from .terrain import drainage_density

    return drainage_density(self.dem(**kwargs) if dem is None else dem,
                            thresholds=thresholds, chosen_km2=chosen_km2)

explore(**kwargs)

Interactive map of the basin. Needs pip install 'basinkit[viz]'.

Source code in basinkit/basin.py

def explore(self, **kwargs):
    """Interactive map of the basin. Needs ``pip install 'basinkit[viz]'``."""
    from .viz import explore

    return explore(self, **kwargs)

export_3d(path, **kwargs)

Write an interactive 3D page for this basin: terrain, imagery, rivers.

One self-contained HTML file with everything embedded, so it opens with no network. See :func:basinkit.viz3d.export_3d for the options.

basin.export_3d("koshi.html")
basin.export_3d("koshi.html", texture=None)   # elevation only, small

This is a way of looking at the layers this package fetches. It makes no claim the other methods do not already make.

Source code in basinkit/basin.py

def export_3d(self, path: str | Path, **kwargs):
    """Write an interactive 3D page for this basin: terrain, imagery, rivers.

    One self-contained HTML file with everything embedded, so it opens with
    no network. See :func:`basinkit.viz3d.export_3d` for the options.

        basin.export_3d("koshi.html")
        basin.export_3d("koshi.html", texture=None)   # elevation only, small

    This is a way of looking at the layers this package fetches. It makes no
    claim the other methods do not already make.
    """
    from .viz3d import export_3d

    return export_3d(self, path, **kwargs)

flow_accumulation(*, dem=None, **kwargs)

Cells draining through each cell, routed over this basin's own DEM.

Source code in basinkit/basin.py

def flow_accumulation(self, *, dem=None, **kwargs):
    """Cells draining through each cell, routed over this basin's own DEM."""
    from .terrain import flow_accumulation

    return flow_accumulation(self.dem(**kwargs) if dem is None else dem)

from_file(path) classmethod

Load a basin from any vector file geopandas can read.

Source code in basinkit/basin.py

@classmethod
def from_file(cls, path: str | Path) -> Basin:
    """Load a basin from any vector file geopandas can read."""
    import geopandas as gpd

    gdf = gpd.read_file(path)
    if gdf.crs and gdf.crs.to_epsg() != 4326:
        gdf = gdf.to_crs("EPSG:4326")
    return cls(gdf.union_all(), {"backend": "file", "path": str(path)})

from_geometry(geometry, provenance=None, crs=None) classmethod

Wrap a polygon you already have (a gauge basin, an official boundary).

crs is the CRS the geometry arrives in. Anything other than EPSG:4326 is reprojected, because every step downstream of here works in lon/lat. Leaving it out keeps the old behaviour and takes the geometry as EPSG:4326 already, so existing callers are unaffected.

Source code in basinkit/basin.py

@classmethod
def from_geometry(cls, geometry, provenance: dict | None = None,
                  crs: str | None = None) -> Basin:
    """Wrap a polygon you already have (a gauge basin, an official boundary).

    ``crs`` is the CRS the geometry arrives in. Anything other than
    EPSG:4326 is reprojected, because every step downstream of here works
    in lon/lat. Leaving it out keeps the old behaviour and takes the
    geometry as EPSG:4326 already, so existing callers are unaffected.
    """
    if crs is not None:
        import geopandas as gpd

        gs = gpd.GeoSeries([geometry], crs=crs)
        if gs.crs is not None and gs.crs.to_epsg() != 4326:
            geometry = gs.to_crs("EPSG:4326").iloc[0]
    return cls(geometry, provenance or {"backend": "user-supplied"})

from_point(lat, lon, *, backend='auto', **kwargs) classmethod

Delineate the basin upstream of an outlet coordinate.

Parameters:

Name	Type	Description	Default
backend	str	'auto' (default), 'hydrobasins', 'dem' or 'api'. See :mod:basinkit.delineate for what each one is good at.	'auto'

Source code in basinkit/basin.py

@classmethod
def from_point(
    cls, lat: float, lon: float, *, backend: str = "auto", **kwargs
) -> Basin:
    """Delineate the basin upstream of an outlet coordinate.

    Parameters
    ----------
    backend
        ``'auto'`` (default), ``'hydrobasins'``, ``'dem'`` or ``'api'``.
        See :mod:`basinkit.delineate` for what each one is good at.
    """
    from .delineate import delineate

    # Coerced before it is compared, because a string reached the
    # comparison and raised "'<=' not supported between str and int" from
    # inside this check -- a message about Python's type system in place
    # of the one written here about coordinates.
    try:
        lat = float(lat)
        lon = float(lon)
    except (TypeError, ValueError):
        raise ValueError(
            f"({lat!r}, {lon!r}) is not a coordinate pair. Pass numbers: "
            "from_point(26.87, 87.15)."
        ) from None
    if not (lat == lat and lon == lon) or lat in (float("inf"), float("-inf")) \
            or lon in (float("inf"), float("-inf")):
        raise ValueError(
            f"({lat}, {lon}) is not a coordinate pair: a latitude or "
            "longitude cannot be NaN or infinite."
        )
    if not -90 <= lat <= 90 or not -180 <= lon <= 180:
        raise ValueError(
            f"({lat}, {lon}) is not a valid lat/lon. Note the order is "
            "(lat, lon) -- swapping them is the usual cause."
        )
    geom, prov = delineate(lat, lon, backend=backend, **kwargs)
    return cls(geom, prov)

hand(*, dem=None, min_area_km2=1.0, **kwargs)

Height above the nearest drainage: the terrain layer flood work wants.

Source code in basinkit/basin.py

def hand(self, *, dem=None, min_area_km2: float = 1.0, **kwargs):
    """Height above the nearest drainage: the terrain layer flood work wants."""
    from .terrain import hand

    return hand(self.dem(**kwargs) if dem is None else dem,
                min_area_km2=min_area_km2)

hillshade(*, dem=None, azimuth=315.0, altitude=45.0, z_factor=1.0, **kwargs)

Shaded relief, 0 to 1, for a figure that reads as terrain.

Source code in basinkit/basin.py

def hillshade(self, *, dem=None, azimuth: float = 315.0,
              altitude: float = 45.0, z_factor: float = 1.0, **kwargs):
    """Shaded relief, 0 to 1, for a figure that reads as terrain."""
    from .terrain import hillshade

    return hillshade(self.dem(**kwargs) if dem is None else dem,
                     azimuth=azimuth, altitude=altitude, z_factor=z_factor)

lakes(min_area_km2=0.0, **kwargs)

HydroLAKES water bodies inside the basin.

Source code in basinkit/basin.py

def lakes(self, min_area_km2: float = 0.0, **kwargs):
    """HydroLAKES water bodies inside the basin."""
    from .sources.vectors import hydrolakes

    return self._acquired(
        hydrolakes(self.geometry, min_area_km2=min_area_km2, **kwargs), "lakes")

landcover(year=None, source='worldcover', **kwargs)

Land cover, as a 2-D array with its class legend in .attrs.

worldcover is ESA WorldCover at 10 m (2020 or 2021); esri is the ESRI / Impact Observatory annual series. Both return the same shape of object, and each carries its own legend -- the two number their classes differently, and code 10 is tree cover in one and cloud in the other.

year=None means 2021 for WorldCover and the latest year published for this location for ESRI.

Source code in basinkit/basin.py

def landcover(self, year: int | None = None, source: str = "worldcover",
              **kwargs):
    """Land cover, as a 2-D array with its class legend in ``.attrs``.

    ``worldcover`` is ESA WorldCover at 10 m (2020 or 2021); ``esri`` is the
    ESRI / Impact Observatory annual series. Both return the same shape of
    object, and each carries its own legend -- the two number their classes
    differently, and code 10 is tree cover in one and cloud in the other.

    ``year=None`` means 2021 for WorldCover and the latest year published
    for this location for ESRI.
    """
    from .sources.landcover import esri_lulc, worldcover

    if source == "worldcover":
        return self._acquired(
            worldcover(self.geometry, year=year or 2021, **kwargs), "landcover")
    if source == "esri":
        return self._acquired(
            esri_lulc(self.geometry, year=year, **kwargs), "landcover")
    raise ValueError(f"Unknown land cover source {source!r}: use 'worldcover' or 'esri'.")

landcover_change(start, end, *, source='esri', **kwargs)

What became what, between two years, in square kilometres.

Only the ESRI annual series carries a year for every year, so that is the default here even though landcover() defaults to WorldCover. The diagonal is the ground that did not change; everything off it is a transition, and the table is the honest form of a deforestation or urban-growth figure because it shows what the loss became.

Source code in basinkit/basin.py

def landcover_change(self, start: int, end: int, *, source: str = "esri",
                     **kwargs):
    """What became what, between two years, in square kilometres.

    Only the ESRI annual series carries a year for every year, so that is
    the default here even though ``landcover()`` defaults to WorldCover.
    The diagonal is the ground that did not change; everything off it is a
    transition, and the table is the honest form of a deforestation or
    urban-growth figure because it shows what the loss became.
    """
    import numpy as np
    import pandas as pd

    from .sources.landcover import _declared_classes
    from .terrain import cell_area_km2

    first = self.landcover(year=start, source=source, **kwargs)
    second = self.landcover(year=end, source=source, **kwargs)
    if first.shape != second.shape:
        second = second.rio.reproject_match(first)

    legend = _declared_classes(first) or {}
    a = np.asarray(first.values, dtype="float64")
    b = np.asarray(second.values, dtype="float64")
    areas = cell_area_km2(first)
    usable = np.isfinite(a) & np.isfinite(b) & (a > 0) & (b > 0)

    rows = []
    for code_from in np.unique(a[usable]):
        for code_to in np.unique(b[usable]):
            picked = usable & (a == code_from) & (b == code_to)
            if not picked.any():
                continue
            rows.append({
                "from": legend.get(int(code_from), int(code_from)),
                "to": legend.get(int(code_to), int(code_to)),
                "area_km2": round(float(areas[picked].sum()), 3),
                "changed": bool(code_from != code_to),
            })
    frame = pd.DataFrame(rows).sort_values("area_km2", ascending=False)
    frame.attrs["years"] = (start, end)
    return frame.reset_index(drop=True)

landform(*, dem=None, window=11, **kwargs)

Six slope-position classes: valley, slopes, flat, ridge.

Source code in basinkit/basin.py

def landform(self, *, dem=None, window: int = 11, **kwargs):
    """Six slope-position classes: valley, slopes, flat, ridge."""
    from .terrain import landform

    return landform(self.dem(**kwargs) if dem is None else dem, window=window)

landsat(start, end, *, cloud_cover=20, bands=None, composite='median', **kwargs)

Landsat Collection 2 Level-2 over the basin (1982 to present).

Source code in basinkit/basin.py

def landsat(self, start: str, end: str, *, cloud_cover: float = 20,
            bands: list[str] | None = None, composite: str | None = "median",
            **kwargs):
    """Landsat Collection 2 Level-2 over the basin (1982 to present)."""
    from .sources.stac import composite as reduce_time
    from .sources.stac import stac_search, stac_stack

    stack_kw = _stack_kwargs(kwargs)
    items = stac_search(
        "landsat", geometry=self.geometry, start=start, end=end,
        cloud_cover=cloud_cover, **kwargs
    )
    ds = stac_stack(items, self.geometry, bands=bands or ["blue", "green", "red", "nir08"],
                    collection="landsat", **stack_kw)
    return reduce_time(ds, composite) if composite else ds

license_report(layers=None)

Attribution and licence text for the layers you used.

Print this into your methods section. Every layer basinkit fetches by default is CC BY 4.0 or more permissive, which means it can be redistributed and used commercially -- but only if it is attributed.

Source code in basinkit/basin.py

def license_report(self, layers: tuple[str, ...] | None = None) -> str:
    """Attribution and licence text for the layers you used.

    Print this into your methods section. Every layer basinkit fetches by
    default is CC BY 4.0 or more permissive, which means it can be
    redistributed and used commercially -- but only if it is attributed.
    """
    keys = layers or catalog.DEFAULT_STACK
    lines = ["Data sources and licences", "=" * 26, ""]
    for key in keys:
        try:
            ds = catalog.get(key)
        except KeyError:
            continue
        lines.append(f"{ds.name}")
        lines.append(f"    Licence : {ds.license}")
        lines.append(f"    Access  : {ds.route}")
        if ds.citation:
            lines.append(f"    Cite    : {ds.citation}")
        if not ds.commercial_ok:
            lines.append("    WARNING : commercial use not permitted")
        if not ds.redistributable:
            lines.append("    WARNING : redistribution not permitted")
        lines.append("")
    return "\n".join(lines)

morphometry(**kwargs)

The classical Horton-Strahler-Schumm morphometric parameters.

m = basin.morphometry()
m["areal"]["drainage_density_km_per_km2"]
m["network"]        # one row per Strahler order

Counts Strahler streams, not the reaches a river dataset splits them into, and measures area, perimeter and every length in one equal-area projection. Both matter: on the Koshi, counting reaches turns the bifurcation ratios into values that are not physically possible.

See :mod:basinkit.morphometry for what each parameter is and for what the numbers can and cannot be compared against.

Source code in basinkit/basin.py

def morphometry(self, **kwargs) -> dict:
    """The classical Horton-Strahler-Schumm morphometric parameters.

        m = basin.morphometry()
        m["areal"]["drainage_density_km_per_km2"]
        m["network"]        # one row per Strahler order

    Counts Strahler *streams*, not the reaches a river dataset splits them
    into, and measures area, perimeter and every length in one equal-area
    projection. Both matter: on the Koshi, counting reaches turns the
    bifurcation ratios into values that are not physically possible.

    See :mod:`basinkit.morphometry` for what each parameter is and for what
    the numbers can and cannot be compared against.
    """
    from .morphometry import morphometry

    return morphometry(self, **kwargs)

plot(**kwargs)

Static matplotlib figure: hypsometry, boundary and river network.

Source code in basinkit/basin.py

def plot(self, **kwargs):
    """Static matplotlib figure: hypsometry, boundary and river network."""
    from .viz import plot

    return plot(self, **kwargs)

precipitation(start=2000, end=None, source='chirps', **kwargs)

Basin-mean rainfall time series. chirps, persiann or terraclimate.

Every source returns the same kind of object: a named xarray.DataArray of precipitation indexed by time. It did not used to -- terraclimate handed back a whole Dataset, so code that worked on one source raised an attribute error on another, and the function's own name promised one quantity while returning thirteen. For the rest of the TerraClimate variables call :meth:water_balance, or sources.climate.terraclimate directly.

Source code in basinkit/basin.py

def precipitation(self, start=2000, end=None, source: str = "chirps", **kwargs):
    """Basin-mean rainfall time series. ``chirps``, ``persiann`` or ``terraclimate``.

    Every source returns the same kind of object: a named
    ``xarray.DataArray`` of precipitation indexed by time. It did not used
    to -- ``terraclimate`` handed back a whole ``Dataset``, so code that
    worked on one source raised an attribute error on another, and the
    function's own name promised one quantity while returning thirteen.
    For the rest of the TerraClimate variables call
    :meth:`water_balance`, or ``sources.climate.terraclimate`` directly.
    """
    from .sources.climate import chirps, persiann, terraclimate

    if source == "chirps":
        return self._acquired(
            chirps(self.geometry, start, end, **kwargs),
            "precipitation", source="chirps", years=f"{start}-{end or 'latest'}")
    if source == "persiann":
        return self._acquired(
            persiann(self.geometry, str(start), end, **kwargs),
            "precipitation", source="persiann", years=f"{start}-{end or 'latest'}")
    if source == "terraclimate":
        ds = terraclimate(self.geometry, ("ppt",), int(start), end, **kwargs)
        out = ds["ppt"]
        out.name = "precipitation"
        out.attrs.setdefault("units", "mm/month")
        # The dataset-level provenance -- which years are missing, whether
        # the clip held -- belongs to this series too, and is the only
        # record that it is incomplete.
        for key, value in ds.attrs.items():
            out.attrs.setdefault(key, value)
        return self._acquired(out, "precipitation",
                              source="terraclimate",
                              years=f"{start}-{end or 'latest'}")
    raise ValueError(
        f"Unknown precipitation source {source!r}: use 'chirps', 'persiann' "
        "or 'terraclimate'."
    )

precipitation_trend(*, start=1985, end=None, alpha=0.05, **kwargs)

Is rainfall over this basin trending, and by how much a year.

Source code in basinkit/basin.py

def precipitation_trend(self, *, start: int = 1985, end: int | None = None,
                        alpha: float = 0.05, **kwargs):
    """Is rainfall over this basin trending, and by how much a year."""
    from .climate import trend

    series = self.precipitation(start=start, end=end, **kwargs)
    s = series.to_series() if hasattr(series, "to_series") else series
    yearly = s.groupby(s.index.year).sum()
    yearly.index = __import__("pandas").to_datetime(
        [f"{int(y)}-01-01" for y in yearly.index]
    )
    return trend(yearly, alpha=alpha)

report(path, *, title=None, dem=None, rivers=None, **kwargs)

Write the eight-page PDF report for this basin.

Cover with the elevation data suitability grade, elevation and hypsometry, slope and aspect, the shape of the ground, the channel network and Horton's laws, the morphometric table with symbols and references, what the answer rests on, and a methods page.

max_pixels= sets the pixel budget for the elevation grid the report is built from; the cover states the resolution it worked out to. The default is :data:basinkit.report.REPORT_MAX_PIXELS and not the package-wide one, because every page derives arrays of its own from this grid.

Source code in basinkit/basin.py

def report(self, path, *, title=None, dem=None, rivers=None, **kwargs):
    """Write the eight-page PDF report for this basin.

    Cover with the elevation data suitability grade, elevation and
    hypsometry, slope and aspect, the shape of the ground, the channel
    network and Horton's laws, the morphometric table with symbols and
    references, what the answer rests on, and a methods page.

    ``max_pixels=`` sets the pixel budget for the elevation grid the
    report is built from; the cover states the resolution it worked out
    to. The default is :data:`basinkit.report.REPORT_MAX_PIXELS` and not
    the package-wide one, because every page derives arrays of its own
    from this grid.
    """
    from .report import REPORT_MAX_PIXELS, report

    # The budget has to be applied here, because this method fetched the
    # grid itself and handed report() a finished one -- so report()'s own
    # max_pixels was never reached, and every caller that goes through
    # this method, which is all of them including the QGIS and ArcGIS
    # algorithms, still read at the package-wide 100 Mpx default. That is
    # what reached 6.1 GB on the Koshi and was killed by the kernel. The
    # earlier fix was real and was in the wrong function.
    if dem is None:
        kwargs.setdefault("max_pixels", REPORT_MAX_PIXELS)
        dem = self.dem(**kwargs)
    return report(self, path, title=title, dem=dem, rivers=rivers)

rivers(min_order=0, **kwargs)

HydroRIVERS reaches inside the basin, with discharge and stream order.

Source code in basinkit/basin.py

def rivers(self, min_order: int = 0, **kwargs):
    """HydroRIVERS reaches inside the basin, with discharge and stream order."""
    from .sources.vectors import hydrorivers

    return self._acquired(
        hydrorivers(self.geometry, min_order=min_order, **kwargs), "rivers")

roughness(*, dem=None, **kwargs)

Local relief within each cell's 3x3 neighbourhood.

Source code in basinkit/basin.py

def roughness(self, *, dem=None, **kwargs):
    """Local relief within each cell's 3x3 neighbourhood."""
    from .terrain import roughness

    return roughness(self.dem(**kwargs) if dem is None else dem)

sentinel1(start, end, *, bands=None, composite='median', **kwargs)

Sentinel-1 RTC: terrain-corrected radar, so it works through cloud.

Source code in basinkit/basin.py

def sentinel1(self, start: str, end: str, *, bands: list[str] | None = None,
              composite: str | None = "median", **kwargs):
    """Sentinel-1 RTC: terrain-corrected radar, so it works through cloud."""
    from .sources.stac import composite as reduce_time
    from .sources.stac import stac_search, stac_stack

    stack_kw = _stack_kwargs(kwargs)
    items = stac_search(
        "sentinel1_rtc", geometry=self.geometry, start=start, end=end, **kwargs
    )
    ds = stac_stack(items, self.geometry, bands=bands or ["vv", "vh"],
                    collection="sentinel1_rtc", **stack_kw)
    return reduce_time(ds, composite) if composite else ds

sentinel2(start, end, *, cloud_cover=20, bands=None, composite='median', **kwargs)

Sentinel-2 L2A over the basin, cloud-filtered and optionally composited.

Source code in basinkit/basin.py

def sentinel2(self, start: str, end: str, *, cloud_cover: float = 20,
              bands: list[str] | None = None, composite: str | None = "median",
              **kwargs):
    """Sentinel-2 L2A over the basin, cloud-filtered and optionally composited."""
    from .sources.stac import composite as reduce_time
    from .sources.stac import stac_search, stac_stack

    stack_kw = _stack_kwargs(kwargs)
    items = stac_search(
        "sentinel2", geometry=self.geometry, start=start, end=end,
        cloud_cover=cloud_cover, **kwargs
    )
    ds = stac_stack(items, self.geometry, bands=bands or ["blue", "green", "red", "nir"],
                    collection="sentinel2", **stack_kw)
    return reduce_time(ds, composite) if composite else ds

slope(*, dem=None, degrees=True, **kwargs)

Steepest descent at every cell, in degrees by default.

Source code in basinkit/basin.py

def slope(self, *, dem=None, degrees: bool = True, **kwargs):
    """Steepest descent at every cell, in degrees by default."""
    from .terrain import slope

    return slope(self.dem(**kwargs) if dem is None else dem, degrees=degrees)

soil(prop='clay', depth='0-5cm', **kwargs)

A SoilGrids property. See :data:basinkit.sources.soil.PROPERTIES.

Source code in basinkit/basin.py

def soil(self, prop: str = "clay", depth: str = "0-5cm", **kwargs):
    """A SoilGrids property. See :data:`basinkit.sources.soil.PROPERTIES`."""
    from .sources.soil import soilgrids

    return self._acquired(
        soilgrids(self.geometry, prop=prop, depth=depth, **kwargs),
        "soil", property=prop, depth=depth)

spi(*, scale=3, start=1985, end=None, **kwargs)

Standardized Precipitation Index over this basin, month by month.

Source code in basinkit/basin.py

def spi(self, *, scale: int = 3, start: int = 1985, end: int | None = None,
        **kwargs):
    """Standardized Precipitation Index over this basin, month by month."""
    from .climate import spi

    series = self.precipitation(start=start, end=end, **kwargs)
    s = series.to_series() if hasattr(series, "to_series") else series
    return spi(s, scale=scale)

streams(*, dem=None, min_area_km2=1.0, **kwargs)

The channel network routed from this basin's own elevation.

Source code in basinkit/basin.py

def streams(self, *, dem=None, min_area_km2: float = 1.0, **kwargs):
    """The channel network routed from this basin's own elevation."""
    from .terrain import streams

    return streams(self.dem(**kwargs) if dem is None else dem,
                   min_area_km2=min_area_km2)

subbasins(**kwargs)

The sub-catchments this basin is assembled from, with their routing.

delineate dissolves the HydroBASINS units it walked into a single polygon. This hands back the pieces instead, each with NEXT_DOWN, which is the routing graph itself: every distributed model wants sub-catchments and the links between them, and nothing here has to be inferred from geometry afterwards.

Available for the hydrobasins backend, which is the default.

Source code in basinkit/basin.py

def subbasins(self, **kwargs):
    """The sub-catchments this basin is assembled from, with their routing.

    ``delineate`` dissolves the HydroBASINS units it walked into a single
    polygon. This hands back the pieces instead, each with ``NEXT_DOWN``,
    which is the routing graph itself: every distributed model wants
    sub-catchments and the links between them, and nothing here has to be
    inferred from geometry afterwards.

    Available for the ``hydrobasins`` backend, which is the default.
    """
    from .delineate.hydrobasins import upstream_units

    prov = self.provenance
    if prov.get("backend") != "hydrobasins":
        raise ValueError(
            "subbasins() reads the units the HydroBASINS traversal walked, "
            f"and this basin came from the {prov.get('backend', 'unknown')!r} "
            "backend. Delineate with backend='hydrobasins' to get them."
        )
    return upstream_units(
        prov["region"], prov["outlet_hybas_id"], prov.get("level", 12), **kwargs
    )

summary(*, terrain=True, landcover=True)

A one-call characterisation of the basin.

Source code in basinkit/basin.py

def summary(self, *, terrain: bool = True, landcover: bool = True) -> dict:
    """A one-call characterisation of the basin."""
    out: dict[str, Any] = {
        "area_km2": round(self.area_km2, 2),
        "centroid_lat_lon": [round(v, 5) for v in self.centroid],
        "bounds": [round(v, 5) for v in self.bounds],
        "bbox_efficiency": round(self.bbox_efficiency, 3),
        "provenance": self.provenance,
    }
    if terrain:
        try:
            out["terrain"] = self.terrain_stats()
        except Exception as exc:
            out["terrain"] = {"error": str(exc)}
    if landcover:
        try:
            from .sources.landcover import class_fractions

            out["landcover_fractions"] = class_fractions(self.landcover())
        except Exception as exc:
            out["landcover_fractions"] = {"error": str(exc)}
    return out

surface_water(layer='occurrence', **kwargs)

JRC Global Surface Water: a pre-reduced 37-year Landsat water stack.

Source code in basinkit/basin.py

def surface_water(self, layer: str = "occurrence", **kwargs):
    """JRC Global Surface Water: a pre-reduced 37-year Landsat water stack."""
    from .sources.water import global_surface_water

    return self._acquired(
        global_surface_water(self.geometry, layer=layer, **kwargs),
        "surface_water", layer=layer)

terrain_stats(*, dem=None, **kwargs)

Elevation, relief and mean slope: the standard morphometry.

Pass dem= the elevation you already hold. Fetching a second one here is not only wasteful: if the first was coarsened to a pixel budget and this one is not, the two disagree about the highest point of the same basin, and a report that prints both looks wrong whichever number the reader believes.

Source code in basinkit/basin.py

def terrain_stats(self, *, dem=None, **kwargs) -> dict:
    """Elevation, relief and mean slope: the standard morphometry.

    Pass ``dem=`` the elevation you already hold. Fetching a second one
    here is not only wasteful: if the first was coarsened to a pixel
    budget and this one is not, the two disagree about the highest point
    of the same basin, and a report that prints both looks wrong whichever
    number the reader believes.
    """
    import numpy as np

    elev = dem if dem is not None else self.dem(**kwargs)
    vals = np.asarray(elev.values, dtype="float64")
    vals = vals[np.isfinite(vals)]
    if vals.size == 0:
        return {}

    # Slope comes from the same routine the slope map uses, which keeps
    # nodata as nodata. Filling the space around the polygon with zeros
    # before differencing puts a plain around the basin and a cliff at its
    # edge, and the mean then answers to those instead: on the Helmand
    # above the Kajaki Dam, which fills 43% of its bounding box, it read 7
    # degrees where the basin's own slopes average 16. That is the
    # bounding-box mistake this package exists to avoid, made against
    # itself.
    from .terrain import slope as _slope_of

    slope = np.asarray(_slope_of(elev).values, dtype="float64")
    finite_slope = slope[np.isfinite(slope)]

    return {
        "area_km2": round(self.area_km2, 2),
        "elev_min_m": round(float(vals.min()), 1),
        "elev_max_m": round(float(vals.max()), 1),
        "elev_mean_m": round(float(vals.mean()), 1),
        "relief_m": round(float(vals.max() - vals.min()), 1),
        "slope_mean_deg": (round(float(finite_slope.mean()), 2)
                           if finite_slope.size else None),
        "bbox_efficiency": round(self.bbox_efficiency, 3),
    }

tpi(*, dem=None, window=11, **kwargs)

Topographic position index: height above the surrounding ground.

Source code in basinkit/basin.py

def tpi(self, *, dem=None, window: int = 11, **kwargs):
    """Topographic position index: height above the surrounding ground."""
    from .terrain import tpi

    return tpi(self.dem(**kwargs) if dem is None else dem, window=window)

tri(*, dem=None, **kwargs)

Terrain ruggedness index: how far a cell sits from its neighbours.

Source code in basinkit/basin.py

def tri(self, *, dem=None, **kwargs):
    """Terrain ruggedness index: how far a cell sits from its neighbours."""
    from .terrain import tri

    return tri(self.dem(**kwargs) if dem is None else dem)

twi(*, dem=None, **kwargs)

Topographic wetness index: where water gathers and ground saturates.

Source code in basinkit/basin.py

def twi(self, *, dem=None, **kwargs):
    """Topographic wetness index: where water gathers and ground saturates."""
    from .terrain import twi

    return twi(self.dem(**kwargs) if dem is None else dem)

validation_envelope(*, suitability=None)

Which part of the validation set this basin belongs to.

The accuracy figures published with basinkit are medians over 2,740 delineations at 2,550 gauges, and they differ by two orders of magnitude between the largest and the smallest band. This reports the error that was measured at this basin's size, and names the conditions the validation never covered, so that a figure which does not apply here is not quoted as though it did.

Reads only the area, so it downloads nothing. Pass a :meth:dem_suitability result to compare the terrain descriptors too.

Source code in basinkit/basin.py

def validation_envelope(self, *, suitability=None):
    """Which part of the validation set this basin belongs to.

    The accuracy figures published with basinkit are medians over 2,740
    delineations at 2,550 gauges, and they differ by two orders of
    magnitude between the largest and the smallest band. This reports the
    error that was measured at *this* basin's size, and names the
    conditions the validation never covered, so that a figure which does
    not apply here is not quoted as though it did.

    Reads only the area, so it downloads nothing. Pass a
    :meth:`dem_suitability` result to compare the terrain descriptors too.
    """
    from .envelope import validation_envelope

    return validation_envelope(self, suitability=suitability)

water_balance(start=2000, end=None)

Monthly P / AET / PET / Q / soil-moisture balance from TerraClimate.

Source code in basinkit/basin.py

def water_balance(self, start: int = 2000, end: int | None = None):
    """Monthly P / AET / PET / Q / soil-moisture balance from TerraClimate."""
    from .sources.climate import water_balance

    return water_balance(self.geometry, start, end)

zonal(values, zones=None, *, bins=None, labels=None, **kwargs)

Summarise one layer inside the classes of another.

basin.zonal(basin.precipitation_grid, basin.landcover()) answers how much rain falls on each land-cover class; passing bins cuts a continuous layer such as elevation into bands instead. Area is summed from the true size of every cell, so a basin spanning several degrees of latitude is not weighted towards its southern edge.

Source code in basinkit/basin.py

def zonal(self, values, zones=None, *, bins=None, labels=None, **kwargs):
    """Summarise one layer inside the classes of another.

    ``basin.zonal(basin.precipitation_grid, basin.landcover())`` answers
    how much rain falls on each land-cover class; passing ``bins`` cuts a
    continuous layer such as elevation into bands instead. Area is summed
    from the true size of every cell, so a basin spanning several degrees
    of latitude is not weighted towards its southern edge.
    """
    from .zonal import zonal

    if zones is None:
        zones = self.landcover(**kwargs)
    return zonal(values, zones, labels=labels, bins=bins)

Upstream basin delineation backends.

Three backends, because no single one is right for every basin:

hydrobasins Graph traversal over HydroBASINS level-12 sub-basins. Global, CC BY 4.0, offline once cached, and fast at any basin size because upstream aggregation is a walk over the NEXT_DOWN field rather than a raster fill.

**Base grid: 15 arc-seconds, about 460 m, from February 2000 SRTM.**
HydroBASINS is extracted from the HydroSHEDS core layers at that
resolution, so this backend inherits it. Two consequences worth stating
plainly: the resolution floor is the level-12 unit (~130 km2 median), and
the flow network is a quarter-century-old DEM. For most basins that is
fine. For a small or heavily modified catchment it is not, and the other
two backends exist for exactly that case.

dem D8 flow routing with pyflwdir over a freshly downloaded Copernicus DEM window. Base grid: 1 arc-second, about 30 m, from 2011-2015 radar. Fifteen times finer than the default and a decade newer, so it is the right instrument for small catchments; its working range reaches about 10,000 km2, since cost grows with basin area and the window has to contain the whole basin.

api The public Global Watersheds service, backed by MERIT-Hydro. Base grid: 3 arc-seconds, about 90 m, multi-error-removed. Five times finer than the default and hydrologically conditioned rather than raw SRTM. No download at all, so it is the fastest first look -- but it is one research group's server and MERIT-Hydro's licence is non-commercial, so basinkit never makes it the default and records both facts in provenance.

auto picks between them from the drainage area implied by the outlet.

Base resolution, side by side:

=============== ================= ================== ==================== backend grid source conditioned =============== ================= ================== ==================== hydrobasins 15 arc-sec, 460 m SRTM, Feb 2000 HydroSHEDS api 3 arc-sec, 90 m MERIT-Hydro yes, error-removed dem 1 arc-sec, 30 m Copernicus, 2011-15 no, routed on the fly =============== ================= ================== ====================

delineate(lat, lon, backend='auto', **kwargs)

Delineate the upstream basin of (lat, lon).

Returns:

Type	Description
(geometry, dict)	The basin polygon in EPSG:4326 and a provenance dict recording which backend and dataset version produced it. The provenance travels with the basin so a result is always attributable.

Source code in basinkit/delineate/__init__.py

def delineate(lat: float, lon: float, backend: str = "auto", **kwargs):
    """Delineate the upstream basin of ``(lat, lon)``.

    Returns
    -------
    (shapely.geometry, dict)
        The basin polygon in EPSG:4326 and a provenance dict recording which
        backend and dataset version produced it. The provenance travels with
        the basin so a result is always attributable.
    """
    if backend == "auto":
        return _auto(lat, lon, **kwargs)
    try:
        fn = _BACKENDS[backend]
    except KeyError:
        raise ValueError(
            f"Unknown backend {backend!r}. Choose from: "
            f"{', '.join(_BACKENDS)} or 'auto'."
        ) from None
    geom, prov = fn(lat, lon, **kwargs)
    _weigh(lat, lon, prov, verify=kwargs.get("verify", True))
    return geom, prov

delineate_api(lat, lon, *, precision='high', **_)

Fetch an upstream basin polygon from the Global Watersheds service.

No local data is downloaded, which makes this the quickest way to look at a basin. The service auto-downgrades to low precision above ~50,000 km2.

This is a courtesy endpoint run by one research group, not managed infrastructure. basinkit therefore never selects it automatically, and stamps backend='api' into provenance so a downstream reader can tell that the geometry did not come from a versioned dataset.

Source code in basinkit/delineate/api.py

def delineate_api(lat: float, lon: float, *, precision: str = "high", **_):
    """Fetch an upstream basin polygon from the Global Watersheds service.

    No local data is downloaded, which makes this the quickest way to look at a
    basin. The service auto-downgrades to low precision above ~50,000 km2.

    This is a courtesy endpoint run by one research group, not managed
    infrastructure. basinkit therefore never selects it automatically, and
    stamps ``backend='api'`` into provenance so a downstream reader can tell
    that the geometry did not come from a versioned dataset.
    """
    from shapely.geometry import shape
    from shapely.ops import unary_union

    data = get_json(
        _ENDPOINT, params={"lat": lat, "lng": lon, "precision": precision}, timeout=120
    )
    feats = data.get("features") or []
    if not feats:
        raise DelineationError(
            f"The Global Watersheds service returned no basin for ({lat}, {lon}). "
            "The point is probably off the mapped river network, or outside its "
            "60S-85N coverage. Try backend='dem'."
        )

    geom = unary_union([shape(f["geometry"]) for f in feats])
    from ..clip import basin_area_km2

    return geom, {
        "backend": "api",
        "service": "mghydro Global Watersheds",
        "source_dataset": "MERIT-Hydro (~90 m)",
        "precision": precision,
        "outlet": (lat, lon),
        "area_km2": round(basin_area_km2(geom), 2),
        "license_note": (
            "Geometry derived from MERIT-Hydro (CC BY-NC 4.0 / ODbL). Not "
            "redistributable under a permissive licence -- use the "
            "'hydrobasins' backend (CC BY 4.0) for anything you intend to publish."
        ),
    }

delineate_dem(lat, lon, *, window_deg=0.5, product='cop30', snap_px=12, min_uparea_km2=1.0, max_window_deg=4.0, max_memory_gb=DEFAULT_MAX_MEMORY_GB, streams=None, burn_depth_m=20.0, progress=True, **_)

Delineate the upstream basin by D8 routing on a Copernicus DEM window.

Parameters:

Name	Type	Description	Default
window_deg	float	Half-width of the initial DEM window in degrees. Grown automatically if the basin reaches the edge.	0.5
snap_px	int	Radius, in pixels, of the search for the true channel cell.	12
streams		Lines to force the flow along, in EPSG:4326: a vector file path, a GeoDataFrame, or shapely geometries. Use it where the elevation model cannot see the channel -- a culverted urban stream, a canal across a divide, a valley the model dams with a road embankment. The flow path is then the one supplied rather than the one found, which is recorded in the provenance because it changes what the answer means.	None
burn_depth_m	float	How far the burned cells are lowered. Deeper than the relief it has to beat, shallower than a trench the flow cannot leave.	20.0
max_window_deg	float	Stop growing at this half-width and raise instead of silently downloading the continent.	4.0
max_memory_gb	float	Refuse, with the figure named, rather than attempt a routing window whose estimated peak memory exceeds this. Above roughly 10,000 km2 the window needed is large enough for that estimate to matter; below it, this never comes up.	DEFAULT_MAX_MEMORY_GB

Source code in basinkit/delineate/dem.py

@_quiet_matmul
def delineate_dem(
    lat: float,
    lon: float,
    *,
    window_deg: float = 0.5,
    product: str = "cop30",
    snap_px: int = 12,
    min_uparea_km2: float = 1.0,
    max_window_deg: float = 4.0,
    max_memory_gb: float = DEFAULT_MAX_MEMORY_GB,
    streams=None,
    burn_depth_m: float = 20.0,
    progress: bool = True,
    **_,
):
    """Delineate the upstream basin by D8 routing on a Copernicus DEM window.

    Parameters
    ----------
    window_deg
        Half-width of the initial DEM window in degrees. Grown automatically if
        the basin reaches the edge.
    snap_px
        Radius, in pixels, of the search for the true channel cell.
    streams
        Lines to force the flow along, in EPSG:4326: a vector file path, a
        GeoDataFrame, or shapely geometries. Use it where the elevation model
        cannot see the channel -- a culverted urban stream, a canal across a
        divide, a valley the model dams with a road embankment. The flow path
        is then the one supplied rather than the one found, which is recorded
        in the provenance because it changes what the answer means.
    burn_depth_m
        How far the burned cells are lowered. Deeper than the relief it has to
        beat, shallower than a trench the flow cannot leave.
    max_window_deg
        Stop growing at this half-width and raise instead of silently
        downloading the continent.
    max_memory_gb
        Refuse, with the figure named, rather than attempt a routing window
        whose estimated peak memory exceeds this. Above roughly 10,000 km2 the
        window needed is large enough for that estimate to matter; below it,
        this never comes up.
    """
    pyflwdir = _require_pyflwdir()
    from shapely.geometry import shape
    from shapely.ops import unary_union

    from ..sources.dem import dem as fetch_dem

    half = window_deg
    while True:
        # Checked before the download, not after: a window too large to route
        # is also too large to be worth fetching, and the old order spent the
        # transfer first and was killed second.
        planned = int(round(2 * half * 3600)) ** 2
        if _estimated_gb(planned) > max_memory_gb:
            raise _refuse(planned, half, max_memory_gb)

        bounds = (lon - half, lat - half, lon + half, lat + half)
        elev = fetch_dem(bounds=bounds, product=product, clip=False, progress=progress)
        elev = elev.squeeze()

        arr = np.asarray(elev.values, dtype="float32")
        if not np.isfinite(arr).any():
            raise DelineationError(
                f"DEM window around ({lat}, {lon}) is entirely nodata -- an "
                "ocean-only extent, or outside the product's coverage."
            )
        arr = np.where(np.isfinite(arr), arr, -9999.0)

        transform = elev.rio.transform()

        burn_record = {"burned": False}
        if streams is not None:
            from ..condition import burn_streams

            arr, burn_record = burn_streams(
                arr, transform, streams, depth_m=burn_depth_m)
        # outlets='edge', not 'min'. With outlets='min' pyflwdir routes the
        # whole window toward its single lowest cell, which on a cropped DEM
        # drags the network away from the real channels -- a river with
        # thousands of km2 upstream can end up with a fraction of a km2 of
        # accumulation. 'edge' lets flow leave wherever it reaches the boundary,
        # which is what a window cut out of a larger landscape actually does.
        flw = pyflwdir.from_dem(
            data=arr, nodata=-9999.0, transform=transform, latlon=True,
            outlets="edge",
        )
        uparea = flw.upstream_area(unit="km2")

        xs = np.asarray(elev.x.values)
        ys = np.asarray(elev.y.values)
        col0 = int(np.abs(xs - lon).argmin())
        row0 = int(np.abs(ys - lat).argmin())

        row, col, snapped_area = _snap_to_stream(
            flw, uparea, row0, col0, search_px=snap_px, min_uparea_km2=min_uparea_km2
        )
        snap_px_moved = int(max(abs(row - row0), abs(col - col0)))

        mask = flw.basins(idxs=np.array([row * flw.shape[1] + col]))
        mask = (mask > 0).astype("uint8")

        if mask.sum() == 0:
            raise DelineationError("Flow routing produced an empty basin.")

        touches_edge = bool(
            mask[0, :].any() or mask[-1, :].any() or mask[:, 0].any() or mask[:, -1].any()
        )
        if touches_edge and half < max_window_deg:
            half *= 2
            nxt = int(round(2 * half * 3600)) ** 2
            warnings.warn(
                f"The basin reaches the edge of the window, so it is being "
                f"doubled to {2 * half:.1f} degrees -- "
                f"{nxt / 1e6:,.0f} megapixels, about "
                f"{_estimated_gb(nxt):.1f} GB. A basin that needs this much "
                "window is at the upper end of what the DEM backend is for; "
                "backend='hydrobasins' costs neither.",
                stacklevel=2)
            continue
        break

    if touches_edge:
        raise DelineationError(
            f"The basin still reaches the edge of a {2 * half:.1f} degree DEM window. "
            "It is too large for the DEM backend -- use backend='hydrobasins', "
            "which handles any size."
        )

    from rasterio.features import shapes as rio_shapes

    geoms = [
        shape(geom)
        for geom, val in rio_shapes(mask, mask=mask.astype(bool), transform=transform)
        if val == 1
    ]
    if not geoms:
        raise DelineationError("Could not vectorise the delineated basin mask.")

    geom = unary_union(geoms).buffer(0)

    from ..clip import basin_area_km2

    return geom, {
        "backend": "dem",
        "source_dataset": f"{product} via D8 routing (pyflwdir)",
        "outlet": (lat, lon),
        "conditioning": burn_record,
        "snapped_outlet": (float(ys[row]), float(xs[col])),
        "snap_distance_px": snap_px_moved,
        "flow_accum_at_outlet_km2": round(float(snapped_area), 3),
        "window_deg": half * 2,
        "window_megapixels": round(arr.size / 1e6, 1),
        "estimated_peak_gb": round(_estimated_gb(arr.size), 2),
        "area_km2": round(basin_area_km2(geom), 3),
        "license": "Copernicus DEM free-and-open licence",
    }

delineate_hydrobasins(lat, lon, *, level=12, snap_km=5.0, river_snap_km=1.0, river_snap_ratio=10.0, progress=True, **_)

Delineate the upstream basin from HydroBASINS level-12 units.

Parameters:

Name	Type	Description	Default
level	int	Pfafstetter level of the source file. 12 is finest and the default.	12
snap_km	float	If the outlet is not inside any unit at all, snap to the nearest one within this distance. Coastal outlets often sit just offshore.	5.0
river_snap_km	float	Guard against the bank-of-a-big-river problem: if a unit within river_snap_km drains at least river_snap_ratio times more than the unit containing the point, snap to it and warn. Set river_snap_ratio=None to always take the containing unit.	1.0
river_snap_ratio	float	Guard against the bank-of-a-big-river problem: if a unit within river_snap_km drains at least river_snap_ratio times more than the unit containing the point, snap to it and warn. Set river_snap_ratio=None to always take the containing unit.	1.0

Source code in basinkit/delineate/hydrobasins.py

def delineate_hydrobasins(
    lat: float,
    lon: float,
    *,
    level: int = 12,
    snap_km: float = 5.0,
    river_snap_km: float = 1.0,
    river_snap_ratio: float | None = 10.0,
    progress: bool = True,
    **_,
):
    """Delineate the upstream basin from HydroBASINS level-12 units.

    Parameters
    ----------
    level
        Pfafstetter level of the source file. 12 is finest and the default.
    snap_km
        If the outlet is not inside any unit at all, snap to the nearest one
        within this distance. Coastal outlets often sit just offshore.
    river_snap_km, river_snap_ratio
        Guard against the bank-of-a-big-river problem: if a unit within
        ``river_snap_km`` drains at least ``river_snap_ratio`` times more than
        the unit containing the point, snap to it and warn. Set
        ``river_snap_ratio=None`` to always take the containing unit.
    """
    errors = []
    for region in candidate_regions(lat, lon):
        shp = fetch_region(region, level, progress=progress)
        seed, snap_info = _find_outlet_unit(
            shp, lat, lon, snap_km, river_snap_km, river_snap_ratio or 0
        )
        if seed is None:
            errors.append(region)
            continue

        ids = _upstream_ids(shp, int(seed["HYBAS_ID"]))
        id_list = sorted(ids)

        # Read back only the geometries we need. For a large basin this is
        # still tens of thousands of polygons, so chunk the SQL IN clause --
        # OGR will refuse a single statement with 200k literals.
        units = _read_units(shp, id_list, ["HYBAS_ID", "SUB_AREA"])

        geom = units.union_all()
        # Dissolving thousands of adjacent polygons leaves hairline slivers on
        # shared edges. A tiny buffer out-and-back welds them without moving
        # the outer boundary perceptibly (1e-5 deg is about a metre).
        geom = geom.buffer(1e-5).buffer(-1e-5)
        if geom.geom_type == "GeometryCollection":
            from shapely.geometry import MultiPolygon

            geom = MultiPolygon(
                [g for g in geom.geoms if g.geom_type in ("Polygon", "MultiPolygon")]
            )

        from ..clip import basin_area_km2

        return geom, {
            "backend": "hydrobasins",
            "source_dataset": f"HydroBASINS v1c level {level} ({REGION_NAMES[region]})",
            "region": region,
            "level": level,
            "n_units": len(units),
            "outlet": (lat, lon),
            "outlet_hybas_id": int(seed["HYBAS_ID"]),
            "reported_up_area_km2": float(seed.get("UP_AREA", 0) or 0),
            **_place(seed),
            **snap_info,
            "area_km2": round(basin_area_km2(geom), 2),
            "license": "CC BY 4.0",
            "citation": "Lehner, B. & Grill, G. (2013). Hydrological Processes 27(15), 2171-2186.",
        }

    raise OutletSnapError(
        f"No HydroBASINS level-{level} unit within {snap_km} km of ({lat}, {lon}); "
        f"searched region(s): {', '.join(errors) or 'none'}.\n"
        "Usually this means the coordinate is offshore or in a closed basin. "
        "Try a larger snap_km, move the point onto the river, or use backend='dem'."
    )

delineate_tdx(lat, lon, *, snap_km=2.0, progress=True, **_)

Delineate the upstream basin from TDX-Hydro unit catchments.

Parameters:

Name	Type	Description	Default
snap_km	float	How far the outlet may sit from the nearest mapped reach. The network is dense at 12 m, so the default is tighter than the HydroBASINS one: a point two kilometres from any reach is more likely to be off the river than to be a reach that is missing.	2.0

Source code in basinkit/delineate/tdx.py

def delineate_tdx(lat: float, lon: float, *, snap_km: float = 2.0,
                  progress: bool = True, **_):
    """Delineate the upstream basin from TDX-Hydro unit catchments.

    Parameters
    ----------
    snap_km
        How far the outlet may sit from the nearest mapped reach. The network
        is dense at 12 m, so the default is tighter than the HydroBASINS one:
        a point two kilometres from any reach is more likely to be off the
        river than to be a reach that is missing.
    """
    from ..clip import basin_area_km2

    # The index gives the processing unit reliably -- units are continental --
    # and the catchment polygons give the reach. Distance only decides which
    # unit to open.
    near_link, vpu, snapped_km = _nearest_reach(lat, lon, snap_km, progress)
    link = _reach_containing(lat, lon, vpu, progress)
    chosen_by = "containing catchment"
    if link is None:
        link, chosen_by = near_link, "nearest recorded point"
    wanted = _upstream(link, vpu, progress)

    units = _catchments(vpu, wanted, progress)

    if units.empty:
        raise DelineationError(
            f"Reach {link} has {len(wanted)} upstream reaches but unit {vpu} "
            "holds no catchment polygon for any of them. One GEOGLOWS unit "
            "(302) is known to be missing most of its catchments; use "
            "backend='hydrobasins' there."
        )

    geom = units.union_all()
    # Dissolving adjacent polygons leaves hairline slivers on shared edges; a
    # buffer out and back welds them without moving the outer boundary
    # perceptibly. 1e-6 degrees is about a tenth of a metre, an order finer
    # than the HydroBASINS pass because these edges are ten times sharper.
    geom = geom.buffer(1e-6).buffer(-1e-6)
    if geom.geom_type == "GeometryCollection":
        from shapely.geometry import MultiPolygon

        geom = MultiPolygon(
            [g for g in geom.geoms if g.geom_type in ("Polygon", "MultiPolygon")]
        )

    provenance = {
        "backend": "tdx",
        "source_dataset": "TDX-Hydro (TanDEM-X 12 m), GEOGLOWS v2 hydrography",
        "grid": "12 m TanDEM-X",
        "vpu": vpu,
        "outlet": (lat, lon),
        "outlet_link": link,
        "n_reaches": len(wanted),
        "n_units": len(units),
        "snapped_km": round(snapped_km, 4),
        "outlet_reach_chosen_by": chosen_by,
        "area_km2": round(basin_area_km2(geom), 2),
        "license": LICENSE,
        "citation": CITATION,
    }
    # Every other backend's answer can be weighed against the river network
    # at the outlet; this one's could not, so a reach picked badly by distance
    # came back looking exactly like a reach picked well. The check is the same
    # one `auto` uses, and it is skipped rather than started if HydroRIVERS is
    # not already downloaded, so asking for this backend never silently pulls
    # a few hundred megabytes.
    try:
        from ..verify import check_outlet

        check = check_outlet(lat, lon, provenance["area_km2"],
                             allow_download=False)
        provenance["outlet_check"] = check.as_dict() if hasattr(check, "as_dict") \
            else dict(vars(check))
    except Exception as exc:                                # noqa: BLE001
        provenance["outlet_check"] = {
            "ok": True, "reason": "check-failed", "detail": str(exc)[:200],
        }

    if chosen_by == "nearest recorded point":
        provenance["note"] = (
            "No unit catchment contains this point, so the reach was taken as "
            "the nearest recorded one. That is a weaker choice: the index "
            "holds one point per reach, so the nearest point is not always on "
            "the nearest channel. Check the area against another backend."
        )

    if len(units) < len(wanted):
        provenance["warning"] = (
            f"{len(wanted) - len(units)} of {len(wanted)} upstream reaches have "
            "no catchment polygon in this GEOGLOWS unit, so the basin is "
            "incomplete by that much."
        )
    return geom, provenance

Whether the elevation model supports the terrain analysis asked of it.

Every terrain product in basinkit -- slope, aspect, curvature, the wetness index, height above drainage, the hypsometric integral, every relief parameter in :func:~basinkit.morphometry.morphometry -- is computed from one elevation raster. Whether those numbers mean anything depends on whether the raster can resolve the terrain in this particular basin, which is a different question in a Himalayan gorge and on a coastal plain.

The accuracy tables published with basinkit answer the first half of that: across thousands of gauges, at this catchment size, this is the error to expect. That is a population figure. It cannot say whether this basin is one of the good cases. This module answers the second half by measuring the basin in front of it.

Five tests, each a measurement against a stated threshold:

filling How much of the surface depression filling had to invent. A basin where a fifth of the cells were raised is one where the drainage network is partly the conditioning algorithm's own work.

slope The noise floor of the model, arctan(vertical error / cell size), and how many of the basin's slopes fall below it. Below that angle a slope value is the model's error expressed as a gradient.

relief Total relief against the model's vertical error. A 40 m basin measured with a 4 m error carries a tenth of its own signal as uncertainty.

water The largest level surface inside the basin. Standing water has no gradient, so slope, aspect, curvature and the wetness index over it are not measurements of anything.

coverage Cells inside the basin with no elevation at all. Global models have withheld tiles and voids, and a missing fifth of a basin biases every statistic computed over it.

The grade is HIGH when every test passes, MODERATE when one is close to its threshold, and LIMITED when one is past it. It is deliberately about the terrain products and not about the basin boundary: the boundary comes from the delineation backend, which is measured separately in the verification page.

suitability(basin, *, dem=None, product=None, support_map=False)

Grade how well the elevation model supports terrain analysis here.

Parameters:

Name	Type	Description	Default
basin		A :class:~basinkit.basin.Basin, or any elevation array carrying its own georeferencing.	required
dem		An already-fetched elevation raster, to avoid downloading it twice.	None
product	str | None	The elevation product, used to pick the vertical error. Taken from the raster's own attributes when it is not given, and assumed to be the least accurate of the four when it cannot be determined.	None
support_map	bool	Also return support: 1 where a cell was neither raised by filling nor below the slope noise floor, 0 where it was, and nodata outside the basin.	False

Returns:

Type	Description
dict	grade, the five tests with their measurements and thresholds, unmet and advisory lists, statement, and the constants the thresholds were applied against.

Source code in basinkit/suitability.py

def suitability(basin, *, dem=None, product: str | None = None,
                support_map: bool = False) -> dict[str, Any]:
    """Grade how well the elevation model supports terrain analysis here.

    Parameters
    ----------
    basin
        A :class:`~basinkit.basin.Basin`, or any elevation array carrying its
        own georeferencing.
    dem
        An already-fetched elevation raster, to avoid downloading it twice.
    product
        The elevation product, used to pick the vertical error. Taken from the
        raster's own attributes when it is not given, and assumed to be the
        least accurate of the four when it cannot be determined.
    support_map
        Also return ``support``: 1 where a cell was neither raised by filling
        nor below the slope noise floor, 0 where it was, and nodata outside
        the basin.

    Returns
    -------
    dict
        ``grade``, the five ``tests`` with their measurements and thresholds,
        ``unmet`` and ``advisory`` lists, ``statement``, and the constants the
        thresholds were applied against.
    """
    elevation = dem
    if elevation is None:
        elevation = basin.dem() if hasattr(basin, "dem") else basin
    values = _as_array(elevation)

    if product is None:
        product = getattr(elevation, "attrs", {}).get("basinkit_product")
    key = str(product).lower() if product else ""
    error_m = VERTICAL_ERROR_M.get(key, DEFAULT_VERTICAL_ERROR_M)

    dy_m, dx_m = _spacing(elevation, values)
    cell_m = float(np.mean([dy_m, float(np.mean(dx_m))]))
    cell_km2 = float(np.mean(dx_m)) * dy_m / 1e6

    inside = np.isfinite(values)
    known = int(inside.sum())
    polygon = _inside_polygon(basin, elevation, values)
    covered = polygon if polygon is not None else np.ones(values.shape, dtype=bool)
    total_cells = int(covered.sum())

    tests: list[dict[str, Any]] = []
    unmet: list[str] = []
    advisory: list[str] = []

    def record(name, verdict, measured, unit, advisory_at, unmet_at, statement,
               evaluated=True, **extra):
        entry = {"test": name, "verdict": verdict, "measured": measured,
                 "unit": unit, "advisory_at": advisory_at, "unmet_at": unmet_at,
                 "evaluated": evaluated, "statement": statement, **extra}
        tests.append(entry)
        if verdict == "unmet":
            unmet.append(name)
        elif verdict == "advisory":
            advisory.append(name)

    # --- filling -----------------------------------------------------------
    depth = _fill_depth(elevation, values)
    raised_frac = None
    if depth is None:
        record("filling", "pass", None, "percent of cells", 5.0, 20.0,
               "Not evaluated: depression filling needs pyflwdir, which is in "
               'the delineate extra. Install it with pip install '
               '"basinkit[delineate]" to include this test.',
               evaluated=False)
    else:
        d = depth[np.isfinite(depth)]
        raised = (d > 0.01)
        raised_pct = float(raised.sum() / d.size * 100) if d.size else 0.0
        max_fill = float(d.max()) if d.size else 0.0
        raised_frac = raised_pct
        record("filling", _verdict(raised_pct, 5.0, 20.0), round(raised_pct, 2),
               "percent of cells", 5.0, 20.0,
               f"Depression filling raised {raised_pct:.1f}% of cells, the "
               f"deepest by {max_fill:.1f} m. Drainage across raised ground is "
               "the conditioning algorithm's construction, not mapped relief.",
               max_fill_m=round(max_fill, 2))

    # --- slope -------------------------------------------------------------
    noise_floor_deg = float(math.degrees(math.atan(error_m / cell_m)))
    slope_deg = np.asarray(getattr(_slope(elevation), "values", None), dtype="float64")
    s = slope_deg[np.isfinite(slope_deg)]
    below_pct = float((s < noise_floor_deg).sum() / s.size * 100) if s.size else 0.0
    record("slope", _verdict(below_pct, 25.0, 60.0), round(below_pct, 2),
           "percent of cells", 25.0, 60.0,
           f"At {error_m:.1f} m vertical error over a {cell_m:.0f} m cell the "
           f"noise floor is {noise_floor_deg:.1f} degrees, and {below_pct:.0f}% "
           "of this basin's slopes fall below it. Those values are the model's "
           "own error expressed as a gradient.",
           noise_floor_deg=round(noise_floor_deg, 3))

    # --- relief ------------------------------------------------------------
    z = values[inside]
    relief_m = float(z.max() - z.min()) if z.size else 0.0
    ratio = relief_m / error_m if error_m else float("inf")
    record("relief", _verdict(ratio, 100.0, 20.0, lower_is_better=False),
           round(ratio, 1), "relief divided by vertical error", 100.0, 20.0,
           f"Relief is {relief_m:.0f} m against a vertical error of "
           f"{error_m:.1f} m, a ratio of {ratio:.0f}. Every elevation-derived "
           "parameter carries that ratio as its relative uncertainty.",
           relief_m=round(relief_m, 1))

    # --- water -------------------------------------------------------------
    body_km2, body_elev = _largest_level_surface(values, cell_km2)
    basin_km2 = known * cell_km2
    body_pct = (body_km2 / basin_km2 * 100) if basin_km2 else 0.0
    # Judged as a share of the basin, not as an absolute area. A 24 km2 lake
    # is most of a small catchment and a rounding error in the Koshi, and an
    # absolute threshold grades the second one as though it were the first.
    if body_km2 >= _MIN_BODY_KM2:
        water_verdict = _verdict(body_pct, 1.0, 5.0)
        water_statement = (
            f"A level surface of {body_km2:.2f} km2 at {body_elev:.0f} m covers "
            f"{body_pct:.2f}% of the basin: standing water, or permanent snow "
            "and ice, which the sensor returns the same way. Slope, aspect, "
            "curvature and the wetness index over it describe a flat return "
            "rather than ground. The basin boundary is unaffected, since "
            "delineation here does not route on this raster."
        )
    else:
        water_verdict = "pass"
        water_statement = (
            "No level surface large enough to distort the terrain statistics."
        )
    record("water", water_verdict, round(body_pct, 3),
           "percent of basin", 1.0, 5.0, water_statement,
           surface_km2=round(body_km2, 3),
           surface_elev_m=None if body_elev is None else round(body_elev, 1))

    # --- coverage ----------------------------------------------------------
    absent = int((covered & ~inside).sum())
    missing_pct = float(absent / total_cells * 100) if total_cells else 0.0
    record("coverage", _verdict(missing_pct, 0.5, 2.0), round(missing_pct, 3),
           "percent of cells", 0.5, 2.0,
           f"{missing_pct:.2f}% of the cells inside this basin carry no "
           "elevation. Global models have withheld tiles and unfilled voids, "
           "and every statistic computed over the basin is biased by whatever "
           "is missing.",
           measured_inside="polygon" if polygon is not None else "whole raster")

    grade = "LIMITED" if unmet else ("MODERATE" if advisory else "HIGH")

    out: dict[str, Any] = {
        "grade": grade,
        "product": product or "unknown",
        "vertical_error_m": error_m,
        "cell_size_m": round(cell_m, 2),
        "basin_km2": round(basin_km2, 3),
        "tests": tests,
        "unmet": unmet,
        "advisory": advisory,
        "statement": _statement(grade, tests),
        "note": "A grade for the terrain products computed from this raster. "
                "The basin boundary is delineated separately and its accuracy "
                "is reported in the verification page, not here.",
    }

    if support_map:
        supported = np.ones(values.shape, dtype="float64")
        if depth is not None:
            supported[np.nan_to_num(depth, nan=0.0) > 0.01] = 0.0
        supported[np.nan_to_num(slope_deg, nan=noise_floor_deg) < noise_floor_deg] = 0.0
        supported[~inside] = np.nan
        out["support"] = _wrap_support(elevation, supported)
        out["support_fraction"] = (
            round(float(np.nansum(supported) / known), 4) if known else None
        )
    elif raised_frac is not None:
        out["support_fraction"] = None

    return out

Whether this basin is one the accuracy figures were measured on.

Every accuracy number published with basinkit is a population figure. The blind validation drew 2,740 delineations at 2,550 gauges in 99 countries and reported the median area error per size band, and that table is the honest answer to "how accurate is it" -- for the population. It is not an answer about the basin in front of you.

The difference is not academic. The headline figures, the ones worth quoting, belong to large catchments: 0.3% above 100,000 km2, 1.3% from 10,000 to 100,000, 1.7% from 2,000 to 10,000. Below 500 km2 the same validation measured 30%, and below 100 km2 it measured 181%. A basin of 60 km2 is inside the validated population and nowhere near the quoted accuracy, and until this module existed nothing in the output said so. The user saw an area, a boundary and a verdict on the outlet, and had no way to know that at this size the measurement behind the software is four figures worse than the one on the front page.

So this module answers a narrower question than "is it accurate". It answers: of the basins that were measured, which ones does this one resemble, and what was the error there -- and, where this basin resembles nothing that was measured, it says that instead of quoting a number that does not apply.

Three kinds of answer come out of it:

inside This basin sits where the validation measured a median area error at or below :data:INSIDE_MEDIAN_PCT. The published figure applies to it.

edge Measured, but where the error was several times the headline figure, or with a descriptor out in the tail of what was validated. The number applies, and it is not the number most people have in mind.

outside Either the measured error here was above :data:EDGE_MEDIAN_PCT, or a descriptor falls beyond anything in the validation set. No published figure describes this case, and the envelope says so rather than interpolating one.

It is deliberately about the basin boundary, which is what the validation measured. :func:~basinkit.suitability.suitability grades the elevation model's support for the terrain products, which is the other half and a different question. The two share their measurements: the four terrain axes below are the same four numbers suitability already computes, compared against the range they covered in the validation rather than against a threshold, so passing a suitability result in costs nothing.

"Envelope" here is the statistical sense -- the region of descriptor space the measurements cover -- and not the geometric one. A basin's bounding box is :attr:~basinkit.basin.Basin.bounds. The function is called :func:validation_envelope so that the two cannot be confused at a call site.

Not everything is characterised. Aridity, flow regulation and karst are each known to matter and none of them was recorded in the validation set, so this module names them as uncharacterised instead of implying coverage it does not have. :data:NOT_CHARACTERISED carries the list and how to measure each one.

validation_envelope(basin=None, *, area_km2=None, suitability=None)

Say which part of the validation set this basin belongs to.

Parameters:

Name	Type	Description	Default
basin		A :class:~basinkit.basin.Basin. Only its area is read, so this downloads nothing and costs nothing.	None
area_km2	float | None	The basin area, when there is no Basin to hand. Given alongside basin it wins, which is what lets a caller ask about a size it has not delineated yet.	None
suitability	dict[str, Any] | None	A :func:~basinkit.suitability.suitability result. Its five measurements are compared against the range they covered in the terrain study, adding those axes to the answer. Without it the envelope rests on size alone, which is the axis the full validation measured.	None

Returns:

Type	Description
dict	position (inside, edge or outside), basin_km2, size with the band and the error measured in it, axes with one entry per descriptor compared, not_characterised, statement and source.

Examples:

>>> validation_envelope(area_km2=61.0)["position"]
'outside'
>>> validation_envelope(area_km2=61.0)["size"]["median_error_pct"]
181.0

Source code in basinkit/envelope.py

def validation_envelope(basin=None, *, area_km2: float | None = None,
             suitability: dict[str, Any] | None = None) -> dict[str, Any]:
    """Say which part of the validation set this basin belongs to.

    Parameters
    ----------
    basin
        A :class:`~basinkit.basin.Basin`. Only its area is read, so this
        downloads nothing and costs nothing.
    area_km2
        The basin area, when there is no ``Basin`` to hand. Given alongside
        ``basin`` it wins, which is what lets a caller ask about a size it has
        not delineated yet.
    suitability
        A :func:`~basinkit.suitability.suitability` result. Its five
        measurements are compared against the range they covered in the
        terrain study, adding those axes to the answer. Without it the
        envelope rests on size alone, which is the axis the full validation
        measured.

    Returns
    -------
    dict
        ``position`` (``inside``, ``edge`` or ``outside``), ``basin_km2``,
        ``size`` with the band and the error measured in it, ``axes`` with one
        entry per descriptor compared, ``not_characterised``,
        ``statement`` and ``source``.

    Examples
    --------
    >>> validation_envelope(area_km2=61.0)["position"]
    'outside'
    >>> validation_envelope(area_km2=61.0)["size"]["median_error_pct"]
    181.0
    """
    if area_km2 is None:
        if basin is None:
            raise TypeError("validation_envelope() needs a basin or an area_km2.")
        area_km2 = float(basin.area_km2)
    area_km2 = float(area_km2)

    band = _band(area_km2)
    size_position = _size_position(band["median_error_pct"])
    axes: list[dict[str, Any]] = [{
        "axis": "size",
        "measured": round(area_km2, 3),
        "unit": "km2",
        "position": size_position,
        "validated_from": (
            f"{VALIDATION_DELINEATIONS:,} delineations at "
            f"{VALIDATION_GAUGES:,} gauges in {VALIDATION_COUNTRIES} countries"
        ),
        "validated_range": band["band"] + " km2",
        "statement": (
            f"At {area_km2:,.1f} km2 this basin falls in the {band['band']} km2 "
            f"band, where the blind validation measured a median area error of "
            f"{band['median_error_pct']:g}%, with {band['within_5_pct']}% of "
            f"delineations within 5% and {band['within_20_pct']}% within 20%."
        ),
    }]

    comparable = TERRAIN_STUDY_KM2[0] <= area_km2 <= TERRAIN_STUDY_KM2[1]
    if suitability:
        for test in suitability.get("tests", ()):
            ref = TERRAIN_STUDY.get(test.get("test"))
            measured = test.get("measured")
            if ref is None or measured is None or not test.get("evaluated", True):
                continue
            measured = float(measured)
            position = _terrain_position(measured, ref)
            axes.append({
                "axis": test["test"],
                "measured": measured,
                "unit": ref["unit"],
                "position": position,
                "validated_from": (
                    f"{TERRAIN_STUDY_GAUGES} of those gauges, with the "
                    f"suitability measurements recorded"
                ),
                "validated_range": f"{ref['min']:g} to {ref['max']:g}",
                "statement": _terrain_statement(test["test"], measured, ref,
                                                position, comparable),
            })

    position = max((a["position"] for a in axes), key=lambda p: _ORDER[p])
    out = {
        "position": position,
        "basin_km2": round(area_km2, 3),
        "size": {
            "band": band["band"] + " km2",
            "median_error_pct": band["median_error_pct"],
            "within_5_pct": band["within_5_pct"],
            "within_20_pct": band["within_20_pct"],
        },
        "axes": axes,
        "outside": [a["axis"] for a in axes if a["position"] == "outside"],
        "edge": [a["axis"] for a in axes if a["position"] == "edge"],
        "not_characterised": [dict(entry) for entry in NOT_CHARACTERISED],
        "statement": _statement(position, axes, suitability is not None),
        "source": "docs/verification.md, from "
                  "verify/benchmarks/accuracy_by_size.csv and "
                  "verify/benchmarks/grade_validation.csv",
        "note": "Where this basin sits in the set the accuracy figures were "
                "measured on. It says nothing about whether this particular "
                "boundary is right: check_outlet tests that against the river "
                "network, and suitability grades the terrain products.",
    }
    return out

A quality indicator for every layer, not only for the elevation model.

:func:~basinkit.suitability.suitability grades the elevation raster, because every terrain product is derived from it. The other layers arrive with no such statement: a land cover map, a soil raster and a rainfall series all come back looking equally authoritative, and nothing in the output says which of them can carry the weight a study puts on it.

This module measures each layer against something independent of it and reports what it found:

elevation The suitability grade, and the same measurement split by terrain class, so that a basin which is part plain and part mountain is not reported as one number. The classes are cut on local relief, which is the quantity the model's vertical error competes with. Both the relief and the noise floor are per-cell quantities, so the split is reported together with the cell size it was measured at.

landcover Two independently produced maps of the same year -- ESA WorldCover and the ESRI annual series -- reduced to the classes they both define, and the share of the basin where they agree. Where two mapping teams disagree, the class under a given pixel is not settled.

soil SoilGrids publishes a 5th and a 95th percentile beside its mean. The width of that interval, in the units of the property, is the model's own statement about how well it knows this ground.

precipitation CHIRPS against TerraClimate over the same years: two products built from different inputs. Their correlation and the difference between their means bound how much of a rainfall figure is the choice of product.

surface_water How much of the water in the basin is permanent, how much seasonal, and how much of the basin the Landsat record never observed cleanly.

delineation The boundary is not measured here -- it comes from the backend -- so this reports the checks that do exist: the river-network consistency check and the accuracy regime for a basin of this size.

Every threshold below is a choice, stated in the output, not a measurement. The measurements are the numbers beside them.

data_quality(basin, *, layers=('elevation', 'landcover', 'soil', 'precipitation', 'surface_water', 'delineation'), dem=None, year=2021, progress=True)

A grade for every layer this basin can return, each against something independent.

Every entry carries the measurement, the threshold it was judged against and a sentence saying what was compared with what. A layer that could not be fetched is reported as an error rather than dropped, so the absence is visible.

Source code in basinkit/quality.py

def data_quality(basin, *, layers: tuple[str, ...] = ("elevation", "landcover", "soil",
                                                      "precipitation", "surface_water",
                                                      "delineation"),
                 dem=None, year: int = 2021, progress: bool = True) -> dict[str, Any]:
    """A grade for every layer this basin can return, each against something independent.

    Every entry carries the measurement, the threshold it was judged against and
    a sentence saying what was compared with what. A layer that could not be
    fetched is reported as an error rather than dropped, so the absence is
    visible.
    """
    runners = {
        "elevation": lambda: elevation_quality(basin, dem=dem, progress=progress),
        "landcover": lambda: landcover_quality(basin, year=year, progress=progress),
        "soil": lambda: soil_quality(basin, progress=progress),
        "precipitation": lambda: precipitation_quality(basin, progress=progress),
        "surface_water": lambda: surface_water_quality(basin, progress=progress),
        "delineation": lambda: delineation_quality(basin),
    }
    unknown = [name for name in layers if name not in runners]
    if unknown:
        raise ValueError(f"Unknown layer(s) {unknown}: choose from {sorted(runners)}")

    rows = []
    for name in layers:
        try:
            rows.append(runners[name]())
        except Exception as exc:
            rows.append({"layer": name, "grade": None, "value": None,
                         "indicator": None, "error": f"{type(exc).__name__}: {exc}",
                         "statement": f"{name} could not be checked: {exc}"})
    graded = [r["grade"] for r in rows if r.get("grade")]
    return {"basin_km2": round(float(basin.area_km2), 2),
            "overall": _worst(graded) if graded else None,
            "layers": rows,
            "note": ("Each layer is judged against something produced independently of it. "
                     "The thresholds are choices and are printed with every result; the "
                     "measurements beside them are what was found in this basin.")}

delineation_quality(basin)

What is actually known about the boundary: the checks that ran.

Source code in basinkit/quality.py

def delineation_quality(basin) -> dict[str, Any]:
    """What is actually known about the boundary: the checks that ran."""
    prov = dict(basin.provenance or {})
    check = prov.get("outlet_check") or {}
    area = float(basin.area_km2)
    regimes = ((100_000, 0.3), (10_000, 1.3), (2_000, 1.7))
    expected = next((err for floor, err in regimes if area >= floor), None)
    if expected is None:
        grade = "MODERATE"
        band = ("Below about 2,000 km2 the default backend's median error grows and the DEM "
                "backend is the better instrument.")
    else:
        grade = "HIGH"
        band = (f"At {area:,.0f} km2 the blind validation across 2,550 gauges in 99 countries "
                f"puts the median area error at {expected}%.")
    if check and check.get("ok") is False:
        grade = "LIMITED"
    statement = band
    if check:
        statement += (" The river-network check "
                      + ("agrees with the result" if check.get("ok") else
                         f"disagrees: {check.get('reason', 'see provenance')}") + ".")
    return {"layer": "delineation", "indicator": "network check and size regime",
            "grade": grade, "value": area, "backend": prov.get("backend"),
            "network_check": check or None, "statement": statement,
            "source": "basinkit validation table; the boundary itself is not remeasured here"}

elevation_quality(basin, dem=None, **kwargs)

The suitability grade, plus the noise-floor test per terrain class.

Source code in basinkit/quality.py

def elevation_quality(basin, dem=None, **kwargs) -> dict[str, Any]:
    """The suitability grade, plus the noise-floor test per terrain class."""
    from .suitability import _range_3x3, suitability
    from .terrain import _as_array
    from .terrain import slope as _slope

    if dem is None:
        kwargs.setdefault("max_pixels", QUALITY_DEM_PIXELS)
        dem = basin.dem(**kwargs)
    graded = suitability(basin, dem=dem)
    values = _as_array(dem)
    step = max(1, int(np.ceil(np.sqrt(values.size / CLASS_SAMPLE_PIXELS))))
    sampled = dem[::step, ::step] if step > 1 else dem
    relief = _range_3x3(_as_array(sampled))
    slopes = np.asarray(_slope(sampled).values, dtype="float64")
    floor = next((t for t in graded["tests"] if t["test"] == "slope"), {}).get("noise_floor_deg")
    cell_m = float(graded["cell_size_m"]) * step
    if step > 1 and graded.get("vertical_error_m"):
        # The floor is atan(vertical error / cell size): on a decimated copy the
        # cells are wider, so the angle the model cannot resolve is smaller.
        floor = round(_noise_floor_deg(float(graded["vertical_error_m"]), cell_m), 3)

    classes = []
    usable = np.isfinite(relief) & np.isfinite(slopes)
    total = int(usable.sum())
    for low, high, name in RELIEF_CLASSES:
        picked = usable & (relief >= low) & (relief < high)
        n = int(picked.sum())
        if n == 0:
            continue
        row = {"terrain": name, "share_of_basin": round(n / total, 4),
               "mean_slope_deg": round(float(np.mean(slopes[picked])), 2)}
        if floor is not None:
            row["below_noise_floor"] = round(float(np.mean(slopes[picked] < floor)), 4)
        classes.append(row)

    statement = graded["statement"]
    if floor is not None and classes:
        worst = max(classes, key=lambda c: c.get("below_noise_floor", 0.0))
        statement += (
            f" Split by terrain, {worst['share_of_basin']:.0%} of the basin is "
            f"{worst['terrain']} ground, where {worst['below_noise_floor']:.0%} of "
            f"slopes fall below the {floor:.1f} degree noise floor of a "
            f"{cell_m:.0f} m cell."
        )
    return {"layer": "elevation", "indicator": "suitability grade, by terrain class",
            "class_sample_step": step, "class_cell_size_m": round(cell_m, 2),
            "class_noise_floor_deg": floor,
            "grade": graded["grade"], "value": graded["support_fraction"],
            "terrain_classes": classes, "tests": graded["tests"],
            "statement": statement,
            "source": f"{graded['product']}, cell {graded['cell_size_m']:.0f} m"}

landcover_quality(basin, year=2021, max_pixels=None, **kwargs)

Agreement between ESA WorldCover and the ESRI annual map for one year.

Source code in basinkit/quality.py

def landcover_quality(basin, year: int = 2021, max_pixels: int | None = None,
                      **kwargs) -> dict[str, Any]:
    """Agreement between ESA WorldCover and the ESRI annual map for one year."""
    budget = LANDCOVER_PIXELS if max_pixels is None else max_pixels
    first = basin.landcover(year=year, source="worldcover", max_pixels=budget, **kwargs)
    # The ESRI reader takes no progress flag, and reads through a different
    # path, so its budget is passed on its own terms.
    second = basin.landcover(year=year, source="esri", max_pixels=budget)
    if second.shape != first.shape:
        second = second.rio.reproject_match(first)
    a = np.asarray(first.values, dtype="float64").squeeze()
    b = np.asarray(second.values, dtype="float64").squeeze()

    common = np.zeros(a.shape, dtype="int16")
    other = np.zeros(a.shape, dtype="int16")
    for i, (_name, (wc_codes, esri_codes)) in enumerate(_COMMON.items(), start=1):
        common[np.isin(a, list(wc_codes))] = i
        other[np.isin(b, list(esri_codes))] = i
    usable = (common > 0) & (other > 0)
    n = int(usable.sum())
    if n == 0:
        return {"layer": "landcover", "grade": "LIMITED", "value": None,
                "indicator": "agreement between two independent maps",
                "statement": "No cell fell in a class both maps define, so no comparison was possible.",
                "source": f"ESA WorldCover {year} against ESRI {year}"}
    agree = float(np.mean(common[usable] == other[usable]))
    per_class, counts = {}, {}
    for i, name in enumerate(_COMMON, start=1):
        picked = usable & (common == i)
        if picked.any():
            per_class[name] = round(float(np.mean(other[picked] == i)), 4)
            counts[name] = int(picked.sum())
    good, bad = THRESHOLDS["landcover_agreement"]
    # A class covering a handful of cells can disagree completely without
    # saying anything about the map, so the weakest class is only named when
    # it holds at least one per cent of the ground compared.
    substantial = {k: v for k, v in per_class.items() if counts[k] >= 0.01 * n}
    weakest = min(substantial, key=substantial.get) if substantial else None
    statement = (f"Two independently produced maps of {year} agree on "
                 f"{agree:.0%} of the basin, over the classes both of them define "
                 f"({n:,} cells).")
    if weakest is not None:
        statement += (f" Of the classes covering more than one per cent of that ground, they "
                      f"agree least on {weakest.replace('_', ' ')}, at {per_class[weakest]:.0%}.")
    return {"layer": "landcover", "indicator": "agreement between two independent maps",
            "grade": _grade(agree, good, bad, higher_is_better=True), "value": round(agree, 4),
            "per_class_agreement": per_class, "class_cells": counts, "compared_cells": n,
            "threshold": {"high_at_or_above": good, "limited_below": bad},
            "statement": statement,
            "source": f"ESA WorldCover {year} against ESRI Annual LULC {year}"}

precipitation_quality(basin, start=2000, end=None, **kwargs)

CHIRPS against TerraClimate: two products, the same basin and years.

Source code in basinkit/quality.py

def precipitation_quality(basin, start: int = 2000, end: int | None = None, **kwargs) -> dict[str, Any]:
    """CHIRPS against TerraClimate: two products, the same basin and years."""
    import pandas as pd

    chirps = basin.precipitation(start=start, end=end, source="chirps", **kwargs)
    terra = basin.precipitation(start=start, end=end, source="terraclimate", **kwargs)
    def _annual(obj, column=None):
        # A Dataset is picked from by name; anything with to_series() becomes
        # one directly. The old code converted to a DataFrame and fell back to
        # column zero when the name it wanted was absent -- and on a
        # DataArray that column is `spatial_ref`, a constant. Correlating
        # against a constant is NaN, which was then graded LIMITED. Never
        # guess a column: a series that cannot be identified is an error.
        if hasattr(obj, "data_vars"):
            name = column if column in obj.data_vars else next(iter(obj.data_vars))
            obj = obj[name]
        s = obj.to_series() if hasattr(obj, "to_series") else obj
        if isinstance(s, pd.DataFrame):
            if column and column in s:
                s = s[column]
            else:
                usable = [c for c in s.columns
                          if c not in ("spatial_ref", "crs", "band")]
                if not usable:
                    raise ValueError(
                        "the rainfall series carries no data column, only "
                        f"{list(s.columns)}")
                s = s[usable[0]]
        s = pd.Series(s)
        index = pd.to_datetime(pd.Index(s.index).get_level_values(-1))
        s = pd.Series(np.asarray(s, dtype="float64"), index=index)
        s = s[~s.index.duplicated()]
        return s.groupby(s.index.year).sum()
    a = _annual(chirps)
    b = _annual(terra, column="ppt")
    years = sorted(set(a.index) & set(b.index))
    years = [y for y in years if y != max(years)] or years   # drop a partial final year
    a, b = a.loc[years], b.loc[years]
    if len(years) < 3:
        return {"layer": "precipitation", "grade": None, "value": None,
                "years": len(years),
                "not_graded_because":
                    (f"only {len(years)} year(s) are covered by both products, and a "
                     "correlation over that is not a measurement."),
                "statement": "Too few overlapping years to compare the two products.",
                "source": "CHIRPS v3.0 against TerraClimate"}
    r = float(np.corrcoef(a.values, b.values)[0, 1])
    bias = float(a.mean() / b.mean() - 1.0) * 100.0
    good, bad = THRESHOLDS["precipitation_r"]
    if r != r:
        return {"layer": "precipitation", "grade": None, "value": None,
                "years": len(years),
                "chirps_mm_per_year": round(float(a.mean()), 1),
                "terraclimate_mm_per_year": round(float(b.mean()), 1),
                "not_graded_because":
                    ("the correlation came out undefined, which happens when one "
                     "of the two series has no variance over these years. The "
                     "two means are reported so the comparison is not lost."),
                "statement": "The two products could not be correlated over these years.",
                "source": "CHIRPS v3.0 against TerraClimate"}
    return {"layer": "precipitation", "indicator": "two independent products, annual totals",
            "grade": _grade(r, good, bad, higher_is_better=True),
            "value": round(r, 3), "bias_pct": round(bias, 1), "years": len(years),
            "chirps_mm_per_year": round(float(a.mean()), 1),
            "terraclimate_mm_per_year": round(float(b.mean()), 1),
            "threshold": {"high_at_or_above_r": good, "limited_below_r": bad},
            "statement": (f"Over {len(years)} years CHIRPS and TerraClimate agree year to year at "
                          f"r = {r:.2f}, with CHIRPS {bias:+.0f}% against TerraClimate "
                          f"({a.mean():,.0f} against {b.mean():,.0f} mm a year). The gap between "
                          "them is the part of a rainfall figure that is the choice of product."),
            "source": "CHIRPS v3.0 against TerraClimate"}

soil_quality(basin, prop='clay', depth='0-5cm', **kwargs)

SoilGrids' own 90 per cent interval, as a share of its mean.

Source code in basinkit/quality.py

def soil_quality(basin, prop: str = "clay", depth: str = "0-5cm", **kwargs) -> dict[str, Any]:
    """SoilGrids' own 90 per cent interval, as a share of its mean."""
    mean = basin.soil(prop, depth, stat="mean", **kwargs)
    low = basin.soil(prop, depth, stat="Q0.05", **kwargs)
    high = basin.soil(prop, depth, stat="Q0.95", **kwargs)
    m = np.asarray(mean.values, dtype="float64").squeeze()
    lo = np.asarray(low.values, dtype="float64").squeeze()
    hi = np.asarray(high.values, dtype="float64").squeeze()
    ok = np.isfinite(m) & np.isfinite(lo) & np.isfinite(hi) & (m > 0)
    if not ok.any():
        return {"layer": "soil", "grade": "LIMITED", "value": None,
                "indicator": "SoilGrids 90% interval width",
                "statement": "SoilGrids returned no usable values for this basin.",
                "source": f"SoilGrids {prop} {depth}"}
    width = float(np.median(hi[ok] - lo[ok]))
    relative = float(np.median((hi[ok] - lo[ok]) / m[ok])) * 100.0
    factor = 10.0 if prop in ("clay", "sand", "silt") else 1.0   # SoilGrids maps g/kg for these
    return {"layer": "soil", "indicator": "SoilGrids 90% interval width, relative to its mean",
            "grade": None, "not_graded_because":
                ("SoilGrids' published quantiles are wide in every basin tested, so any "
                 "pass/fail threshold here would be invented. The width is reported so that "
                 "basins, depths and properties can be compared with each other."),
            "value": round(relative, 1),
            "interval_width_percent_points": round(width / factor, 1),
            "median_mean_percent": round(float(np.median(m[ok])) / factor, 1),
            "statement": (f"At the median cell SoilGrids reports {np.median(m[ok])/factor:.0f}% {prop} "
                          f"with a 5th to 95th percentile range {width/factor:,.0f} points wide, "
                          f"{relative:.0f}% of that mean. Read the map for pattern, not for a field."),
            "source": f"SoilGrids {prop} {depth}, mean against its own Q0.05 and Q0.95"}

surface_water_quality(basin, **kwargs)

Permanent against seasonal water, and how much was never observed.

Source code in basinkit/quality.py

def surface_water_quality(basin, **kwargs) -> dict[str, Any]:
    """Permanent against seasonal water, and how much was never observed."""
    occurrence = basin.surface_water(**kwargs)
    v = np.asarray(occurrence.values, dtype="float64").squeeze()
    inside = np.isfinite(v)
    total = int(inside.sum())
    if total == 0:
        return {"layer": "surface_water", "grade": "LIMITED", "value": None,
                "indicator": "observed share of the basin",
                "statement": "The surface water record covers none of this basin.",
                "source": "JRC Global Surface Water occurrence"}
    permanent = float(np.mean(v[inside] >= 75))
    seasonal = float(np.mean((v[inside] > 5) & (v[inside] < 75)))
    return {"layer": "surface_water", "indicator": "permanent against seasonal water",
            "grade": None, "not_graded_because":
                ("The JRC record carries no per-pixel confidence, so there is nothing here to "
                 "grade against. The split between permanent and seasonal water is reported "
                 "instead, because a basin whose water is mostly seasonal is one where a single "
                 "date would have misled."),
            "value": round(permanent, 5), "seasonal": round(seasonal, 5),
            "statement": (f"{permanent:.2%} of the basin holds water more than three months in four "
                          f"and {seasonal:.2%} holds it seasonally, over 1984 to 2021. "
                          "Water hidden under cloud or canopy is not in the record, so the seasonal "
                          "figure is a floor."),
            "source": "JRC Global Surface Water occurrence, 1984 to 2021"}

An eight-page report a thesis or a manuscript can use directly.

Everything basinkit computes about a basin, laid out on A4 with the numbers carrying their symbols, units and original references, and with the suitability grade on the cover rather than in a footnote. The point is that a reader can see what the result rests on without being handed the code.

The pages, in order:

1. Cover: where, how big, which backend, and the suitability grade
2. Elevation, hillshade, the hypsometric curve and the elevation distribution
3. Slope and aspect, with their distributions
4. Curvature, position index, ruggedness and the slope-position classes
5. The channel network and Horton's laws
6. The morphometric parameters, with symbols and references
7. The suitability scorecard, the support map and the drainage-density curve
8. Methods: sources, licences, versions and what to cite

Pages 5 and 6 need the river network, which is a separate download. Without it they are replaced by a page saying so rather than quietly omitted.

estimated_peak_gb(n_pixels)

What the report will need, in gigabytes, for a grid of this size.

Source code in basinkit/report.py

def estimated_peak_gb(n_pixels: int) -> float:
    """What the report will need, in gigabytes, for a grid of this size."""
    return FIXED_OVERHEAD_GB + n_pixels * BYTES_PER_PIXEL / 1e9

report(basin, path, *, title=None, dem=None, rivers=None, morphometry=None, max_pixels=REPORT_MAX_PIXELS, progress=True)

Write the eight-page report for basin to path.

Parameters:

Name	Type	Description	Default
path		Where to write the PDF.	required
title	str | None	Printed on the cover. Defaults to the basin's coordinates.	None
dem		Already-computed layers, to avoid fetching or recomputing them. A grid passed in is used as it stands, budget included -- the caller chose it.	None
rivers		Already-computed layers, to avoid fetching or recomputing them. A grid passed in is used as it stands, budget included -- the caller chose it.	None
morphometry		Already-computed layers, to avoid fetching or recomputing them. A grid passed in is used as it stands, budget included -- the caller chose it.	None
max_pixels	int | None	Pixel budget for the elevation grid this report is built from. The cover names the resolution that budget worked out to.	REPORT_MAX_PIXELS

Returns:

Type	Description
str	The path written.

Source code in basinkit/report.py

def report(basin, path, *, title: str | None = None, dem=None, rivers=None,
           morphometry=None, max_pixels: int | None = REPORT_MAX_PIXELS,
           progress: bool = True) -> str:
    """Write the eight-page report for ``basin`` to ``path``.

    Parameters
    ----------
    path
        Where to write the PDF.
    title
        Printed on the cover. Defaults to the basin's coordinates.
    dem, rivers, morphometry
        Already-computed layers, to avoid fetching or recomputing them. A grid
        passed in is used as it stands, budget included -- the caller chose it.
    max_pixels
        Pixel budget for the elevation grid this report is built from. The
        cover names the resolution that budget worked out to.

    Returns
    -------
    str
        The path written.
    """
    _, plt = _matplotlib()
    from matplotlib.backends.backend_pdf import PdfPages

    from .suitability import suitability
    from .terrain import drainage_density

    elevation = basin.dem(max_pixels=max_pixels) if dem is None else dem

    # Said before the work, not discovered during it. The figure is an
    # estimate from a measurement, so this warns and carries on rather than
    # refusing: a report is cheap to retry at a smaller budget and the caller
    # may know more about the machine than this does.
    pixels = int(getattr(elevation, "size", 0) or 0)
    if pixels:
        want = estimated_peak_gb(pixels)
        have = _available_gb()
        if have is not None and want > have:
            warnings.warn(
                f"This report is built from a {pixels / 1e6:.1f} megapixel "
                f"grid, which needs about {want:.1f} GB -- more than the "
                f"{have:.1f} GB this machine reports as available. Every page "
                "derives arrays of its own from the grid, so the cost is "
                f"roughly {BYTES_PER_PIXEL} bytes per pixel and not the grid "
                "alone. If the process is killed, pass a smaller max_pixels= "
                f"(the default is {REPORT_MAX_PIXELS / 1e6:.0f} Mpx) or "
                "report a sub-basin.",
                stacklevel=2)
    lat, lon = basin.centroid
    title = title or f"Catchment at {lat:.4f}, {lon:.4f}"

    grade = suitability(basin, dem=elevation, support_map=True)

    try:
        density: dict[str, Any] | None = drainage_density(elevation)
    except (MissingDependency, ValueError):
        density = None

    network_reason = None
    if rivers is None:
        try:
            rivers = basin.rivers()
        except Exception as exc:
            network_reason = f"{type(exc).__name__}: {exc}"
            rivers = None
    if morphometry is None and rivers is not None:
        try:
            from .morphometry import morphometry as compute

            morphometry = compute(basin, dem=elevation, rivers=rivers)
        except Exception as exc:
            network_reason = f"{type(exc).__name__}: {exc}"
            morphometry = None

    with PdfPages(str(path)) as pdf:
        _cover(plt, pdf, basin, title, grade, elevation, 1)
        _elevation_page(plt, pdf, basin, elevation, 2)
        _slope_page(plt, pdf, basin, elevation, 3)
        _landform_page(plt, pdf, basin, elevation, 4)
        if morphometry is not None:
            _network_page(plt, pdf, basin, rivers, morphometry, 5)
            _parameter_page(plt, pdf, morphometry, 6)
        else:
            _missing_network_page(plt, pdf, network_reason or "not requested", 5)
            _parameter_page(plt, pdf, {}, 6)
        _suitability_page(plt, pdf, basin, elevation, grade, density, 7)
        _methods_page(plt, pdf, basin, grade, density, 8)

        # Author and Keywords go into the PDF's own metadata, so the
        # attribution survives the file being renamed, forwarded or attached.
        from . import __author__, __doi__, __version__

        info = pdf.infodict()
        info["Title"] = title
        info["Author"] = __author__
        info["Creator"] = f"basinkit {__version__}"
        info["Subject"] = (
            f"Catchment analysis, elevation data suitability {grade['grade']}"
        )
        info["Keywords"] = (
            f"basinkit; hydrology; catchment; cite https://doi.org/{__doi__}"
        )

    del progress
    return str(path)

Clip and mask rasters to a basin polygon.

This is the part every generic downloader skips. eodag, earthaccess, pystac-client and friends all take a bounding box and give you whole scenes or tiles. For a river basin that is the wrong shape by a wide margin: a long dendritic catchment can occupy under a third of its own bbox, so two thirds of what you download, store and average over belongs to a neighbouring basin. Everything here works on the polygon.

basin_area_km2(geometry, crs='EPSG:4326')

Area of a lon/lat polygon in km2, via an equal-area projection.

Computing area in degrees is a common and badly wrong shortcut. This reprojects to a Lambert azimuthal equal-area centred on the basin itself, which is accurate to well under a percent at catchment scale.

Source code in basinkit/clip.py

def basin_area_km2(geometry, crs: str = "EPSG:4326") -> float:
    """Area of a lon/lat polygon in km2, via an equal-area projection.

    Computing area in degrees is a common and badly wrong shortcut. This
    reprojects to a Lambert azimuthal equal-area centred on the basin itself,
    which is accurate to well under a percent at catchment scale.
    """
    import geopandas as gpd
    from shapely.geometry import shape

    geoms = [shape(g) for g in _as_geoms(geometry)]
    gdf = gpd.GeoDataFrame(geometry=geoms, crs=crs)
    if gdf.crs is not None and gdf.crs.to_epsg() != 4326:
        gdf = gdf.to_crs("EPSG:4326")
    cx, cy = gdf.union_all().centroid.x, gdf.union_all().centroid.y
    return float(gdf.to_crs(laea_crs(cy, cx)).area.sum() / 1e6)

clip_raster(path_or_ds, geometry, *, crs='EPSG:4326', all_touched=False, nodata=None, drop_empty=True)

Open a raster, clip it to geometry and mask everything outside.

Parameters:

Name	Type	Description	Default
path_or_ds	Any	File path, URL, or an already-open xarray object.	required
geometry	Any	Basin polygon in crs.	required
all_touched	bool	Include pixels merely touched by the boundary. Set True for small basins on coarse grids, where strict centroid-in-polygon can return an empty array.	False
drop_empty	bool	Raise a helpful error instead of returning an all-nodata array.	True

Returns:

Type	Description
DataArray	Clipped to the polygon envelope and masked outside the polygon itself.

Source code in basinkit/clip.py

def clip_raster(
    path_or_ds: Any,
    geometry: Any,
    *,
    crs: str = "EPSG:4326",
    all_touched: bool = False,
    nodata: float | None = None,
    drop_empty: bool = True,
):
    """Open a raster, clip it to ``geometry`` and mask everything outside.

    Parameters
    ----------
    path_or_ds
        File path, URL, or an already-open ``xarray`` object.
    geometry
        Basin polygon in ``crs``.
    all_touched
        Include pixels merely touched by the boundary. Set ``True`` for small
        basins on coarse grids, where strict centroid-in-polygon can return an
        empty array.
    drop_empty
        Raise a helpful error instead of returning an all-nodata array.

    Returns
    -------
    xarray.DataArray
        Clipped to the polygon envelope and masked outside the polygon itself.
    """
    import rioxarray  # noqa: F401  (registers the .rio accessor)
    import xarray as xr

    geoms = _as_geoms(geometry)

    if isinstance(path_or_ds, (xr.DataArray, xr.Dataset)):
        da = path_or_ds
    else:
        da = rioxarray.open_rasterio(str(path_or_ds), masked=True, chunks="auto")

    if da.rio.crs is None:
        da = da.rio.write_crs(crs)

    if str(da.rio.crs).upper() != str(crs).upper():
        import geopandas as gpd
        from shapely.geometry import shape

        gdf = gpd.GeoDataFrame(
            geometry=[shape(g) for g in geoms], crs=crs
        ).to_crs(da.rio.crs)
        from shapely.geometry import mapping

        geoms = [mapping(g) for g in gdf.geometry]

    with warnings.catch_warnings():
        warnings.simplefilter("ignore")
        clipped = da.rio.clip(
            geoms, crs=da.rio.crs, all_touched=all_touched, drop=True, from_disk=True
        )

    if nodata is not None:
        clipped = clipped.rio.write_nodata(nodata)

    if drop_empty and not all_touched and _is_worth_checking(clipped):
        if _all_nodata(clipped):
            # Retry once with all_touched before giving up: this is the classic
            # small-basin-on-a-coarse-grid failure, not a real absence of data.
            return clip_raster(
                path_or_ds, geometry, crs=crs, all_touched=True,
                nodata=nodata, drop_empty=False,
            )
    return clipped

clip_stack(da, geometry, *, crs='EPSG:4326', all_touched=False)

Clip an already-loaded (possibly multi-temporal) DataArray to a polygon.

Source code in basinkit/clip.py

def clip_stack(da, geometry, *, crs: str = "EPSG:4326", all_touched: bool = False):
    """Clip an already-loaded (possibly multi-temporal) DataArray to a polygon."""
    return clip_raster(da, geometry, crs=crs, all_touched=all_touched)

laea_crs(lat, lon)

A Lambert azimuthal equal-area CRS centred on a lon/lat point.

Two things this does that a bare proj4 string did not.

It refuses a centre that is not in degrees. A polygon in a projected CRS has a centroid in metres, and passing a northing as lat_0 produces an internal PROJ error that says nothing about the real mistake.

It builds the CRS from a dict. The proj4-string path needs a PROJ able to parse it, and inside QGIS that is not always the PROJ the rest of the stack was built against.

Source code in basinkit/clip.py

def laea_crs(lat: float, lon: float):
    """A Lambert azimuthal equal-area CRS centred on a lon/lat point.

    Two things this does that a bare proj4 string did not.

    It refuses a centre that is not in degrees. A polygon in a projected CRS
    has a centroid in metres, and passing a northing as ``lat_0`` produces an
    internal PROJ error that says nothing about the real mistake.

    It builds the CRS from a dict. The proj4-string path needs a PROJ able to
    parse it, and inside QGIS that is not always the PROJ the rest of the
    stack was built against.
    """
    from pyproj import CRS

    if not (-90.0 <= lat <= 90.0 and -180.0 <= lon <= 180.0):
        raise ValueError(
            "An equal-area projection has to be centred on a lon/lat point, "
            f"but this geometry's centroid is lat={lat:,.3f}, lon={lon:,.3f}. "
            "Those are metres, not degrees, so the geometry is in a projected "
            "CRS. Reproject it to EPSG:4326, or say which CRS it is in: "
            "Basin.from_geometry(geom, crs='EPSG:32645')."
        )
    return CRS.from_dict({"proj": "laea", "lat_0": float(lat), "lon_0": float(lon),
                          "x_0": 0, "y_0": 0, "datum": "WGS84", "units": "m",
                          "no_defs": True})

zonal_mean(da, geometry=None, *, dims=('y', 'x'))

Area-weighted-ish spatial mean over a clipped array.

On a geographic grid, pixel area shrinks with cos(latitude). Ignoring that biases a basin mean toward its poleward end. For a small basin the error is negligible; for a Nile- or Ob-sized one it is not, so this weights by cos(lat) whenever a latitude coordinate is present.

Source code in basinkit/clip.py

def zonal_mean(da, geometry=None, *, dims: tuple[str, ...] = ("y", "x")):
    """Area-weighted-ish spatial mean over a clipped array.

    On a geographic grid, pixel area shrinks with ``cos(latitude)``. Ignoring
    that biases a basin mean toward its poleward end. For a small basin the
    error is negligible; for a Nile- or Ob-sized one it is not, so this
    weights by ``cos(lat)`` whenever a latitude coordinate is present.
    """
    if geometry is not None:
        da = clip_raster(da, geometry)

    lat_name = next((n for n in ("y", "lat", "latitude") if n in da.coords), None)
    if lat_name is not None:
        weights = np.cos(np.deg2rad(da[lat_name]))
        weights.name = "weights"
        present = [d for d in dims if d in da.dims]
        return da.weighted(weights.fillna(0)).mean(dim=present, skipna=True)
    return da.mean(dim=[d for d in dims if d in da.dims], skipna=True)

Tile mosaicking with a hard ceiling on how much comes into memory.

The naive version of this -- merge every tile at native resolution, then clip -- is fine for a 200 km2 catchment and fatal for a 50,000 km2 one. ESA WorldCover over the Koshi basin is nine tiles and a five-degree extent; at 10 m that is about 2.5 billion pixels, which no laptop will hold.

So the mosaic is budgeted. basinkit estimates the output size first and, if it exceeds max_pixels, coarsens by an integer factor and records that in the array's attributes. Nearest-neighbour resampling is used for categorical layers so class codes are never averaged into meaningless intermediates.

merge_tiles(paths, bounds, *, max_pixels=DEFAULT_MAX_PIXELS, categorical=False, nodata_below=None, src_nodata=None)

Merge tiles into one DataArray, coarsening if the result would be huge.

Returns:

Type	Description
(DataArray, dict)	The mosaic, and a dict describing whether and by how much it was coarsened. The caller writes that into .attrs so the decision is visible in the output rather than silent.

Source code in basinkit/mosaic.py

def merge_tiles(
    paths: list[Path],
    bounds: tuple[float, float, float, float],
    *,
    max_pixels: int = DEFAULT_MAX_PIXELS,
    categorical: bool = False,
    nodata_below: float | None = None,
    src_nodata: float | None = None,
):
    """Merge tiles into one DataArray, coarsening if the result would be huge.

    Returns
    -------
    (xarray.DataArray, dict)
        The mosaic, and a dict describing whether and by how much it was
        coarsened. The caller writes that into ``.attrs`` so the decision is
        visible in the output rather than silent.
    """
    import rasterio
    import rioxarray  # noqa: F401
    import xarray as xr
    from rasterio.enums import Resampling
    from rasterio.merge import merge

    srcs = [rasterio.open(p) for p in paths]
    crs_assumed = False
    try:
        native_res = abs(srcs[0].transform.a)
        want = estimate_pixels(bounds, native_res)

        factor = 1
        if want > max_pixels:
            factor = int(math.ceil(math.sqrt(want / max_pixels)))
        res = native_res * factor

        if factor > 1:
            warnings.warn(
                f"This extent is {want / 1e6:.0f} megapixels at native "
                f"{native_res * 111_320:.0f} m resolution, above the "
                f"{max_pixels / 1e6:.0f} Mpx budget. Coarsening {factor}x to "
                f"~{res * 111_320:.0f} m. Pass max_pixels= to change this, or "
                "work on a sub-basin to keep native resolution.",
                stacklevel=3,
            )

        # Several products use a sentinel value (255 in JRC surface water, 0 in
        # WorldCover) without declaring it as nodata in the file header.
        # rasterio only masks what a dataset calls nodata, so an averaging
        # resample happily blends the sentinel into real values -- a
        # water-occurrence grid that must top out at 100% comes back with
        # pixels at 101, and every basin mean is slightly wrong. Nearest-
        # neighbour never mixes two values, so it is used whenever the sentinel
        # is undeclared. Averaging is kept where nodata is declared properly
        # (the DEMs), because there it is the better downsampler.
        undeclared = src_nodata is not None and any(s_.nodata is None for s_ in srcs)
        resampling = (
            Resampling.nearest
            if (categorical or undeclared)
            else Resampling.average
        )

        merge_kwargs = {}
        if src_nodata is not None:
            merge_kwargs["nodata"] = src_nodata

        arr, transform = merge(
            srcs, bounds=bounds, res=(res, res), resampling=resampling,
            **merge_kwargs,
        )
        crs = srcs[0].crs
        if crs is None:
            # rasterio returns None when it cannot read the file's CRS at all,
            # which a PROJ installation missing its data directory will do --
            # QGIS's bundled stack among them. Every tile source here is
            # lon/lat, so say so rather than let an array with no CRS travel
            # on and fail later inside rioxarray, where the message names the
            # symptom and not the cause.
            crs = "EPSG:4326"
            crs_assumed = True
    finally:
        for handle in srcs:
            handle.close()

    ny, nx = arr.shape[-2], arr.shape[-1]
    xs = transform.c + transform.a * (np.arange(nx) + 0.5)
    ys = transform.f + transform.e * (np.arange(ny) + 0.5)

    da = xr.DataArray(
        arr if categorical else arr.astype("float32"),
        dims=("band", "y", "x"),
        coords={"band": np.arange(1, arr.shape[0] + 1), "y": ys, "x": xs},
    ).rio.write_crs(crs)

    if nodata_below is not None:
        da = da.where(da > nodata_below)
    if src_nodata is not None and not categorical:
        da = da.where(da != src_nodata)

    return da, {
        "basinkit_resampling": resampling.name,
        "basinkit_native_res_m": round(native_res * 111_320, 1),
        "basinkit_output_res_m": round(res * 111_320, 1),
        "basinkit_coarsen_factor": factor,
        "basinkit_tiles_merged": len(paths),
        "basinkit_crs_assumed": crs_assumed,
    }

Machine-readable catalogue of the open datasets basinkit can fetch.

Every entry records not just where the data lives but the two facts that decide whether a pipeline is reproducible: whether an account is needed, and what the licence permits. basinkit.catalog.table() prints it; the CLI and Basin.license_report() read from the same dict, so the documentation can never drift from the code.

Verified live 2026-08-25.

Dataset dataclass

One fetchable open dataset.

Source code in basinkit/catalog.py

@dataclass(frozen=True)
class Dataset:
    """One fetchable open dataset."""

    key: str
    name: str
    category: str
    resolution: str
    temporal: str
    coverage: str
    license: str
    auth: Auth
    route: str
    commercial_ok: bool
    redistributable: bool
    notes: str = ""
    citation: str = ""
    extras: tuple[str, ...] = field(default_factory=tuple)
    #: Whether basinkit can actually fetch this, as opposed to merely knowing
    #: about it. A catalogue that lists datasets it cannot deliver is worse
    #: than a shorter one: it reads as a feature list. Entries with
    #: ``implemented=False`` are documented pointers, and asking for one raises
    #: an error that says exactly how to obtain the data by hand.
    implemented: bool = True

    @property
    def clean(self) -> bool:
        """True when the layer is safe to use commercially *and* redistribute."""
        return self.commercial_ok and self.redistributable and self.auth == "none"

    @property
    def how_to_get(self) -> str:
        """Human-readable instructions for a dataset basinkit cannot fetch."""
        lines = [f"{self.name} is catalogued but basinkit cannot fetch it yet.",
                 f"    Access : {self.route}",
                 f"    Licence: {self.license}"]
        if self.auth != "none":
            lines.append(f"    Auth   : {self.auth} required")
        if self.notes:
            lines.append(f"    Note   : {self.notes}")
        lines.append(
            "    Once you have it locally, pass the geometry to "
            "basinkit.clip.clip_raster() to cut it to the basin."
        )
        return "\n".join(lines)

clean property

True when the layer is safe to use commercially and redistribute.

how_to_get property

Human-readable instructions for a dataset basinkit cannot fetch.

anonymous()

Datasets needing no account of any kind.

Source code in basinkit/catalog.py

def anonymous() -> list[Dataset]:
    """Datasets needing no account of any kind."""
    return [d for d in DATASETS.values() if d.auth == "none"]

implemented()

Datasets basinkit can actually fetch today.

Source code in basinkit/catalog.py

def implemented() -> list[Dataset]:
    """Datasets basinkit can actually fetch today."""
    return [d for d in DATASETS.values() if d.implemented]

require(key)

Return a dataset, raising with instructions if it cannot be fetched.

Source code in basinkit/catalog.py

def require(key: str) -> Dataset:
    """Return a dataset, raising with instructions if it cannot be fetched."""
    from .exceptions import NotImplementedSource

    ds = get(key)
    if not ds.implemented:
        raise NotImplementedSource(ds.how_to_get)
    return ds

table(datasets=None)

Render the catalogue as a plain-text table.

Source code in basinkit/catalog.py

def table(datasets: list[Dataset] | None = None) -> str:
    """Render the catalogue as a plain-text table."""
    rows = datasets if datasets is not None else list(DATASETS.values())
    rows = sorted(rows, key=lambda d: (d.category, d.key))
    head = f"{'key':<14} {'category':<12} {'auth':<8} {'comm':<5} {'fetch':<6} {'name'}"
    lines = [head, "-" * len(head)]
    for d in rows:
        lines.append(
            f"{d.key:<14} {d.category:<12} {d.auth:<8} "
            f"{'yes' if d.commercial_ok else 'NO':<5} "
            f"{'yes' if d.implemented else 'DOC':<6} {d.name}"
        )
    lines.append("")
    lines.append("fetch=DOC means basinkit documents the dataset but cannot download it;")
    lines.append("ask for it and you get instructions rather than a stack trace.")
    return "\n".join(lines)

unimplemented()

Datasets basinkit documents but cannot fetch.

Source code in basinkit/catalog.py

def unimplemented() -> list[Dataset]:
    """Datasets basinkit documents but cannot fetch."""
    return [d for d in DATASETS.values() if not d.implemented]