Skip to content

Post-Processing

Polygonisation, merging, area filtering, smoothing, simplification, regularisation and the LULC crop filter.

postprocess

Post-processing utilities for field boundary polygons.

Provides shared functionality for polygonization, simplification, regularization, area filtering, and cross-tile merging.

filter_polygons 

filter_polygons(gdf: GeoDataFrame, min_area_m2: float = 2500.0, max_area_m2: float | None = None, remove_holes_below_m2: float | None = None, lulc_mask_path: str | None = None, lulc_agricultural_classes: list[int] | None = None, lulc_min_fraction: float = 0.5) -> gpd.GeoDataFrame

Filter field boundary polygons by area and an optional LULC raster.

Steps, in order:

  1. Area filter: keep polygons with min_area_m2 <= area (skipped when min_area_m2 <= 0) and area <= max_area_m2 (when given). Areas are computed in EPSG:6933 and exclude holes. Rows with a missing geometry have no area and are removed by an active area filter.
  2. Hole removal: interior rings whose own area is below remove_holes_below_m2 are filled.
  3. LULC mask (only when both lulc_mask_path and lulc_agricultural_classes are given): keep polygons for which more than lulc_min_fraction of the valid (non-nodata) raster pixels whose centres fall inside the polygon belong to lulc_agricultural_classes. The pipeline's LULC crop filter is :func:agribound.postprocess.lulc_filter.filter_by_lulc, not this.

No area columns are added to the output; the pipeline writes the final metrics:area/metrics:perimeter after all geometry edits.

Parameters:

Name Type Description Default
gdf GeoDataFrame

Input polygons (must have a CRS).

 required
min_area_m2 float

Minimum polygon area in m² (default 2500).

 2500.0
max_area_m2 float or None

Maximum polygon area in m². None means no upper limit.

 None
remove_holes_below_m2 float or None

Fill interior rings (holes) smaller than this area in m².

 None
lulc_mask_path str or None

Path to a classified LULC raster.

 None
lulc_agricultural_classes list[int] or None

LULC class values considered agricultural.

 None
lulc_min_fraction float

Minimum agricultural pixel fraction for the LULC mask (exclusive, default 0.5).

 0.5

Returns:

Type Description
GeoDataFrame

Filtered polygons in the input CRS (index reset, attributes kept).

Raises:

Type Description
ValueError

If gdf has no CRS (areas in m² cannot be computed).
Source code in agribound/postprocess/filter.py 
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
def filter_polygons(
    gdf: gpd.GeoDataFrame,
    min_area_m2: float = 2500.0,
    max_area_m2: float | None = None,
    remove_holes_below_m2: float | None = None,
    lulc_mask_path: str | None = None,
    lulc_agricultural_classes: list[int] | None = None,
    lulc_min_fraction: float = 0.5,
) -> gpd.GeoDataFrame:
    """Filter field boundary polygons by area and an optional LULC raster.

    Steps, in order:

    1. Area filter: keep polygons with ``min_area_m2 <= area`` (skipped when
       ``min_area_m2 <= 0``) and ``area <= max_area_m2`` (when given). Areas
       are computed in EPSG:6933 and exclude holes. Rows with a missing
       geometry have no area and are removed by an active area filter.
    2. Hole removal: interior rings whose own area is below
       ``remove_holes_below_m2`` are filled.
    3. LULC mask (only when both *lulc_mask_path* and
       *lulc_agricultural_classes* are given): keep polygons for which more
       than *lulc_min_fraction* of the valid (non-nodata) raster pixels whose
       centres fall inside the polygon belong to *lulc_agricultural_classes*.
       The pipeline's LULC crop filter is
       :func:`agribound.postprocess.lulc_filter.filter_by_lulc`, not this.

    No area columns are added to the output; the pipeline writes the final
    ``metrics:area``/``metrics:perimeter`` after all geometry edits.

    Parameters
    ----------
    gdf : geopandas.GeoDataFrame
        Input polygons (must have a CRS).
    min_area_m2 : float
        Minimum polygon area in m² (default 2500).
    max_area_m2 : float or None
        Maximum polygon area in m². *None* means no upper limit.
    remove_holes_below_m2 : float or None
        Fill interior rings (holes) smaller than this area in m².
    lulc_mask_path : str or None
        Path to a classified LULC raster.
    lulc_agricultural_classes : list[int] or None
        LULC class values considered agricultural.
    lulc_min_fraction : float
        Minimum agricultural pixel fraction for the LULC mask (exclusive,
        default 0.5).

    Returns
    -------
    geopandas.GeoDataFrame
        Filtered polygons in the input CRS (index reset, attributes kept).

    Raises
    ------
    ValueError
        If *gdf* has no CRS (areas in m² cannot be computed).
    """
    if len(gdf) == 0:
        return gdf
    if gdf.crs is None:
        raise ValueError("filter_polygons needs a GeoDataFrame with a CRS to compute areas in m²")

    result = gdf.copy()
    ea_crs = get_equal_area_crs()
    initial_count = len(result)

    # 1. Area filtering ------------------------------------------------------
    if min_area_m2 > 0 or max_area_m2 is not None:
        area = result.geometry.to_crs(ea_crs).area
        keep = np.ones(len(result), dtype=bool)
        if min_area_m2 > 0:
            keep &= (area >= min_area_m2).to_numpy()
        if max_area_m2 is not None:
            keep &= (area <= max_area_m2).to_numpy()
        result = result[keep]

    # 2. Remove small holes ----------------------------------------------------
    if remove_holes_below_m2 is not None and len(result) > 0:
        geoms_ea = result.geometry.to_crs(ea_crs)
        result[result.geometry.name] = [
            _remove_small_holes(g, g_ea, remove_holes_below_m2)
            for g, g_ea in zip(result.geometry, geoms_ea, strict=True)
        ]

    # 3. LULC filtering -----------------------------------------------------------
    if lulc_mask_path is not None and lulc_agricultural_classes is not None and len(result) > 0:
        result = _filter_by_lulc(
            result, lulc_mask_path, lulc_agricultural_classes, lulc_min_fraction
        )

    removed = initial_count - len(result)
    if removed > 0:
        logger.info(
            "Filtered %d polygons (min=%.0f m2, max=%s m2): %d remaining",
            removed,
            min_area_m2,
            max_area_m2 or "inf",
            len(result),
        )

    return result.reset_index(drop=True)

merge_polygons 

merge_polygons(gdf: GeoDataFrame, iou_threshold: float = 0.3, containment_threshold: float = 0.8, *, return_groups: bool = False) -> gpd.GeoDataFrame | tuple[gpd.GeoDataFrame, list[list[int]]]

Merge polygons that overlap strongly.

Two polygons a and b whose interiors or boundaries intersect are linked when IoU = area(a ∩ b) / area(a ∪ b) >= iou_threshold or when area(a ∩ b) / min(area(a), area(b)) >= containment_threshold. Areas are measured in the frame's CRS units (the ratios are unit-free). Links are transitive (union-find): if a is linked to b and b to c, all three are merged, even when a and c do not overlap. Each group of linked polygons is replaced by the union of its (repaired) geometries.

Attribute rule for merged rows: all non-geometry columns are copied from the group's largest input polygon (by area); values are not aggregated. If a group's union is a MultiPolygon it is split into one row per part, each carrying those attributes. Polygons that are not linked to any other are returned unchanged, attributes included. Output rows are ordered by the position of each group's first member in the input. gdf.attrs is preserved and attrs["merge_stats"] records the counts. A frame with fewer than two rows is returned as it is (the same object, without merge_stats).

Parameters:

