I/O Utilities

Raster and vector I/O, study-area parsing (files, GEE assets, bbox: strings, WKT), CRS/UTM helpers and value-scale conversions.

io

I/O utilities for raster and vector data.

Provides functions for reading/writing GeoTIFF files, vector formats (GeoJSON, GeoPackage, GeoParquet), study-area parsing (files, GEE assets, bbox: strings, WKT), CRS/UTM helpers and radiometric value-scale conversions.

get_equal_area_crs

get_equal_area_crs() -> pyproj.CRS

Return the NSIDC EASE-Grid 2.0 equal-area CRS (EPSG:6933).

Areas in this CRS are exact; lengths and distances are not (they are stretched away from the standard parallels at 30°N/S, e.g. by about 9 % at 55° latitude), so use a local UTM CRS or geodesic lengths for perimeters.

Returns:

Type	Description
CRS	EPSG:6933 equal-area cylindrical projection.

Source code in agribound/io/crs.py

def get_equal_area_crs() -> pyproj.CRS:
    """Return the NSIDC EASE-Grid 2.0 equal-area CRS (EPSG:6933).

    Areas in this CRS are exact; lengths and distances are not (they are
    stretched away from the standard parallels at 30°N/S, e.g. by about 9 % at
    55° latitude), so use a local UTM CRS or geodesic lengths for perimeters.

    Returns
    -------
    pyproj.CRS
        EPSG:6933 equal-area cylindrical projection.
    """
    return pyproj.CRS.from_epsg(6933)

get_utm_crs

get_utm_crs(lon: float, lat: float) -> pyproj.CRS

Determine the UTM CRS for a given longitude/latitude.

Parameters:

Name	Type	Description	Default
lon	float	Longitude in degrees.	required
lat	float	Latitude in degrees (lat >= 0 selects the northern zone).	required

Returns:

Type	Description
CRS	WGS 84 / UTM CRS (EPSG:326xx or EPSG:327xx).

Source code in agribound/io/crs.py

def get_utm_crs(lon: float, lat: float) -> pyproj.CRS:
    """Determine the UTM CRS for a given longitude/latitude.

    Parameters
    ----------
    lon : float
        Longitude in degrees.
    lat : float
        Latitude in degrees (``lat >= 0`` selects the northern zone).

    Returns
    -------
    pyproj.CRS
        WGS 84 / UTM CRS (EPSG:326xx or EPSG:327xx).
    """
    return pyproj.CRS.from_epsg(utm_epsg(utm_zone_for_lon(lon), south=lat < 0))

reproject_raster

reproject_raster(src_path: str | Path, dst_path: str | Path, dst_crs: Any, resolution: float | None = None, resampling: str = 'nearest') -> str

Reproject a raster to a different CRS.

Parameters:

Name	Type	Description	Default
src_path	str or Path	Source raster file.	required
dst_path	str or Path	Destination raster file.	required
dst_crs	CRS or str	Target coordinate reference system.	required
resolution	float or None	Output pixel resolution. If None, computed automatically.	None
resampling	str	Resampling method: "nearest", "bilinear", "cubic".	'nearest'

Returns:

Type	Description
str	Path to the reprojected raster.

Source code in agribound/io/crs.py

def reproject_raster(
    src_path: str | Path,
    dst_path: str | Path,
    dst_crs: Any,
    resolution: float | None = None,
    resampling: str = "nearest",
) -> str:
    """Reproject a raster to a different CRS.

    Parameters
    ----------
    src_path : str or Path
        Source raster file.
    dst_path : str or Path
        Destination raster file.
    dst_crs : CRS or str
        Target coordinate reference system.
    resolution : float or None
        Output pixel resolution. If *None*, computed automatically.
    resampling : str
        Resampling method: ``"nearest"``, ``"bilinear"``, ``"cubic"``.

    Returns
    -------
    str
        Path to the reprojected raster.
    """
    resample_map = {
        "nearest": Resampling.nearest,
        "bilinear": Resampling.bilinear,
        "cubic": Resampling.cubic,
    }
    resample_method = resample_map.get(resampling, Resampling.nearest)

    src_path = Path(src_path)
    dst_path = Path(dst_path)
    dst_path.parent.mkdir(parents=True, exist_ok=True)

    with rasterio.open(src_path) as src:
        kwargs = {}
        if resolution is not None:
            kwargs["resolution"] = resolution

        transform, width, height = calculate_default_transform(
            src.crs, dst_crs, src.width, src.height, *src.bounds, **kwargs
        )

        meta = src.meta.copy()
        meta.update(
            {
                "crs": dst_crs,
                "transform": transform,
                "width": width,
                "height": height,
                "compress": "lzw",
            }
        )

        with rasterio.open(dst_path, "w", **meta) as dst:
            for band_idx in range(1, src.count + 1):
                reproject(
                    source=rasterio.band(src, band_idx),
                    destination=rasterio.band(dst, band_idx),
                    src_transform=src.transform,
                    src_crs=src.crs,
                    dst_transform=transform,
                    dst_crs=dst_crs,
                    resampling=resample_method,
                )

    return str(dst_path)

