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 {}
        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

        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) -> Basin:
        """Wrap a polygon you already have (a gauge basin, an official boundary)."""
        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}>"
        )

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

        return dem(self.geometry, product=product, **kwargs)

    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 worldcover(self.geometry, year=year or 2021, **kwargs)
        if source == "esri":
            return esri_lulc(self.geometry, year=year, **kwargs)
        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 soilgrids(self.geometry, prop=prop, depth=depth, **kwargs)

    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``."""
        from .sources.climate import chirps, persiann, terraclimate

        if source == "chirps":
            return chirps(self.geometry, start, end, **kwargs)
        if source == "persiann":
            return persiann(self.geometry, str(start), end, **kwargs)
        if source == "terraclimate":
            return terraclimate(self.geometry, ("ppt",), int(start), end, **kwargs)
        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 global_surface_water(self.geometry, layer=layer, **kwargs)

    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")
        if hybas_id is None:
            raise ValueError(
                "BasinATLAS is keyed by HydroBASINS id, which only the "
                "'hydrobasins' backend records. Re-delineate with "
                "backend='hydrobasins', or pass a geometry to "
                "basinkit.sources.attributes.hydroatlas() directly."
            )
        gdf = hydroatlas(hybas_id=hybas_id, prefixes=prefixes, **kwargs)
        return describe(gdf.iloc[0])

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

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

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

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

    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

        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")
        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

        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")
        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

        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")
        return reduce_time(ds, composite) if composite else ds

    # -- summaries ---------------------------------------------------------
    def terrain_stats(self) -> dict:
        """Elevation, relief and mean slope: the standard morphometry."""
        import numpy as np

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

        res = abs(float(elev.rio.resolution()[0]))
        lat = self.centroid[0]
        cell_m = res * 111_320 * np.cos(np.deg2rad(lat))
        gy, gx = np.gradient(np.nan_to_num(np.asarray(elev.values, dtype="float64")))
        slope = np.degrees(np.arctan(np.hypot(gx, gy) / max(cell_m, 1e-6)))

        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(np.nanmean(slope)), 2),
            "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}"

        (outdir / "LICENSES.txt").write_text(self.license_report())
        manifest["layers"]["licenses"] = "LICENSES.txt"
        (outdir / "manifest.json").write_text(json.dumps(manifest, indent=2, default=str))
        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.

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")
    if hybas_id is None:
        raise ValueError(
            "BasinATLAS is keyed by HydroBASINS id, which only the "
            "'hydrobasins' backend records. Re-delineate with "
            "backend='hydrobasins', or pass a geometry to "
            "basinkit.sources.attributes.hydroatlas() directly."
        )
    gdf = hydroatlas(hybas_id=hybas_id, prefixes=prefixes, **kwargs)
    return describe(gdf.iloc[0])

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}"
        )

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

    return dem(self.geometry, product=product, **kwargs)

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}"

    (outdir / "LICENSES.txt").write_text(self.license_report())
    manifest["layers"]["licenses"] = "LICENSES.txt"
    (outdir / "manifest.json").write_text(json.dumps(manifest, indent=2, default=str))
    return manifest

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)

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) classmethod

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

Source code in basinkit/basin.py

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

    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)

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 hydrolakes(self.geometry, min_area_km2=min_area_km2, **kwargs)

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 worldcover(self.geometry, year=year or 2021, **kwargs)
    if source == "esri":
        return esri_lulc(self.geometry, year=year, **kwargs)
    raise ValueError(f"Unknown land cover source {source!r}: use 'worldcover' or 'esri'.")

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

    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")
    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.

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``."""
    from .sources.climate import chirps, persiann, terraclimate

    if source == "chirps":
        return chirps(self.geometry, start, end, **kwargs)
    if source == "persiann":
        return persiann(self.geometry, str(start), end, **kwargs)
    if source == "terraclimate":
        return terraclimate(self.geometry, ("ppt",), int(start), end, **kwargs)
    raise ValueError(
        f"Unknown precipitation source {source!r}: use 'chirps', 'persiann' "
        "or 'terraclimate'."
    )

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 hydrorivers(self.geometry, min_order=min_order, **kwargs)

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

    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")
    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

    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")
    return reduce_time(ds, composite) if composite else ds

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 soilgrids(self.geometry, prop=prop, depth=depth, **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 global_surface_water(self.geometry, layer=layer, **kwargs)

terrain_stats()

Elevation, relief and mean slope: the standard morphometry.

Source code in basinkit/basin.py

def terrain_stats(self) -> dict:
    """Elevation, relief and mean slope: the standard morphometry."""
    import numpy as np

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

    res = abs(float(elev.rio.resolution()[0]))
    lat = self.centroid[0]
    cell_m = res * 111_320 * np.cos(np.deg2rad(lat))
    gy, gx = np.gradient(np.nan_to_num(np.asarray(elev.values, dtype="float64")))
    slope = np.degrees(np.arctan(np.hypot(gx, gy) / max(cell_m, 1e-6)))

    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(np.nanmean(slope)), 2),
        "bbox_efficiency": round(self.bbox_efficiency, 3),
    }

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)

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 answer for small catchments, and the wrong one for large ones: cost grows with basin area and the window must 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
    return fn(lat, lon, **kwargs)

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, 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
max_window_deg	float	Stop growing at this half-width and raise instead of silently downloading the continent.	4.0

Source code in basinkit/delineate/dem.py

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,
    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.
    max_window_deg
        Stop growing at this half-width and raise instead of silently
        downloading the continent.
    """
    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:
        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()
        # 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
            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),
        "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,
        "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.
    """
    import geopandas as gpd

    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.
        frames = []
        for i in range(0, len(id_list), 2000):
            chunk = id_list[i : i + 2000]
            where = "HYBAS_ID IN (" + ",".join(str(x) for x in chunk) + ")"
            frames.append(
                gpd.read_file(
                    shp, where=where, columns=["HYBAS_ID", "SUB_AREA"], engine="pyogrio"
                )
            )
        units = gpd.GeoDataFrame(
            __import__("pandas").concat(frames, ignore_index=True), crs=frames[0].crs
        )

        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),
            **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'."
    )

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)
    cx, cy = gdf.union_all().centroid.x, gdf.union_all().centroid.y
    aea = f"+proj=laea +lat_0={cy} +lon_0={cx} +datum=WGS84 +units=m +no_defs"
    return float(gdf.to_crs(aea).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)

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]
    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
    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),
    }

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]