Name Type Description Default
gdf GeoDataFrame

Input polygons (potentially with duplicates).

 required
iou_threshold float

IoU at or above which two polygons are merged (default 0.3).

 0.3
containment_threshold float

Fraction of the smaller polygon covered by the intersection at or above which two polygons are merged (default 0.8).

 0.8
return_groups bool

Also return, for every output row, the input positions (0-based) that were merged into it.

 False

Returns:

Type Description
GeoDataFrame or tuple

Polygons with all input columns (index reset); with return_groups, (frame, groups) where groups[i] lists the input positions of output row i.
Source code in agribound/postprocess/merge.py 
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
def merge_polygons(
    gdf: gpd.GeoDataFrame,
    iou_threshold: float = 0.3,
    containment_threshold: float = 0.8,
    *,
    return_groups: bool = False,
) -> gpd.GeoDataFrame | tuple[gpd.GeoDataFrame, list[list[int]]]:
    """Merge polygons that overlap strongly.

    Two polygons *a* and *b* whose interiors or boundaries intersect are
    linked when ``IoU = area(a ∩ b) / area(a ∪ b) >= iou_threshold`` or when
    ``area(a ∩ b) / min(area(a), area(b)) >= containment_threshold``. Areas
    are measured in the frame's CRS units (the ratios are unit-free). Links
    are transitive (union-find): if *a* is linked to *b* and *b* to *c*, all
    three are merged, even when *a* and *c* do not overlap. Each group of
    linked polygons is replaced by the union of its (repaired) geometries.

    Attribute rule for merged rows: all non-geometry columns are copied from
    the group's largest input polygon (by area); values are not aggregated.
    If a group's union is a MultiPolygon it is split into one row per part,
    each carrying those attributes. Polygons that are not linked to any other
    are returned unchanged, attributes included. Output rows are ordered by
    the position of each group's first member in the input. ``gdf.attrs`` is
    preserved and ``attrs["merge_stats"]`` records the counts. A frame with
    fewer than two rows is returned as it is (the same object, without
    ``merge_stats``).

    Parameters
    ----------
    gdf : geopandas.GeoDataFrame
        Input polygons (potentially with duplicates).
    iou_threshold : float
        IoU at or above which two polygons are merged (default 0.3).
    containment_threshold : float
        Fraction of the smaller polygon covered by the intersection at or
        above which two polygons are merged (default 0.8).
    return_groups : bool
        Also return, for every output row, the input positions (0-based)
        that were merged into it.

    Returns
    -------
    geopandas.GeoDataFrame or tuple
        Polygons with all input columns (index reset); with *return_groups*,
        ``(frame, groups)`` where ``groups[i]`` lists the input positions of
        output row ``i``.
    """
    if len(gdf) <= 1:
        return (gdf, [[i] for i in range(len(gdf))]) if return_groups else gdf

    geom_col = gdf.geometry.name
    geoms = np.asarray(gdf.geometry.array, dtype=object)
    n = len(geoms)
    present = np.array([g is not None and not g.is_empty for g in geoms], dtype=bool)

    # Repaired copies used for the overlap tests and unions (invalid rings would
    # make GEOS set operations fail).
    from agribound.postprocess.simplify import make_polygonal

    work = np.array(
        [
            make_polygonal(g) if ok and not g.is_valid else g
            for g, ok in zip(geoms, present, strict=True)
        ],
        dtype=object,
    )
    areas = np.array([g.area if ok else 0.0 for g, ok in zip(work, present, strict=True)])

    parent = np.arange(n)

    def find(i: int) -> int:
        while parent[i] != i:
            parent[i] = parent[parent[i]]
            i = parent[i]
        return int(i)

    idx = np.flatnonzero(present)
    if len(idx) > 1:
        tree = shapely.STRtree(work[idx])
        left, right = tree.query(work[idx], predicate="intersects")
        pairs = left < right
        a = idx[left[pairs]]
        b = idx[right[pairs]]
        if len(a):
            inter = shapely.area(shapely.intersection(work[a], work[b]))
            union_area = areas[a] + areas[b] - inter
            min_area = np.minimum(areas[a], areas[b])
            with np.errstate(divide="ignore", invalid="ignore"):
                iou = np.where(union_area > 0, inter / union_area, 0.0)
                contain = np.where(min_area > 0, inter / min_area, 0.0)
            link = (iou >= iou_threshold) | (contain >= containment_threshold)
            for i, j in zip(a[link], b[link], strict=True):
                ri, rj = find(int(i)), find(int(j))
                if ri != rj:
                    parent[ri] = rj

    groups: dict[int, list[int]] = {}
    for i in range(n):
        groups.setdefault(find(i), []).append(i)

    row_positions: list[int] = []
    out_geoms: list = []
    out_groups: list[list[int]] = []
    n_groups_merged = 0
    for members in groups.values():
        if len(members) == 1:
            row_positions.append(members[0])
            out_geoms.append(geoms[members[0]])
            out_groups.append(list(members))
            continue
        n_groups_merged += 1
        representative = max(members, key=lambda m: areas[m])
        merged = make_polygonal(shapely.union_all([work[m] for m in members]))
        parts = list(merged.geoms) if isinstance(merged, MultiPolygon) else [merged]
        for part in parts:
            row_positions.append(representative)
            out_geoms.append(part)
            out_groups.append(list(members))

    result = gdf.iloc[row_positions].copy()
    result[geom_col] = gpd.GeoSeries(out_geoms, index=result.index, crs=gdf.crs)
    result = result.reset_index(drop=True)
    result.attrs = dict(gdf.attrs)
    result.attrs["merge_stats"] = {
        "n_in": int(n),
        "n_out": int(len(result)),
        "n_groups_merged": int(n_groups_merged),
        "iou_threshold": float(iou_threshold),
        "containment_threshold": float(containment_threshold),
    }
    if n_groups_merged:
        logger.info(
            "Merged overlapping polygons: %d -> %d (%d groups)", n, len(result), n_groups_merged
        )
    return (result, out_groups) if return_groups else result

polygonize_mask 

polygonize_mask(mask_path: str, min_area_m2: float = 2500.0, band: int = 1, connectivity: int = 4, field_value: int | None = None) -> gpd.GeoDataFrame

Convert a raster segmentation mask to field boundary polygons.

Every connected region (under connectivity) of equal non-zero value becomes one polygon, with holes where other values are enclosed. Non-finite pixels and pixels equal to the raster's (finite) nodata value are treated as background (0). Float rasters are accepted only when all finite values are whole numbers (e.g. labels stored as float32); probability rasters must be thresholded first. Polygons that rasterio.features.shapes returns invalid are repaired with :func:agribound.postprocess.simplify.make_polygonal rather than dropped.

Parameters:

Name Type Description Default
mask_path str

Path to the segmentation mask GeoTIFF. Non-zero values are treated as field pixels (or only field_value if specified).

 required
min_area_m2 float

Minimum polygon area in m² to keep (default 2500; 0 keeps all), applied with :func:agribound.postprocess.filter.filter_polygons.

 2500.0
band int

Band index to polygonize (1-based, default 1).

 1
connectivity int

Pixel connectivity for grouping, 4 or 8 (default 4).

 4
field_value int or None

If set, only pixels equal to this value are polygonized. If None, all non-zero pixels are included.

 None

Returns:

Type Description
GeoDataFrame

Polygons with class_value (the pixel value of the region) and geometry columns, in the raster's CRS.

Raises:

Type Description
ValueError