utm_crs_for_geometry

utm_crs_for_geometry(geom_4326: Any) -> pyproj.CRS

Return the UTM CRS of the centroid of a geometry in EPSG:4326.

Parameters:

Name	Type	Description	Default
geom_4326	shapely geometry, GeoSeries or GeoDataFrame	Geometry in EPSG:4326. GeoSeries/GeoDataFrames are reprojected to EPSG:4326 if needed and unioned.	required

Returns:

Type	Description
CRS	UTM CRS of the centroid (see :func:get_utm_crs).

Raises:

Type	Description
ValueError	If the geometry is empty.

Source code in agribound/io/crs.py

def utm_crs_for_geometry(geom_4326: Any) -> pyproj.CRS:
    """Return the UTM CRS of the centroid of a geometry in EPSG:4326.

    Parameters
    ----------
    geom_4326 : shapely geometry, GeoSeries or GeoDataFrame
        Geometry in EPSG:4326. GeoSeries/GeoDataFrames are reprojected to
        EPSG:4326 if needed and unioned.

    Returns
    -------
    pyproj.CRS
        UTM CRS of the centroid (see :func:`get_utm_crs`).

    Raises
    ------
    ValueError
        If the geometry is empty.
    """
    geom = geom_4326
    if hasattr(geom, "geometry") and hasattr(geom, "crs"):
        series = geom.geometry
        if series.crs is not None and not series.crs.equals("EPSG:4326"):
            series = series.to_crs("EPSG:4326")
        geom = series.union_all()
    if geom is None or geom.is_empty:
        raise ValueError("Cannot determine a UTM zone for an empty geometry")
    centroid = geom.centroid
    return get_utm_crs(centroid.x, centroid.y)

utm_zones_for_bounds

utm_zones_for_bounds(bounds: tuple[float, float, float, float]) -> list[int]

Return the UTM zone numbers spanned by an EPSG:4326 bounding box.

Parameters:

Name	Type	Description	Default
bounds	tuple[float, float, float, float]	(min_lon, min_lat, max_lon, max_lat). A box crossing the antimeridian is given with min_lon > max_lon.	required

Returns:

Type	Description
list[int]	Zone numbers from west to east. A box whose eastern edge lies exactly on a zone boundary does not include the next zone.

Source code in agribound/io/crs.py

def utm_zones_for_bounds(bounds: tuple[float, float, float, float]) -> list[int]:
    """Return the UTM zone numbers spanned by an EPSG:4326 bounding box.

    Parameters
    ----------
    bounds : tuple[float, float, float, float]
        ``(min_lon, min_lat, max_lon, max_lat)``. A box crossing the
        antimeridian is given with ``min_lon > max_lon``.

    Returns
    -------
    list[int]
        Zone numbers from west to east. A box whose eastern edge lies exactly
        on a zone boundary does not include the next zone.
    """
    min_lon, _min_lat, max_lon, _max_lat = (float(v) for v in bounds)
    first = utm_zone_for_lon(min_lon)
    # Step 1e-9 degrees (~0.1 mm) west so an edge on a zone boundary stays in the zone.
    east = max_lon - 1e-9 if max_lon != min_lon else max_lon
    last = utm_zone_for_lon(east)
    if min_lon <= max_lon and first <= last:
        return list(range(first, last + 1))
    # Crossing the antimeridian (or wrapping at +/-180).
    return list(range(first, 61)) + list(range(1, last + 1))

get_raster_info

get_raster_info(path: str | Path) -> RasterInfo

Read metadata from a raster file without loading pixel data.

Parameters:

Name	Type	Description	Default
path	str or Path	Path to the raster file (GeoTIFF).	required

Returns:

Type	Description
RasterInfo	Raster metadata.

Raises:

Type	Description
FileNotFoundError	If the file does not exist.

Source code in agribound/io/raster.py

