HPC Tiling

agribound.hpc backs the agribound tiles commands. See HPC and large areas.

hpc

Batch execution of large study areas (HPC job arrays).

• :mod:agribound.hpc.tiles: cut a study area into tiles with halos, write a manifest with one configuration per tile, run tiles idempotently, report their status and merge the results.
• :mod:agribound.hpc.regions: read the region definitions in examples/regions/*.yaml.
• :mod:agribound.hpc.cli: the agribound tiles command group.

Submodules are imported on first attribute access, so importing this package (which agribound does to register the CLI group) stays cheap.

Tiles

tiles

Tiling of large study areas for batch (HPC) runs.

A region that is too large for one composite download or one GPU job is cut into square tiles. Each tile is an ordinary Agribound run with its own YAML configuration, output file and provenance record, so tiles can be executed as independent Slurm array tasks and restarted at any time.

Workflow

1. :func:make_tiles cuts the study area into a grid of core cells and adds a halo around each core.
2. :func:write_tile_manifest writes tiles.gpkg, tiles.txt, one configuration per tile (tiles/<tile_id>/config.yaml) and manifest.json.
3. :func:run_tile runs one tile: stage="composite" (downloads only: composite or embeddings, the LULC raster for lulc_mode="raster" and the FTW window composites of two-window FTW models, also as ensemble members), stage="delineate" (delineation of an already staged tile) or stage="all" (the composite stage, then the delineation), using :func:agribound.pipeline.build_composite and :func:agribound.pipeline.delineate.
4. :func:merge_tiles combines the tile outputs into one file and writes a merged provenance summary; :func:tile_status reports progress.

No-data tiles

A rectangular grid over a real region contains tiles without input data: open water, or areas outside a source's coverage (NAIP outside the US, missing TESSERA tiles). The composite builders raise :class:agribound.composites.base.NoDataError (a :class:ValueError) for them: no image intersects the study-area extent, no valid pixel inside the study area, no TESSERA / Google embedding data, no USGS NAIP Plus imagery, a local raster that does not overlap. :func:run_tile recognises these errors by their type (:func:no_data_reason; a :class:ValueError whose message matches :data:NO_DATA_PATTERNS is accepted as a fallback) and records the tile as "no-data" instead of failing: it writes a content-addressed nodata_<key>.json marker in the tile's cache directory (plus no_data.json in the tile directory) with the error, and returns status="no-data". The marker is keyed like the stage markers, so it is reused by later runs with the same inputs and ignored after a configuration change; overwrite=True (--overwrite) retries the download. Any other error still fails the tile. :func:merge_tiles merges the other tiles and lists the no-data tiles, their reasons and their core area in the summary; it raises if no tile produced output, which usually means that the source has no data for the year.

Grids

crs="utm" (default) The study area is split at UTM zone boundaries (6 degree longitude bands, zone = floor((lon + 180) / 6) + 1) and at the equator. Each part is tiled on its own WGS 84 / UTM grid (EPSG:326xx north, 327xx south), whose origin is the lower-left corner of that part rounded down to whole metres. Tile IDs are "<zone><N|S>_<col>_<row>", e.g. "55S_003_012". crs="equal-area" One Lambert azimuthal equal-area grid centred on the study-area centroid (+proj=laea, WGS 84). Tile IDs are "ea_<col>_<row>". Distances and tile sizes are close to true within about 1,000 km of the centre.

In both cases every tile's composite is exported in the UTM zone of the tile (export_crs="EPSG:<utm_epsg>") unless the base configuration sets an explicit EPSG: code.

Cores and halos

The core of a tile is its grid cell; with clip=True (default) it is the cell intersected with the study area, otherwise the full cell (always cut at UTM zone and equator boundaries for crs="utm"). The halo is the core's bounding box in the grid CRS expanded by halo_m on every side. The tile's study_area is the EPSG:4326 bounding box of the halo (a "bbox:..." string), so the downloaded composite covers at least the halo.

Halo rule

A field is delineated whole only if it lies completely inside the halo of the tile that owns it. Every field that crosses a core boundary extends at most its own size past that boundary, so the halo must exceed the largest expected field dimension: fields larger than the halo that cross a core boundary can be truncated. Centre pivots are about 800 m across, so the default halo_m=1000 is a minimum; use more for regions with larger fields. :func:merge_tiles counts kept polygons that reach the edge of their halo (n_reaching_halo_edge), which flags possible truncation.

Merge rule

A polygon from tile T is kept only if T owns its representative point (:func:shapely.point_on_surface, computed in EPSG:4326). Ownership is decided arithmetically from the grid definition (zone/hemisphere of the point, then floor of its grid coordinates), which assigns every point to exactly one grid cell; with clip=True the point must also lie in the study area. Identical polygons delineated by two neighbouring tiles are therefore kept exactly once. Two different detections of the same field (one per tile) have different representative points and can, rarely, both be kept or both be dropped; :func:merge_tiles reports overlapping polygons from different tiles (n_cross_tile_overlap_pairs) so this can be checked. Study areas that cross the antimeridian are not supported.

With clip=True this representative-point test is the only selection at the region's study-area outline; with clip=False there is none. The base configuration's aoi_selection is applied in each tile run to the tile's study area, which is the tile's halo box, not to the region's outline. Reference polygons given to :func:merge_tiles are selected with the matching rule (representative point with clip=True, else intersecting the study area).

make_tiles

make_tiles(study_area: Any, tile_size_m: float = 20000.0, halo_m: float = 1000.0, crs: str = 'utm', clip: bool = True, config: Any = None) -> gpd.GeoDataFrame

Cut a study area into square tiles with halos.

Parameters:

Name	Type	Description	Default
study_area	str, GeoDataFrame, GeoSeries or shapely geometry	Anything :func:agribound.io.vector.read_study_area accepts (vector file, "bbox:minx,miny,maxx,maxy", WKT, GEE asset ID), or an in-memory geometry (a bare shapely geometry is taken as EPSG:4326).	required
tile_size_m	float	Edge length of the core cells in metres of the grid CRS (default 20,000).	20000.0
halo_m	float	Margin added around each core's bounding box (default 1,000 m). See the halo rule in the module docstring.	1000.0
crs	str	"utm" (per-zone UTM grids, default) or "equal-area" (one Lambert azimuthal equal-area grid centred on the study area).	'utm'
clip	bool	Intersect the cores with the study area (default True). With False the cores are whole grid cells (still cut at UTM zone and equator boundaries for crs="utm").	True
config	AgriboundConfig or None	Base configuration; used only to read a GEE-asset study area with the configured Earth Engine project and credentials (:func:agribound.auth.ensure_gee).	None

Returns:

Type

Description

GeoDataFrame

One row per tile in EPSG:4326, ordered by grid system, row and column, with columns index, tile_id, system, grid_crs, zone and hemisphere (UTM grids; None otherwise), col, row, utm_epsg (EPSG code of the tile's UTM zone, used as the export CRS), core_area_km2 (in the grid CRS), halo_minx, halo_miny, halo_maxx, halo_maxy (halo box in the grid CRS), bounds (EPSG:4326 bounds of the halo), study_area (the same bounds as a "bbox:..." string), the core as the active geometry and the halo as a second geometry column halo. attrs["grid"] holds the grid definition used by :func:assign_tile_ids, attrs["study_area_geometry"] the unioned study area and attrs["study_area_label"] its description.

Raises:

Type	Description
ValueError	For invalid sizes, an unknown crs, or an empty study area.

Source code in agribound/hpc/tiles.py

def make_tiles(
    study_area: Any,
    tile_size_m: float = 20000.0,
    halo_m: float = 1000.0,
    crs: str = "utm",
    clip: bool = True,
    config: Any = None,
) -> gpd.GeoDataFrame:
    """Cut a study area into square tiles with halos.

    Parameters
    ----------
    study_area : str, GeoDataFrame, GeoSeries or shapely geometry
        Anything :func:`agribound.io.vector.read_study_area` accepts (vector
        file, ``"bbox:minx,miny,maxx,maxy"``, WKT, GEE asset ID), or an
        in-memory geometry (a bare shapely geometry is taken as EPSG:4326).
    tile_size_m : float
        Edge length of the core cells in metres of the grid CRS (default
        20,000).
    halo_m : float
        Margin added around each core's bounding box (default 1,000 m). See
        the halo rule in the module docstring.
    crs : str
        ``"utm"`` (per-zone UTM grids, default) or ``"equal-area"`` (one
        Lambert azimuthal equal-area grid centred on the study area).
    clip : bool
        Intersect the cores with the study area (default *True*). With
        *False* the cores are whole grid cells (still cut at UTM zone and
        equator boundaries for ``crs="utm"``).
    config : AgriboundConfig or None
        Base configuration; used only to read a GEE-asset study area with the
        configured Earth Engine project and credentials
        (:func:`agribound.auth.ensure_gee`).

    Returns
    -------
    geopandas.GeoDataFrame
        One row per tile in EPSG:4326, ordered by grid system, row and column,
        with columns ``index``, ``tile_id``, ``system``, ``grid_crs``,
        ``zone`` and ``hemisphere`` (UTM grids; *None* otherwise), ``col``,
        ``row``, ``utm_epsg`` (EPSG code of the tile's UTM zone, used as the
        export CRS), ``core_area_km2`` (in the grid CRS), ``halo_minx``,
        ``halo_miny``, ``halo_maxx``, ``halo_maxy`` (halo box in the grid
        CRS), ``bounds`` (EPSG:4326 bounds of the halo), ``study_area`` (the
        same bounds as a ``"bbox:..."`` string), the core as the active
        ``geometry`` and the halo as a second geometry column ``halo``.
        ``attrs["grid"]`` holds the grid definition used by
        :func:`assign_tile_ids`, ``attrs["study_area_geometry"]`` the unioned
        study area and ``attrs["study_area_label"]`` its description.

    Raises
    ------
    ValueError
        For invalid sizes, an unknown *crs*, or an empty study area.
    """
    tile_size_m = float(tile_size_m)
    halo_m = float(halo_m)
    if not math.isfinite(tile_size_m) or tile_size_m <= 0:
        raise ValueError(f"tile_size_m must be > 0, got {tile_size_m}")
    if not math.isfinite(halo_m) or halo_m < 0:
        raise ValueError(f"halo_m must be >= 0, got {halo_m}")
    crs = str(crs).lower().strip()
    if crs not in GRID_KINDS:
        raise ValueError(f"crs must be one of {GRID_KINDS}, got {crs!r}")
    if halo_m < 1000:
        logger.warning(
            "halo_m=%.0f m: fields larger than the halo that cross a core boundary can be "
            "truncated (centre pivots are ~800 m across).",
            halo_m,
        )
    if halo_m > tile_size_m:
        logger.warning(
            "halo_m=%.0f m exceeds tile_size_m=%.0f m: each composite covers >9x its core.",
            halo_m,
            tile_size_m,
        )

    from agribound.io.crs import get_utm_crs

    aoi, label = _study_area_geometry(study_area, config=config)

    if crs == "utm":
        parts = _utm_parts(aoi)
    else:
        parts = [
            {
                "system": "ea",
                "zone": None,
                "hemisphere": None,
                "crs_string": _laea_definition(aoi),
                "part": aoi,
            }
        ]

    systems: dict[str, dict[str, Any]] = {}
    records: list[dict[str, Any]] = []
    halos_4326: list[Any] = []
    seg_m = max(10.0, tile_size_m / 20.0)

    for part in parts:
        crs_string = f"EPSG:{part['epsg']}" if crs == "utm" else part["crs_string"]
        grid_crs = pyproj.CRS.from_user_input(crs_string)
        to_grid = pyproj.Transformer.from_crs("EPSG:4326", grid_crs, always_xy=True)
        to_4326 = pyproj.Transformer.from_crs(grid_crs, "EPSG:4326", always_xy=True)
        part_dense = shapely.segmentize(part["part"], _AOI_SEGMENT_DEG)
        part_grid = shapely.make_valid(_project(part_dense, to_grid))
        # Whole-metre origin; the 1e-6 m slack absorbs reprojection round-off so that
        # an AOI starting at 600000 m does not get its origin at 599999 m.
        origin = (
            math.floor(part_grid.bounds[0] + 1e-6),
            math.floor(part_grid.bounds[1] + 1e-6),
        )

        if clip:
            region = part_grid
        elif crs == "utm":
            # Whole cells, cut only at the zone/equator boundaries.
            minx, miny, maxx, maxy = part["part"].bounds
            lat_abs = min(89.0, max(abs(miny), abs(maxy)))
            margin = 2.0 * tile_size_m / (111_320.0 * max(math.cos(math.radians(lat_abs)), 0.05))
            band = part["band"]
            window = box(
                max(band.bounds[0], minx - margin),
                max(band.bounds[1], miny - margin),
                min(band.bounds[2], maxx + margin),
                min(band.bounds[3], maxy + margin),
            )
            region = shapely.make_valid(
                _project(shapely.segmentize(window, _AOI_SEGMENT_DEG), to_grid)
            )
        else:
            region = None

        cols, rows, cells = _grid_cells(part_grid, origin, tile_size_m)
        if region is None:
            cores = cells
        else:
            cores = np.array([_polygonal(g) for g in shapely.intersection(cells, region)])
        systems[part["system"]] = {
            "crs": crs_string,
            "zone": part["zone"],
            "hemisphere": part["hemisphere"],
            "origin": [float(origin[0]), float(origin[1])],
        }

        for col, row, core in zip(cols, rows, cores, strict=True):
            if core.is_empty or core.area <= 0:
                continue
            cx0, cy0, cx1, cy1 = core.bounds
            halo = box(cx0 - halo_m, cy0 - halo_m, cx1 + halo_m, cy1 + halo_m)
            core_4326 = _polygonal(
                shapely.make_valid(_project(shapely.segmentize(core, seg_m), to_4326))
            )
            halo_4326 = _project(shapely.segmentize(halo, seg_m), to_4326)
            if crs == "utm":
                utm_epsg = int(part["epsg"])
            else:
                c = core_4326.centroid
                utm_epsg = int(get_utm_crs(c.x, c.y).to_epsg())
            hb = tuple(float(v) for v in halo_4326.bounds)
            records.append(
                {
                    "tile_id": _tile_id(part["system"], col, row),
                    "system": part["system"],
                    "grid_crs": systems[part["system"]]["crs"],
                    "zone": part["zone"],
                    "hemisphere": part["hemisphere"],
                    "col": int(col),
                    "row": int(row),
                    "utm_epsg": utm_epsg,
                    "core_area_km2": float(core.area) / 1e6,
                    "halo_minx": float(halo.bounds[0]),
                    "halo_miny": float(halo.bounds[1]),
                    "halo_maxx": float(halo.bounds[2]),
                    "halo_maxy": float(halo.bounds[3]),
                    "bounds": hb,
                    "study_area": _fmt_bbox(hb),
                    "geometry": core_4326,
                }
            )
            halos_4326.append(halo_4326)

    if not records:
        raise ValueError(f"No tiles intersect the study area {label}")

    order = sorted(
        range(len(records)),
        key=lambda i: (records[i]["system"], records[i]["row"], records[i]["col"]),
    )
    records = [records[i] for i in order]
    halos_4326 = [halos_4326[i] for i in order]
    for i, rec in enumerate(records):
        rec["index"] = i

    columns = ["index", "tile_id"] + [c for c in records[0] if c not in ("index", "tile_id")]
    frame = pd.DataFrame.from_records(records, columns=columns)
    tiles = gpd.GeoDataFrame(frame, geometry="geometry", crs="EPSG:4326")
    tiles["halo"] = gpd.GeoSeries(halos_4326, crs="EPSG:4326", index=tiles.index)
    tiles.attrs["grid"] = {
        "kind": crs,
        "tile_size_m": tile_size_m,
        "halo_m": halo_m,
        "clip": bool(clip),
        "systems": systems,
    }
    tiles.attrs["study_area_geometry"] = aoi
    tiles.attrs["study_area_label"] = label
    logger.info(
        "make_tiles: %d tiles (%s grid, %.0f m cores, %.0f m halo) over %s",
        len(tiles),
        crs,
        tile_size_m,
        halo_m,
        label,
    )
    return tiles

write_tile_manifest

write_tile_manifest(tiles: GeoDataFrame, base_config: Any, out_dir: str | Path, *, cache_root: str | Path | None = None, overwrite: bool = False, keep_reference: bool = False, allow_fine_tune_per_tile: bool = False) -> Path

Write the tile manifest, per-tile configurations and tiles.gpkg.

Files written under out_dir:

• tiles.gpkg with layers tiles (cores), halos and study_area (EPSG:4326);
• tiles.txt: one tile ID per line, in index order (line i + 1 is tile i, i.e. Slurm array index i);
• tiles/<tile_id>/config.yaml: the base configuration with study_area set to the tile's halo bbox, output_path set to tiles/<tile_id>/fields.<ext>, cache_dir set to <cache_root>/<tile_id>, export_crs set to the tile's UTM EPSG code when the base uses "utm", overwrite=False and provenance=True (both needed for idempotent restarts; a base value that differs is replaced with a WARNING), and relative paths of the base configuration made absolute against the current directory, so tile jobs can start in any directory: local_tif_path, gee_service_account_key, reference_boundaries, sam_model and the engine parameters checkpoint_path, weights_path, model_path when they name an existing file or directory (a relative local_tif_path, key or reference that does not exist is kept, with a WARNING), and embedding_cache_dir always;
• manifest.json (paths relative to out_dir).

Writing is idempotent: an existing manifest with the same signature (output directory, grid, tile entries, base configuration and cache root) is kept, and a different one is refused unless overwrite.

Parameters:

Name	Type	Description	Default
tiles	GeoDataFrame	Output of :func:make_tiles.	required
base_config	(AgriboundConfig, dict or path)	Configuration shared by all tiles (a YAML path is loaded with :meth:AgriboundConfig.from_yaml). Its study_area and output_path are replaced per tile.	required
out_dir	str or Path	Directory for the manifest and all tile outputs.	required
cache_root	(str, Path or None)	Parent of the per-tile cache directories. Default: the base configuration's cache_dir if set, else <out_dir>/tiles (i.e. tiles/<tile_id>/cache). Tile IDs depend only on the study area and tiling parameters, so manifests for several engines over the same source and year can share one cache_root and hence one download per tile -- provided only one job writes a given tile's cache at a time (see examples/hpc/README.md).	None
overwrite	bool	Replace an existing manifest that differs from the new one. Tile outputs are never deleted; outputs whose configuration changed are protected by their provenance hash.	False
keep_reference	bool	Keep reference_boundaries in the tile configurations. By default it is removed (with an INFO message), because per-tile evaluation against the reference polygons that intersect each halo double-counts fields; evaluate the merged output instead (merge_tiles(..., reference=...)).	False
allow_fine_tune_per_tile	bool	Allow fine_tune=True, which fine-tunes a separate model in every tile on the reference polygons inside that tile. Refused by default: fine-tune once, then pass the checkpoint with engine_params["checkpoint_path"].	False

Returns:

Type	Description
Path	Path of manifest.json.

Raises:

Type	Description
FileExistsError	If a different manifest exists in out_dir and overwrite is False.
ValueError	If tiles did not come from :func:make_tiles, or the base configuration fine-tunes and allow_fine_tune_per_tile is False.

Source code in agribound/hpc/tiles.py

def write_tile_manifest(
    tiles: gpd.GeoDataFrame,
    base_config: Any,
    out_dir: str | Path,
    *,
    cache_root: str | Path | None = None,
    overwrite: bool = False,
    keep_reference: bool = False,
    allow_fine_tune_per_tile: bool = False,
) -> Path:
    """Write the tile manifest, per-tile configurations and ``tiles.gpkg``.

    Files written under *out_dir*:

    - ``tiles.gpkg`` with layers ``tiles`` (cores), ``halos`` and
      ``study_area`` (EPSG:4326);
    - ``tiles.txt``: one tile ID per line, in index order (line ``i + 1`` is
      tile ``i``, i.e. Slurm array index ``i``);
    - ``tiles/<tile_id>/config.yaml``: the base configuration with
      ``study_area`` set to the tile's halo bbox, ``output_path`` set to
      ``tiles/<tile_id>/fields.<ext>``, ``cache_dir`` set to
      ``<cache_root>/<tile_id>``, ``export_crs`` set to the tile's UTM EPSG
      code when the base uses ``"utm"``, ``overwrite=False`` and
      ``provenance=True`` (both needed for idempotent restarts; a base value
      that differs is replaced with a WARNING), and relative paths of the
      base configuration made absolute against the current directory, so
      tile jobs can start in any directory: ``local_tif_path``,
      ``gee_service_account_key``, ``reference_boundaries``, ``sam_model``
      and the engine parameters ``checkpoint_path``, ``weights_path``,
      ``model_path`` when they name an existing file or directory (a
      relative ``local_tif_path``, key or reference that does not exist is
      kept, with a WARNING), and ``embedding_cache_dir`` always;
    - ``manifest.json`` (paths relative to *out_dir*).

    Writing is idempotent: an existing manifest with the same signature
    (output directory, grid, tile entries, base configuration and cache
    root) is kept, and a different one is refused unless *overwrite*.

    Parameters
    ----------
    tiles : geopandas.GeoDataFrame
        Output of :func:`make_tiles`.
    base_config : AgriboundConfig, dict or path
        Configuration shared by all tiles (a YAML path is loaded with
        :meth:`AgriboundConfig.from_yaml`). Its ``study_area`` and
        ``output_path`` are replaced per tile.
    out_dir : str or Path
        Directory for the manifest and all tile outputs.
    cache_root : str, Path or None
        Parent of the per-tile cache directories. Default: the base
        configuration's ``cache_dir`` if set, else ``<out_dir>/tiles``
        (i.e. ``tiles/<tile_id>/cache``). Tile IDs depend only on the study
        area and tiling parameters, so manifests for several engines over the
        same source and year can share one *cache_root* and hence one
        download per tile -- provided only one job writes a given tile's
        cache at a time (see ``examples/hpc/README.md``).
    overwrite : bool
        Replace an existing manifest that differs from the new one. Tile
        outputs are never deleted; outputs whose configuration changed are
        protected by their provenance hash.
    keep_reference : bool
        Keep ``reference_boundaries`` in the tile configurations. By default
        it is removed (with an INFO message), because per-tile evaluation
        against the reference polygons that intersect each halo double-counts
        fields; evaluate the merged output instead (``merge_tiles(...,
        reference=...)``).
    allow_fine_tune_per_tile : bool
        Allow ``fine_tune=True``, which fine-tunes a separate model in every
        tile on the reference polygons inside that tile. Refused by default:
        fine-tune once, then pass the checkpoint with
        ``engine_params["checkpoint_path"]``.

    Returns
    -------
    pathlib.Path
        Path of ``manifest.json``.

    Raises
    ------
    FileExistsError
        If a different manifest exists in *out_dir* and *overwrite* is False.
    ValueError
        If *tiles* did not come from :func:`make_tiles`, or the base
        configuration fine-tunes and *allow_fine_tune_per_tile* is False.
    """
    from agribound._version import __version__
    from agribound.provenance import config_hash, to_jsonable

    grid = tiles.attrs.get("grid")
    if not grid or len(tiles) == 0:
        raise ValueError("tiles must be a non-empty GeoDataFrame returned by make_tiles()")
    config = _as_config(base_config)
    path_overrides = _absolute_path_overrides(config)
    if path_overrides:
        config = config.merged(**path_overrides)
    if config.fine_tune and not allow_fine_tune_per_tile:
        raise ValueError(
            "The base configuration has fine_tune=True, which would fine-tune a separate model "
            "in every tile. Fine-tune once (e.g. 'agribound delineate --fine-tune --reference "
            "REF' over the reference area), then pass the checkpoint to the tiles with "
            "engine_params checkpoint_path=<path> (CLI: --engine-param checkpoint_path=...). "
            "Pass allow_fine_tune_per_tile=True (--allow-fine-tune-per-tile) to override."
        )

    out_dir = Path(out_dir).expanduser().resolve()
    tiles_dir = out_dir / "tiles"
    if cache_root is None:
        cache_root_path = (
            Path(config.cache_dir).expanduser().resolve() if config.cache_dir else None
        )
    else:
        cache_root_path = Path(cache_root).expanduser().resolve()

    if config.overwrite:
        logger.warning(
            "Tile configurations use overwrite=False (base had overwrite=True); "
            "use 'agribound tiles run --overwrite' to recompute tiles."
        )
    if not config.provenance:
        logger.warning(
            "Tile configurations use provenance=True (base had provenance=False): "
            "provenance records are required to skip finished tiles on restart."
        )
    drop_reference = (
        bool(config.reference_boundaries) and not keep_reference and not config.fine_tune
    )
    if drop_reference:
        logger.info(
            "reference_boundaries removed from the tile configurations; evaluate the merged "
            "output instead ('agribound tiles merge --reference %s').",
            config.reference_boundaries,
        )

    ext = config.get_output_extension()
    tile_entries: list[dict[str, Any]] = []
    tile_configs: list[Any] = []
    for row in tiles.itertuples(index=False):
        tile_dir = tiles_dir / row.tile_id
        cache_dir = (cache_root_path / row.tile_id) if cache_root_path else (tile_dir / "cache")
        overrides: dict[str, Any] = {
            "study_area": row.study_area,
            "output_path": str(tile_dir / f"fields{ext}"),
            "cache_dir": str(cache_dir),
            "overwrite": False,
            "provenance": True,
        }
        if config.export_crs == "utm":
            overrides["export_crs"] = f"EPSG:{int(row.utm_epsg)}"
        if drop_reference:
            overrides["reference_boundaries"] = None
        tile_config = config.merged(**overrides)
        tile_configs.append(tile_config)
        rel_dir = tile_dir.relative_to(out_dir)
        tile_entries.append(
            {
                "index": int(row.index),
                "tile_id": row.tile_id,
                "system": row.system,
                "col": int(row.col),
                "row": int(row.row),
                "utm_epsg": int(row.utm_epsg),
                "core_area_km2": round(float(row.core_area_km2), 6),
                "halo_grid_bounds": [
                    float(row.halo_minx),
                    float(row.halo_miny),
                    float(row.halo_maxx),
                    float(row.halo_maxy),
                ],
                "study_area": row.study_area,
                "dir": str(rel_dir),
                "config": str(rel_dir / TILE_CONFIG_FILENAME),
                "output": str(rel_dir / f"fields{ext}"),
                "config_hash": config_hash(tile_config),
            }
        )

    base_dict = to_jsonable(config.to_dict())
    # cache_dir is excluded from the configuration hash, so the cache root is
    # part of the signature explicitly (a different root must not be ignored).
    signature = _sha1_json(
        {
            "out_dir": str(out_dir),
            "grid": grid,
            "tiles": tile_entries,
            "base": base_dict,
            "cache_root": str(cache_root_path) if cache_root_path else None,
        }
    )
    manifest = {
        "schema_version": MANIFEST_SCHEMA_VERSION,
        "agribound_version": __version__,
        "created_utc": _utc_now(),
        "signature": signature,
        "study_area": tiles.attrs.get("study_area_label"),
        "grid": grid,
        "base_config": base_dict,
        "base_config_hash": config_hash(config),
        "output_format": config.output_format,
        "cache_root": str(cache_root_path) if cache_root_path else None,
        "reference_boundaries": config.reference_boundaries,
        "approx_tile_raster_mb": _estimate_tile_raster_mb(config, tiles),
        "n_tiles": len(tile_entries),
        "files": {
            "tiles_gpkg": TILES_GPKG_FILENAME,
            "tile_list": TILE_LIST_FILENAME,
        },
        "tiles": tile_entries,
    }

    manifest_path = out_dir / MANIFEST_FILENAME
    existing = _read_json(manifest_path) if manifest_path.exists() else None
    if existing is not None and existing.get("signature") == signature:
        logger.info("Manifest %s is up to date (%d tiles)", manifest_path, len(tile_entries))
        _write_tile_configs(out_dir, tile_entries, tile_configs, only_missing=True)
        if not (out_dir / TILES_GPKG_FILENAME).exists():
            _write_tiles_gpkg(tiles, out_dir / TILES_GPKG_FILENAME)
        if not (out_dir / TILE_LIST_FILENAME).exists():
            _write_tile_list(out_dir, tile_entries)
        return manifest_path
    if manifest_path.exists() and not overwrite:
        raise FileExistsError(
            f"{manifest_path} exists and describes different tiles, a different base "
            "configuration or a different cache root. Use another --out-dir, or "
            "overwrite=True (--overwrite) to replace the manifest and tile configurations "
            "(tile outputs are kept)."
        )

    out_dir.mkdir(parents=True, exist_ok=True)
    _write_tile_configs(out_dir, tile_entries, tile_configs, only_missing=False)
    _write_tiles_gpkg(tiles, out_dir / TILES_GPKG_FILENAME)
    _write_tile_list(out_dir, tile_entries)
    _write_json_atomic(manifest_path, manifest)
    logger.info("Wrote %s (%d tiles)", manifest_path, len(tile_entries))
    return manifest_path

load_manifest

load_manifest(manifest: str | Path | dict) -> dict[str, Any]

Load manifest.json (or a directory containing it).

Returns:

Type	Description
dict	The manifest with an extra key "_root" (absolute directory of the manifest, against which its relative paths are resolved).

Raises:

Type	Description
FileNotFoundError	If the manifest does not exist.
ValueError	If the file is not a manifest of a supported schema version.

Source code in agribound/hpc/tiles.py

def load_manifest(manifest: str | Path | dict) -> dict[str, Any]:
    """Load ``manifest.json`` (or a directory containing it).

    Returns
    -------
    dict
        The manifest with an extra key ``"_root"`` (absolute directory of the
        manifest, against which its relative paths are resolved).

    Raises
    ------
    FileNotFoundError
        If the manifest does not exist.
    ValueError
        If the file is not a manifest of a supported schema version.
    """
    if isinstance(manifest, dict):
        if "_root" not in manifest:
            raise ValueError("A manifest dict must come from load_manifest() (missing '_root')")
        return manifest
    path = Path(manifest).expanduser()
    if path.is_dir():
        path = path / MANIFEST_FILENAME
    if not path.exists():
        raise FileNotFoundError(f"Tile manifest not found: {path}")
    data = _read_json(path)
    if data is None or "tiles" not in data or "grid" not in data:
        raise ValueError(f"{path} is not an agribound tile manifest")
    if str(data.get("schema_version")) != MANIFEST_SCHEMA_VERSION:
        raise ValueError(
            f"{path} has manifest schema {data.get('schema_version')!r}; this agribound reads "
            f"schema {MANIFEST_SCHEMA_VERSION!r}. Re-create it with 'agribound tiles make'."
        )
    data["_root"] = str(path.parent.resolve())
    return data

load_tile_config

load_tile_config(manifest: str | Path | dict, tile: int | str) -> Any

Load the :class:~agribound.config.AgriboundConfig of one tile.

Source code in agribound/hpc/tiles.py

def load_tile_config(manifest: str | Path | dict, tile: int | str) -> Any:
    """Load the :class:`~agribound.config.AgriboundConfig` of one tile."""
    from agribound.config import AgriboundConfig

    m = load_manifest(manifest)
    return AgriboundConfig.from_yaml(_tile_paths(m, _tile_entry(m, tile))["config"])

run_tile

run_tile(manifest: str | Path | dict, index: int | str, stage: str = 'all', *, overwrite: bool = False) -> dict[str, Any]

Run one tile of a manifest (idempotent).

Parameters:

Name	Type	Description	Default
manifest	(str, Path or dict)	manifest.json, its directory, or a loaded manifest.	required
index	int or str	Tile index (Slurm array index, 0-based) or tile ID.	required
stage	str	"composite": build the composite / embeddings with :func:agribound.pipeline.build_composite, download the LULC raster when lulc_mode="raster" (a failed LULC download fails the stage, whatever lulc_on_error says, because the delineation stage would need network access), and for engine="ftw" with a two-window model on an Earth Engine source build the two seasonal window composites with the FTW engine's own window builder; then write the stage markers. "delineate": run :func:agribound.pipeline.delineate on a tile whose composite stage is done (raises otherwise), so the composite, LULC raster and FTW windows are read from the cache. "all": the composite stage without the LULC download, then :func:agribound.pipeline.delineate, which downloads the LULC raster itself and applies lulc_on_error.	'all'
overwrite	bool	Re-run even if the stage is already done, and retry tiles recorded as no-data (the delineation output is replaced; builders still reuse cached rasters -- delete the tile's cache directory to force new downloads). A re-run delineation gets a new run ID, so a merged output made before it no longer matches the tile outputs: :func:merge_tiles then raises :class:FileExistsError unless it is also given overwrite=True (agribound tiles merge --overwrite).	False

Returns:

Type	Description
dict	{"tile_id", "index", "stage", "status"} where status is "done" (ran now), "skipped" (already done) or "no-data" (no input data for the tile, see the module docstring; with "reason"), plus raster_path (composite stage) or output, n_output, run_id (delineation), and wall_s.

Raises:

Type	Description
RuntimeError	For stage="delineate" on a tile that has not been staged, or a failed LULC download in the composite stage.
Exception	Any other pipeline error is re-raised after <stage>.failed.json (error and traceback) is written to the tile directory.

Notes

Completion is decided from content, not from markers alone: a delineation is done when the tile output exists and its provenance record reports success with the current configuration hash and results versions (the pipeline's reuse test, :func:agribound.provenance.reuse_mismatch); a composite stage is done when the content-addressed stage markers in the tile's cache directory point to existing rasters (composite; LULC raster when needed; FTW window rasters when needed).

stage="delineate" guarantees only that these inputs are cached. Other network access during delineation is not prevented: model weights (agribound tiles prefetch), and Earth Engine for lulc_mode="server". Ensemble tiles stage the window composites of their two-window FTW members; an ensemble member whose staging failed with on_member_error="skip" builds its inputs again (online) during delineation or is skipped.

Source code in agribound/hpc/tiles.py

def run_tile(
    manifest: str | Path | dict,
    index: int | str,
    stage: str = "all",
    *,
    overwrite: bool = False,
) -> dict[str, Any]:
    """Run one tile of a manifest (idempotent).

    Parameters
    ----------
    manifest : str, Path or dict
        ``manifest.json``, its directory, or a loaded manifest.
    index : int or str
        Tile index (Slurm array index, 0-based) or tile ID.
    stage : str
        ``"composite"``: build the composite / embeddings with
        :func:`agribound.pipeline.build_composite`, download the LULC raster
        when ``lulc_mode="raster"`` (a failed LULC download fails the stage,
        whatever ``lulc_on_error`` says, because the delineation stage would
        need network access), and for ``engine="ftw"`` with a two-window
        model on an Earth Engine source build the two seasonal window
        composites with the FTW engine's own window builder; then write the
        stage markers.
        ``"delineate"``: run :func:`agribound.pipeline.delineate` on a tile
        whose composite stage is done (raises otherwise), so the composite,
        LULC raster and FTW windows are read from the cache.
        ``"all"``: the composite stage without the LULC download, then
        :func:`agribound.pipeline.delineate`, which downloads the LULC raster
        itself and applies ``lulc_on_error``.
    overwrite : bool
        Re-run even if the stage is already done, and retry tiles recorded as
        no-data (the delineation output is replaced; builders still reuse
        cached rasters -- delete the tile's cache directory to force new
        downloads). A re-run delineation gets a new run ID, so a merged
        output made before it no longer matches the tile outputs:
        :func:`merge_tiles` then raises :class:`FileExistsError` unless it is
        also given ``overwrite=True`` (``agribound tiles merge --overwrite``).

    Returns
    -------
    dict
        ``{"tile_id", "index", "stage", "status"}`` where *status* is
        ``"done"`` (ran now), ``"skipped"`` (already done) or ``"no-data"``
        (no input data for the tile, see the module docstring; with
        ``"reason"``), plus ``raster_path`` (composite stage) or ``output``,
        ``n_output``, ``run_id`` (delineation), and ``wall_s``.

    Raises
    ------
    RuntimeError
        For ``stage="delineate"`` on a tile that has not been staged, or a
        failed LULC download in the composite stage.
    Exception
        Any other pipeline error is re-raised after ``<stage>.failed.json``
        (error and traceback) is written to the tile directory.

    Notes
    -----
    Completion is decided from content, not from markers alone: a
    delineation is done when the tile output exists and its provenance
    record reports success with the current configuration hash and results
    versions (the pipeline's reuse test,
    :func:`agribound.provenance.reuse_mismatch`); a composite
    stage is done when the content-addressed stage markers in the tile's
    cache directory point to existing rasters (composite; LULC raster when
    needed; FTW window rasters when needed).

    ``stage="delineate"`` guarantees only that these inputs are cached.
    Other network access during delineation is not prevented: model weights
    (``agribound tiles prefetch``), and Earth Engine for
    ``lulc_mode="server"``. Ensemble tiles stage the window composites of
    their two-window FTW members; an ensemble member whose staging failed
    with ``on_member_error="skip"`` builds its inputs again (online) during
    delineation or is skipped.
    """
    stage = str(stage).lower().strip()
    if stage not in STAGES:
        raise ValueError(f"stage must be one of {STAGES}, got {stage!r}")
    m = load_manifest(manifest)
    entry = _tile_entry(m, index)
    paths = _tile_paths(m, entry)
    from agribound.config import AgriboundConfig

    config = AgriboundConfig.from_yaml(paths["config"])
    t0 = time.perf_counter()
    base = {"tile_id": entry["tile_id"], "index": entry["index"], "stage": stage}

    if stage == "composite":
        return _run_composite_stage(config, entry, paths, base, t0, overwrite, require_lulc=True)
    if not overwrite:
        reused = _reuse_current_output(config, entry, paths, base)
        if reused is not None:
            return reused
    if stage == "all":
        staged = _run_composite_stage(config, entry, paths, base, t0, overwrite, require_lulc=False)
        if staged["status"] == _NO_DATA:
            return staged
    return _run_delineate_stage(config, entry, paths, base, t0, overwrite, stage)

tile_status

tile_status(manifest: str | Path | dict) -> pd.DataFrame

Report the state of every tile.

Returns:

Type

Description

DataFrame

One row per tile with columns index, tile_id, composite and delineate ("done", "failed", "pending" or "no-data" -- no input data, a final state, see the module docstring; the delineation can also be "stale" -- an output exists but was made with a different configuration, or by an agribound release whose results for it differ (:mod:agribound._results) -- and either column is "error" when the tile configuration cannot be loaded), n_output, run_id, error (first line of the latest failure, if any) and no_data_reason.

Source code in agribound/hpc/tiles.py

def tile_status(manifest: str | Path | dict) -> pd.DataFrame:
    """Report the state of every tile.

    Returns
    -------
    pandas.DataFrame
        One row per tile with columns ``index``, ``tile_id``, ``composite``
        and ``delineate`` (``"done"``, ``"failed"``, ``"pending"`` or
        ``"no-data"`` -- no input data, a final state, see the module
        docstring; the delineation can also be ``"stale"`` -- an output
        exists but was made with a different configuration, or by an
        agribound release whose results for it differ
        (:mod:`agribound._results`) -- and either
        column is ``"error"`` when the tile configuration cannot be loaded),
        ``n_output``, ``run_id``, ``error`` (first line of the latest
        failure, if any) and ``no_data_reason``.
    """
    from agribound.config import AgriboundConfig
    from agribound.provenance import read_provenance, reuse_mismatch

    m = load_manifest(manifest)
    rows = []
    for entry in m["tiles"]:
        paths = _tile_paths(m, entry)
        row: dict[str, Any] = {
            "index": entry["index"],
            "tile_id": entry["tile_id"],
            "composite": _PENDING,
            "delineate": _PENDING,
            "n_output": None,
            "run_id": None,
            "error": None,
            "no_data_reason": None,
        }
        try:
            config = AgriboundConfig.from_yaml(paths["config"])
        except Exception as exc:
            row.update(composite="error", delineate="error", error=f"{type(exc).__name__}: {exc}")
            rows.append(row)
            continue

        no_data = _no_data_marker(config)
        if _staged_marker(config) is not None:
            row["composite"] = _DONE
        elif no_data is not None:
            row["composite"] = _NO_DATA
        elif paths["composite_failed"].exists():
            row["composite"] = _FAILED
            row["error"] = _short_error(paths["composite_failed"])

        output = paths["output"]
        record = read_provenance(output) if output.exists() else None
        if record is not None and record.get("status") == "success":
            if reuse_mismatch(record, config, output) is None:
                row["delineate"] = _DONE
                row["n_output"] = (record.get("facts") or {}).get("n_output")
                row["run_id"] = record.get("run_id")
                row["error"] = None
            else:
                row["delineate"] = _STALE
        elif no_data is not None:
            row["delineate"] = _NO_DATA
            row["no_data_reason"] = no_data.get("reason")
            row["error"] = None
        elif paths["delineate_failed"].exists():
            row["delineate"] = _FAILED
            row["error"] = _short_error(paths["delineate_failed"])
        elif record is not None and record.get("status") == "failed":
            row["delineate"] = _FAILED
            row["error"] = record.get("error")
        rows.append(row)
    return pd.DataFrame(
        rows,
        columns=[
            "index",
            "tile_id",
            "composite",
            "delineate",
            "n_output",
            "run_id",
            "error",
            "no_data_reason",
        ],
    )

merge_tiles

merge_tiles(manifest: str | Path | dict, output: str | Path | None = None, *, crs: str = 'EPSG:4326', allow_missing: bool = False, overwrite: bool = False, reference: str | Path | None = None, check_overlaps: bool = True) -> gpd.GeoDataFrame

Merge finished tile outputs into one vector file.

Each tile keeps only the polygons whose representative point it owns (see the merge rule in the module docstring); the kept polygons get an agribound:tile_id column, are reprojected to crs and concatenated. A summary is written to <output>.provenance.json.

Parameters:

Name	Type	Description	Default
manifest	(str, Path or dict)	Tile manifest.	required
output	(str, Path or None)	Output file (format from the extension). Default :func:default_merge_output.	None
crs	str	CRS of the merged output (default EPSG:4326, valid for any region).	'EPSG:4326'
allow_missing	bool	Merge even if some tiles are not done (they are listed in the summary and a WARNING is logged). Default: raise. Tiles recorded as "no-data" (see the module docstring) are never "missing": they are merged as empty and listed under no_data_tiles with their reasons, with a WARNING.	False
overwrite	bool	Recompute and replace an existing merged output. Without it, an output made from exactly the same tile runs is returned without recomputation, and one made from other tile runs raises :class:FileExistsError.	False
reference	(str, Path or None)	Reference boundaries; when given, they are restricted to the study area (the study_area layer of tiles.gpkg) with the rule the merged polygons follow: by representative point with clip=True, else those intersecting the study area. They are compared with the merged output using :func:agribound.evaluate.evaluate; the metrics are stored in the summary (evaluation, with the reference counts and rule in evaluation_reference) and in gdf.attrs["evaluation_metrics"]. Reference polygons in no-data tiles count as false negatives.	None
check_overlaps	bool	Count pairs of kept polygons from different tiles whose overlap exceeds half of the smaller polygon's area (default True).	True

Returns:

Type	Description
GeoDataFrame	The merged polygons. attrs["merge_summary"] holds the summary.

Raises:

Type	Description
RuntimeError	If tiles are missing and allow_missing is False, or if no tile is done and some are no-data (nothing to merge; usually the source has no data for the year).
FileExistsError	If output exists, was made from different tile runs, and overwrite is False.

Source code in agribound/hpc/tiles.py

def merge_tiles(
    manifest: str | Path | dict,
    output: str | Path | None = None,
    *,
    crs: str = "EPSG:4326",
    allow_missing: bool = False,
    overwrite: bool = False,
    reference: str | Path | None = None,
    check_overlaps: bool = True,
) -> gpd.GeoDataFrame:
    """Merge finished tile outputs into one vector file.

    Each tile keeps only the polygons whose representative point it owns
    (see the merge rule in the module docstring); the kept polygons get an
    ``agribound:tile_id`` column, are reprojected to *crs* and concatenated.
    A summary is written to ``<output>.provenance.json``.

    Parameters
    ----------
    manifest : str, Path or dict
        Tile manifest.
    output : str, Path or None
        Output file (format from the extension). Default
        :func:`default_merge_output`.
    crs : str
        CRS of the merged output (default EPSG:4326, valid for any region).
    allow_missing : bool
        Merge even if some tiles are not done (they are listed in the summary
        and a WARNING is logged). Default: raise. Tiles recorded as
        ``"no-data"`` (see the module docstring) are never "missing": they
        are merged as empty and listed under ``no_data_tiles`` with their
        reasons, with a WARNING.
    overwrite : bool
        Recompute and replace an existing merged output. Without it, an
        output made from exactly the same tile runs is returned without
        recomputation, and one made from other tile runs raises
        :class:`FileExistsError`.
    reference : str, Path or None
        Reference boundaries; when given, they are restricted to the study
        area (the ``study_area`` layer of ``tiles.gpkg``) with the rule the
        merged polygons follow: by representative point with ``clip=True``,
        else those intersecting the study area. They are compared with the
        merged output using :func:`agribound.evaluate.evaluate`; the metrics
        are stored in the summary (``evaluation``, with the reference counts
        and rule in ``evaluation_reference``) and in
        ``gdf.attrs["evaluation_metrics"]``. Reference polygons in no-data
        tiles count as false negatives.
    check_overlaps : bool
        Count pairs of kept polygons from different tiles whose overlap
        exceeds half of the smaller polygon's area (default *True*).

    Returns
    -------
    geopandas.GeoDataFrame
        The merged polygons. ``attrs["merge_summary"]`` holds the summary.

    Raises
    ------
    RuntimeError
        If tiles are missing and *allow_missing* is False, or if no tile is
        done and some are no-data (nothing to merge; usually the source has
        no data for the year).
    FileExistsError
        If *output* exists, was made from different tile runs, and
        *overwrite* is False.
    """
    from agribound._repro import collect_versions
    from agribound._version import __version__
    from agribound.io.vector import read_vector, write_vector
    from agribound.provenance import provenance_path, read_provenance, write_provenance

    m = load_manifest(manifest)
    grid = m["grid"]
    output = Path(output).expanduser() if output else default_merge_output(m)
    status = tile_status(m)
    done_mask = status["delineate"] == _DONE
    no_data_mask = status["delineate"] == _NO_DATA
    missing = status.loc[~(done_mask | no_data_mask), ["index", "tile_id", "delineate", "error"]]
    no_data = status.loc[no_data_mask, ["index", "tile_id", "no_data_reason"]].rename(
        columns={"no_data_reason": "reason"}
    )
    if len(missing) and not allow_missing:
        preview = ", ".join(missing["tile_id"].head(20))
        raise RuntimeError(
            f"{len(missing)} of {len(status)} tiles are not done ({preview}"
            f"{', ...' if len(missing) > 20 else ''}). Indices: "
            f"{format_index_ranges(missing['index'])}. Re-run them ('agribound tiles run') or "
            "pass allow_missing=True (--allow-missing)."
        )
    if not done_mask.any() and len(no_data):
        raise RuntimeError(
            f"No tile produced output: {len(no_data)} of {len(status)} tiles have no input data "
            f"(e.g. {no_data['tile_id'].iloc[0]}: {no_data['reason'].iloc[0]}). Check that the "
            "source covers the study area in this year ('agribound tiles status')."
        )

    signature = _sha1_json(
        {
            "tiles": sorted(
                zip(status.loc[done_mask, "tile_id"], status.loc[done_mask, "run_id"], strict=True)
            ),
            "crs": crs,
            "missing": sorted(missing["tile_id"]),
            "no_data": sorted(no_data["tile_id"]),
            "reference": str(reference) if reference else None,
            "manifest": m.get("signature"),
        }
    )
    if output.exists():
        previous = read_provenance(output)
        if not overwrite:
            if previous is not None and previous.get("inputs_signature") == signature:
                logger.info("Merged output %s is up to date", output)
                gdf = read_vector(output)
                gdf.attrs["merge_summary"] = previous
                gdf.attrs["reused"] = True
                return gdf
            raise FileExistsError(
                f"{output} exists and was not made from the current tile outputs. Pass "
                "overwrite=True (--overwrite) or choose another output."
            )

    study_area = _read_study_area_layer(m)
    shapely.prepare(study_area)
    aoi = study_area if grid.get("clip") else None
    entries = {e["tile_id"]: e for e in m["tiles"]}
    parts: list[gpd.GeoDataFrame] = []
    per_tile: list[dict[str, Any]] = []
    tile_records: list[dict[str, Any]] = []
    warnings: list[str] = []
    candidates: list[np.ndarray] = []

    for tile_id in status.loc[done_mask, "tile_id"]:
        entry = entries[tile_id]
        paths = _tile_paths(m, entry)
        record = read_provenance(paths["output"]) or {}
        tile_records.append(record)
        gdf = read_vector(paths["output"])
        info = {
            "tile_id": tile_id,
            "index": entry["index"],
            "run_id": record.get("run_id"),
            "config_hash": record.get("config_hash"),
            "n_in": len(gdf),
            "n_kept": 0,
            "n_reaching_halo_edge": 0,
            "wall_s": record.get("wall_s"),
        }
        if len(gdf) == 0:
            per_tile.append(info)
            continue
        if gdf.crs is None:
            raise ValueError(f"Tile output {paths['output']} has no CRS")
        points = representative_points(gdf.geometry)
        owner = assign_tile_ids(points, grid)
        keep = owner == tile_id
        if aoi is not None:
            keep &= shapely.covers(aoi, points)
        kept = gdf.loc[keep].copy()
        info["n_kept"] = len(kept)
        if len(kept):
            system = entry["system"]
            grid_crs = grid["systems"][system]["crs"]
            geoms_grid = np.asarray(kept.geometry.to_crs(grid_crs).values, dtype=object)
            halo_box = box(*entry["halo_grid_bounds"])
            reaching = ~shapely.within(geoms_grid, halo_box)
            info["n_reaching_halo_edge"] = int(reaching.sum())
            cell = cell_box(grid, system, entry["col"], entry["row"])
            candidates.append(~shapely.within(geoms_grid, cell))
            kept["agribound:tile_id"] = tile_id
            parts.append(kept.to_crs(crs))
        per_tile.append(info)

    if parts:
        merged = gpd.GeoDataFrame(pd.concat(parts, ignore_index=True), geometry="geometry", crs=crs)
    else:
        merged = gpd.GeoDataFrame(
            {"agribound:tile_id": []}, geometry=gpd.GeoSeries([], crs=crs), crs=crs
        )

    n_edge = int(sum(t["n_reaching_halo_edge"] for t in per_tile))
    if n_edge:
        msg = (
            f"{n_edge} kept polygons reach the edge of their tile's halo and may be truncated; "
            f"consider a halo larger than the largest field (current halo_m={grid['halo_m']})."
        )
        logger.warning(msg)
        warnings.append(msg)
    n_overlap = None
    if check_overlaps and len(merged):
        n_overlap = _cross_tile_overlap_pairs(merged, np.concatenate(candidates))
        if n_overlap:
            msg = (
                f"{n_overlap} pairs of polygons from different tiles overlap by more than half "
                "of the smaller polygon (the same field delineated by two tiles)."
            )
            logger.warning(msg)
            warnings.append(msg)
    if len(missing):
        msg = f"Merged with {len(missing)} tiles missing: {format_index_ranges(missing['index'])}"
        logger.warning(msg)
        warnings.append(msg)
    if len(no_data):
        msg = (
            f"{len(no_data)} tiles have no input data and contribute no polygons: "
            f"{format_index_ranges(no_data['index'])} (reasons in no_data_tiles)"
        )
        logger.warning(msg)
        warnings.append(msg)

    metrics = None
    reference_selection = None
    if reference:
        metrics, reference_selection = _evaluate_merged(
            merged, reference, study_area, bool(grid.get("clip"))
        )

    write_vector(merged, output)
    summary = _merge_summary(
        m,
        per_tile,
        tile_records,
        missing,
        no_data,
        merged,
        crs,
        n_edge,
        n_overlap,
        metrics,
        warnings,
        signature,
        output,
    )
    summary["evaluation_reference"] = reference_selection
    summary["agribound_version"] = __version__
    summary["versions"] = collect_versions()
    write_provenance(output, summary)
    merged.attrs["merge_summary"] = summary
    if metrics is not None:
        merged.attrs["evaluation_metrics"] = metrics
    logger.info(
        "Merged %d polygons from %d tiles -> %s (summary %s)",
        len(merged),
        int(done_mask.sum()),
        output,
        provenance_path(output).name,
    )
    return merged

representative_points

representative_points(geometries: GeoSeries) -> np.ndarray

Return :func:shapely.point_on_surface of each geometry, computed in EPSG:4326.

Source code in agribound/hpc/tiles.py

def representative_points(geometries: gpd.GeoSeries) -> np.ndarray:
    """Return :func:`shapely.point_on_surface` of each geometry, computed in EPSG:4326."""
    series = geometries
    if series.crs is not None and not series.crs.equals("EPSG:4326"):
        series = series.to_crs("EPSG:4326")
    return shapely.point_on_surface(np.asarray(series.values, dtype=object))

assign_tile_ids

assign_tile_ids(points_4326: Any, grid: dict[str, Any]) -> np.ndarray

Return the ID of the grid cell that owns each EPSG:4326 point.

The cell is found arithmetically: for UTM grids the point's zone (floor((lon + 180) / 6) + 1) and hemisphere (lat < 0 is south) select the grid system, then col = floor((x - x0) / size) and row = floor((y - y0) / size) in that system's CRS. Every point is assigned to exactly one cell; cells on a shared edge belong to the cell on their east/north side.

Parameters:

Name	Type	Description	Default
points_4326	array-like of shapely Points	Points in EPSG:4326 (e.g. from :func:representative_points).	required
grid	dict	Grid definition (tiles.attrs["grid"] or manifest["grid"]).	required

Returns:

Type	Description
ndarray	Object array of tile IDs; None for empty points or points outside every grid system. An ID is returned whether or not a tile with that ID exists; callers compare it with their tile's ID.

Source code in agribound/hpc/tiles.py

def assign_tile_ids(points_4326: Any, grid: dict[str, Any]) -> np.ndarray:
    """Return the ID of the grid cell that owns each EPSG:4326 point.

    The cell is found arithmetically: for UTM grids the point's zone
    (``floor((lon + 180) / 6) + 1``) and hemisphere (``lat < 0`` is south)
    select the grid system, then ``col = floor((x - x0) / size)`` and
    ``row = floor((y - y0) / size)`` in that system's CRS. Every point is
    assigned to exactly one cell; cells on a shared edge belong to the cell
    on their east/north side.

    Parameters
    ----------
    points_4326 : array-like of shapely Points
        Points in EPSG:4326 (e.g. from :func:`representative_points`).
    grid : dict
        Grid definition (``tiles.attrs["grid"]`` or ``manifest["grid"]``).

    Returns
    -------
    numpy.ndarray
        Object array of tile IDs; *None* for empty points or points outside
        every grid system. An ID is returned whether or not a tile with that
        ID exists; callers compare it with their tile's ID.
    """
    points = np.asarray(points_4326, dtype=object)
    result = np.full(points.shape, None, dtype=object)
    if points.size == 0:
        return result
    valid = ~shapely.is_empty(points) & ~shapely.is_missing(points)
    xs = np.full(points.shape, np.nan)
    ys = np.full(points.shape, np.nan)
    xs[valid] = shapely.get_x(points[valid])
    ys[valid] = shapely.get_y(points[valid])

    systems = grid["systems"]
    keys = np.full(points.shape, None, dtype=object)
    if grid["kind"] == "utm":
        # Same formula as agribound.io.crs.utm_zone_for_lon, vectorised.
        wrapped = np.mod(xs[valid] + 180.0, 360.0) - 180.0
        zones = np.mod(np.floor((wrapped + 180.0) / 6.0).astype(np.int64), 60) + 1
        hemis = np.where(ys[valid] < 0, "S", "N")
        keys[valid] = [f"{z:02d}{h}" for z, h in zip(zones, hemis, strict=True)]
    else:
        keys[valid] = "ea"

    size = float(grid["tile_size_m"])
    for system, info in systems.items():
        mask = keys == system
        if not mask.any():
            continue
        transformer = pyproj.Transformer.from_crs("EPSG:4326", info["crs"], always_xy=True)
        gx, gy = _transform_xy(transformer, xs[mask], ys[mask])
        x0, y0 = info["origin"]
        cols = np.floor((gx - x0) / size).astype(np.int64)
        rows = np.floor((gy - y0) / size).astype(np.int64)
        result[mask] = [_tile_id(system, c, r) for c, r in zip(cols, rows, strict=True)]
    return result

no_data_reason

no_data_reason(exc: BaseException) -> str | None

Return the no-data message in exc's chain, or None.

Walks exc, its __cause__ and __context__ (e.g. the FTW engine's :class:RuntimeError wrapping a window composite's :class:~agribound.composites.base.NoDataError) and returns "<ExceptionType>: <message>" of the first :class:~agribound.composites.base.NoDataError in the chain; if there is none, of the first :class:ValueError whose message matches :data:NO_DATA_PATTERNS (fallback).

Source code in agribound/hpc/tiles.py

def no_data_reason(exc: BaseException) -> str | None:
    """Return the no-data message in *exc*'s chain, or *None*.

    Walks ``exc``, its ``__cause__`` and ``__context__`` (e.g. the FTW
    engine's :class:`RuntimeError` wrapping a window composite's
    :class:`~agribound.composites.base.NoDataError`) and returns
    ``"<ExceptionType>: <message>"`` of the first
    :class:`~agribound.composites.base.NoDataError` in the chain; if there is
    none, of the first :class:`ValueError` whose message matches
    :data:`NO_DATA_PATTERNS` (fallback).
    """
    from agribound.composites.base import NoDataError

    chain: list[BaseException] = []
    seen: set[int] = set()
    current: BaseException | None = exc
    while current is not None and id(current) not in seen:
        seen.add(id(current))
        chain.append(current)
        current = current.__cause__ or current.__context__
    for item in chain:
        if isinstance(item, NoDataError):
            return f"{type(item).__name__}: {item}"
    for item in chain:
        if isinstance(item, ValueError):
            message = str(item)
            if any(p.search(message) for p in NO_DATA_PATTERNS):
                return f"{type(item).__name__}: {message}"
    return None

format_index_ranges

format_index_ranges(indices: Iterable[int]) -> str

Compress integers into a Slurm --array style list, e.g. "0-3,7,9-10".

Parameters:

Name	Type	Description	Default
indices	iterable of int	Indices (duplicates and order are ignored).	required

Returns:

Type	Description
str	Comma-separated ranges; empty string for no indices.

Source code in agribound/hpc/tiles.py

def format_index_ranges(indices: Iterable[int]) -> str:
    """Compress integers into a Slurm ``--array`` style list, e.g. ``"0-3,7,9-10"``.

    Parameters
    ----------
    indices : iterable of int
        Indices (duplicates and order are ignored).

    Returns
    -------
    str
        Comma-separated ranges; empty string for no indices.
    """
    values = sorted({int(i) for i in indices})
    if not values:
        return ""
    ranges: list[str] = []
    start = prev = values[0]
    for value in values[1:]:
        if value == prev + 1:
            prev = value
            continue
        ranges.append(f"{start}-{prev}" if prev > start else f"{start}")
        start = prev = value
    ranges.append(f"{start}-{prev}" if prev > start else f"{start}")
    return ",".join(ranges)

Regions

regions

Region definitions for large-area runs (examples/regions/*.yaml).

A region file describes one agricultural region: a verified EPSG:4326 bounding box, a small test box, recommended years, sources and engines, data availability notes and tiling parameters. The run block holds the machine-read defaults used by examples/run_region_delineation.sh::

name: iowa_corn_belt_us
title: ...
bbox: [minx, miny, maxx, maxy]        # EPSG:4326
test_bbox: [minx, miny, maxx, maxy]   # optional small box for quick checks
study_area: null                      # optional vector file (relative to this file)
run:
  years: [2023, 2024]
  sources: [sentinel2, landsat]
  engines: [delineate-anything, ftw]
  tile_size_km: 20
  halo_m: 1500
  reference: null                     # optional, relative to this file
  tessera_version: v1                 # optional
  lulc_dataset: auto                  # optional (--lulc-dataset)

All other keys are free-form documentation (verification results, source availability, reference data, FTW notes). :func:plan_runs expands the years x sources x engines matrix with the registry rules.

load_region

load_region(name_or_path: str | Path) -> dict[str, Any]

Load and validate a region file.

Returns:

Type	Description
dict	The YAML content plus "_path" (absolute file path), "study_area" resolved to an absolute path or a "bbox:..." string, "test_study_area" ("bbox:..." or None) and run["reference"] resolved to an absolute path (or None).

Raises:

Type	Description
ValueError	If required keys are missing or the boxes are invalid.

Source code in agribound/hpc/regions.py

def load_region(name_or_path: str | Path) -> dict[str, Any]:
    """Load and validate a region file.

    Returns
    -------
    dict
        The YAML content plus ``"_path"`` (absolute file path),
        ``"study_area"`` resolved to an absolute path or a ``"bbox:..."``
        string, ``"test_study_area"`` (``"bbox:..."`` or *None*) and
        ``run["reference"]`` resolved to an absolute path (or *None*).

    Raises
    ------
    ValueError
        If required keys are missing or the boxes are invalid.
    """
    path = find_region_file(name_or_path)
    with open(path) as f:
        data = yaml.safe_load(f) or {}
    if not isinstance(data, dict):
        raise ValueError(f"Region file {path} must contain a mapping")
    missing = [k for k in _REQUIRED if k not in data]
    if missing:
        raise ValueError(f"Region file {path} is missing keys {missing}")
    name = str(data["name"])
    run = data["run"] or {}
    missing = [k for k in _RUN_REQUIRED if k not in run]
    if missing:
        raise ValueError(f"Region file {path}: run block is missing keys {missing}")
    bbox = _check_bbox(data["bbox"], "bbox", name)
    test_bbox = data.get("test_bbox")
    if test_bbox is not None:
        test_bbox = _check_bbox(test_bbox, "test_bbox", name)

    base = path.parent
    study_area = data.get("study_area")
    if study_area:
        resolved = (base / study_area).resolve()
        if not resolved.exists():
            raise ValueError(f"Region {name!r}: study_area file {resolved} does not exist")
        data["study_area"] = str(resolved)
    else:
        data["study_area"] = bbox_string(bbox)
    data["test_study_area"] = bbox_string(test_bbox) if test_bbox else None
    reference = run.get("reference")
    run["reference"] = str((base / reference).resolve()) if reference else None
    for key in ("years", "sources", "engines"):
        value = run[key]
        run[key] = [value] if isinstance(value, str | int) else list(value)
    data["run"] = run
    data["_path"] = str(path)
    return data

find_region_file

find_region_file(name_or_path: str | Path) -> Path

Resolve a region name (e.g. "punjab_in") or a path to a YAML file.

Names are looked up as <name>.yaml in $AGRIBOUND_REGIONS_DIR (os.pathsep-separated), ./examples/regions and the examples/regions directory of the agribound source checkout.

Raises:

Type	Description
FileNotFoundError	If no region file is found (the message lists the searched places).

Source code in agribound/hpc/regions.py

def find_region_file(name_or_path: str | Path) -> Path:
    """Resolve a region name (e.g. ``"punjab_in"``) or a path to a YAML file.

    Names are looked up as ``<name>.yaml`` in ``$AGRIBOUND_REGIONS_DIR``
    (``os.pathsep``-separated), ``./examples/regions`` and the
    ``examples/regions`` directory of the agribound source checkout.

    Raises
    ------
    FileNotFoundError
        If no region file is found (the message lists the searched places).
    """
    path = Path(name_or_path).expanduser()
    if path.suffix in (".yaml", ".yml") or path.exists():
        if path.exists():
            return path.resolve()
        raise FileNotFoundError(f"Region file not found: {path}")
    searched = []
    for directory in _search_dirs():
        candidate = directory / f"{name_or_path}.yaml"
        searched.append(str(candidate))
        if candidate.exists():
            return candidate.resolve()
    raise FileNotFoundError(
        f"Region {name_or_path!r} not found. Searched: {searched}. Pass a YAML path or set "
        f"{REGIONS_DIR_ENV}."
    )

plan_runs

plan_runs(years: list[int], sources: list[str], engines: list[str], *, tessera_version: str | None = None, fine_tune: bool = False, has_checkpoint: bool = False, include_restricted: bool = False) -> list[dict[str, Any]]

Expand years x sources x engines into runs and skipped combinations.

A combination is skipped (with a reason) when the engine does not support the source (:func:agribound.registry.engine_supports_source), the year is outside the source's range (:func:agribound.registry.source_year_range, TESSERA by version), the source is restricted (SPOT) and include_restricted is False, or the engine is not label-free and neither fine_tune (the engine must be fine-tunable) nor has_checkpoint is set.

Returns:

Type	Description
list of dict	action ("run" or "skip"), year, source, engine, fine_tune (bool: fine-tune this engine on the reference), note (comma-separated flags: "gfm-env" for engines that need the GFM environment, "cpu" for engines whose registry entry has gpu_recommended=False; empty otherwise) and reason (for skips).

Raises:

Type	Description
ValueError	For unknown sources or engines.

Source code in agribound/hpc/regions.py

def plan_runs(
    years: list[int],
    sources: list[str],
    engines: list[str],
    *,
    tessera_version: str | None = None,
    fine_tune: bool = False,
    has_checkpoint: bool = False,
    include_restricted: bool = False,
) -> list[dict[str, Any]]:
    """Expand years x sources x engines into runs and skipped combinations.

    A combination is skipped (with a reason) when the engine does not support
    the source (:func:`agribound.registry.engine_supports_source`), the year
    is outside the source's range
    (:func:`agribound.registry.source_year_range`, TESSERA by version), the
    source is restricted (SPOT) and *include_restricted* is False, or the
    engine is not label-free and neither *fine_tune* (the engine must be
    fine-tunable) nor *has_checkpoint* is set.

    Returns
    -------
    list of dict
        ``action`` (``"run"`` or ``"skip"``), ``year``, ``source``,
        ``engine``, ``fine_tune`` (bool: fine-tune this engine on the
        reference), ``note`` (comma-separated flags: ``"gfm-env"`` for
        engines that need the GFM environment, ``"cpu"`` for engines whose
        registry entry has ``gpu_recommended=False``; empty otherwise) and
        ``reason`` (for skips).

    Raises
    ------
    ValueError
        For unknown sources or engines.
    """
    from agribound.registry import (
        ENGINE_REGISTRY,
        SOURCE_REGISTRY,
        engine_supports_source,
        source_year_range,
    )

    unknown = [s for s in sources if s not in SOURCE_REGISTRY]
    unknown += [e for e in engines if e not in ENGINE_REGISTRY]
    if unknown:
        raise ValueError(
            f"Unknown source/engine names {unknown}. Sources: {sorted(SOURCE_REGISTRY)}; "
            f"engines: {sorted(ENGINE_REGISTRY)}"
        )
    plan = []
    for year in years:
        for source in sources:
            version = tessera_version if source == "tessera-embedding" else None
            years_ok = source_year_range(source, tessera_version=version)
            for engine in engines:
                info = ENGINE_REGISTRY[engine]
                notes = ["gfm-env"] if engine in GFM_ENV_ENGINES else []
                if not info["gpu_recommended"]:
                    notes.append("cpu")
                row: dict[str, Any] = {
                    "action": "run",
                    "year": int(year),
                    "source": source,
                    "engine": engine,
                    "fine_tune": False,
                    "note": ",".join(notes),
                    "reason": "",
                }
                if not engine_supports_source(engine, source):
                    row.update(action="skip", reason=f"{engine} does not support {source}")
                elif years_ok is not None and not (
                    years_ok[0] <= int(year) <= (years_ok[1] or 9999)
                ):
                    last = years_ok[1] if years_ok[1] is not None else "present"
                    label = f" {version}" if version else ""
                    row.update(
                        action="skip",
                        reason=f"no {source}{label} data for {year} ({years_ok[0]}-{last})",
                    )
                elif SOURCE_REGISTRY[source].get("restricted") and not include_restricted:
                    row.update(action="skip", reason=f"{source} is restricted (--include-spot)")
                elif not info["label_free"]:
                    if fine_tune and info["fine_tunable"]:
                        row["fine_tune"] = True
                    elif not has_checkpoint:
                        row.update(
                            action="skip",
                            reason=(
                                f"{engine} has no label-free weights: pass --fine-tune (needs a "
                                "reference) or --engine-param checkpoint_path=..."
                            ),
                        )
                plan.append(row)
    return plan