If connectivity is not 4 or 8, or a float raster holds non-integer values.
Source code in agribound/postprocess/polygonize.py 
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
def polygonize_mask(
    mask_path: str,
    min_area_m2: float = 2500.0,
    band: int = 1,
    connectivity: int = 4,
    field_value: int | None = None,
) -> gpd.GeoDataFrame:
    """Convert a raster segmentation mask to field boundary polygons.

    Every connected region (under *connectivity*) of equal non-zero value
    becomes one polygon, with holes where other values are enclosed.
    Non-finite pixels and pixels equal to the raster's (finite) nodata value
    are treated as background (0). Float rasters are accepted only when all
    finite values are whole numbers (e.g. labels stored as float32);
    probability rasters must be thresholded first. Polygons that
    ``rasterio.features.shapes`` returns invalid are repaired with
    :func:`agribound.postprocess.simplify.make_polygonal` rather than dropped.

    Parameters
    ----------
    mask_path : str
        Path to the segmentation mask GeoTIFF. Non-zero values are
        treated as field pixels (or only *field_value* if specified).
    min_area_m2 : float
        Minimum polygon area in m² to keep (default 2500; ``0`` keeps all),
        applied with :func:`agribound.postprocess.filter.filter_polygons`.
    band : int
        Band index to polygonize (1-based, default 1).
    connectivity : int
        Pixel connectivity for grouping, 4 or 8 (default 4).
    field_value : int or None
        If set, only pixels equal to this value are polygonized.
        If *None*, all non-zero pixels are included.

    Returns
    -------
    geopandas.GeoDataFrame
        Polygons with ``class_value`` (the pixel value of the region) and
        ``geometry`` columns, in the raster's CRS.

    Raises
    ------
    ValueError
        If *connectivity* is not 4 or 8, or a float raster holds non-integer
        values.
    """
    if connectivity not in (4, 8):
        raise ValueError(f"connectivity must be 4 or 8, got {connectivity}")

    with warnings.catch_warnings():
        # Some engines write nodata values the dtype cannot represent (e.g. -inf for uint8).
        warnings.filterwarnings("ignore", message=".*nodata.*")
        with rasterio.open(mask_path) as src:
            mask = src.read(band)
            transform = src.transform
            crs = src.crs
            nodata = src.nodata

    background = ~np.isfinite(mask) if np.issubdtype(mask.dtype, np.floating) else None
    if nodata is not None and np.isfinite(nodata):
        is_nodata = mask == nodata
        background = is_nodata if background is None else (background | is_nodata)

    if np.issubdtype(mask.dtype, np.floating):
        finite = mask[np.isfinite(mask)]
        if finite.size and not np.all(finite == np.round(finite)):
            raise ValueError(
                f"{mask_path} band {band} holds non-integer values; polygonize_mask expects "
                "an integer label/mask raster (threshold probability rasters first)"
            )
        mask = np.where(np.isfinite(mask), mask, 0).astype(np.int32)
    elif mask.dtype not in (np.uint8, np.int16, np.uint16, np.int32, np.float32):
        # rasterio.features.shapes supports int16, int32, uint8, uint16 and float32 only.
        mask = mask.astype(np.int32)
    if background is not None and background.any():
        mask = np.where(background, 0, mask).astype(mask.dtype)

    from agribound.postprocess.simplify import make_polygonal

    polygons = []
    values = []
    n_repaired = 0
    for geom, val in rio_shapes(mask, connectivity=connectivity, transform=transform):
        if val == 0:
            continue
        if field_value is not None and int(val) != field_value:
            continue
        poly = shapely_shape(geom)
        if not poly.is_valid:
            poly = make_polygonal(poly)
            n_repaired += 1
        if poly is not None and not poly.is_empty:
            polygons.append(poly)
            values.append(int(val))
    if n_repaired:
        logger.info("Repaired %d invalid polygons from %s", n_repaired, mask_path)

    if not polygons:
        logger.warning("No polygons extracted from mask %s", mask_path)
        return gpd.GeoDataFrame({"class_value": np.array([], dtype=np.int64)}, geometry=[], crs=crs)

    gdf = gpd.GeoDataFrame({"class_value": values, "geometry": polygons}, crs=crs)

    if min_area_m2 > 0:
        from agribound.postprocess.filter import filter_polygons

        gdf = filter_polygons(gdf, min_area_m2=min_area_m2)

    logger.info("Polygonized %d features from %s", len(gdf), mask_path)
    return gdf

regularize_polygons 

regularize_polygons(gdf: GeoDataFrame, method: str = 'adaptive', *, simplify_tolerance_m: float = 0.5, area_threshold: float = 0.9, parallel_threshold_m: float = 1.0, allow_45_degree: bool = True, allow_circles: bool = True, circle_threshold: float = 0.9, num_cores: int = 1) -> gpd.GeoDataFrame

Regularize field boundary polygon geometry.

Parameters:

Name Type Description Default
gdf GeoDataFrame

Input polygons (with a CRS).

 required
method str

"adaptive" (default), "orthogonal" or "none" (returns gdf unchanged). See the module docstring for what each does.

 'adaptive'
simplify_tolerance_m float

Douglas-Peucker tolerance in metres used by both methods (geoai default 0.5).

 0.5
area_threshold float

"adaptive" only: a rectangle replacement is rejected when rect_area / area is below this value or above its inverse (geoai default 0.9).

 0.9
parallel_threshold_m float

"orthogonal" only: distance in metres below which nearly parallel adjacent edges are merged (default 1.0).

 1.0
allow_45_degree (bool, bool, float)

"orthogonal" only: allow 45-degree edges; replace polygons whose IoU with a fitted circle is at least circle_threshold by that circle (defaults True, True, 0.9).

 True
allow_circles (bool, bool, float)

"orthogonal" only: allow 45-degree edges; replace polygons whose IoU with a fitted circle is at least circle_threshold by that circle (defaults True, True, 0.9).

 True
circle_threshold (bool, bool, float)

"orthogonal" only: allow 45-degree edges; replace polygons whose IoU with a fitted circle is at least circle_threshold by that circle (defaults True, True, 0.9).

 True
num_cores int

"orthogonal" only: worker processes (default 1, sequential).

 1

Returns:

Type Description
GeoDataFrame

Regularized polygons in the input CRS with the input's columns (index reset). "adaptive" keeps one row per input row. "orthogonal" processes each polygon part separately, so a MultiPolygon row (or a polygon that buildingregulariser splits) yields one row per part with the same attributes; rows for which it returns no geometry keep their original geometry (the count is logged as a warning).

Raises:

Type Description
ValueError

If method is unknown or gdf has no CRS.
ImportError