def get_raster_info(path: str | Path) -> RasterInfo:
    """Read metadata from a raster file without loading pixel data.

    Parameters
    ----------
    path : str or Path
        Path to the raster file (GeoTIFF).

    Returns
    -------
    RasterInfo
        Raster metadata.

    Raises
    ------
    FileNotFoundError
        If the file does not exist.
    """
    path = Path(path)
    if not path.exists():
        raise FileNotFoundError(f"Raster file not found: {path}")

    with rasterio.open(path) as src:
        return RasterInfo(
            path=str(path),
            width=src.width,
            height=src.height,
            count=src.count,
            crs=src.crs,
            transform=src.transform,
            bounds=src.bounds,
            dtype=str(src.dtypes[0]),
            nodata=src.nodata,
            res=src.res,
        )

percentile_stretch_uint8

percentile_stretch_uint8(arr: ndarray, nodata: float | None = None, low: float = 1.0, high: float = 99.0, per_band: bool = True, *, return_bounds: bool = False, return_bounds_only: bool = False, max_sample_side: int = MAX_SAMPLE_SIDE)

Stretch an array to uint8 between per-band percentiles.

Mirrors Delineate-Anything's DataAnalyser/DataLoaderCached normalisation: uint8 input is returned unchanged; otherwise the low and high percentiles are computed from valid, finite, strictly positive pixels on a nearest-neighbour sampling grid of at most max_sample_side pixels per side, and values are mapped with clip(255 * (v - lo) / (hi - lo), 0, 255) and truncated to uint8. Invalid pixels (non-finite in any band, or equal to nodata in all bands) are set to 0.

Parameters:

Name	Type	Description	Default
arr	ndarray	(bands, height, width) or (height, width) array.	required
nodata	float or None	Nodata value; a pixel is nodata when all bands equal it.	None
low	float	Percentiles (default 1 and 99).	1.0
high	float	Percentiles (default 1 and 99).	1.0
per_band	bool	Compute bounds per band (default) or pooled over all bands.	True
return_bounds	bool	Also return the (lows, highs) lists used.	False
return_bounds_only	bool	Return only (lows, highs), without building the uint8 output. The bounds are identical to those of the full call. Only the pixels of the sampling grid are examined (all pixels when the array is at most max_sample_side pixels per side). With per_band at most one float64 copy of one band's valid, positive sampled pixels exists at a time; pooled bounds concatenate the copies of all bands. On a 4000 x 4000 x 3 float32 array (192 MB; the sample is the whole array) the tracemalloc peak of this call is 207 MB, against 320 MB for return_bounds=True and 784 MB for return_bounds=True in the implementation before 2026-09-27 (measured 2026-09-27; excludes the input array).	False
max_sample_side	int	Maximum side of the percentile sampling grid (default 4096).	MAX_SAMPLE_SIDE

Returns:

Type	Description
ndarray or tuple	uint8 array of the input shape; (array, lows, highs) when return_bounds is True; (lows, highs) when return_bounds_only is True.

Raises:

Type	Description
ValueError	If a band has no valid positive pixels.

Source code in agribound/io/raster.py

