Composites

Composite builders create the stage-A raster: Earth Engine composites, USGS NAIP Plus exports, local GeoTIFFs and embedding rasters. See Satellite sources.

composites

Composite builders for satellite imagery acquisition.

Handles creating annual composites from various satellite sources via Google Earth Engine or loading local GeoTIFF files.

CompositeBuilder

Bases: ABC

Abstract base class for satellite composite builders.

Subclasses must implement :meth:build and :meth:get_band_mapping.

Source code in agribound/composites/base.py

class CompositeBuilder(ABC):
    """Abstract base class for satellite composite builders.

    Subclasses must implement :meth:`build` and :meth:`get_band_mapping`.
    """

    @abstractmethod
    def build(self, config: AgriboundConfig) -> str:
        """Build a composite (or fetch embeddings) and write it as a local GeoTIFF.

        Parameters
        ----------
        config : AgriboundConfig
            Pipeline configuration.

        Returns
        -------
        str
            Path to the composite GeoTIFF.
        """

    @abstractmethod
    def get_band_mapping(self, source: str) -> dict[str, str]:
        """Return the band name mapping for a satellite source.

        Parameters
        ----------
        source : str
            Satellite source name.

        Returns
        -------
        dict[str, str]
            Mapping of canonical names (R, G, B, NIR, ...) to source band names.
        """

    def get_resolution(self, source: str) -> float | None:
        """Return the default export resolution in metres for a source.

        Parameters
        ----------
        source : str
            Satellite source name.

        Returns
        -------
        float or None
            Resolution in metres, or *None* when it depends on the input.
        """
        return SOURCE_REGISTRY.get(source, {}).get("resolution_m")

build abstractmethod

build(config: AgriboundConfig) -> str

Build a composite (or fetch embeddings) and write it as a local GeoTIFF.

Parameters:

Name	Type	Description	Default
config	AgriboundConfig	Pipeline configuration.	required

Returns:

Type	Description
str	Path to the composite GeoTIFF.

Source code in agribound/composites/base.py

@abstractmethod
def build(self, config: AgriboundConfig) -> str:
    """Build a composite (or fetch embeddings) and write it as a local GeoTIFF.

    Parameters
    ----------
    config : AgriboundConfig
        Pipeline configuration.

    Returns
    -------
    str
        Path to the composite GeoTIFF.
    """

get_band_mapping abstractmethod

get_band_mapping(source: str) -> dict[str, str]

Return the band name mapping for a satellite source.

Parameters:

Name	Type	Description	Default
source	str	Satellite source name.	required

Returns:

Type	Description
dict[str, str]	Mapping of canonical names (R, G, B, NIR, ...) to source band names.

Source code in agribound/composites/base.py

@abstractmethod
def get_band_mapping(self, source: str) -> dict[str, str]:
    """Return the band name mapping for a satellite source.

    Parameters
    ----------
    source : str
        Satellite source name.

    Returns
    -------
    dict[str, str]
        Mapping of canonical names (R, G, B, NIR, ...) to source band names.
    """

get_resolution

get_resolution(source: str) -> float | None

Return the default export resolution in metres for a source.

Parameters:

Name	Type	Description	Default
source	str	Satellite source name.	required

Returns:

Type	Description
float or None	Resolution in metres, or None when it depends on the input.

Source code in agribound/composites/base.py

def get_resolution(self, source: str) -> float | None:
    """Return the default export resolution in metres for a source.

    Parameters
    ----------
    source : str
        Satellite source name.

    Returns
    -------
    float or None
        Resolution in metres, or *None* when it depends on the input.
    """
    return SOURCE_REGISTRY.get(source, {}).get("resolution_m")

NoDataError

Bases: ValueError

No input imagery or embeddings for this study area and year.

Raised by the composite builders when the source has no data for the study area in the requested period: no image intersects the study-area extent (the message lists the years that have images there), the downloaded composite or embedding has no valid pixel inside the study area, TESSERA or the Source Cooperative mirror returns no tile, USGS NAIP Plus has no imagery there, or a local raster does not overlap the study area. Invalid configurations and service errors raise other exceptions. A subclass of :class:ValueError, so except ValueError handlers keep catching it; :func:agribound.hpc.tiles.no_data_reason uses it to record such tiles as "no-data" instead of failed.

Source code in agribound/composites/base.py

class NoDataError(ValueError):
    """No input imagery or embeddings for this study area and year.

    Raised by the composite builders when the source has no data for the
    study area in the requested period: no image intersects the study-area
    extent (the message lists the years that have images there), the
    downloaded composite or embedding has no valid pixel inside the study
    area, TESSERA or the Source Cooperative mirror returns no tile, USGS NAIP
    Plus has no imagery there, or a local raster does not overlap the study
    area. Invalid configurations and service errors raise other exceptions.
    A subclass of :class:`ValueError`, so ``except ValueError`` handlers
    keep catching it; :func:`agribound.hpc.tiles.no_data_reason` uses it to
    record such tiles as "no-data" instead of failed.
    """

get_composite_builder

get_composite_builder(source: str) -> CompositeBuilder

Factory function to get the composite builder for a source.

Parameters:

Name	Type	Description	Default
source	str	Satellite source name.	required

Returns:

Type	Description
CompositeBuilder	Builder instance for the given source.

Raises:

Type	Description
ValueError	If the source is not recognised.

Source code in agribound/composites/base.py

def get_composite_builder(source: str) -> CompositeBuilder:
    """Factory function to get the composite builder for a source.

    Parameters
    ----------
    source : str
        Satellite source name.

    Returns
    -------
    CompositeBuilder
        Builder instance for the given source.

    Raises
    ------
    ValueError
        If the source is not recognised.
    """
    import importlib

    if source not in SOURCE_REGISTRY:
        raise ValueError(f"Unknown source {source!r}. Available: {list(SOURCE_REGISTRY.keys())}")

    target = _BUILDERS.get(source, _DEFAULT_BUILDER)
    module_name, _, class_name = target.partition(":")
    module = importlib.import_module(module_name)
    return getattr(module, class_name)()

list_sources

list_sources() -> dict[str, dict[str, Any]]

List all satellite sources and their metadata.

Returns:

Type	Description
dict[str, dict]	Deep copy of :data:SOURCE_REGISTRY (safe to mutate).

Examples:

>>> from agribound import list_sources
>>> for name, info in list_sources().items():
...     print(name, info["resolution_m"], info["value_scale"])

Source code in agribound/registry.py

def list_sources() -> dict[str, dict[str, Any]]:
    """List all satellite sources and their metadata.

    Returns
    -------
    dict[str, dict]
        Deep copy of :data:`SOURCE_REGISTRY` (safe to mutate).

    Examples
    --------
    >>> from agribound import list_sources
    >>> for name, info in list_sources().items():
    ...     print(name, info["resolution_m"], info["value_scale"])
    """
    return copy.deepcopy(SOURCE_REGISTRY)

Earth Engine

gee

Google Earth Engine composite builder and Earth Engine download helpers.

:class:GEECompositeBuilder builds a composite for landsat, landsat-pan, sentinel2, hls, naip, spot and spot-pan on Google Earth Engine and downloads it as one local GeoTIFF with geedim 2 (the ee.Image.gd accessor).

Radiometry of the written composites (value_scale in :mod:agribound.registry):

• Sentinel-2 (COPERNICUS/S2_SR_HARMONIZED): surface reflectance x 10000 as stored in the collection (12 bands B1-B12 without B10).
• Landsat 5/7/8/9 Collection 2 Level-2: every image is converted to surface reflectance (DN * 2.75e-5 - 0.2), clipped at 0 and multiplied by 10000 before compositing. Landsat 5/7 bands are renamed to the Landsat 8/9 names (SR_B2 ... SR_B7 = blue, green, red, NIR, SWIR1, SWIR2).
• HLS v2.0: Earth Engine stores HLS as 0-1 reflectance; values are multiplied by 10000. HLSS30 bands B1, B2, B3, B4, B8A, B11, B12 are renamed to the HLSL30 names B1 ... B7 (so B5 is NIR narrow and B6/B7 are SWIR1/SWIR2 for both sensors).
• Landsat PAN: native 15 m B8 TOA reflectance from Landsat 7/8/9, exported as float32 without scaling (unit). The two PAN bandpasses (Landsat 7 0.52-0.90 um, Landsat 8/9 0.50-0.68 um) are not mixed unless landsat_pan_missions names both (see :func:landsat_pan_missions_for_window); AGRIBOUND_SENSORS and AGRIBOUND_SENSOR_IMAGES record the missions whose images were used.
• NAIP (USDA/NAIP/DOQQ): 8-bit digital numbers (R, G, B, N), mosaicked (not composited), exported as uint8 with nodata 0 at config.naip_resolution_m. Only 4-band images are used (some early years are RGB only). Without date_range, images from year - 1 to year + 1 are mosaicked with the exact-year images on top and the newest image on top within each group.
• SPOT 6/7 (AIRBUS/SPOT6_7, restricted): per-band medians (or greenest-pixel selections) of the scenes' raw integer digital numbers, exported as float32 (a median of an even number of scenes can be a half-integer). The radiometry of these DN has not been verified, so they are treated as uncalibrated (value scale "dn").

Float composites are written as float32 with NaN nodata. Masked pixels (clouds, cloud shadows, no valid observation) are NaN. After download, the share of valid pixels inside the study-area polygons is written to the AGRIBOUND_VALID_FRACTION tag; the builder raises when it is 0 and logs a WARNING when it is below :data:LOW_VALID_FRACTION.

Cloud masks: Landsat QA_PIXEL bits 0-4 (fill, dilated cloud, cirrus, cloud, cloud shadow); HLS Fmask bits 1-3 (cloud, adjacent to cloud/shadow, cloud shadow); Sentinel-2 either the scene classification (SCL classes 3 cloud shadow, 8 cloud medium probability, 9 cloud high probability, 10 thin cirrus) or Cloud Score+ (cs_cdf >= config.cloud_score_threshold). Scenes are pre-filtered on their scene cloud-cover property (<= config.cloud_cover_max).