If geoai-py (or, for "orthogonal", buildingregulariser) is not installed.
Source code in agribound/postprocess/regularize.py 
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
def regularize_polygons(
    gdf: gpd.GeoDataFrame,
    method: str = "adaptive",
    *,
    simplify_tolerance_m: float = 0.5,
    area_threshold: float = 0.9,
    parallel_threshold_m: float = 1.0,
    allow_45_degree: bool = True,
    allow_circles: bool = True,
    circle_threshold: float = 0.9,
    num_cores: int = 1,
) -> gpd.GeoDataFrame:
    """Regularize field boundary polygon geometry.

    Parameters
    ----------
    gdf : geopandas.GeoDataFrame
        Input polygons (with a CRS).
    method : str
        ``"adaptive"`` (default), ``"orthogonal"`` or ``"none"`` (returns
        *gdf* unchanged). See the module docstring for what each does.
    simplify_tolerance_m : float
        Douglas-Peucker tolerance in metres used by both methods (geoai
        default 0.5).
    area_threshold : float
        ``"adaptive"`` only: a rectangle replacement is rejected when
        ``rect_area / area`` is below this value or above its inverse
        (geoai default 0.9).
    parallel_threshold_m : float
        ``"orthogonal"`` only: distance in metres below which nearly parallel
        adjacent edges are merged (default 1.0).
    allow_45_degree, allow_circles, circle_threshold : bool, bool, float
        ``"orthogonal"`` only: allow 45-degree edges; replace polygons whose
        IoU with a fitted circle is at least *circle_threshold* by that
        circle (defaults True, True, 0.9).
    num_cores : int
        ``"orthogonal"`` only: worker processes (default 1, sequential).

    Returns
    -------
    geopandas.GeoDataFrame
        Regularized polygons in the input CRS with the input's columns (index
        reset). ``"adaptive"`` keeps one row per input row. ``"orthogonal"``
        processes each polygon part separately, so a MultiPolygon row (or a
        polygon that buildingregulariser splits) yields one row per part with
        the same attributes; rows for which it returns no geometry keep their
        original geometry (the count is logged as a warning).

    Raises
    ------
    ValueError
        If *method* is unknown or *gdf* has no CRS.
    ImportError
        If geoai-py (or, for ``"orthogonal"``, buildingregulariser) is not
        installed.
    """
    method = str(method).lower().strip()
    if method not in VALID_METHODS:
        raise ValueError(f"Unknown regularization method {method!r}. Choose from {VALID_METHODS}")
    if len(gdf) == 0 or method == "none":
        return gdf
    if gdf.crs is None:
        raise ValueError("regularize_polygons needs a GeoDataFrame with a CRS (metre parameters)")
    if not (~gdf.geometry.isna() & ~gdf.geometry.is_empty).any():
        return gdf  # nothing to regularize (geoai.regularize rejects empty input)

    try:
        from geoai.utils import geometry as geoai_geometry
    except ImportError as exc:
        raise ImportError(
            "Polygon regularization uses geoai-py (>= 0.43.1). Install it with: "
            "pip install 'agribound[geoai]'"
        ) from exc

    from agribound.postprocess.simplify import _to_original, make_polygonal, to_metric_crs

    metric, original_crs = to_metric_crs(gdf)
    geom_col = metric.geometry.name

    if method == "adaptive":
        present = ~metric.geometry.isna() & ~metric.geometry.is_empty
        subset = metric.loc[present, [geom_col]]
        out = geoai_geometry.adaptive_regularization(
            subset,
            simplify_tolerance=simplify_tolerance_m,
            area_threshold=area_threshold,
            preserve_shape=True,
        )
        new_geoms = [make_polygonal(g) for g in out.geometry]
        result = metric.copy()
        result.loc[present, geom_col] = gpd.GeoSeries(new_geoms, index=subset.index, crs=metric.crs)
        n_changed = sum(not a.equals(b) for a, b in zip(subset.geometry, new_geoms, strict=True))
        detail = f"{n_changed} geometries changed"
    else:
        if importlib.util.find_spec("buildingregulariser") is None:
            raise ImportError(
                "method='orthogonal' uses buildingregulariser through "
                "geoai.utils.geometry.regularize. Install it with: pip install buildingregulariser"
            )
        work = metric.copy()
        work[_ROW_ID] = np.arange(len(work))
        present = ~work.geometry.isna() & ~work.geometry.is_empty
        subset = work.loc[present, [_ROW_ID, geom_col]]
        if geom_col != "geometry":
            # buildingregulariser's regularize_geodataframe hard-codes the column name
            # "geometry" (its regularisation loop and its cleanup step).
            subset = subset.rename_geometry("geometry")
        out = geoai_geometry.regularize(
            subset,
            parallel_threshold=parallel_threshold_m,
            target_crs=None,
            simplify=True,
            simplify_tolerance=simplify_tolerance_m,
            allow_45_degree=allow_45_degree,
            allow_circles=allow_circles,
            circle_threshold=circle_threshold,
            num_cores=num_cores,
            include_metadata=False,
        )
        out = out[~out.geometry.isna() & ~out.geometry.is_empty]
        returned = set(out[_ROW_ID].astype(int))
        missing = [int(r) for r in subset[_ROW_ID] if int(r) not in returned]
        if missing:
            logger.warning(
                "Orthogonal regularization returned no geometry for %d polygons; their "
                "original geometry is kept",
                len(missing),
            )
        attrs = work.drop(columns=[geom_col]).set_index(_ROW_ID, drop=False)
        pieces_rows = list(out[_ROW_ID].astype(int))
        pieces_geoms = [make_polygonal(g) for g in out.geometry]
        untouched = [int(r) for r in work.loc[~present, _ROW_ID]] + missing
        pieces_rows += untouched
        pieces_geoms += [work.geometry.iloc[r] for r in untouched]
        order = np.argsort(np.asarray(pieces_rows), kind="stable")
        rows = [pieces_rows[k] for k in order]
        result = gpd.GeoDataFrame(
            attrs.loc[rows].reset_index(drop=True),
            geometry=gpd.GeoSeries([pieces_geoms[k] for k in order], crs=metric.crs),
            crs=metric.crs,
        )
        if result.geometry.name != geom_col:
            result = result.rename_geometry(geom_col)
        result = result.drop(columns=[_ROW_ID])
        result = result[list(gdf.columns)]
        detail = f"{len(out)} regularized parts, {len(missing)} kept unchanged"

    result = _to_original(result, original_crs).reset_index(drop=True)
    result.attrs = dict(gdf.attrs)
    logger.info(
        "Regularized polygons (method=%s): %d rows in, %d rows out (%s)",
        method,
        len(gdf),
        len(result),
        detail,
    )
    return result

simplify_polygons 

simplify_polygons(gdf: GeoDataFrame, tolerance: float = 2.0, preserve_topology: bool = True) -> gpd.GeoDataFrame

Simplify field boundary polygons with a tolerance in metres.

Unless the frame is already in a UTM zone it is projected to the UTM zone estimated from its extent (:func:to_metric_crs), simplified with :meth:geopandas.GeoSeries.simplify, and projected back, so the tolerance is in ground metres for any input CRS. A frame without a CRS is simplified in its native coordinate units (logged as a warning).

Geometries that become invalid are repaired (:func:make_polygonal); geometries that collapse to empty are removed (the number removed is logged). Attribute columns are preserved.

Applied after :func:smooth_polygons, as in the pipeline's default post-processing, the 2 m default removed area from 95-98 % of the polygons of four real engine outputs (median -0.2 % to -1.9 % of the smoothed area); see :func:agribound.pipeline._postprocess.

Parameters:

Name Type Description Default
gdf GeoDataFrame

Input polygons.

 required
tolerance float

Douglas-Peucker tolerance in metres (default 2.0). <= 0 returns gdf unchanged.

 2.0
preserve_topology bool

Use the topology-preserving variant (default True), which avoids collapsed or self-intersecting rings.

 True

Returns:

Type Description
GeoDataFrame