def percentile_stretch_uint8(
    arr: np.ndarray,
    nodata: float | None = None,
    low: float = 1.0,
    high: float = 99.0,
    per_band: bool = True,
    *,
    return_bounds: bool = False,
    return_bounds_only: bool = False,
    max_sample_side: int = MAX_SAMPLE_SIDE,
):
    """Stretch an array to uint8 between per-band percentiles.

    Mirrors Delineate-Anything's ``DataAnalyser``/``DataLoaderCached``
    normalisation: uint8 input is returned unchanged; otherwise the *low* and
    *high* percentiles are computed from valid, finite, strictly positive
    pixels on a nearest-neighbour sampling grid of at most
    ``max_sample_side`` pixels per side, and values are mapped with
    ``clip(255 * (v - lo) / (hi - lo), 0, 255)`` and truncated to uint8.
    Invalid pixels (non-finite in any band, or equal to *nodata* in all
    bands) are set to 0.

    Parameters
    ----------
    arr : numpy.ndarray
        ``(bands, height, width)`` or ``(height, width)`` array.
    nodata : float or None
        Nodata value; a pixel is nodata when all bands equal it.
    low, high : float
        Percentiles (default 1 and 99).
    per_band : bool
        Compute bounds per band (default) or pooled over all bands.
    return_bounds : bool
        Also return the ``(lows, highs)`` lists used.
    return_bounds_only : bool
        Return only ``(lows, highs)``, without building the uint8 output.
        The bounds are identical to those of the full call. Only the pixels
        of the sampling grid are examined (all pixels when the array is at
        most ``max_sample_side`` pixels per side). With *per_band* at most
        one float64 copy of one band's valid, positive sampled pixels exists
        at a time; pooled bounds concatenate the copies of all bands. On a
        4000 x 4000 x 3 float32 array (192 MB; the sample is the whole
        array) the tracemalloc peak of this call is 207 MB, against 320 MB
        for ``return_bounds=True`` and 784 MB for ``return_bounds=True`` in
        the implementation before 2026-09-27 (measured 2026-09-27; excludes
        the input array).
    max_sample_side : int
        Maximum side of the percentile sampling grid (default 4096).

    Returns
    -------
    numpy.ndarray or tuple
        uint8 array of the input shape; ``(array, lows, highs)`` when
        *return_bounds* is *True*; ``(lows, highs)`` when
        *return_bounds_only* is *True*.

    Raises
    ------
    ValueError
        If a band has no valid positive pixels.
    """
    data = np.asarray(arr)
    squeeze = data.ndim == 2
    if squeeze:
        data = data[np.newaxis, ...]
    if data.ndim != 3:
        raise ValueError(f"Expected a 2-D or 3-D array, got shape {np.asarray(arr).shape}")

    if data.dtype == np.uint8:
        lows = [0.0] * data.shape[0]
        highs = [255.0] * data.shape[0]
        if return_bounds_only:
            return lows, highs
        out = data.copy()
    else:
        lows, highs = _stretch_bounds(data, nodata, low, high, per_band, max_sample_side)
        if return_bounds_only:
            return lows, highs
        valid = _valid_pixels(data, nodata)
        out = np.zeros(data.shape, dtype=np.uint8)
        for i, (lo, hi) in enumerate(zip(lows, highs, strict=True)):
            span = hi - lo if hi > lo else 1e-12
            # 255 * ((v - lo) / span), evaluated in place in the same order.
            band = data[i].astype(np.float64)
            band -= lo
            band /= span
            band *= 255.0
            np.clip(band, 0, 255, out=band)
            band[~(valid & np.isfinite(band))] = 0
            out[i] = band.astype(np.uint8)

    if squeeze:
        out = out[0]
    if return_bounds:
        return out, lows, highs
    return out

read_raster

read_raster(path: str | Path, bands: list[int] | None = None, window: Window | None = None) -> tuple[np.ndarray, dict[str, Any]]

Read a raster file into a NumPy array.

Parameters:

Name	Type	Description	Default
path	str or Path	Path to the raster file.	required
bands	list[int] or None	1-based band indices to read. None reads all bands.	None
window	Window or None	Spatial sub-window to read. None reads the full extent.	None

Returns:

Name	Type	Description
data	ndarray	Pixel data with shape (bands, height, width).
meta	dict	Rasterio metadata dictionary (crs, transform, width, height, etc.).

Source code in agribound/io/raster.py

def read_raster(
    path: str | Path,
    bands: list[int] | None = None,
    window: rasterio.windows.Window | None = None,
) -> tuple[np.ndarray, dict[str, Any]]:
    """Read a raster file into a NumPy array.

    Parameters
    ----------
    path : str or Path
        Path to the raster file.
    bands : list[int] or None
        1-based band indices to read. *None* reads all bands.
    window : rasterio.windows.Window or None
        Spatial sub-window to read. *None* reads the full extent.

    Returns
    -------
    data : numpy.ndarray
        Pixel data with shape ``(bands, height, width)``.
    meta : dict
        Rasterio metadata dictionary (crs, transform, width, height, etc.).
    """
    path = Path(path)

    with rasterio.open(path) as src:
        if bands is None:
            bands = list(range(1, src.count + 1))

        data = src.read(bands, window=window)
        meta = src.meta.copy()

        if window is not None:
            meta.update(
                {
                    "width": window.width,
                    "height": window.height,
                    "transform": src.window_transform(window),
                }
            )

        meta["count"] = len(bands)
        return data, meta

to_s2_dn

to_s2_dn(arr: ndarray, source: str, allow_unknown: bool = False, *, value_scale: str | None = None) -> np.ndarray

Convert composite pixel values to Sentinel-2 L2A units (reflectance x 10000).

COPERNICUS/S2_SR_HARMONIZED stores surface reflectance x 10000 with the processing-baseline offset removed, so this is to_unit_reflectance(arr) * 10000. Landsat/HLS composites are already on this scale after the 1.0 harmonisation; uint8 imagery is only rescaled (0-255 -> 0-10000), which does not make it radiometrically comparable.

Parameters:

Name	Type	Description	Default
arr	ndarray	See :func:to_unit_reflectance.	required
source	ndarray	See :func:to_unit_reflectance.	required
allow_unknown	ndarray	See :func:to_unit_reflectance.	required
value_scale	ndarray	See :func:to_unit_reflectance.	required

Returns:

Type	Description
ndarray	float32 array of the same shape.

Source code in agribound/io/raster.py

def to_s2_dn(
    arr: np.ndarray,
    source: str,
    allow_unknown: bool = False,
    *,
    value_scale: str | None = None,
) -> np.ndarray:
    """Convert composite pixel values to Sentinel-2 L2A units (reflectance x 10000).

    ``COPERNICUS/S2_SR_HARMONIZED`` stores surface reflectance x 10000 with
    the processing-baseline offset removed, so this is
    ``to_unit_reflectance(arr) * 10000``. Landsat/HLS composites are already
    on this scale after the 1.0 harmonisation; uint8 imagery is only rescaled
    (0-255 -> 0-10000), which does not make it radiometrically comparable.

    Parameters
    ----------
    arr, source, allow_unknown, value_scale
        See :func:`to_unit_reflectance`.

    Returns
    -------
    numpy.ndarray
        float32 array of the same shape.
    """
    unit = to_unit_reflectance(arr, source, allow_unknown, value_scale=value_scale)
    return (unit * np.float32(10000.0)).astype(np.float32)

to_unit_reflectance

to_unit_reflectance(arr: ndarray, source: str, allow_unknown: bool = False, *, value_scale: str | None = None) -> np.ndarray

Convert composite pixel values to 0-1 reflectance-like floats.

"reflectance_x10000" sources are divided by 10000 and "uint8" sources by 255 (NAIP digital numbers are not calibrated reflectance; the result is only a 0-1 rescaling). NaN is preserved.

Parameters:

Name	Type	Description	Default
arr	ndarray	Pixel values as written by the source's composite builder.	required
source	str	Source name (its value_scale is looked up in the registry).	required
allow_unknown	bool	For "dn"/"unknown" sources (SPOT, local): infer the scale with :func:infer_value_scale and log a warning instead of raising.	False
value_scale	str or None	Explicit scale ("reflectance_x10000", "uint8" or "unit") overriding the registry, e.g. for a local file of known radiometry.	None

Returns:

Type	Description
ndarray	float32 array of the same shape.

Raises:

Type	Description
ValueError	For embedding sources, and for "dn"/"unknown" sources unless allow_unknown or value_scale is given.

Source code in agribound/io/raster.py

def to_unit_reflectance(
    arr: np.ndarray,
    source: str,
    allow_unknown: bool = False,
    *,
    value_scale: str | None = None,
) -> np.ndarray:
    """Convert composite pixel values to 0-1 reflectance-like floats.

    ``"reflectance_x10000"`` sources are divided by 10000 and ``"uint8"``
    sources by 255 (NAIP digital numbers are not calibrated reflectance; the
    result is only a 0-1 rescaling). NaN is preserved.

    Parameters
    ----------
    arr : numpy.ndarray
        Pixel values as written by the source's composite builder.
    source : str
        Source name (its ``value_scale`` is looked up in the registry).
    allow_unknown : bool
        For ``"dn"``/``"unknown"`` sources (SPOT, local): infer the scale with
        :func:`infer_value_scale` and log a warning instead of raising.
    value_scale : str or None
        Explicit scale (``"reflectance_x10000"``, ``"uint8"`` or ``"unit"``)
        overriding the registry, e.g. for a local file of known radiometry.

    Returns
    -------
    numpy.ndarray
        float32 array of the same shape.

    Raises
    ------
    ValueError
        For embedding sources, and for ``"dn"``/``"unknown"`` sources unless
        *allow_unknown* or *value_scale* is given.
    """
    divisor = _scale_divisor(source, arr, allow_unknown, value_scale)
    return (np.asarray(arr, dtype=np.float32) / np.float32(divisor)).astype(np.float32)

write_raster

write_raster(path: str | Path, data: ndarray, crs: Any, transform: Any, nodata: float | None = None, dtype: str | None = None, compress: str = 'lzw') -> str

Write a NumPy array as a GeoTIFF.

Parameters:

Name	Type	Description	Default
path	str or Path	Destination file path.	required
data	ndarray	Pixel data with shape (bands, height, width) or (height, width).	required
crs	CRS or str	Coordinate reference system.	required
transform	Affine	Affine transform.	required
nodata	float or None	Nodata value to encode in the file.	None
dtype	str or None	Output data type. Defaults to the array dtype.	None
compress	str	Compression method (default "lzw").	'lzw'

Returns:

Type	Description
str	Path to the written file.

Source code in agribound/io/raster.py