Extent: the export pixel grid covers the bounding box of the study area in config.export_crs ("utm" = WGS 84 / UTM zone of the study-area centroid, used for the whole study area), with pixel edges on multiples of the resolution, so composites of the same study area, CRS and resolution (e.g. FTW's two date windows) share one grid. Images are selected with filterBounds on the outline of that grid. Pixels are not masked to the study-area polygons: when the study area is a set of field polygons (for example the reference boundaries), masking would imprint their outlines on the engine input.

Download: component images are resampled to the grid by Earth Engine (nearest neighbour, the Earth Engine default). Grids larger than config.tile_size pixels per side are split into aligned tiles that are downloaded one after another (each tile by geedim with up to config.gee_max_requests concurrent requests), cached, and assembled into one GeoTIFF.

GEECompositeBuilder

Bases: CompositeBuilder

Composite builder for the Earth Engine imagery sources.

Handles landsat, landsat-pan, sentinel2, hls, naip, spot and spot-pan (see the module docstring for masks, radiometry and the export grid). The result is cached under cache_path(config, f"{source}_composite", ".tif", ...), whose key includes the study area, year, date range, compositing and export settings, so different windows (e.g. FTW's two windows) get different files.

Attributes:

Name	Type	Description
last_metadata	dict	Facts about the most recent :meth:build (also written as AGRIBOUND_* GeoTIFF tags).

Source code in agribound/composites/gee.py

class GEECompositeBuilder(CompositeBuilder):
    """Composite builder for the Earth Engine imagery sources.

    Handles ``landsat``, ``landsat-pan``, ``sentinel2``, ``hls``, ``naip``,
    ``spot`` and ``spot-pan`` (see the module docstring for masks, radiometry
    and the export grid). The result is cached under
    ``cache_path(config, f"{source}_composite", ".tif", ...)``, whose key
    includes the study area, year, date range, compositing and export
    settings, so different windows (e.g. FTW's two windows) get different
    files.

    Attributes
    ----------
    last_metadata : dict
        Facts about the most recent :meth:`build` (also written as
        ``AGRIBOUND_*`` GeoTIFF tags).
    """

    def __init__(self) -> None:
        self.last_metadata: dict[str, Any] = {}

    def build(self, config: AgriboundConfig) -> str:
        """Build and download the composite for *config*.

        Parameters
        ----------
        config : AgriboundConfig
            Pipeline configuration.

        Returns
        -------
        str
            Path to the composite GeoTIFF.

        Raises
        ------
        NoDataError
            (A :class:`ValueError`.) If no image matches the filters (the
            message lists the years that have images over the study-area
            extent), or the downloaded composite has no valid pixel inside the
            study area (the file is then deleted).
        ValueError
            If the composite method does not apply to the source.
        ExportTaskStartedError
            For ``export_method="gdrive"``/``"gcs"``: after starting the batch
            export task, or without starting one when the task recorded by an
            earlier run (``<composite stem>_task.json``) is queued, running or
            completed.

        Notes
        -----
        After a local download, the share of pixels inside the study-area
        polygons that are valid in every band is written to the
        ``AGRIBOUND_VALID_FRACTION`` tag (:func:`raster_valid_fraction`); a
        WARNING is logged when it is below :data:`LOW_VALID_FRACTION`.
        """
        from agribound._cache import cache_key, cache_path
        from agribound.auth import ensure_gee

        if config.source not in _COLLECTION_BUILDERS:
            raise ValueError(f"No Earth Engine collection builder for source {config.source!r}")
        if config.source == "naip" and config.composite_method != "median":
            raise ValueError(
                "NAIP is mosaicked (one image per location), not composited, so "
                f"composite_method={config.composite_method!r} does not apply. Use the default "
                "composite_method='median' (ignored for NAIP)."
            )

        stem = f"{config.source}_composite"
        out_path = cache_path(config, stem, ".tif", COMPOSITE_RECIPE_VERSION)
        if out_path.exists():
            logger.info("Using cached composite: %s", out_path)
            self.last_metadata = read_composite_tags(out_path)
            cached_fraction = self.last_metadata.get("AGRIBOUND_VALID_FRACTION")
            if cached_fraction is not None:
                _warn_low_coverage(
                    float(cached_fraction),
                    f"The cached {config.source} composite",
                    _coverage_advice(config.source),
                )
            return str(out_path)

        ensure_gee(config)
        geom = study_area_geometry_4326(config)
        crs = resolve_export_crs(config.export_crs, geom)
        grid = compute_export_grid(geom, crs, export_resolution_m(config))
        # Select images over the whole export grid (the study-area bounding box in
        # the export CRS); the composite is not masked to the study-area polygons.
        region = ee_geometry(grid_footprint_4326(grid))
        spec = _COLLECTION_BUILDERS[config.source](config, region)
        window = date_window(config)

        label = f"{config.source} {config.year}"
        summary = collection_summary(
            spec.raw, context=f"{label} image count", count_property=spec.count_property
        )
        if summary["n_images"] == 0:
            years = available_years(spec.history, context=f"{label} available years")
            if config.source == "landsat-pan":
                # The history covers the missions searched, so the years are labelled with them.
                searched = spec.notes["missions_selected"]
                hint = _landsat7_hint(config, region, window, searched.split(","))
                error = NoDataError(
                    f"{_years_error(config, years, window, of=searched)} Missions searched: "
                    f"{searched} (landsat_pan_missions={spec.notes['landsat_pan_missions']!r})."
                    f"{hint}"
                )
            else:
                error = _years_error(config, years, window)
            if config.source == "naip":
                all_years = available_years(spec.history_all_bands, context="naip available years")
                rgb_only = sorted(set(all_years) - set(years))
                if rgb_only:
                    error = NoDataError(f"{error} Years with RGB-only NAIP: {rgb_only}.")
            raise error
        logger.info(
            "%s: %d image(s) in %s..%s from years %s",
            label,
            summary["n_images"],
            window[0],
            window[1],
            summary["years"],
        )
        if config.source == "naip":
            year = int(config.year)
            other = [y for y in summary["years"] if y != year]
            if config.date_range is None and year not in summary["years"]:
                logger.warning(
                    "No 4-band NAIP image from %d intersects the study-area extent; the mosaic "
                    "uses images from %s (newest on top).",
                    year,
                    other,
                )
            elif config.date_range is None and other:
                logger.warning(
                    "NAIP mosaic for %d: %d images are on top; images from %s fill gaps.",
                    year,
                    year,
                    other,
                )
            composite = spec.prepared.mosaic()
            method = "mosaic"
        else:
            method = config.composite_method
            composite = apply_composite_method(spec.prepared, method, config.source)
            if method == "max_ndvi":
                method = "greenest"  # documented alias
        from agribound._version import __version__

        tags = {
            "AGRIBOUND_VERSION": __version__,
            "AGRIBOUND_SOURCE": config.source,
            "AGRIBOUND_VALUE_SCALE": SOURCE_REGISTRY[config.source]["value_scale"],
            "AGRIBOUND_YEAR": config.year,
            "AGRIBOUND_DATE_START": window[0],
            "AGRIBOUND_DATE_END_EXCLUSIVE": window[1],
            "AGRIBOUND_COMPOSITE_METHOD": method,
            "AGRIBOUND_N_IMAGES": summary["n_images"],
            "AGRIBOUND_IMAGE_YEARS": ",".join(str(y) for y in summary["years"]),
            "AGRIBOUND_COLLECTIONS": ",".join(spec.collections),
            "AGRIBOUND_CLOUD_COVER_MAX": config.cloud_cover_max,
            "AGRIBOUND_EXPORT_CRS": crs,
            "AGRIBOUND_RESOLUTION_M": spec.resolution_m,
            "AGRIBOUND_CACHE_KEY": cache_key(config, COMPOSITE_RECIPE_VERSION),
            "AGRIBOUND_RECIPE_VERSION": COMPOSITE_RECIPE_VERSION,
        }
        for key, value in spec.notes.items():
            if isinstance(value, str | int | float):
                tags[f"AGRIBOUND_{key.upper()}"] = value
        if spec.count_property:
            tags.update(_sensor_tags(config, label, summary.get("counts", {}), spec.count_labels))
        self.last_metadata = dict(tags)

        if config.export_method != "local":
            self._batch_export(config, composite, spec, grid, out_path)

        export_ee_image(
            composite,
            out_path,
            grid=grid,
            dtype=spec.dtype,
            band_names=spec.bands,
            max_requests=config.gee_max_requests,
            tile_size=config.tile_size,
            tags=tags,
            label=f"{label} composite",
        )
        require = "all" if np.issubdtype(np.dtype(spec.dtype), np.floating) else "any"
        fraction = raster_valid_fraction(out_path, geom, require=require)
        if fraction <= 0.0:
            out_path.unlink(missing_ok=True)
            advice = _coverage_advice(config.source)
            if config.source == "sentinel2" and config.s2_cloud_mask == "cloud_score_plus":
                advice += " A lower cloud_score_threshold keeps more pixels."
            raise NoDataError(
                f"The {label} composite has no valid pixel inside the study area although "
                f"{summary['n_images']} image(s) matched the filters. {advice}"
            )
        import rasterio

        with rasterio.open(out_path, "r+") as dst:
            dst.update_tags(AGRIBOUND_VALID_FRACTION=str(round(fraction, 6)))
        tags["AGRIBOUND_VALID_FRACTION"] = round(fraction, 6)
        self.last_metadata = dict(tags)
        _warn_low_coverage(fraction, f"The {label} composite", _coverage_advice(config.source))
        return str(out_path)

    def _batch_export(
        self,
        config: AgriboundConfig,
        composite: Any,
        spec: CollectionSpec,
        grid: ExportGrid,
        out_path: Path,
    ) -> None:
        """Start (or find) the Drive/GCS batch export and raise ExportTaskStartedError."""
        import json

        marker = export_task_marker(out_path)
        previous = recorded_export_task(marker)
        if previous is not None:
            state = previous["state"]
            if state == "COMPLETED":
                what = (
                    f"Earth Engine export task {previous['id']} (started by an earlier run) has "
                    f"completed -> {previous.get('destination')}. Download the GeoTIFF and run "
                    "with source='local' (local_tif_path=...)"
                )
            else:
                what = (
                    f"Earth Engine export task {previous['id']} (started by an earlier run) is "
                    f"{state} -> {previous.get('destination')}; no new task was started. When it "
                    "has finished, download the GeoTIFF and run with source='local' "
                    "(local_tif_path=...)"
                )
            raise ExportTaskStartedError(
                f"{what}, or use export_method='local' to download directly. Delete {marker} "
                "to start a new export task.",
                [previous],
            )
        task = start_batch_export(
            composite,
            grid=grid,
            dtype=spec.dtype,
            band_names=spec.bands,
            method=config.export_method,
            description=out_path.stem,
            gcs_bucket=config.gcs_bucket,
        )
        record = {
            **task,
            "started_utc": _dt.datetime.now(_dt.UTC).isoformat(timespec="seconds"),
            "crs": grid.crs,
            "crs_transform": grid.crs_transform,
            "width": grid.width,
            "height": grid.height,
        }
        marker.write_text(json.dumps(record, indent=2), encoding="utf-8")
        raise ExportTaskStartedError(
            f"Started Earth Engine export task {task['id']} -> {task['destination']} "
            f"({grid.width} x {grid.height} px, {grid.crs}). The composite is not "
            "available locally until the task finishes; download the GeoTIFF and run "
            "with source='local' (local_tif_path=...), or use export_method='local' to "
            f"download directly. The task is recorded in {marker}, so rerunning does not "
            "start a second export.",
            [task],
        )

    def get_band_mapping(self, source: str) -> dict[str, str]:
        """Return the canonical band mapping (canonical name -> band name) of a source."""
        info = SOURCE_REGISTRY.get(source, {})
        return dict(info.get("canonical_bands") or {})

build

build(config: AgriboundConfig) -> str

Build and download the composite for config.

Parameters:

Name	Type	Description	Default
config	AgriboundConfig	Pipeline configuration.	required

Returns:

Type	Description
str	Path to the composite GeoTIFF.

Raises:

Type	Description
NoDataError	(A :class:ValueError.) If no image matches the filters (the message lists the years that have images over the study-area extent), or the downloaded composite has no valid pixel inside the study area (the file is then deleted).
ValueError	If the composite method does not apply to the source.
ExportTaskStartedError	For export_method="gdrive"/"gcs": after starting the batch export task, or without starting one when the task recorded by an earlier run (<composite stem>_task.json) is queued, running or completed.

Notes

After a local download, the share of pixels inside the study-area polygons that are valid in every band is written to the AGRIBOUND_VALID_FRACTION tag (:func:raster_valid_fraction); a WARNING is logged when it is below :data:LOW_VALID_FRACTION.

Source code in agribound/composites/gee.py

def build(self, config: AgriboundConfig) -> str:
    """Build and download the composite for *config*.

    Parameters
    ----------
    config : AgriboundConfig
        Pipeline configuration.

    Returns
    -------
    str
        Path to the composite GeoTIFF.

    Raises
    ------
    NoDataError
        (A :class:`ValueError`.) If no image matches the filters (the
        message lists the years that have images over the study-area
        extent), or the downloaded composite has no valid pixel inside the
        study area (the file is then deleted).
    ValueError
        If the composite method does not apply to the source.
    ExportTaskStartedError
        For ``export_method="gdrive"``/``"gcs"``: after starting the batch
        export task, or without starting one when the task recorded by an
        earlier run (``<composite stem>_task.json``) is queued, running or
        completed.

    Notes
    -----
    After a local download, the share of pixels inside the study-area
    polygons that are valid in every band is written to the
    ``AGRIBOUND_VALID_FRACTION`` tag (:func:`raster_valid_fraction`); a
    WARNING is logged when it is below :data:`LOW_VALID_FRACTION`.
    """
    from agribound._cache import cache_key, cache_path
    from agribound.auth import ensure_gee

    if config.source not in _COLLECTION_BUILDERS:
        raise ValueError(f"No Earth Engine collection builder for source {config.source!r}")
    if config.source == "naip" and config.composite_method != "median":
        raise ValueError(
            "NAIP is mosaicked (one image per location), not composited, so "
            f"composite_method={config.composite_method!r} does not apply. Use the default "
            "composite_method='median' (ignored for NAIP)."
        )

    stem = f"{config.source}_composite"
    out_path = cache_path(config, stem, ".tif", COMPOSITE_RECIPE_VERSION)
    if out_path.exists():
        logger.info("Using cached composite: %s", out_path)
        self.last_metadata = read_composite_tags(out_path)
        cached_fraction = self.last_metadata.get("AGRIBOUND_VALID_FRACTION")
        if cached_fraction is not None:
            _warn_low_coverage(
                float(cached_fraction),
                f"The cached {config.source} composite",
                _coverage_advice(config.source),
            )
        return str(out_path)

    ensure_gee(config)
    geom = study_area_geometry_4326(config)
    crs = resolve_export_crs(config.export_crs, geom)
    grid = compute_export_grid(geom, crs, export_resolution_m(config))
    # Select images over the whole export grid (the study-area bounding box in
    # the export CRS); the composite is not masked to the study-area polygons.
    region = ee_geometry(grid_footprint_4326(grid))
    spec = _COLLECTION_BUILDERS[config.source](config, region)
    window = date_window(config)

    label = f"{config.source} {config.year}"
    summary = collection_summary(
        spec.raw, context=f"{label} image count", count_property=spec.count_property
    )
    if summary["n_images"] == 0:
        years = available_years(spec.history, context=f"{label} available years")
        if config.source == "landsat-pan":
            # The history covers the missions searched, so the years are labelled with them.
            searched = spec.notes["missions_selected"]
            hint = _landsat7_hint(config, region, window, searched.split(","))
            error = NoDataError(
                f"{_years_error(config, years, window, of=searched)} Missions searched: "
                f"{searched} (landsat_pan_missions={spec.notes['landsat_pan_missions']!r})."
                f"{hint}"
            )
        else:
            error = _years_error(config, years, window)
        if config.source == "naip":
            all_years = available_years(spec.history_all_bands, context="naip available years")
            rgb_only = sorted(set(all_years) - set(years))
            if rgb_only:
                error = NoDataError(f"{error} Years with RGB-only NAIP: {rgb_only}.")
        raise error
    logger.info(
        "%s: %d image(s) in %s..%s from years %s",
        label,
        summary["n_images"],
        window[0],
        window[1],
        summary["years"],
    )
    if config.source == "naip":
        year = int(config.year)
        other = [y for y in summary["years"] if y != year]
        if config.date_range is None and year not in summary["years"]:
            logger.warning(
                "No 4-band NAIP image from %d intersects the study-area extent; the mosaic "
                "uses images from %s (newest on top).",
                year,
                other,
            )
        elif config.date_range is None and other:
            logger.warning(
                "NAIP mosaic for %d: %d images are on top; images from %s fill gaps.",
                year,
                year,
                other,
            )
        composite = spec.prepared.mosaic()
        method = "mosaic"
    else:
        method = config.composite_method
        composite = apply_composite_method(spec.prepared, method, config.source)
        if method == "max_ndvi":
            method = "greenest"  # documented alias
    from agribound._version import __version__

    tags = {
        "AGRIBOUND_VERSION": __version__,
        "AGRIBOUND_SOURCE": config.source,
        "AGRIBOUND_VALUE_SCALE": SOURCE_REGISTRY[config.source]["value_scale"],
        "AGRIBOUND_YEAR": config.year,
        "AGRIBOUND_DATE_START": window[0],
        "AGRIBOUND_DATE_END_EXCLUSIVE": window[1],
        "AGRIBOUND_COMPOSITE_METHOD": method,
        "AGRIBOUND_N_IMAGES": summary["n_images"],
        "AGRIBOUND_IMAGE_YEARS": ",".join(str(y) for y in summary["years"]),
        "AGRIBOUND_COLLECTIONS": ",".join(spec.collections),
        "AGRIBOUND_CLOUD_COVER_MAX": config.cloud_cover_max,
        "AGRIBOUND_EXPORT_CRS": crs,
        "AGRIBOUND_RESOLUTION_M": spec.resolution_m,
        "AGRIBOUND_CACHE_KEY": cache_key(config, COMPOSITE_RECIPE_VERSION),
        "AGRIBOUND_RECIPE_VERSION": COMPOSITE_RECIPE_VERSION,
    }
    for key, value in spec.notes.items():
        if isinstance(value, str | int | float):
            tags[f"AGRIBOUND_{key.upper()}"] = value
    if spec.count_property:
        tags.update(_sensor_tags(config, label, summary.get("counts", {}), spec.count_labels))
    self.last_metadata = dict(tags)

    if config.export_method != "local":
        self._batch_export(config, composite, spec, grid, out_path)

    export_ee_image(
        composite,
        out_path,
        grid=grid,
        dtype=spec.dtype,
        band_names=spec.bands,
        max_requests=config.gee_max_requests,
        tile_size=config.tile_size,
        tags=tags,
        label=f"{label} composite",
    )
    require = "all" if np.issubdtype(np.dtype(spec.dtype), np.floating) else "any"
    fraction = raster_valid_fraction(out_path, geom, require=require)
    if fraction <= 0.0:
        out_path.unlink(missing_ok=True)
        advice = _coverage_advice(config.source)
        if config.source == "sentinel2" and config.s2_cloud_mask == "cloud_score_plus":
            advice += " A lower cloud_score_threshold keeps more pixels."
        raise NoDataError(
            f"The {label} composite has no valid pixel inside the study area although "
            f"{summary['n_images']} image(s) matched the filters. {advice}"
        )
    import rasterio

    with rasterio.open(out_path, "r+") as dst:
        dst.update_tags(AGRIBOUND_VALID_FRACTION=str(round(fraction, 6)))
    tags["AGRIBOUND_VALID_FRACTION"] = round(fraction, 6)
    self.last_metadata = dict(tags)
    _warn_low_coverage(fraction, f"The {label} composite", _coverage_advice(config.source))
    return str(out_path)

get_band_mapping

get_band_mapping(source: str) -> dict[str, str]

Return the canonical band mapping (canonical name -> band name) of a source.

Source code in agribound/composites/gee.py

def get_band_mapping(self, source: str) -> dict[str, str]:
    """Return the canonical band mapping (canonical name -> band name) of a source."""
    info = SOURCE_REGISTRY.get(source, {})
    return dict(info.get("canonical_bands") or {})

ExportTaskStartedError

Bases: RuntimeError

Raised when the composite is exported by an Earth Engine batch task (Drive/GCS).

Raised after a task was started, or when a task recorded by an earlier run is still queued, running, or has completed (no second task is started then). The composite is not available locally until the task finishes, so the pipeline cannot continue. tasks lists {"id", "description", "destination"} (plus "state" for a recorded task) for each task.

Source code in agribound/composites/gee.py

class ExportTaskStartedError(RuntimeError):
    """Raised when the composite is exported by an Earth Engine batch task (Drive/GCS).

    Raised after a task was started, or when a task recorded by an earlier
    run is still queued, running, or has completed (no second task is
    started then). The composite is not available locally until the task
    finishes, so the pipeline cannot continue. ``tasks`` lists ``{"id",
    "description", "destination"}`` (plus ``"state"`` for a recorded task)
    for each task.
    """

    def __init__(self, message: str, tasks: list[dict[str, str]]):
        super().__init__(message)
        self.tasks = tasks

apply_composite_method

apply_composite_method(collection: Any, method: str, source: str) -> Any

Reduce a prepared collection to one image.

Parameters:

Name	Type	Description	Default
collection	ImageCollection	Masked, scaled images with the source's output bands.	required
method	str	"median": per-band median of the unmasked values. "greenest": per pixel, the image with the highest NDVI (qualityMosaic) computed from the canonical NIR and R bands. "max_ndvi" is an alias of "greenest".	required
source	str	Source name (to look up the NIR/red bands).	required

Raises:

Type	Description
ValueError	For an unknown method, or greenest/max_ndvi on a source without NIR and red bands.

Source code in agribound/composites/gee.py

def apply_composite_method(collection: Any, method: str, source: str) -> Any:
    """Reduce a prepared collection to one image.

    Parameters
    ----------
    collection : ee.ImageCollection
        Masked, scaled images with the source's output bands.
    method : str
        ``"median"``: per-band median of the unmasked values.
        ``"greenest"``: per pixel, the image with the highest NDVI
        (``qualityMosaic``) computed from the canonical ``NIR`` and ``R``
        bands. ``"max_ndvi"`` is an alias of ``"greenest"``.
    source : str
        Source name (to look up the NIR/red bands).

    Raises
    ------
    ValueError
        For an unknown method, or ``greenest``/``max_ndvi`` on a source
        without NIR and red bands.
    """
    if method == "median":
        return collection.median()
    if method in ("greenest", "max_ndvi"):
        canonical = SOURCE_REGISTRY[source].get("canonical_bands") or {}
        nir, red = canonical.get("NIR"), canonical.get("R")
        if not nir or not red or nir == red:
            raise ValueError(
                f"composite_method={method!r} needs NIR and red bands, which source "
                f"{source!r} does not have. Use composite_method='median'."
            )
        bands = list(SOURCE_REGISTRY[source]["all_bands"])

        def add_ndvi(image):
            return image.addBands(image.normalizedDifference([nir, red]).rename("NDVI"))

        return collection.map(add_ndvi).qualityMosaic("NDVI").select(bands)
    raise ValueError(f"Unknown composite_method {method!r}")

date_window

date_window(config: AgriboundConfig) -> tuple[str, str]

Return (start, end_exclusive) dates for ee.ImageCollection.filterDate.

config.date_range is inclusive of its end date (as in the rest of Agribound), so one day is added to it. Without a date range the window is the calendar year [year-01-01, (year+1)-01-01).

Source code in agribound/composites/gee.py

def date_window(config: AgriboundConfig) -> tuple[str, str]:
    """Return ``(start, end_exclusive)`` dates for ``ee.ImageCollection.filterDate``.

    ``config.date_range`` is inclusive of its end date (as in the rest of
    Agribound), so one day is added to it. Without a date range the window is
    the calendar year ``[year-01-01, (year+1)-01-01)``.
    """
    if config.date_range is not None:
        start = _dt.date.fromisoformat(str(config.date_range[0]))
        end = _dt.date.fromisoformat(str(config.date_range[1])) + _dt.timedelta(days=1)
        return start.isoformat(), end.isoformat()
    return f"{int(config.year)}-01-01", f"{int(config.year) + 1}-01-01"

landsat_pan_missions_for_window

landsat_pan_missions_for_window(setting: str | tuple[str, ...], window: tuple[str, str]) -> tuple[str, ...]

Return the missions a landsat-pan composite of window uses.

Parameters:

Name	Type	Description	Default
setting	str or tuple[str, ...]	config.landsat_pan_missions: "auto" or mission IDs.	required
window	tuple[str, str]	(start, end_exclusive) dates (:func:date_window).	required

Returns:

Type	Description
tuple[str, ...]	With "auto", the OLI missions (LC08, LC09) whose acquisition period overlaps the window, else LE07 if its period does, so the two PAN bandpasses are never mixed. With mission IDs, those of them whose period overlaps the window. Empty when no mission applies.

Source code in agribound/composites/gee.py

def landsat_pan_missions_for_window(
    setting: str | tuple[str, ...], window: tuple[str, str]
) -> tuple[str, ...]:
    """Return the missions a ``landsat-pan`` composite of *window* uses.

    Parameters
    ----------
    setting : str or tuple[str, ...]
        ``config.landsat_pan_missions``: ``"auto"`` or mission IDs.
    window : tuple[str, str]
        ``(start, end_exclusive)`` dates (:func:`date_window`).

    Returns
    -------
    tuple[str, ...]
        With ``"auto"``, the OLI missions (``LC08``, ``LC09``) whose
        acquisition period overlaps the window, else ``LE07`` if its period
        does, so the two PAN bandpasses are never mixed. With mission IDs,
        those of them whose period overlaps the window. Empty when no mission
        applies.
    """
    overlapping = [
        mission
        for mission, (_, first, last) in LANDSAT_PAN_COLLECTIONS.items()
        if _overlaps(window, first, last)
    ]
    if setting == "auto":
        oli = [m for m in overlapping if m in LANDSAT_PAN_OLI_MISSIONS]
        return tuple(oli or [m for m in overlapping if m == "LE07"])
    return tuple(m for m in overlapping if m in setting)

USGS NAIP Plus

usgs

USGS NAIP Plus ImageServer composite builder.

The USGS NAIP Plus ImageServer serves only the latest NAIP/HRO vintage of each state (years 2012-2023 as of 2026-09, most states 2019-2023), not the historical NAIP archive, so most years are unavailable for a given state. Use source="naip" (Earth Engine, 2002-2023) for other years. The service's pixel size is 0.3 m; state vintages are 0.3-0.6 m.

The output covers the export grid: the bounding box of the study area in config.export_crs ("utm": the UTM zone of the study-area centroid), as for the other builders; pixels are not masked to the study-area polygons.

:class:USGSNAIPPlusCompositeBuilder:

1. queries catalogue items (Category = 1, Year = config.year, optional State) intersecting the export-grid outline; when none match it raises :class:ValueError listing the years the service has for the area;
2. selects footprints greedily (closest year, 4 bands, largest overlap, finest resolution) until the grid outline is covered or the service's maxMosaicImageCount is reached, and locks the mosaic to them (esriMosaicLockRaster); a WARNING is logged when the selected footprints cover less than 99 % of the grid outline;
3. exports tiles of at most min(config.tile_size, maxImageWidth, maxImageHeight) pixels in EPSG:3857 at the finest ground resolution of the selected footprints (the service resamples with bilinear interpolation), with the Web Mercator pixel size scaled by 1 / cos(latitude) of the study-area centroid so that the ground sampling matches;
4. merges the tiles and resamples them (bilinear) onto the export grid;
5. writes a uint8 GeoTIFF (nodata 0; pixels no selected footprint covers are 0) with a JSON manifest of the query and selection, and records the share of study-area pixels that have imagery (AGRIBOUND_VALID_FRACTION; a WARNING is logged below 99 %, and :class:ValueError is raised at 0).

The band order is assumed to be R, G, B, N (the service reports generic band names Band_1-Band_4 and sensor type CNIR).

USGSNAIPPlusCompositeBuilder

Bases: CompositeBuilder

Build a local uint8 raster for the study area from the USGS NAIP Plus ImageServer.

Attributes:

Name	Type	Description
SELECTION_VERSION	int	Version of the selection/export recipe (part of the cache key).
last_metadata	dict	Facts about the most recent :meth:build.

Source code in agribound/composites/usgs.py

class USGSNAIPPlusCompositeBuilder(CompositeBuilder):
    """Build a local uint8 raster for the study area from the USGS NAIP Plus ImageServer.

    Attributes
    ----------
    SELECTION_VERSION : int
        Version of the selection/export recipe (part of the cache key).
    last_metadata : dict
        Facts about the most recent :meth:`build`.
    """

    SELECTION_VERSION = 3

    def __init__(self) -> None:
        self.last_metadata: dict[str, Any] = {}

    def build(self, config: AgriboundConfig) -> str:
        """Download and mosaic a USGS NAIP Plus raster covering the study-area extent.

        Returns
        -------
        str
            Path to the uint8 GeoTIFF (cached with
            :func:`agribound._cache.cache_path`).

        Raises
        ------
        NoDataError
            (A :class:`ValueError`.) If the service has no imagery for the
            requested year (the message lists the years it has for the area),
            or no pixel inside the study area has imagery.
        ValueError
            If no study area is set or the study area is empty.
        """
        from agribound._cache import cache_path

        if not config.study_area:
            raise ValueError("study_area is required when source='usgs-naip-plus'")

        from agribound.composites.gee import (
            compute_export_grid,
            grid_footprint_4326,
            raster_valid_fraction,
            resolve_export_crs,
            split_grid,
            study_area_geometry_4326,
        )

        geom_4326 = study_area_geometry_4326(config)

        final_path = cache_path(config, "usgs_naip_plus", ".tif", self.SELECTION_VERSION)
        manifest_path = final_path.with_suffix(".json")
        if final_path.exists() and manifest_path.exists():
            logger.info("Using cached USGS NAIP Plus raster: %s", final_path)
            return str(final_path)

        client = USGSNAIPPlusClient(
            service_url=getattr(config, "usgs_service_url", DEFAULT_USGS_NAIPPLUS_URL),
            timeout_s=getattr(config, "usgs_timeout_s", 120),
            retries=getattr(config, "usgs_retries", 3),
        )

        service_metadata = client.get_service_metadata()
        # The query/selection target is the outline of a provisional export grid
        # (the study-area bounding box in the export CRS), in EPSG:3857.
        export_crs = resolve_export_crs(config.export_crs, geom_4326)
        provisional = compute_export_grid(geom_4326, export_crs, _DEFAULT_RESOLUTION_M)
        aoi_3857 = _to_3857(grid_footprint_4326(provisional))
        where = self._build_where_clause(config)
        candidates = client.query_candidates(aoi_3857.bounds, where, out_sr=3857)
        candidates = self._filter_candidates(candidates, aoi_3857)

        if not candidates and getattr(config, "usgs_allow_year_fallback", False):
            where = self._build_where_clause(config, allow_year_fallback=True)
            candidates = client.query_candidates(aoi_3857.bounds, where, out_sr=3857)
            candidates = self._filter_candidates(candidates, aoi_3857)

        if not candidates:
            raise self._no_imagery_error(client, config, aoi_3857)

        max_ids = int(service_metadata.get("maxMosaicImageCount", 50))
        selected_ids, ranked_candidates = self._select_lock_raster_ids(
            candidates,
            aoi_3857,
            target_year=int(config.year),
            max_ids=max_ids,
        )
        if not selected_ids:
            raise RuntimeError("Candidate selection returned no LockRaster IDs")
        selected = [c for c in ranked_candidates if c.object_id in set(selected_ids)]
        footprint_share = self._footprint_coverage(selected, aoi_3857)
        if footprint_share < LOW_COVERAGE_FRACTION:
            reason = (
                f"the selection reached the service's maxMosaicImageCount ({max_ids})"
                if len(selected_ids) >= max_ids
                else "no other candidate footprint adds coverage"
            )
            logger.warning(
                "The %d selected USGS NAIP Plus footprint(s) cover %.1f%% of the study-area "
                "extent (%s); uncovered pixels are 0 (nodata).",
                len(selected_ids),
                100.0 * footprint_share,
                reason,
            )

        export_max_px = min(
            int(getattr(config, "tile_size", 10_000)),
            int(service_metadata.get("maxImageWidth", 4000)),
            int(service_metadata.get("maxImageHeight", 4000)),
        )
        resolution_m = self._estimate_resolution_m(selected)
        centroid_lat = float(geom_4326.centroid.y)
        resolution_3857 = resolution_m / math.cos(math.radians(centroid_lat))

        grid = compute_export_grid(geom_4326, export_crs, resolution_m)
        grid_3857 = compute_export_grid(grid_footprint_4326(grid), "EPSG:3857", resolution_3857)
        windows = split_grid(grid_3857, export_max_px)
        tile_dir = final_path.with_name(final_path.stem + "_tiles")
        tile_dir.mkdir(parents=True, exist_ok=True)

        raw_tile_paths: list[str] = []
        for idx, (row_off, col_off, height, width) in enumerate(windows):
            tile_grid = grid_3857.window(row_off, col_off, height, width)
            t = tile_grid.transform
            bounds = (t.c, t.f + height * t.e, t.c + width * t.a, t.f)
            tile_path = tile_dir / f"tile_{idx:03d}.tif"
            if not tile_path.exists():
                logger.info(
                    "Exporting USGS NAIP Plus tile %d/%d (%d x %d px) with %d LockRaster IDs",
                    idx + 1,
                    len(windows),
                    width,
                    height,
                    len(selected_ids),
                )
                client.export_image(
                    bbox_3857=bounds,
                    width=width,
                    height=height,
                    lock_raster_ids=selected_ids,
                    output_path=tile_path,
                )
            self._validate_export(tile_path)
            raw_tile_paths.append(str(tile_path))

        band_count = self._write_final(raw_tile_paths, final_path, grid)
        if band_count < 4:
            logger.warning(
                "USGS NAIP Plus export has %d bands; the NIR band (N) is missing", band_count
            )
        # Uncovered pixels are 0 in every band; a dark pixel can be 0 in one band.
        valid_fraction = raster_valid_fraction(final_path, geom_4326, require="any")
        if valid_fraction <= 0.0:
            final_path.unlink(missing_ok=True)
            raise NoDataError(
                f"The USGS NAIP Plus export for {config.year} has no imagery inside the study "
                f"area (LockRaster IDs {selected_ids}). Use source='naip' (Earth Engine)."
            )
        if valid_fraction < LOW_COVERAGE_FRACTION:
            logger.warning(
                "USGS NAIP Plus imagery covers %.1f%% of the study area; the other pixels are 0 "
                "(nodata).",
                100.0 * valid_fraction,
            )

        self.last_metadata = {
            "AGRIBOUND_SOURCE": "usgs-naip-plus",
            "AGRIBOUND_VALUE_SCALE": SOURCE_REGISTRY["usgs-naip-plus"]["value_scale"],
            "AGRIBOUND_YEAR": config.year,
            "AGRIBOUND_EXPORT_CRS": grid.crs,
            "AGRIBOUND_RESOLUTION_M": resolution_m,
            "AGRIBOUND_IMAGE_YEARS": ",".join(
                sorted({str(c.year) for c in selected if c.year is not None})
            ),
            "AGRIBOUND_LOCK_RASTER_IDS": ",".join(str(i) for i in selected_ids),
            "AGRIBOUND_FOOTPRINT_COVERAGE": round(footprint_share, 6),
            "AGRIBOUND_VALID_FRACTION": round(valid_fraction, 6),
        }
        with rasterio.open(final_path, "r+") as dst:
            dst.update_tags(**{k: str(v) for k, v in self.last_metadata.items()})

        self._write_manifest(
            manifest_path=manifest_path,
            config=config,
            service_metadata=service_metadata,
            where=where,
            selected_ids=selected_ids,
            ranked_candidates=ranked_candidates,
            resolution_m=resolution_m,
            resolution_3857=resolution_3857,
            export_grid=grid,
            raw_tile_paths=raw_tile_paths,
            final_path=final_path,
            aoi_3857=aoi_3857,
            coverage={
                "footprint_coverage_of_extent": round(footprint_share, 6),
                "valid_fraction_in_study_area": round(valid_fraction, 6),
                "max_mosaic_image_count_reached": len(selected_ids) >= max_ids,
            },
        )
        import shutil

        shutil.rmtree(tile_dir, ignore_errors=True)
        logger.info("USGS NAIP Plus raster written to: %s", final_path)
        return str(final_path)

    def get_band_mapping(self, source: str) -> dict[str, str]:
        """Return the canonical band mapping (R, G, B, NIR -> R, G, B, N)."""
        info = SOURCE_REGISTRY.get(source, {})
        return dict(info.get("canonical_bands") or {})

    # ------------------------------------------------------------------
    # Query and selection
    # ------------------------------------------------------------------

    def _build_where_clause(
        self,
        config: AgriboundConfig,
        *,
        allow_year_fallback: bool = False,
        include_year: bool = True,
    ) -> str:
        parts = ["Category = 1"]

        state = getattr(config, "usgs_state", None)
        if state:
            state_clean = str(state).strip().upper().replace("'", "''")
            parts.append(f"State = '{state_clean}'")

        if not include_year:
            pass
        elif allow_year_fallback:
            parts.append(f"Year >= {int(config.year) - 1}")
            parts.append(f"Year <= {int(config.year) + 1}")
        else:
            parts.append(f"Year = {int(config.year)}")

        return " AND ".join(parts)

    def _no_imagery_error(
        self, client: USGSNAIPPlusClient, config: AgriboundConfig, aoi_3857: BaseGeometry
    ) -> NoDataError:
        """Build the 'no imagery' error, listing the years the service has for the area."""
        state = getattr(config, "usgs_state", None)
        where_any_year = self._build_where_clause(config, include_year=False)
        try:
            counts = client.query_year_counts(aoi_3857.bounds, where_any_year)
            years = ", ".join(f"{y} ({n} items)" for y, n in counts.items()) or "none"
        except Exception as exc:  # the listing is only diagnostic
            logger.warning("Could not list the available USGS NAIP Plus years: %s", exc)
            years = "unknown (year query failed)"
        requested = (
            f"{int(config.year) - 1}-{int(config.year) + 1}"
            if getattr(config, "usgs_allow_year_fallback", False)
            else str(config.year)
        )
        where = f" in state {state!r}" if state else ""
        return NoDataError(
            f"USGS NAIP Plus has no imagery for {requested}{where} over the study-area extent. "
            f"Years available for this area: {years}. The service holds only the latest "
            "vintage of each state; use source='naip' (Earth Engine, 2002-2023) for other years."
        )

    def _filter_candidates(
        self,
        candidates: list[USGSRasterCandidate],
        aoi_3857: BaseGeometry,
    ) -> list[USGSRasterCandidate]:
        filtered: list[USGSRasterCandidate] = []
        for candidate in candidates:
            if candidate.geometry is None or candidate.geometry.is_empty:
                continue
            if not candidate.geometry.intersects(aoi_3857):
                continue
            filtered.append(candidate)
        return filtered

    def _select_lock_raster_ids(
        self,
        candidates: list[USGSRasterCandidate],
        aoi_3857: BaseGeometry,
        *,
        target_year: int,
        max_ids: int,
    ) -> tuple[list[int], list[USGSRasterCandidate]]:
        ranked = sorted(
            candidates,
            key=lambda candidate: self._candidate_sort_key(candidate, aoi_3857, target_year),
        )

        selected: list[USGSRasterCandidate] = []
        covered = None

        for candidate in ranked:
            if len(selected) >= max_ids:
                break

            footprint = candidate.geometry.intersection(aoi_3857)
            if footprint.is_empty:
                continue

            if covered is None:
                selected.append(candidate)
                covered = footprint
                if covered.area >= aoi_3857.area * 0.999:
                    break
                continue

            additional = footprint.difference(covered)
            if additional.is_empty:
                continue

            selected.append(candidate)
            covered = covered.union(footprint)
            if covered.area >= aoi_3857.area * 0.999:
                break

        if not selected and ranked:
            selected = [ranked[0]]

        selected_ids = [candidate.object_id for candidate in selected]
        return selected_ids, ranked

    @staticmethod
    def _footprint_coverage(selected: list[USGSRasterCandidate], aoi_3857: BaseGeometry) -> float:
        """Share of *aoi_3857* (the export-grid outline) covered by the selected footprints."""
        from shapely.ops import unary_union

        if aoi_3857.is_empty or aoi_3857.area <= 0:
            return 0.0
        parts = [
            c.geometry.intersection(aoi_3857)
            for c in selected
            if c.geometry is not None and not c.geometry.is_empty
        ]
        if not parts:
            return 0.0
        return float(min(1.0, unary_union(parts).area / aoi_3857.area))

    def _candidate_sort_key(
        self,
        candidate: USGSRasterCandidate,
        aoi_3857: BaseGeometry,
        target_year: int,
    ) -> tuple[Any, ...]:
        intersection_area = 0.0
        if candidate.geometry is not None and not candidate.geometry.is_empty:
            intersection_area = candidate.geometry.intersection(aoi_3857).area

        year_distance = abs(candidate.year - target_year) if candidate.year is not None else 9999

        return (
            0 if candidate.category == 1 else 1,
            year_distance,
            0 if candidate.band_count and candidate.band_count >= 4 else 1,
            -intersection_area,
            candidate.resolution_value if candidate.resolution_value is not None else float("inf"),
            candidate.acquisition_date or "",
            candidate.object_id,
        )

    def _estimate_resolution_m(self, candidates: list[USGSRasterCandidate]) -> float:
        """Finest ground resolution (m) among *candidates* (the selected footprints).

        ``resolution_value`` is interpreted with ``resolution_units`` (metres
        or feet); candidates with other or missing units are ignored. Returns
        1.0 m, with a WARNING, when no candidate has a usable resolution.
        """
        values = []
        for candidate in candidates:
            value = candidate.resolution_value
            if value is None or value <= 0:
                continue
            units = (candidate.resolution_units or "").strip().lower().replace("-", "_")
            if units in _METRE_UNITS:
                values.append(float(value))
            elif units in _FOOT_UNITS:
                values.append(float(value) * 0.3048)
            else:
                logger.warning(
                    "Ignoring resolution %s %r of LockRaster %d (unknown units)",
                    value,
                    candidate.resolution_units,
                    candidate.object_id,
                )
        if values:
            return float(min(values))
        logger.warning(
            "No usable resolution in the selected USGS footprints; exporting at %.1f m",
            _DEFAULT_RESOLUTION_M,
        )
        return _DEFAULT_RESOLUTION_M

    # ------------------------------------------------------------------
    # Raster handling
    # ------------------------------------------------------------------

    def _validate_export(self, path: str | Path) -> None:
        path = Path(path)
        if not path.exists():
            raise FileNotFoundError(f"Expected export not found: {path}")

        with rasterio.open(path) as src:
            if src.width <= 0 or src.height <= 0:
                raise RuntimeError(f"Invalid raster dimensions in {path}")
            if src.count < 3:
                raise RuntimeError(f"Expected at least 3 bands in {path}, found {src.count}")
            if src.crs is None:
                raise RuntimeError(f"Raster CRS is missing in {path}")

    @staticmethod
    def _write_final(tile_paths: list[str], out_path: Path, grid: Any) -> int:
        """Merge the exported tiles and resample them (bilinear) onto *grid*.

        Tiles are merged in their own CRS (``rasterio.merge.merge``), then read
        through a ``WarpedVRT`` on *grid* block by block; pixels that no tile
        covers are set to 0 (nodata). Pixels are not masked to the study-area
        polygons.

        Returns
        -------
        int
            Number of bands written.
        """
        from rasterio.enums import Resampling
        from rasterio.merge import merge
        from rasterio.vrt import WarpedVRT
        from rasterio.windows import Window

        mosaic_path = out_path.with_name(out_path.stem + "_mosaic.tif")
        part = out_path.with_name(out_path.name + ".part")
        try:
            if len(tile_paths) == 1:
                src_path = tile_paths[0]
            else:
                merge(
                    tile_paths,
                    dst_path=str(mosaic_path),
                    dst_kwds={"compress": "deflate", "tiled": True, "BIGTIFF": "IF_SAFER"},
                )
                src_path = str(mosaic_path)

            with rasterio.open(src_path) as src:
                count = src.count
                profile = {
                    "driver": "GTiff",
                    "dtype": "uint8",
                    "count": count,
                    "width": grid.width,
                    "height": grid.height,
                    "crs": grid.crs,
                    "transform": grid.transform,
                    "nodata": 0,
                    "tiled": True,
                    "blockxsize": 512,
                    "blockysize": 512,
                    "compress": "deflate",
                    "predictor": 2,
                    "BIGTIFF": "IF_SAFER",
                }
                with (
                    WarpedVRT(
                        src,
                        crs=grid.crs,
                        transform=grid.transform,
                        width=grid.width,
                        height=grid.height,
                        resampling=Resampling.bilinear,
                    ) as vrt,
                    rasterio.open(part, "w", **profile) as dst,
                ):
                    names = ["R", "G", "B", "N"][:count]
                    for i, name in enumerate(names, start=1):
                        dst.set_band_description(i, name)
                    rows = 1024
                    for r0 in range(0, grid.height, rows):
                        window = Window(0, r0, grid.width, min(rows, grid.height - r0))
                        data = vrt.read(window=window)
                        covered = vrt.read_masks(1, window=window) > 0
                        data[:, ~covered] = 0
                        dst.write(data.astype(np.uint8, copy=False), window=window)
            os.replace(part, out_path)
        finally:
            if part.exists():
                part.unlink()
            if mosaic_path.exists():
                mosaic_path.unlink()
        return count

    def _write_manifest(
        self,
        *,
        manifest_path: Path,
        config: AgriboundConfig,
        service_metadata: dict[str, Any],
        where: str,
        selected_ids: list[int],
        ranked_candidates: list[USGSRasterCandidate],
        resolution_m: float,
        resolution_3857: float,
        export_grid: Any,
        raw_tile_paths: list[str],
        final_path: Path,
        aoi_3857: BaseGeometry,
        coverage: dict[str, Any] | None = None,
    ) -> None:
        manifest = {
            "source": "usgs-naip-plus",
            "selection_version": self.SELECTION_VERSION,
            "service_url": getattr(config, "usgs_service_url", DEFAULT_USGS_NAIPPLUS_URL),
            "year": config.year,
            "usgs_state": getattr(config, "usgs_state", None),
            "allow_year_fallback": getattr(config, "usgs_allow_year_fallback", False),
            "where": where,
            "selected_lock_raster_ids": selected_ids,
            "resolution_m": resolution_m,
            "export_resolution_3857_units": resolution_3857,
            "export_crs": export_grid.crs,
            "export_shape": [export_grid.height, export_grid.width],
            "export_transform": list(export_grid.crs_transform),
            "query_bounds_3857": [float(v) for v in aoi_3857.bounds],
            # The tiles are deleted once the final raster is written.
            "n_export_tiles": len(raw_tile_paths),
            "coverage": dict(coverage or {}),
            "final_raster_path": str(final_path),
            "service_limits": {
                "maxImageWidth": service_metadata.get("maxImageWidth"),
                "maxImageHeight": service_metadata.get("maxImageHeight"),
                "maxMosaicImageCount": service_metadata.get("maxMosaicImageCount"),
            },
            "ranked_candidates": [
                {
                    "object_id": candidate.object_id,
                    "year": candidate.year,
                    "state": candidate.state,
                    "acquisition_date": candidate.acquisition_date,
                    "resolution_value": candidate.resolution_value,
                    "resolution_units": candidate.resolution_units,
                    "band_count": candidate.band_count,
                    "category": candidate.category,
                    "name": candidate.name,
                    "download_url": candidate.download_url,
                    "bounds_3857": (
                        list(candidate.geometry.bounds)
                        if candidate.geometry is not None and not candidate.geometry.is_empty
                        else None
                    ),
                }
                for candidate in ranked_candidates
            ],
        }
        manifest_path.write_text(json.dumps(manifest, indent=2), encoding="utf-8")

build

build(config: AgriboundConfig) -> str

Download and mosaic a USGS NAIP Plus raster covering the study-area extent.

Returns:

Type	Description
str	Path to the uint8 GeoTIFF (cached with :func:agribound._cache.cache_path).

Raises:

Type	Description
NoDataError	(A :class:ValueError.) If the service has no imagery for the requested year (the message lists the years it has for the area), or no pixel inside the study area has imagery.
ValueError	If no study area is set or the study area is empty.

Source code in agribound/composites/usgs.py

def build(self, config: AgriboundConfig) -> str:
    """Download and mosaic a USGS NAIP Plus raster covering the study-area extent.

    Returns
    -------
    str
        Path to the uint8 GeoTIFF (cached with
        :func:`agribound._cache.cache_path`).

    Raises
    ------
    NoDataError
        (A :class:`ValueError`.) If the service has no imagery for the
        requested year (the message lists the years it has for the area),
        or no pixel inside the study area has imagery.
    ValueError
        If no study area is set or the study area is empty.
    """
    from agribound._cache import cache_path

    if not config.study_area:
        raise ValueError("study_area is required when source='usgs-naip-plus'")

    from agribound.composites.gee import (
        compute_export_grid,
        grid_footprint_4326,
        raster_valid_fraction,
        resolve_export_crs,
        split_grid,
        study_area_geometry_4326,
    )

    geom_4326 = study_area_geometry_4326(config)

    final_path = cache_path(config, "usgs_naip_plus", ".tif", self.SELECTION_VERSION)
    manifest_path = final_path.with_suffix(".json")
    if final_path.exists() and manifest_path.exists():
        logger.info("Using cached USGS NAIP Plus raster: %s", final_path)
        return str(final_path)

    client = USGSNAIPPlusClient(
        service_url=getattr(config, "usgs_service_url", DEFAULT_USGS_NAIPPLUS_URL),
        timeout_s=getattr(config, "usgs_timeout_s", 120),
        retries=getattr(config, "usgs_retries", 3),
    )

    service_metadata = client.get_service_metadata()
    # The query/selection target is the outline of a provisional export grid
    # (the study-area bounding box in the export CRS), in EPSG:3857.
    export_crs = resolve_export_crs(config.export_crs, geom_4326)
    provisional = compute_export_grid(geom_4326, export_crs, _DEFAULT_RESOLUTION_M)
    aoi_3857 = _to_3857(grid_footprint_4326(provisional))
    where = self._build_where_clause(config)
    candidates = client.query_candidates(aoi_3857.bounds, where, out_sr=3857)
    candidates = self._filter_candidates(candidates, aoi_3857)

    if not candidates and getattr(config, "usgs_allow_year_fallback", False):
        where = self._build_where_clause(config, allow_year_fallback=True)
        candidates = client.query_candidates(aoi_3857.bounds, where, out_sr=3857)
        candidates = self._filter_candidates(candidates, aoi_3857)

    if not candidates:
        raise self._no_imagery_error(client, config, aoi_3857)

    max_ids = int(service_metadata.get("maxMosaicImageCount", 50))
    selected_ids, ranked_candidates = self._select_lock_raster_ids(
        candidates,
        aoi_3857,
        target_year=int(config.year),
        max_ids=max_ids,
    )
    if not selected_ids:
        raise RuntimeError("Candidate selection returned no LockRaster IDs")
    selected = [c for c in ranked_candidates if c.object_id in set(selected_ids)]
    footprint_share = self._footprint_coverage(selected, aoi_3857)
    if footprint_share < LOW_COVERAGE_FRACTION:
        reason = (
            f"the selection reached the service's maxMosaicImageCount ({max_ids})"
            if len(selected_ids) >= max_ids
            else "no other candidate footprint adds coverage"
        )
        logger.warning(
            "The %d selected USGS NAIP Plus footprint(s) cover %.1f%% of the study-area "
            "extent (%s); uncovered pixels are 0 (nodata).",
            len(selected_ids),
            100.0 * footprint_share,
            reason,
        )

    export_max_px = min(
        int(getattr(config, "tile_size", 10_000)),
        int(service_metadata.get("maxImageWidth", 4000)),
        int(service_metadata.get("maxImageHeight", 4000)),
    )
    resolution_m = self._estimate_resolution_m(selected)
    centroid_lat = float(geom_4326.centroid.y)
    resolution_3857 = resolution_m / math.cos(math.radians(centroid_lat))

    grid = compute_export_grid(geom_4326, export_crs, resolution_m)
    grid_3857 = compute_export_grid(grid_footprint_4326(grid), "EPSG:3857", resolution_3857)
    windows = split_grid(grid_3857, export_max_px)
    tile_dir = final_path.with_name(final_path.stem + "_tiles")
    tile_dir.mkdir(parents=True, exist_ok=True)

    raw_tile_paths: list[str] = []
    for idx, (row_off, col_off, height, width) in enumerate(windows):
        tile_grid = grid_3857.window(row_off, col_off, height, width)
        t = tile_grid.transform
        bounds = (t.c, t.f + height * t.e, t.c + width * t.a, t.f)
        tile_path = tile_dir / f"tile_{idx:03d}.tif"
        if not tile_path.exists():
            logger.info(
                "Exporting USGS NAIP Plus tile %d/%d (%d x %d px) with %d LockRaster IDs",
                idx + 1,
                len(windows),
                width,
                height,
                len(selected_ids),
            )
            client.export_image(
                bbox_3857=bounds,
                width=width,
                height=height,
                lock_raster_ids=selected_ids,
                output_path=tile_path,
            )
        self._validate_export(tile_path)
        raw_tile_paths.append(str(tile_path))

    band_count = self._write_final(raw_tile_paths, final_path, grid)
    if band_count < 4:
        logger.warning(
            "USGS NAIP Plus export has %d bands; the NIR band (N) is missing", band_count
        )
    # Uncovered pixels are 0 in every band; a dark pixel can be 0 in one band.
    valid_fraction = raster_valid_fraction(final_path, geom_4326, require="any")
    if valid_fraction <= 0.0:
        final_path.unlink(missing_ok=True)
        raise NoDataError(
            f"The USGS NAIP Plus export for {config.year} has no imagery inside the study "
            f"area (LockRaster IDs {selected_ids}). Use source='naip' (Earth Engine)."
        )
    if valid_fraction < LOW_COVERAGE_FRACTION:
        logger.warning(
            "USGS NAIP Plus imagery covers %.1f%% of the study area; the other pixels are 0 "
            "(nodata).",
            100.0 * valid_fraction,
        )

    self.last_metadata = {
        "AGRIBOUND_SOURCE": "usgs-naip-plus",
        "AGRIBOUND_VALUE_SCALE": SOURCE_REGISTRY["usgs-naip-plus"]["value_scale"],
        "AGRIBOUND_YEAR": config.year,
        "AGRIBOUND_EXPORT_CRS": grid.crs,
        "AGRIBOUND_RESOLUTION_M": resolution_m,
        "AGRIBOUND_IMAGE_YEARS": ",".join(
            sorted({str(c.year) for c in selected if c.year is not None})
        ),
        "AGRIBOUND_LOCK_RASTER_IDS": ",".join(str(i) for i in selected_ids),
        "AGRIBOUND_FOOTPRINT_COVERAGE": round(footprint_share, 6),
        "AGRIBOUND_VALID_FRACTION": round(valid_fraction, 6),
    }
    with rasterio.open(final_path, "r+") as dst:
        dst.update_tags(**{k: str(v) for k, v in self.last_metadata.items()})

    self._write_manifest(
        manifest_path=manifest_path,
        config=config,
        service_metadata=service_metadata,
        where=where,
        selected_ids=selected_ids,
        ranked_candidates=ranked_candidates,
        resolution_m=resolution_m,
        resolution_3857=resolution_3857,
        export_grid=grid,
        raw_tile_paths=raw_tile_paths,
        final_path=final_path,
        aoi_3857=aoi_3857,
        coverage={
            "footprint_coverage_of_extent": round(footprint_share, 6),
            "valid_fraction_in_study_area": round(valid_fraction, 6),
            "max_mosaic_image_count_reached": len(selected_ids) >= max_ids,
        },
    )
    import shutil

    shutil.rmtree(tile_dir, ignore_errors=True)
    logger.info("USGS NAIP Plus raster written to: %s", final_path)
    return str(final_path)

get_band_mapping

get_band_mapping(source: str) -> dict[str, str]

Return the canonical band mapping (R, G, B, NIR -> R, G, B, N).

Source code in agribound/composites/usgs.py

def get_band_mapping(self, source: str) -> dict[str, str]:
    """Return the canonical band mapping (R, G, B, NIR -> R, G, B, N)."""
    info = SOURCE_REGISTRY.get(source, {})
    return dict(info.get("canonical_bands") or {})

Local GeoTIFFs and embeddings

local

Local GeoTIFF and embedding builders.

• :class:LocalCompositeBuilder (source="local"): validates a user GeoTIFF and, when a study area is configured, crops it to the study area's bounding box (pixel values unchanged). Its radiometry is unknown to Agribound (value scale "unknown").

• :class:EmbeddingCompositeBuilder: pre-computed embeddings.

• tessera-embedding: TESSERA embeddings streamed from the public Zarr stores with :class:geotessera.GeoTesseraZarr (geotessera >= 0.10), for config.tessera_version / config.tessera_variant. The study area is split by UTM zone, because GeoTesseraZarr.read_region only reads the zone that holds the centre of the requested box; each zone is read on its native 10 m grid and resampled (nearest neighbour) onto one grid in config.export_crs (for "utm", the zone of the study-area centroid, whose pixels are copied without resampling).

• google-embedding: Google Satellite Embedding V1 (AlphaEarth Foundations, 64 bands A00-A63, 2017-2025) from Earth Engine (google_embedding_backend="gee", default) or from the Source Cooperative COG mirror (google_embedding_backend="source_coop", no Earth Engine compute), whose int8 tiles are read window by window with rasterio using the mirror's tile index and de-quantised as geoai-py does.

Embedding rasters are float32 with NaN nodata (NaN marks pixels without an embedding: water, no data or missing coverage). As for the imagery builders, they cover the bounding box of the study area in the export CRS and are not masked to the study-area polygons. The TESSERA and Source Cooperative paths hold the whole area in memory while it is assembled (512 bytes per pixel for TESSERA, 256 for Google embeddings, about twice that at peak while the zone reads and the output coexist), so split very large regions into tiles (see agribound.hpc).

:func:tessera_coverage reports TESSERA tile coverage of a bounding box from the dataset manifest without downloading embeddings.

LocalCompositeBuilder

Bases: CompositeBuilder

Builder for user-provided GeoTIFFs (source="local").

Validates the file (:func:validate_local_raster). When a study area is configured, the raster is cropped to the study area's bounding box (:func:crop_to_extent; pixel values are not changed or masked) into the cache (cache_path(config, "local_crop_<name>", ".tif", ...), keyed by the study area and the file's path, size and modification time). Without a study area, or when the bounding box contains the whole raster, the file itself is returned.

Attributes:

Name	Type	Description
last_metadata	dict	Facts about the most recent :meth:build (value scale "unknown").

Source code in agribound/composites/local.py

class LocalCompositeBuilder(CompositeBuilder):
    """Builder for user-provided GeoTIFFs (``source="local"``).

    Validates the file (:func:`validate_local_raster`). When a study area is
    configured, the raster is cropped to the study area's bounding box
    (:func:`crop_to_extent`; pixel values are not changed or masked) into the
    cache (``cache_path(config, "local_crop_<name>", ".tif", ...)``, keyed by
    the study area and the file's path, size and modification time). Without
    a study area, or when the bounding box contains the whole raster, the file
    itself is returned.

    Attributes
    ----------
    last_metadata : dict
        Facts about the most recent :meth:`build` (value scale ``"unknown"``).
    """

    #: Version of the crop recipe (part of the cache key).
    CROP_VERSION = 2

    def __init__(self) -> None:
        self.last_metadata: dict[str, Any] = {}

    def build(self, config: AgriboundConfig) -> str:
        """Validate the local GeoTIFF and crop it to the study-area extent.

        Raises
        ------
        NoDataError
            (A :class:`ValueError`.) If the study-area extent does not overlap
            the raster.
        ValueError
            If ``local_tif_path`` is unset or the raster is unusable.
        FileNotFoundError
            If ``local_tif_path`` does not exist.
        """
        from agribound._cache import cache_path

        if config.local_tif_path is None:
            raise ValueError("local_tif_path must be set when source='local'")
        src_path = Path(config.local_tif_path).expanduser()
        info = validate_local_raster(src_path, config.bands)
        self.last_metadata = {
            "AGRIBOUND_SOURCE": "local",
            "AGRIBOUND_VALUE_SCALE": SOURCE_REGISTRY["local"]["value_scale"],
            "AGRIBOUND_LOCAL_SOURCE_PATH": str(src_path.resolve()),
            "AGRIBOUND_LOCAL_CRS": info["crs"],
            "AGRIBOUND_LOCAL_BAND_COUNT": info["count"],
            "AGRIBOUND_LOCAL_DTYPE": info["dtype"],
        }
        logger.info(
            "Local GeoTIFF %s: %d band(s), %s, %s, nodata=%s (value scale unknown)",
            src_path,
            info["count"],
            info["dtype"],
            info["crs"],
            info["nodata"],
        )
        if not config.study_area:
            return str(src_path)

        from agribound.composites.gee import study_area_geometry_4326

        cropped = cache_path(config, f"local_crop_{src_path.stem}", ".tif", self.CROP_VERSION)
        if cropped.exists():
            logger.info("Using cached cropped raster: %s", cropped)
            return str(cropped)

        geometry = study_area_geometry_4326(config)
        logger.info("Cropping %s to the study-area extent", src_path)
        try:
            result = crop_to_extent(src_path, cropped, geometry, tags=self.last_metadata)
        except ValueError as exc:
            # A study area outside the raster is a no-data condition (NoDataError).
            error_type = NoDataError if isinstance(exc, NoDataError) else ValueError
            raise error_type(
                f"Could not crop {src_path} to the study area {config.study_area!r}: {exc}. "
                "Check that they overlap, or leave study_area empty to use the whole raster."
            ) from exc
        if result is None:
            logger.info("The study-area extent contains the whole raster; using %s", src_path)
            return str(src_path)
        return result

    def get_band_mapping(self, source: str) -> dict[str, str]:
        """Return the positional band mapping assumed for local files.

        Engines read ``R``, ``G``, ``B`` (and ``NIR``) from bands 1, 2, 3 (4)
        unless ``AgriboundConfig.bands`` maps canonical names to other band
        indices (see :func:`agribound.engines.base.get_canonical_band_indices`).
        """
        return {"R": "1", "G": "2", "B": "3", "NIR": "4"}

build

build(config: AgriboundConfig) -> str

Validate the local GeoTIFF and crop it to the study-area extent.

Raises:

Type	Description
NoDataError	(A :class:ValueError.) If the study-area extent does not overlap the raster.
ValueError	If local_tif_path is unset or the raster is unusable.
FileNotFoundError	If local_tif_path does not exist.

Source code in agribound/composites/local.py

def build(self, config: AgriboundConfig) -> str:
    """Validate the local GeoTIFF and crop it to the study-area extent.

    Raises
    ------
    NoDataError
        (A :class:`ValueError`.) If the study-area extent does not overlap
        the raster.
    ValueError
        If ``local_tif_path`` is unset or the raster is unusable.
    FileNotFoundError
        If ``local_tif_path`` does not exist.
    """
    from agribound._cache import cache_path

    if config.local_tif_path is None:
        raise ValueError("local_tif_path must be set when source='local'")
    src_path = Path(config.local_tif_path).expanduser()
    info = validate_local_raster(src_path, config.bands)
    self.last_metadata = {
        "AGRIBOUND_SOURCE": "local",
        "AGRIBOUND_VALUE_SCALE": SOURCE_REGISTRY["local"]["value_scale"],
        "AGRIBOUND_LOCAL_SOURCE_PATH": str(src_path.resolve()),
        "AGRIBOUND_LOCAL_CRS": info["crs"],
        "AGRIBOUND_LOCAL_BAND_COUNT": info["count"],
        "AGRIBOUND_LOCAL_DTYPE": info["dtype"],
    }
    logger.info(
        "Local GeoTIFF %s: %d band(s), %s, %s, nodata=%s (value scale unknown)",
        src_path,
        info["count"],
        info["dtype"],
        info["crs"],
        info["nodata"],
    )
    if not config.study_area:
        return str(src_path)

    from agribound.composites.gee import study_area_geometry_4326

    cropped = cache_path(config, f"local_crop_{src_path.stem}", ".tif", self.CROP_VERSION)
    if cropped.exists():
        logger.info("Using cached cropped raster: %s", cropped)
        return str(cropped)

    geometry = study_area_geometry_4326(config)
    logger.info("Cropping %s to the study-area extent", src_path)
    try:
        result = crop_to_extent(src_path, cropped, geometry, tags=self.last_metadata)
    except ValueError as exc:
        # A study area outside the raster is a no-data condition (NoDataError).
        error_type = NoDataError if isinstance(exc, NoDataError) else ValueError
        raise error_type(
            f"Could not crop {src_path} to the study area {config.study_area!r}: {exc}. "
            "Check that they overlap, or leave study_area empty to use the whole raster."
        ) from exc
    if result is None:
        logger.info("The study-area extent contains the whole raster; using %s", src_path)
        return str(src_path)
    return result

get_band_mapping

get_band_mapping(source: str) -> dict[str, str]

Return the positional band mapping assumed for local files.

Engines read R, G, B (and NIR) from bands 1, 2, 3 (4) unless AgriboundConfig.bands maps canonical names to other band indices (see :func:agribound.engines.base.get_canonical_band_indices).

Source code in agribound/composites/local.py

def get_band_mapping(self, source: str) -> dict[str, str]:
    """Return the positional band mapping assumed for local files.

    Engines read ``R``, ``G``, ``B`` (and ``NIR``) from bands 1, 2, 3 (4)
    unless ``AgriboundConfig.bands`` maps canonical names to other band
    indices (see :func:`agribound.engines.base.get_canonical_band_indices`).
    """
    return {"R": "1", "G": "2", "B": "3", "NIR": "4"}

EmbeddingCompositeBuilder

Bases: CompositeBuilder

Builder for the pre-computed embedding sources.

tessera-embedding uses :func:build_tessera_embedding; google-embedding uses :func:build_google_embedding_gee or, with google_embedding_backend="source_coop", :func:build_google_embedding_source_coop. Outputs are cached with :func:agribound._cache.cache_path (keyed by study area, year, export CRS, dataset version/variant and backend).

Attributes:

Name	Type	Description
last_metadata	dict	Tags of the most recent output.

Source code in agribound/composites/local.py

class EmbeddingCompositeBuilder(CompositeBuilder):
    """Builder for the pre-computed embedding sources.

    ``tessera-embedding`` uses :func:`build_tessera_embedding`;
    ``google-embedding`` uses :func:`build_google_embedding_gee` or, with
    ``google_embedding_backend="source_coop"``,
    :func:`build_google_embedding_source_coop`. Outputs are cached with
    :func:`agribound._cache.cache_path` (keyed by study area, year, export
    CRS, dataset version/variant and backend).

    Attributes
    ----------
    last_metadata : dict
        Tags of the most recent output.
    """

    def __init__(self) -> None:
        self.last_metadata: dict[str, Any] = {}

    def build(self, config: AgriboundConfig) -> str:
        """Fetch the embeddings for the study area and return the GeoTIFF path."""
        if not config.study_area:
            raise ValueError(f"study_area is required for source={config.source!r}")
        if config.source == "tessera-embedding":
            path, meta = build_tessera_embedding(config)
        elif config.source == "google-embedding":
            backend = getattr(config, "google_embedding_backend", "gee") or "gee"
            if backend == "gee":
                path, meta = build_google_embedding_gee(config)
            elif backend == "source_coop":
                path, meta = build_google_embedding_source_coop(config)
            else:
                raise ValueError(
                    f"Unknown google_embedding_backend {backend!r}; use 'gee' or 'source_coop'"
                )
        else:
            raise ValueError(f"Unknown embedding source: {config.source}")
        self.last_metadata = dict(meta)
        return path

    def get_band_mapping(self, source: str) -> dict[str, str]:
        """Embedding sources have no canonical optical bands; returns an empty mapping."""
        return dict(SOURCE_REGISTRY.get(source, {}).get("canonical_bands") or {})

build

build(config: AgriboundConfig) -> str

Fetch the embeddings for the study area and return the GeoTIFF path.

Source code in agribound/composites/local.py

def build(self, config: AgriboundConfig) -> str:
    """Fetch the embeddings for the study area and return the GeoTIFF path."""
    if not config.study_area:
        raise ValueError(f"study_area is required for source={config.source!r}")
    if config.source == "tessera-embedding":
        path, meta = build_tessera_embedding(config)
    elif config.source == "google-embedding":
        backend = getattr(config, "google_embedding_backend", "gee") or "gee"
        if backend == "gee":
            path, meta = build_google_embedding_gee(config)
        elif backend == "source_coop":
            path, meta = build_google_embedding_source_coop(config)
        else:
            raise ValueError(
                f"Unknown google_embedding_backend {backend!r}; use 'gee' or 'source_coop'"
            )
    else:
        raise ValueError(f"Unknown embedding source: {config.source}")
    self.last_metadata = dict(meta)
    return path

get_band_mapping

get_band_mapping(source: str) -> dict[str, str]

Embedding sources have no canonical optical bands; returns an empty mapping.

Source code in agribound/composites/local.py

def get_band_mapping(self, source: str) -> dict[str, str]:
    """Embedding sources have no canonical optical bands; returns an empty mapping."""
    return dict(SOURCE_REGISTRY.get(source, {}).get("canonical_bands") or {})

validate_local_raster

validate_local_raster(path: str | Path, bands: dict[str, int] | None = None) -> dict[str, Any]

Check that a local GeoTIFF can be used as pipeline input.

Parameters:

Name	Type	Description	Default
path	str or Path	Raster file.	required
bands	dict or None	AgriboundConfig.bands (canonical name -> 1-based band index).	None

Returns:

Type	Description
dict	crs, count, dtype, nodata, width, height, res.

Raises:

Type	Description
FileNotFoundError	If the file does not exist.
ValueError	If the raster has no CRS, no bands, or bands holds an index outside 1 .. band count (indices are 1-based).

Notes

A raster without a nodata value is accepted with a WARNING (zeros are then treated as valid data by the engines); fewer than 3 bands without a bands mapping is accepted with a WARNING (RGB engines need 3 bands).

Source code in agribound/composites/local.py

def validate_local_raster(path: str | Path, bands: dict[str, int] | None = None) -> dict[str, Any]:
    """Check that a local GeoTIFF can be used as pipeline input.

    Parameters
    ----------
    path : str or Path
        Raster file.
    bands : dict or None
        ``AgriboundConfig.bands`` (canonical name -> 1-based band index).

    Returns
    -------
    dict
        ``crs``, ``count``, ``dtype``, ``nodata``, ``width``, ``height``,
        ``res``.

    Raises
    ------
    FileNotFoundError
        If the file does not exist.
    ValueError
        If the raster has no CRS, no bands, or *bands* holds an index outside
        ``1 .. band count`` (indices are 1-based).

    Notes
    -----
    A raster without a nodata value is accepted with a WARNING (zeros are
    then treated as valid data by the engines); fewer than 3 bands without a
    *bands* mapping is accepted with a WARNING (RGB engines need 3 bands).
    """
    import rasterio

    path = Path(path).expanduser()
    if not path.exists():
        raise FileNotFoundError(f"Local TIF not found: {path}")
    with rasterio.open(path) as src:
        info = {
            "crs": src.crs.to_string() if src.crs else None,
            "count": src.count,
            "dtype": str(src.dtypes[0]),
            "nodata": src.nodata,
            "width": src.width,
            "height": src.height,
            "res": tuple(float(v) for v in src.res),
        }
    if info["crs"] is None:
        raise ValueError(
            f"{path} has no CRS. Assign one (e.g. gdal_edit.py -a_srs EPSG:XXXX) before using it."
        )
    if info["count"] < 1 or info["width"] < 1 or info["height"] < 1:
        raise ValueError(f"{path} has no pixels or no bands")
    if bands:
        invalid = {k: v for k, v in bands.items() if not 1 <= int(v) <= info["count"]}
        if invalid:
            raise ValueError(
                f"bands={dict(bands)} has invalid band index(es) {invalid}: indices are 1-based "
                f"and {path} has {info['count']} band(s), so they must be in 1..{info['count']}"
            )
    elif info["count"] < 3:
        logger.warning(
            "%s has %d band(s); engines that need R, G, B read bands 1-3. Set 'bands' to map "
            "canonical names to band indices.",
            path,
            info["count"],
        )
    if info["nodata"] is None:
        logger.warning(
            "%s has no nodata value: every pixel (including zeros) is treated as valid data", path
        )
    return info

tessera_coverage

tessera_coverage(bbox: tuple[float, float, float, float], year: int, version: str = 'v1', variant: str | None = None, cache_dir: str | None = None) -> dict[str, Any]

Report TESSERA tile coverage of a bounding box from the dataset manifest.

Tiles are 0.1 x 0.1 degree cells. The manifest (GeoTessera(...) .registry) is downloaded once into cache_dir (about 212 MB for v1, 59 MB for v1.1); no embeddings are downloaded.

Parameters:

Name	Type	Description	Default
bbox	tuple	(min_lon, min_lat, max_lon, max_lat) in EPSG:4326.	required
year	int	Year to check.	required
version	str	"v1", "v1.1" or "v2".	'v1'
variant	str or None	Dataset variant (None = geotessera default for the version).	None
cache_dir	str or None	geotessera registry cache directory (e.g. shared HPC scratch).	None

Returns:

Type	Description
dict	version, variant, year, bbox, tiles_expected (all 0.1-degree cells intersecting the box, including water, which has no tiles), tiles_available, fraction (available / expected), tiles_by_year (available tiles in the box for every year of the dataset) and dataset_years.

Source code in agribound/composites/local.py

def tessera_coverage(
    bbox: tuple[float, float, float, float],
    year: int,
    version: str = "v1",
    variant: str | None = None,
    cache_dir: str | None = None,
) -> dict[str, Any]:
    """Report TESSERA tile coverage of a bounding box from the dataset manifest.

    Tiles are 0.1 x 0.1 degree cells. The manifest (``GeoTessera(...)
    .registry``) is downloaded once into *cache_dir* (about 212 MB for v1,
    59 MB for v1.1); no embeddings are downloaded.

    Parameters
    ----------
    bbox : tuple
        ``(min_lon, min_lat, max_lon, max_lat)`` in EPSG:4326.
    year : int
        Year to check.
    version : str
        ``"v1"``, ``"v1.1"`` or ``"v2"``.
    variant : str or None
        Dataset variant (*None* = geotessera default for the version).
    cache_dir : str or None
        geotessera registry cache directory (e.g. shared HPC scratch).

    Returns
    -------
    dict
        ``version``, ``variant``, ``year``, ``bbox``, ``tiles_expected``
        (all 0.1-degree cells intersecting the box, including water, which
        has no tiles), ``tiles_available``, ``fraction``
        (available / expected), ``tiles_by_year`` (available tiles in the box
        for every year of the dataset) and ``dataset_years``.
    """
    from geotessera import GeoTessera

    if version not in _TESSERA_NORMALISED:
        raise ValueError(
            f"Unknown tessera version {version!r}. Choose from {tuple(_TESSERA_NORMALISED)}"
        )
    minx, miny, maxx, maxy = (float(v) for v in bbox)
    gt = GeoTessera(dataset_version=version, dataset_variant=variant, cache_dir=cache_dir)

    def _in_box(lon: float, lat: float) -> bool:
        tol = 1e-9
        return (
            lon - 0.05 < maxx - tol
            and lon + 0.05 > minx + tol
            and lat - 0.05 < maxy - tol
            and lat + 0.05 > miny + tol
        )

    def _count(y: int) -> int:
        tiles = gt.registry.load_blocks_for_region(bounds=(minx, miny, maxx, maxy), year=int(y))
        return len({(round(lon, 2), round(lat, 2)) for _y, lon, lat in tiles if _in_box(lon, lat)})

    _, n_lon = _tile_count_1d(minx, maxx)
    _, n_lat = _tile_count_1d(miny, maxy)
    expected = n_lon * n_lat
    dataset_years = [int(y) for y in gt.registry.get_available_years()]
    by_year = {y: _count(y) for y in dataset_years}
    available = by_year.get(int(year), 0) if int(year) in by_year else _count(int(year))
    return {
        "version": version,
        "variant": gt.dataset_variant,
        "year": int(year),
        "bbox": (minx, miny, maxx, maxy),
        "tiles_expected": expected,
        "tiles_available": available,
        "fraction": (available / expected) if expected else 0.0,
        "tiles_by_year": by_year,
        "dataset_years": dataset_years,
    }