Simplified polygons in the input CRS (index reset).
Source code in agribound/postprocess/simplify.py 
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
def simplify_polygons(
    gdf: gpd.GeoDataFrame,
    tolerance: float = 2.0,
    preserve_topology: bool = True,
) -> gpd.GeoDataFrame:
    """Simplify field boundary polygons with a tolerance in metres.

    Unless the frame is already in a UTM zone it is projected to the UTM zone
    estimated from its extent (:func:`to_metric_crs`), simplified with
    :meth:`geopandas.GeoSeries.simplify`, and projected back, so the
    tolerance is in ground metres for any input CRS. A frame without a CRS
    is simplified in its native coordinate units (logged as a warning).

    Geometries that become invalid are repaired (:func:`make_polygonal`);
    geometries that collapse to empty are removed (the number removed is
    logged). Attribute columns are preserved.

    Applied after :func:`smooth_polygons`, as in the pipeline's default
    post-processing, the 2 m default removed area from 95-98 % of the
    polygons of four real engine outputs (median -0.2 % to -1.9 % of the
    smoothed area); see :func:`agribound.pipeline._postprocess`.

    Parameters
    ----------
    gdf : geopandas.GeoDataFrame
        Input polygons.
    tolerance : float
        Douglas-Peucker tolerance in metres (default 2.0). ``<= 0`` returns
        *gdf* unchanged.
    preserve_topology : bool
        Use the topology-preserving variant (default *True*), which avoids
        collapsed or self-intersecting rings.

    Returns
    -------
    geopandas.GeoDataFrame
        Simplified polygons in the input CRS (index reset).
    """
    if len(gdf) == 0 or tolerance <= 0:
        return gdf
    if not (~gdf.geometry.isna() & ~gdf.geometry.is_empty).any():
        # Nothing to simplify and no extent to pick a UTM zone from; as below,
        # empty geometries are removed and missing ones kept.
        return gdf[~gdf.geometry.is_empty].reset_index(drop=True)

    if gdf.crs is None:
        logger.warning(
            "simplify_polygons: input has no CRS; tolerance %.3g is applied in native "
            "coordinate units, not metres",
            tolerance,
        )
    result, original_crs = to_metric_crs(gdf)
    result = result.copy()
    geom_col = result.geometry.name
    result[geom_col] = result.geometry.simplify(tolerance, preserve_topology=preserve_topology)

    invalid = ~result.geometry.is_valid & ~result.geometry.isna()
    n_repaired = int(invalid.sum())
    if n_repaired:
        result.loc[invalid, geom_col] = result.loc[invalid, geom_col].map(make_polygonal)
    keep = ~result.geometry.is_empty
    n_removed = int((~keep).sum())
    result = result[keep]
    if n_repaired or n_removed:
        logger.info(
            "Simplified polygons: %d repaired after simplification, %d removed (collapsed)",
            n_repaired,
            n_removed,
        )

    result = _to_original(result, original_crs)
    return result.reset_index(drop=True)

Smoothing and simplification

simplify

Polygon smoothing and simplification.

  • :func:smooth_polygons applies Chaikin corner cutting to round off the pixel-staircase edges left by raster-to-vector conversion.
  • :func:simplify_polygons applies Douglas-Peucker simplification with a tolerance in metres: unless the data are already in a UTM zone, they are projected to the UTM zone estimated from their extent for the operation and projected back afterwards.

The module also holds two helpers shared by the other post-processing modules: :func:to_metric_crs and :func:make_polygonal.

to_metric_crs 

to_metric_crs(gdf: GeoDataFrame) -> tuple[gpd.GeoDataFrame, Any]

Project gdf to a UTM zone so that distances are in ground metres.

Frames already in a UTM zone (pyproj.CRS.utm_zone is set, e.g. EPSG:326xx/327xx or NAD83 UTM) are returned unchanged. Every other CRS -- geographic, Web Mercator, equal-area, State Plane in feet -- is projected to the UTM zone that :meth:geopandas.GeoDataFrame.estimate_utm_crs picks for the data's extent.

Parameters:

Name Type Description Default
gdf GeoDataFrame

Input frame with a CRS and at least one non-empty geometry (unless it is already in a UTM zone).

 required

Returns:

Type Description
tuple[GeoDataFrame, CRS or None]

(frame_in_utm, original_crs). When the input has no CRS it is returned unchanged together with None.
Source code in agribound/postprocess/simplify.py 
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
def to_metric_crs(gdf: gpd.GeoDataFrame) -> tuple[gpd.GeoDataFrame, Any]:
    """Project *gdf* to a UTM zone so that distances are in ground metres.

    Frames already in a UTM zone (``pyproj.CRS.utm_zone`` is set, e.g.
    EPSG:326xx/327xx or NAD83 UTM) are returned unchanged. Every other CRS --
    geographic, Web Mercator, equal-area, State Plane in feet -- is projected
    to the UTM zone that :meth:`geopandas.GeoDataFrame.estimate_utm_crs`
    picks for the data's extent.

    Parameters
    ----------
    gdf : geopandas.GeoDataFrame
        Input frame with a CRS and at least one non-empty geometry (unless it
        is already in a UTM zone).

    Returns
    -------
    tuple[geopandas.GeoDataFrame, pyproj.CRS or None]
        ``(frame_in_utm, original_crs)``. When the input has no CRS it is
        returned unchanged together with *None*.
    """
    original_crs = gdf.crs
    if original_crs is None:
        return gdf, None
    if original_crs.utm_zone is not None:
        return gdf, original_crs
    return gdf.to_crs(gdf.estimate_utm_crs()), original_crs

make_polygonal 

make_polygonal(geom: Any) -> Polygon | MultiPolygon | None

Return a valid polygonal version of geom.

Invalid geometries are repaired with :func:shapely.validation.make_valid and only their polygonal parts are kept (lines and points produced by the repair are discarded). Returns None for missing input and an empty Polygon when nothing polygonal remains.

Parameters:

Name Type Description Default
geom shapely geometry or None

Input geometry.

 required

Returns:

Type Description
Polygon or MultiPolygon or None
Source code in agribound/postprocess/simplify.py 
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
def make_polygonal(geom: Any) -> Polygon | MultiPolygon | None:
    """Return a valid polygonal version of *geom*.

    Invalid geometries are repaired with :func:`shapely.validation.make_valid`
    and only their polygonal parts are kept (lines and points produced by the
    repair are discarded). Returns *None* for missing input and an empty
    ``Polygon`` when nothing polygonal remains.

    Parameters
    ----------
    geom : shapely geometry or None
        Input geometry.

    Returns
    -------
    shapely.geometry.Polygon or shapely.geometry.MultiPolygon or None
    """
    if geom is None:
        return None
    if geom.is_empty:
        return Polygon()
    fixed = geom if geom.is_valid else make_valid(geom)
    if isinstance(fixed, Polygon | MultiPolygon):
        return fixed
    if isinstance(fixed, GeometryCollection):
        parts: list[Polygon] = []
        for part in fixed.geoms:
            if isinstance(part, Polygon) and not part.is_empty:
                parts.append(part)
            elif isinstance(part, MultiPolygon):
                parts.extend(p for p in part.geoms if not p.is_empty)
        if not parts:
            return Polygon()
        return parts[0] if len(parts) == 1 else MultiPolygon(parts)
    return Polygon()

smooth_polygons 

smooth_polygons(gdf: GeoDataFrame, iterations: int = 3, max_segment_length: float | None = None) -> gpd.GeoDataFrame

Smooth pixel-staircase artefacts with Chaikin's corner-cutting algorithm.

Chaikin's scheme is affine-invariant, so it is applied directly to the vertex coordinates in the frame's CRS (no reprojection). Each iteration doubles the vertex count and cuts every corner at 1/4 and 3/4 of its adjacent segments, so the depth of a cut scales with the segment length and convex outlines shrink. :func:rasterio.features.shapes traces an axis-aligned rectangular pixel region with only its four corners, and three iterations then remove about 16 % of its area (at any scale); staircase outlines of oblique edges have a vertex at every pixel step and lose much less (about 0.02 % for a 60 x 40 pixel rectangle rotated by 30 degrees, about 1 % for a 10 x 6 pixel one). Passing max_segment_length (e.g. the pixel size) first inserts vertices so no segment is longer than that, which confines the rounding to the corners (a 20 x 10 pixel rectangle densified at one pixel loses about 0.1 %).