def write_raster(
    path: str | Path,
    data: np.ndarray,
    crs: Any,
    transform: Any,
    nodata: float | None = None,
    dtype: str | None = None,
    compress: str = "lzw",
) -> str:
    """Write a NumPy array as a GeoTIFF.

    Parameters
    ----------
    path : str or Path
        Destination file path.
    data : numpy.ndarray
        Pixel data with shape ``(bands, height, width)`` or ``(height, width)``.
    crs : rasterio.crs.CRS or str
        Coordinate reference system.
    transform : rasterio.transform.Affine
        Affine transform.
    nodata : float or None
        Nodata value to encode in the file.
    dtype : str or None
        Output data type. Defaults to the array dtype.
    compress : str
        Compression method (default ``"lzw"``).

    Returns
    -------
    str
        Path to the written file.
    """
    path = Path(path)
    path.parent.mkdir(parents=True, exist_ok=True)

    if data.ndim == 2:
        data = data[np.newaxis, ...]

    count, height, width = data.shape
    if dtype is None:
        dtype = str(data.dtype)

    with rasterio.open(
        path,
        "w",
        driver="GTiff",
        height=height,
        width=width,
        count=count,
        dtype=dtype,
        crs=crs,
        transform=transform,
        nodata=nodata,
        compress=compress,
        tiled=True,
        blockxsize=256,
        blockysize=256,
        BIGTIFF="YES",
    ) as dst:
        dst.write(data)

    return str(path)

read_config_study_area

read_config_study_area(config: Any) -> gpd.GeoDataFrame

Read config.study_area, keeping a local copy of a GEE asset.

Vector files, "bbox:..." strings and WKT are read with :func:read_study_area. A GEE asset ID is read from its local copy (:func:study_area_cache_file) when that file exists, without contacting Earth Engine; otherwise it is read from Earth Engine with the configured credentials (read_study_area(..., config=config)) and its geometries are written to the copy. The pipeline stages that read the study area (composite builders, study-area selection, evaluation, FTW window dates, LULC filter) do so through this function, so a GEE-asset study area read once on a node with network access (e.g. by agribound composite) is available to a later agribound delineate with the same cache directory on a node without it.

The copy holds the geometries only (EPSG:4326, no properties). Like the cached composites, it is keyed by the asset ID, not by the asset's content: after changing an asset in place, delete the copy (and the cached composites) or use another cache directory. A copy that cannot be read is replaced; a copy that cannot be written is skipped with a WARNING.

Parameters:

Name	Type	Description	Default
config	AgriboundConfig	Configuration with a study_area.	required

Returns:

Type	Description
GeoDataFrame	Study area geometry (for a GEE asset read from Earth Engine, with the asset's properties; from the copy, without them).

Source code in agribound/io/vector.py

def read_config_study_area(config: Any) -> gpd.GeoDataFrame:
    """Read ``config.study_area``, keeping a local copy of a GEE asset.

    Vector files, ``"bbox:..."`` strings and WKT are read with
    :func:`read_study_area`. A GEE asset ID is read from its local copy
    (:func:`study_area_cache_file`) when that file exists, without contacting
    Earth Engine; otherwise it is read from Earth Engine with the configured
    credentials (``read_study_area(..., config=config)``) and its geometries
    are written to the copy. The pipeline stages that read the study area
    (composite builders, study-area selection, evaluation, FTW window dates,
    LULC filter) do so through this function, so a GEE-asset study area read
    once on a node with network access (e.g. by ``agribound composite``) is
    available to a later ``agribound delineate`` with the same cache
    directory on a node without it.

    The copy holds the geometries only (EPSG:4326, no properties). Like the
    cached composites, it is keyed by the asset ID, not by the asset's
    content: after changing an asset in place, delete the copy (and the
    cached composites) or use another cache directory. A copy that cannot be
    read is replaced; a copy that cannot be written is skipped with a
    WARNING.

    Parameters
    ----------
    config : AgriboundConfig
        Configuration with a ``study_area``.

    Returns
    -------
    geopandas.GeoDataFrame
        Study area geometry (for a GEE asset read from Earth Engine, with the
        asset's properties; from the copy, without them).
    """
    from shapely.geometry import mapping

    path = study_area_cache_file(config)
    if path is None:
        return read_study_area(config.study_area, config=config)
    if path.exists():
        try:
            data = json.loads(path.read_text())
            geometries = [
                shape(f["geometry"]) if f.get("geometry") else None for f in data["features"]
            ]
            logger.debug("Using the local copy of study area %s: %s", config.study_area, path)
            return gpd.GeoDataFrame(geometry=geometries, crs="EPSG:4326")
        except (OSError, ValueError, KeyError, TypeError) as exc:
            logger.warning("Ignoring unreadable study-area copy %s: %s", path, exc)
    gdf = read_study_area(config.study_area, config=config)
    geoms_4326 = gdf.geometry if gdf.crs is None else gdf.geometry.to_crs("EPSG:4326")
    features = [
        {
            "type": "Feature",
            "geometry": None if g is None or g.is_empty else mapping(g),
            "properties": {},
        }
        for g in geoms_4326
    ]
    payload = {
        "type": "FeatureCollection",
        "agribound:asset_id": str(config.study_area).strip(),
        "features": features,
    }
    partial = path.with_name(path.name + ".partial")
    try:
        partial.write_text(json.dumps(payload))
        os.replace(partial, path)
        logger.info("Saved a local copy of study area %s: %s", config.study_area, path)
    except OSError as exc:
        logger.warning("Could not save a local copy of study area %s: %s", config.study_area, exc)
    return gdf

read_study_area

read_study_area(path: str | Path, config: Any = None) -> gpd.GeoDataFrame

Read a study area definition.

Parameters:

Name	Type	Description	Default
path	str or Path	One of: a local vector file (GeoJSON, GeoPackage, Shapefile, GeoParquet, FlatGeobuf); a GEE FeatureCollection asset ID ("projects/..." or "users/..."); a bounding box "bbox:minx,miny,maxx,maxy" in EPSG:4326; a WKT geometry (POLYGON, MULTIPOLYGON, ...) in EPSG:4326, or EWKT with an SRID=<epsg>; prefix.	required
config	AgriboundConfig or None	Used only for GEE asset IDs: Earth Engine is initialised with :func:agribound.auth.ensure_gee, i.e. with the configured project, service-account key, endpoint and workload tag. Without it, an uninitialised Earth Engine client is set up with :func:agribound.auth.setup_gee defaults (GEE_PROJECT, gcloud project or key project).	None

Returns:

Type	Description
GeoDataFrame	Study area geometry.

Raises:

Type	Description
ValueError	If the string is not a recognised format, or the bbox/WKT is invalid.
FileNotFoundError	If a vector file path does not exist.

Source code in agribound/io/vector.py

def read_study_area(path: str | Path, config: Any = None) -> gpd.GeoDataFrame:
    """Read a study area definition.

    Parameters
    ----------
    path : str or Path
        One of:

        - a local vector file (GeoJSON, GeoPackage, Shapefile, GeoParquet,
          FlatGeobuf);
        - a GEE FeatureCollection asset ID (``"projects/..."`` or
          ``"users/..."``);
        - a bounding box ``"bbox:minx,miny,maxx,maxy"`` in EPSG:4326;
        - a WKT geometry (``POLYGON``, ``MULTIPOLYGON``, ...) in EPSG:4326, or
          EWKT with an ``SRID=<epsg>;`` prefix.
    config : AgriboundConfig or None
        Used only for GEE asset IDs: Earth Engine is initialised with
        :func:`agribound.auth.ensure_gee`, i.e. with the configured project,
        service-account key, endpoint and workload tag. Without it, an
        uninitialised Earth Engine client is set up with
        :func:`agribound.auth.setup_gee` defaults (``GEE_PROJECT``, ``gcloud``
        project or key project).

    Returns
    -------
    geopandas.GeoDataFrame
        Study area geometry.

    Raises
    ------
    ValueError
        If the string is not a recognised format, or the bbox/WKT is invalid.
    FileNotFoundError
        If a vector file path does not exist.
    """
    text = str(path).strip()

    # GEE asset ID
    if text.startswith("projects/") or text.startswith("users/"):
        return _read_gee_asset(text, config=config)

    # Bounding box
    if text.lower().startswith("bbox:"):
        minx, miny, maxx, maxy = parse_bbox(text)
        return gpd.GeoDataFrame(
            {"name": ["bbox"]}, geometry=[box(minx, miny, maxx, maxy)], crs="EPSG:4326"
        )

    # WKT geometry
    if _WKT_RE.match(text):
        from shapely import wkt

        crs = "EPSG:4326"
        wkt_text = text
        if text.upper().startswith("SRID="):
            srid, wkt_text = text.split(";", 1)
            crs = f"EPSG:{int(srid.split('=', 1)[1])}"
        try:
            geom = wkt.loads(wkt_text)
        except Exception as exc:
            raise ValueError(f"Invalid WKT study area: {exc}") from exc
        if geom.is_empty:
            raise ValueError("WKT study area is empty")
        return gpd.GeoDataFrame({"name": ["wkt"]}, geometry=[geom], crs=crs)

    return read_vector(text)

read_vector

read_vector(path: str | Path) -> gpd.GeoDataFrame

Read a vector file into a GeoDataFrame.

Supports GeoJSON, GeoPackage, Shapefile, and GeoParquet formats.

Parameters:

Name	Type	Description	Default
path	str or Path	Path to the vector file.	required

Returns:

Type	Description
GeoDataFrame	Loaded vector data.

Raises:

Type	Description
FileNotFoundError	If the file does not exist.
ValueError	If the file format is not supported.

Source code in agribound/io/vector.py

def read_vector(path: str | Path) -> gpd.GeoDataFrame:
    """Read a vector file into a GeoDataFrame.

    Supports GeoJSON, GeoPackage, Shapefile, and GeoParquet formats.

    Parameters
    ----------
    path : str or Path
        Path to the vector file.

    Returns
    -------
    geopandas.GeoDataFrame
        Loaded vector data.

    Raises
    ------
    FileNotFoundError
        If the file does not exist.
    ValueError
        If the file format is not supported.
    """
    path = Path(path)
    if not path.exists():
        raise FileNotFoundError(f"Vector file not found: {path}")

    suffix = path.suffix.lower()
    if suffix in (".parquet", ".geoparquet"):
        return gpd.read_parquet(path)
    elif suffix in (".geojson", ".json", ".gpkg", ".shp", ".fgb"):
        return gpd.read_file(path)
    else:
        raise ValueError(
            f"Unsupported vector format: {suffix!r}. "
            "Supported: .geojson, .json, .gpkg, .shp, .parquet, .geoparquet, .fgb"
        )

write_vector

write_vector(gdf: GeoDataFrame, path: str | Path, format: str | None = None) -> str

Write a GeoDataFrame to a vector file.

When writing to .parquet, the output is fiboa-compliant GeoParquet.

Parameters:

Name	Type	Description	Default
gdf	GeoDataFrame	Vector data to write.	required
path	str or Path	Destination file path.	required
format	str or None	Output format ("gpkg", "geojson", "parquet", "shp", "fgb"). If None, inferred from the file extension. A format that contradicts a recognised extension raises :class:ValueError.	None

Returns:

Type	Description
str	Path to the written file.

Raises:

Type	Description
ValueError	If the format cannot be inferred, is unsupported, or contradicts the file extension.

Source code in agribound/io/vector.py

def write_vector(
    gdf: gpd.GeoDataFrame,
    path: str | Path,
    format: str | None = None,
) -> str:
    """Write a GeoDataFrame to a vector file.

    When writing to ``.parquet``, the output is fiboa-compliant GeoParquet.

    Parameters
    ----------
    gdf : geopandas.GeoDataFrame
        Vector data to write.
    path : str or Path
        Destination file path.
    format : str or None
        Output format (``"gpkg"``, ``"geojson"``, ``"parquet"``, ``"shp"``,
        ``"fgb"``). If *None*, inferred from the file extension. A format
        that contradicts a recognised extension raises :class:`ValueError`.

    Returns
    -------
    str
        Path to the written file.

    Raises
    ------
    ValueError
        If the format cannot be inferred, is unsupported, or contradicts the
        file extension.
    """
    path = Path(path)
    suffix = path.suffix.lower()
    format_map = {
        ".gpkg": "gpkg",
        ".geojson": "geojson",
        ".json": "geojson",
        ".parquet": "parquet",
        ".geoparquet": "parquet",
        ".shp": "shp",
        ".fgb": "fgb",
    }
    implied = format_map.get(suffix)
    if format is None:
        format = implied
        if format is None:
            raise ValueError(f"Cannot infer format from extension {suffix!r}")
    else:
        format = str(format).lower().strip()
        if implied is not None and implied != format:
            raise ValueError(
                f"format={format!r} contradicts the extension of {path.name!r} "
                f"({implied!r}); change the extension or the format."
            )
    path.parent.mkdir(parents=True, exist_ok=True)

    if format == "parquet":
        _write_fiboa_parquet(gdf, path)
    elif format == "geojson":
        # GeoJSON requires EPSG:4326
        if gdf.crs is not None and not gdf.crs.equals("EPSG:4326"):
            gdf = gdf.to_crs("EPSG:4326")
        gdf.to_file(path, driver="GeoJSON")
    elif format == "gpkg":
        if path.exists():
            path.unlink()  # Remove existing to avoid stale layers
        gdf.to_file(path, driver="GPKG", layer="fields")
    elif format == "shp":
        gdf.to_file(path, driver="ESRI Shapefile")
    elif format == "fgb":
        gdf.to_file(path, driver="FlatGeobuf")
    else:
        raise ValueError(f"Unsupported output format: {format!r}")

    return str(path)