The corner cut grows with the segment length, so long straight edges traced with few vertices are rounded over tens of metres: a 400 x 300 m rectangle traced with four corners loses 16.4 % of its area and moves by up to 60 m (Hausdorff distance). On real engine outputs, three iterations changed the polygon area by a median of -0.35 % to -1.6 % (10th percentile -3.3 % to -7.7 %; Hausdorff distance median 6.6-9.3 m, 90th percentile 17-25 m) for Delineate-Anything on 0.9 m NAIP and 10 m Sentinel-2 and FTW on 10 m Sentinel-2; see :func:agribound.pipeline._postprocess for the full measurement.

Holes are smoothed the same way. If a smoothed polygon becomes invalid it is repaired and only its largest part is kept; if nothing remains the original polygon is kept. Attribute columns are preserved.

Parameters:

Name Type Description Default
gdf GeoDataFrame

Input polygons (typically raw output from raster vectorization).

 required
iterations int

Number of Chaikin iterations (default 3). 0 returns gdf unchanged.

 3
max_segment_length float or None

Densify rings with :func:shapely.segmentize to this maximum segment length (CRS units) before smoothing. None (default) smooths the vertices as given.

 None

Returns:

Type Description
GeoDataFrame

Smoothed polygons (empty geometries removed, index reset).
Source code in agribound/postprocess/simplify.py 
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
def smooth_polygons(
    gdf: gpd.GeoDataFrame,
    iterations: int = 3,
    max_segment_length: float | None = None,
) -> gpd.GeoDataFrame:
    """Smooth pixel-staircase artefacts with Chaikin's corner-cutting algorithm.

    Chaikin's scheme is affine-invariant, so it is applied directly to the
    vertex coordinates in the frame's CRS (no reprojection). Each iteration
    doubles the vertex count and cuts every corner at 1/4 and 3/4 of its
    adjacent segments, so the depth of a cut scales with the segment length
    and convex outlines shrink. :func:`rasterio.features.shapes` traces an
    axis-aligned rectangular pixel region with only its four corners, and
    three iterations then remove about 16 % of its area (at any scale);
    staircase outlines of oblique edges have a vertex at every pixel step and
    lose much less (about 0.02 % for a 60 x 40 pixel rectangle rotated by
    30 degrees, about 1 % for a 10 x 6 pixel one). Passing
    *max_segment_length* (e.g. the pixel size) first inserts vertices so no
    segment is longer than that, which confines the rounding to the corners
    (a 20 x 10 pixel rectangle densified at one pixel loses about 0.1 %).

    The corner cut grows with the segment length, so long straight edges
    traced with few vertices are rounded over tens of metres: a 400 x 300 m
    rectangle traced with four corners loses 16.4 % of its area and moves by
    up to 60 m (Hausdorff distance). On real engine outputs, three
    iterations changed the polygon area by a median of -0.35 % to -1.6 %
    (10th percentile -3.3 % to -7.7 %; Hausdorff distance median 6.6-9.3 m,
    90th percentile 17-25 m) for Delineate-Anything on 0.9 m NAIP and 10 m
    Sentinel-2 and FTW on 10 m Sentinel-2; see
    :func:`agribound.pipeline._postprocess` for the full measurement.

    Holes are smoothed the same way. If a smoothed polygon becomes invalid it
    is repaired and only its largest part is kept; if nothing remains the
    original polygon is kept. Attribute columns are preserved.

    Parameters
    ----------
    gdf : geopandas.GeoDataFrame
        Input polygons (typically raw output from raster vectorization).
    iterations : int
        Number of Chaikin iterations (default 3). ``0`` returns *gdf*
        unchanged.
    max_segment_length : float or None
        Densify rings with :func:`shapely.segmentize` to this maximum segment
        length (CRS units) before smoothing. *None* (default) smooths the
        vertices as given.

    Returns
    -------
    geopandas.GeoDataFrame
        Smoothed polygons (empty geometries removed, index reset).
    """
    if len(gdf) == 0 or iterations <= 0:
        return gdf
    if max_segment_length is not None and max_segment_length <= 0:
        raise ValueError(f"max_segment_length must be > 0, got {max_segment_length}")

    result = gdf.copy()

    def _smooth(geom):
        if geom is None or geom.is_empty:
            return geom
        if max_segment_length is not None:
            import shapely

            geom = shapely.segmentize(geom, max_segment_length)
        if isinstance(geom, MultiPolygon):
            return MultiPolygon([_smooth_polygon(p, iterations) for p in geom.geoms])
        if isinstance(geom, Polygon):
            return _smooth_polygon(geom, iterations)
        return geom

    result[result.geometry.name] = result.geometry.map(_smooth)
    result = result[~result.geometry.is_empty]

    logger.info("Smoothed %d polygons (Chaikin x%d)", len(result), iterations)
    return result.reset_index(drop=True)

simplify_polygons 

simplify_polygons(gdf: GeoDataFrame, tolerance: float = 2.0, preserve_topology: bool = True) -> gpd.GeoDataFrame

Simplify field boundary polygons with a tolerance in metres.

Unless the frame is already in a UTM zone it is projected to the UTM zone estimated from its extent (:func:to_metric_crs), simplified with :meth:geopandas.GeoSeries.simplify, and projected back, so the tolerance is in ground metres for any input CRS. A frame without a CRS is simplified in its native coordinate units (logged as a warning).

Geometries that become invalid are repaired (:func:make_polygonal); geometries that collapse to empty are removed (the number removed is logged). Attribute columns are preserved.

Applied after :func:smooth_polygons, as in the pipeline's default post-processing, the 2 m default removed area from 95-98 % of the polygons of four real engine outputs (median -0.2 % to -1.9 % of the smoothed area); see :func:agribound.pipeline._postprocess.

Parameters:

Name Type Description Default
gdf GeoDataFrame

Input polygons.

 required
tolerance float

Douglas-Peucker tolerance in metres (default 2.0). <= 0 returns gdf unchanged.

 2.0
preserve_topology bool

Use the topology-preserving variant (default True), which avoids collapsed or self-intersecting rings.

 True

Returns:

Type Description
GeoDataFrame

Simplified polygons in the input CRS (index reset).
Source code in agribound/postprocess/simplify.py 
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
def simplify_polygons(
    gdf: gpd.GeoDataFrame,
    tolerance: float = 2.0,
    preserve_topology: bool = True,
) -> gpd.GeoDataFrame:
    """Simplify field boundary polygons with a tolerance in metres.

    Unless the frame is already in a UTM zone it is projected to the UTM zone
    estimated from its extent (:func:`to_metric_crs`), simplified with
    :meth:`geopandas.GeoSeries.simplify`, and projected back, so the
    tolerance is in ground metres for any input CRS. A frame without a CRS
    is simplified in its native coordinate units (logged as a warning).

    Geometries that become invalid are repaired (:func:`make_polygonal`);
    geometries that collapse to empty are removed (the number removed is
    logged). Attribute columns are preserved.

    Applied after :func:`smooth_polygons`, as in the pipeline's default
    post-processing, the 2 m default removed area from 95-98 % of the
    polygons of four real engine outputs (median -0.2 % to -1.9 % of the
    smoothed area); see :func:`agribound.pipeline._postprocess`.

    Parameters
    ----------
    gdf : geopandas.GeoDataFrame
        Input polygons.
    tolerance : float
        Douglas-Peucker tolerance in metres (default 2.0). ``<= 0`` returns
        *gdf* unchanged.
    preserve_topology : bool
        Use the topology-preserving variant (default *True*), which avoids
        collapsed or self-intersecting rings.

    Returns
    -------
    geopandas.GeoDataFrame
        Simplified polygons in the input CRS (index reset).
    """
    if len(gdf) == 0 or tolerance <= 0:
        return gdf
    if not (~gdf.geometry.isna() & ~gdf.geometry.is_empty).any():
        # Nothing to simplify and no extent to pick a UTM zone from; as below,
        # empty geometries are removed and missing ones kept.
        return gdf[~gdf.geometry.is_empty].reset_index(drop=True)

    if gdf.crs is None:
        logger.warning(
            "simplify_polygons: input has no CRS; tolerance %.3g is applied in native "
            "coordinate units, not metres",
            tolerance,
        )
    result, original_crs = to_metric_crs(gdf)
    result = result.copy()
    geom_col = result.geometry.name
    result[geom_col] = result.geometry.simplify(tolerance, preserve_topology=preserve_topology)

    invalid = ~result.geometry.is_valid & ~result.geometry.isna()
    n_repaired = int(invalid.sum())
    if n_repaired:
        result.loc[invalid, geom_col] = result.loc[invalid, geom_col].map(make_polygonal)
    keep = ~result.geometry.is_empty
    n_removed = int((~keep).sum())
    result = result[keep]
    if n_repaired or n_removed:
        logger.info(
            "Simplified polygons: %d repaired after simplification, %d removed (collapsed)",
            n_repaired,
            n_removed,
        )

    result = _to_original(result, original_crs)
    return result.reset_index(drop=True)

LULC crop filter

lulc_filter

LULC crop filter.

Removes non-agricultural polygons using a land-use/land-cover (LULC) dataset read from Google Earth Engine. For every polygon the filter computes a crop value in [0, 1] and drops polygons below config.lulc_crop_threshold:

======================== ========================================== ========== ========= lulc_dataset Earth Engine asset / crop rule Years Scale ======================== ========================================== ========== ========= "nlcd" projects/sat-io/open-datasets/USGS/ 1985-2025 30 m ANNUAL_NLCD/LANDCOVER band b1, classes 81 (pasture/hay) and 82 (cultivated crops); fraction of pixels "cdl" USDA/NASS/CDL band cultivated == 2; 2013-2023 30 m fraction of pixels (CONUS) "dynamic_world" GOOGLE/DYNAMICWORLD/V1 annual median of 2016-last 10 m crops; mean probability (not a pixel full year fraction) "c3s" projects/sat-io/open-datasets/ESA/C3S-LC- 2000-2022 300 m L4-LCCS band b1, classes 10, 11, 12, 20, 30; fraction of pixels ======================== ========================================== ========== =========

Tree crops (config.lulc_tree_crops=True): Dynamic World files plantations and orchards under trees, so their crops probability is low and the default rule removes them. With the option, the Dynamic World value is the annual median of the per-image sum of the crops and trees probabilities, and C3S also counts its tree-cover classes (:data:C3S_TREE_CLASSES: 50, 60, 61, 62, 70, 71, 72, 80, 81, 82, 90). With these two datasets the filter then keeps forest as well as tree crops; it still removes water, built-up, bare, grass and shrub land. NLCD and CDL are unchanged, so with them it still removes forest: NLCD class 82 (cultivated crops) includes "perennial woody crops such as orchards and vineyards", and CDL counts orchards as cultivated.

Year ranges are those of the assets on 2026-09-26. When the requested year is outside a dataset's range, the nearest available year is used, recorded in lulc:year and logged as a WARNING. Annual NLCD and C3S come from the community catalogue (projects/sat-io), which has no Google service-level agreement; if the Annual NLCD asset cannot be read, the official USGS/NLCD_RELEASES/2021_REL/NLCD (single year 2021, band landcover) is used instead, with a WARNING. The CDL cultivated band is present in Earth Engine only for 2013-2023.

Dataset routing (lulc_dataset="auto", :func:select_lulc_dataset), for the study area (or, without a study area, the local raster's footprint):

  1. If the area intersects the conterminous-US envelope (-125, 24, -66, 50) and at least 90 % of it has valid Annual NLCD pixels (checked on Earth Engine), use NLCD. Northern Mexico and southern Canada fall inside the envelope but have no NLCD pixels, so they are not routed to NLCD.
  2. Otherwise, for years with a complete Dynamic World year (2016 up to the previous calendar year), use Dynamic World; later years use the last complete year.
  3. Otherwise (before 2016) use C3S (nearest year in 2000-2022).

The GEE catalogue notes that Dynamic World crop probabilities can be comparatively low in the absence of obvious distinguishing features and on high-return surfaces in arid climates, so the default threshold may remove real fields in arid regions.

Output columns: lulc:crop_fraction (float; NaN when the dataset has no valid pixel under the polygon, never 0), lulc:dataset, lulc:year (the dataset year used) and lulc:valid (bool). Polygons with NaN are kept (lulc_nodata_policy="keep") or dropped ("drop").

Modes (config.lulc_mode):

  • "server": per-polygon means with ee.Image.reduceRegions (Earth Engine's pixel-area-weighted mean at the dataset scale), in batches of config.lulc_batch_size polygons, retried on rate limits.
  • "raster": :func:prefetch_lulc_raster downloads a float32 crop raster (0/1 mask, or probability for Dynamic World; NaN nodata; 30 m, 10 m for Dynamic World) during the composite stage; the filter then averages the pixels whose centres fall inside each polygon locally (:func:zonal_mean_from_raster), so the delineation stage can run offline.

filter_by_lulc 

filter_by_lulc(gdf: GeoDataFrame, config: AgriboundConfig) -> gpd.GeoDataFrame

Drop polygons whose crop value is below config.lulc_crop_threshold.

See the module docstring for datasets, routing, modes and columns.

Parameters:

Name Type Description Default
gdf GeoDataFrame

Field polygons (any CRS; the output keeps it).

 required
config AgriboundConfig

Pipeline configuration (lulc_* fields, year, study_area).

 required

Returns:

Type Description
GeoDataFrame

Kept polygons with lulc:crop_fraction, lulc:dataset, lulc:year and lulc:valid columns (index reset). attrs keeps the input attributes and adds lulc_stats = {dataset, year_requested, year_used, n_in, n_kept, n_nan, threshold, mode, nodata_policy, tree_crops, n_below_threshold, n_nan_dropped, asset, band, value, selection_reason, nlcd_valid_fraction}.

Raises:

Type Description
RuntimeError

If the LULC values cannot be computed (Earth Engine or download failure, missing cached raster offline, ...).
Source code in agribound/postprocess/lulc_filter.py 
 947
 948
 949
 950
 951
 952
 953
 954
 955
 956
 957
 958
 959
 960
 961
 962
 963
 964
 965
 966
 967
 968
 969
 970
 971
 972
 973
 974
 975
 976
 977
 978
 979
 980
 981
 982
 983
 984
 985
 986
 987
 988
 989
 990
 991
 992
 993
 994
 995
 996
 997
 998
 999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
def filter_by_lulc(gdf: gpd.GeoDataFrame, config: AgriboundConfig) -> gpd.GeoDataFrame:
    """Drop polygons whose crop value is below ``config.lulc_crop_threshold``.

    See the module docstring for datasets, routing, modes and columns.

    Parameters
    ----------
    gdf : geopandas.GeoDataFrame
        Field polygons (any CRS; the output keeps it).
    config : AgriboundConfig
        Pipeline configuration (``lulc_*`` fields, ``year``, ``study_area``).

    Returns
    -------
    geopandas.GeoDataFrame
        Kept polygons with ``lulc:crop_fraction``, ``lulc:dataset``,
        ``lulc:year`` and ``lulc:valid`` columns (index reset). ``attrs`` keeps
        the input attributes and adds ``lulc_stats`` = {``dataset``,
        ``year_requested``, ``year_used``, ``n_in``, ``n_kept``, ``n_nan``,
        ``threshold``, ``mode``, ``nodata_policy``, ``tree_crops``,
        ``n_below_threshold``, ``n_nan_dropped``, ``asset``, ``band``,
        ``value``, ``selection_reason``, ``nlcd_valid_fraction``}.

    Raises
    ------
    RuntimeError
        If the LULC values cannot be computed (Earth Engine or download
        failure, missing cached raster offline, ...).
    """
    if len(gdf) == 0:
        return _empty_result(gdf, config)
    try:
        values, selection, meta = _compute_values(gdf, config)
    except RuntimeError:
        raise
    except Exception as exc:
        raise RuntimeError(f"LULC filter failed: {type(exc).__name__}: {exc}") from exc

    threshold = float(config.lulc_crop_threshold)
    valid = np.isfinite(values)
    above = valid & (values >= threshold)
    keep_nan = config.lulc_nodata_policy == "keep"
    keep = above | (~valid & keep_nan)

    result = gdf.copy()
    result["lulc:crop_fraction"] = values
    result["lulc:dataset"] = selection.dataset
    result["lulc:year"] = int(meta["year_used"])
    result["lulc:valid"] = valid
    out = result.loc[keep].reset_index(drop=True)
    n_nan = int((~valid).sum())
    stats = {
        "dataset": selection.dataset,
        "year_requested": int(selection.year_requested),
        "year_used": int(meta["year_used"]),
        "n_in": int(len(gdf)),
        "n_kept": int(len(out)),
        "n_nan": n_nan,
        "threshold": threshold,
        "mode": config.lulc_mode,
        "nodata_policy": config.lulc_nodata_policy,
        "tree_crops": bool(config.lulc_tree_crops),
        "n_below_threshold": int((valid & ~above).sum()),
        "n_nan_dropped": 0 if keep_nan else n_nan,
        "asset": meta["asset"],
        "band": meta["band"],
        "value": meta["value"],
        "selection_reason": selection.reason,
        "nlcd_valid_fraction": selection.nlcd_valid_fraction,
    }
    out.attrs = dict(gdf.attrs)
    out.attrs["lulc_stats"] = stats
    if int(meta["year_used"]) != int(selection.year_requested):
        logger.warning(
            "LULC filter used %s %d for requested year %d",
            selection.dataset,
            meta["year_used"],
            selection.year_requested,
        )
    if n_nan:
        logger.warning(
            "LULC filter: %d of %d polygons have no valid %s pixels (%s by lulc_nodata_policy=%r)",
            n_nan,
            len(gdf),
            selection.dataset,
            "kept" if keep_nan else "dropped",
            config.lulc_nodata_policy,
        )
    logger.info(
        "LULC filter (%s %d, %s mode): %d -> %d polygons (threshold=%.2f)",
        selection.dataset,
        meta["year_used"],
        config.lulc_mode,
        len(gdf),
        len(out),
        threshold,
    )
    return out

prefetch_lulc_raster 

prefetch_lulc_raster(config: AgriboundConfig) -> str | None

Download the LULC crop raster for the study area (lulc_mode="raster").

Writes a single-band float32 GeoTIFF (band crop: 1/0 crop mask, or the Dynamic World crop probability; NaN where the dataset has no data) covering the study-area bounding box plus a 3-pixel margin, at 30 m (10 m for Dynamic World), in config.export_crs, to cache_path(config, "lulc_<dataset>", ".tif", ...). An existing file is reused.

Parameters:

Name Type Description Default
config AgriboundConfig

Pipeline configuration.

 required

Returns:

Type Description
str or None

Path to the raster, or None when config.lulc_filter is False.

Raises:

Type Description
RuntimeError

If the dataset cannot be selected or downloaded.
Source code in agribound/postprocess/lulc_filter.py 
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
def prefetch_lulc_raster(config: AgriboundConfig) -> str | None:
    """Download the LULC crop raster for the study area (``lulc_mode="raster"``).

    Writes a single-band float32 GeoTIFF (band ``crop``: 1/0 crop mask, or the
    Dynamic World crop probability; NaN where the dataset has no data) covering
    the study-area bounding box plus a 3-pixel margin, at 30 m (10 m for
    Dynamic World), in ``config.export_crs``, to
    ``cache_path(config, "lulc_<dataset>", ".tif", ...)``. An existing file is
    reused.

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

    Returns
    -------
    str or None
        Path to the raster, or *None* when ``config.lulc_filter`` is *False*.

    Raises
    ------
    RuntimeError
        If the dataset cannot be selected or downloaded.
    """
    if not config.lulc_filter:
        return None
    try:
        return _prefetch(config)
    except RuntimeError:
        raise
    except Exception as exc:
        raise RuntimeError(f"LULC raster prefetch failed: {type(exc).__name__}: {exc}") from exc

select_lulc_dataset 

select_lulc_dataset(config: AgriboundConfig, aoi_gdf: GeoDataFrame | None = None) -> tuple[str, int]

Select the LULC dataset and year for config.

Applies config.lulc_dataset directly, or for "auto" the routing rules in the module docstring to the study area (else the local raster's footprint, else the envelope of aoi_gdf). The NLCD coverage test is an Earth Engine request; it is skipped for areas that do not intersect the conterminous-US envelope. The decision is cached next to the other intermediates (lulc_selection_<key>.json), so later calls with the same configuration (e.g. on an offline node) do not query Earth Engine.

Parameters:

Name Type Description Default
config AgriboundConfig

Pipeline configuration (year, lulc_dataset, study_area).

 required
aoi_gdf GeoDataFrame or None

Polygons used as the area when there is neither a study area nor a local raster.

 None

Returns:

Type Description
tuple[str, int]

(dataset, year_used). year_used is the nearest year the dataset has (for NLCD, the official 2021 release is substituted when the Annual NLCD asset is unreadable; the year actually used is then reported by :func:filter_by_lulc).
Source code in agribound/postprocess/lulc_filter.py 
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
def select_lulc_dataset(
    config: AgriboundConfig, aoi_gdf: gpd.GeoDataFrame | None = None
) -> tuple[str, int]:
    """Select the LULC dataset and year for *config*.

    Applies ``config.lulc_dataset`` directly, or for ``"auto"`` the routing
    rules in the module docstring to the study area (else the local raster's
    footprint, else the envelope of *aoi_gdf*). The NLCD coverage test is an
    Earth Engine request; it is skipped for areas that do not intersect the
    conterminous-US envelope. The decision is cached next to the other
    intermediates (``lulc_selection_<key>.json``), so later calls with the
    same configuration (e.g. on an offline node) do not query Earth Engine.

    Parameters
    ----------
    config : AgriboundConfig
        Pipeline configuration (``year``, ``lulc_dataset``, ``study_area``).
    aoi_gdf : geopandas.GeoDataFrame or None
        Polygons used as the area when there is neither a study area nor a
        local raster.

    Returns
    -------
    tuple[str, int]
        ``(dataset, year_used)``. ``year_used`` is the nearest year the
        dataset has (for NLCD, the official 2021 release is substituted when
        the Annual NLCD asset is unreadable; the year actually used is then
        reported by :func:`filter_by_lulc`).
    """
    selection = _select(config, aoi_gdf)
    return selection.dataset, selection.year_used