Field Boundary Extraction with GeoPandas

Field boundary extraction converts raw drone classification masks, GPS traces, or satellite-derived crop maps into clean, topology-valid polygon layers that drive every downstream precision agriculture operation — producing a GeoPackage ready for variable rate application (VRA) planning, yield map zoning, and regulatory acreage reporting. This guide delivers a tested, production-ready Python workflow for agtech engineers, farm data analysts, and Python GIS developers who need to automate boundary generation at scale without manual digitizing.

This page is part of Ag-GIS Data Fundamentals & Spatial Reference Systems — see there for the full spatial data pipeline context.

Prerequisites

Python packages (pinned to tested versions):

pip install "geopandas>=1.0" "shapely>=2.0" "rasterio>=1.3" "numpy>=1.24" "scipy>=1.11" "scikit-image>=0.22"

Input data requirements:

• A binary or labeled segmentation raster (GeoTIFF, single band) — e.g., crop/non-crop classification from a DJI P4 Multispectral or MicaSense RedEdge-MX flight, or a Sentinel-2 derived mask
• CRS embedded in the raster (read with rasterio.open().crs); a missing CRS must be assigned before any metric operation
• Nodata convention: background pixels as 0, field pixels as 1 (or class integer values for multi-crop masks)
• RTK GPS accuracy of ±2 cm or better if GPS traces are the boundary source; standard GNSS (±3 m) introduces unacceptable boundary drift for sub-field applications

Correct CRS handling is the most critical prerequisite. Metric projections (UTM zone, EPSG:326xx/327xx) are required for accurate area and distance calculations. For foundational guidance on projection selection and datum shifts, read Understanding CRS in Precision Agriculture before scaling this pipeline across multi-farm datasets.

Extraction pipeline overview

1. Concept and algorithm

Field boundary extraction is a raster-to-vector workflow. A segmentation raster encodes crop-classified pixels as integer values; the goal is to trace the edges between foreground (crop) and background (non-crop) regions into closed polygon geometries that accurately represent agronomic parcels.

Why boundary accuracy matters agronomically: Variable rate application controllers (ISOXML-compatible terminals on John Deere, CNH, and AGCO machinery) use field boundary polygons as the spatial envelope within which prescription maps are applied. A boundary displaced by 3 m means headlands receive the wrong rate; a boundary with self-intersecting rings causes some controllers to reject the prescription entirely. For yield mapping, the boundary defines the spatial filter applied to combine GPS positions — missed rows near boundaries inflate or deflate per-hectare yield estimates.

The core algorithm is marching squares contour tracing, implemented in rasterio.features.shapes. It iterates over the raster grid and traces the boundary of each contiguous region sharing the same pixel value. The result is a generator of (geometry, value) tuples in GeoJSON-compatible dict form, each geometry being a Polygon or MultiPolygon. The algorithm runs in O(n) on the number of pixels, making it practical for 10 cm resolution drone tiles covering hundreds of hectares.

Why morphological pre-processing is mandatory: Classification masks from multispectral indices contain mixed-pixel noise at field edges — individual pixels that flip between crop and non-crop due to shadow, BRDF effects, or adjacency to irrigation canals. Without morphological closing and opening, the tracer produces thousands of micro-polygons (1-5 pixel fragments) that bloat the vector file and corrupt downstream spatial joins. A 3×3 structuring element removes artifacts smaller than 9 pixels (roughly 0.09 m² at 10 cm resolution) before a single polygon is created.

2. Step-by-step implementation

Step 1: Ingest the segmentation raster and clean the mask

import rasterio
import numpy as np
from scipy.ndimage import binary_closing, binary_opening
from pathlib import Path


def load_and_clean_mask(
    raster_path: str,
    target_epsg: int = 32618,
    kernel_size: int = 3,
) -> tuple[np.ndarray, object, object]:
    """
    Load a single-band segmentation GeoTIFF, assert CRS, and apply
    morphological cleaning.  Returns (mask, affine_transform, crs).
    """
    with rasterio.open(raster_path) as src:
        mask = src.read(1).astype(np.uint8)
        transform = src.transform
        crs = src.crs

        # Fail fast if CRS is missing or geographic (degrees not metres)
        assert crs is not None, f"Raster has no embedded CRS: {raster_path}"
        assert crs.is_projected, (
            f"CRS {crs.to_epsg()} is geographic — reproject to EPSG:{target_epsg} first. "
            "See: /ag-gis-data-fundamentals-spatial-reference-systems/understanding-crs-in-precision-agriculture/"
        )
        assert crs.to_epsg() == target_epsg, (
            f"Expected EPSG:{target_epsg}, got EPSG:{crs.to_epsg()}. "
            "Transform input before extraction."
        )

    struct = np.ones((kernel_size, kernel_size), dtype=bool)
    # Close small gaps (fills holes inside fields)
    mask = binary_closing(mask, structure=struct).astype(np.uint8)
    # Open isolated noise pixels (removes salt-and-pepper speckling)
    mask = binary_opening(mask, structure=struct).astype(np.uint8)

    return mask, transform, crs

Validate before proceeding:

mask, transform, crs = load_and_clean_mask(
    "farm_a_crop_mask.tif", target_epsg=32618
)
assert mask.ndim == 2, "Expected a single-band 2D array"
assert set(np.unique(mask)).issubset({0, 1}), "Mask should be binary (0/1)"
print(f"Mask shape: {mask.shape}, foreground pixels: {mask.sum()}")

Step 2: Convert raster regions to polygon geometries

import geopandas as gpd
from rasterio.features import shapes
from shapely.geometry import shape


def raster_to_polygons(
    mask: np.ndarray,
    transform,
    crs,
    min_area_ha: float = 0.5,
) -> gpd.GeoDataFrame:
    """
    Trace polygon boundaries for all foreground pixel regions.
    Filters out fragments smaller than min_area_ha hectares.
    """
    geom_records = [
        {"geometry": shape(geom), "class_val": int(val)}
        for geom, val in shapes(mask, mask=(mask > 0), transform=transform)
    ]
    assert len(geom_records) > 0, "No foreground polygons found — check mask values"

    gdf = gpd.GeoDataFrame(geom_records, crs=f"EPSG:{crs.to_epsg()}")

    # Area in hectares requires a projected CRS (metres)
    gdf["area_ha"] = gdf.geometry.area / 10_000.0
    before = len(gdf)
    gdf = gdf[gdf["area_ha"] >= min_area_ha].copy().reset_index(drop=True)
    print(f"Kept {len(gdf)} / {before} polygons >= {min_area_ha} ha")
    return gdf

The mask=(mask > 0) argument to shapes() is critical: it restricts tracing to foreground pixels only, preventing the algorithm from tracing the border of the entire raster extent as a giant polygon.

Validate output:

gdf = raster_to_polygons(mask, transform, crs, min_area_ha=0.5)
assert (gdf.area / 10_000 >= 0.5).all(), "Area filter not applied correctly"
assert gdf.crs.to_epsg() == 32618, "CRS lost during GeoDataFrame construction"
print(gdf[["class_val", "area_ha"]].describe())

Step 3: Geometric cleaning and topology repair

Raw vectorized boundaries carry two artefact types: jagged stair-step edges from the pixel grid, and self-intersecting rings from classification noise at corners. Both must be resolved before export.

def clean_boundaries(
    gdf: gpd.GeoDataFrame,
    simplify_tolerance: float = 2.0,
    sliver_threshold_m2: float = 100.0,
) -> gpd.GeoDataFrame:
    """
    Repair and simplify field boundary geometries.

    simplify_tolerance: Douglas-Peucker tolerance in CRS units (metres for UTM).
    sliver_threshold_m2: polygons smaller than this (m²) are dropped as artefacts.
    """
    # 1. Fix self-intersections from classification noise
    gdf["geometry"] = gdf.geometry.apply(
        lambda g: g.buffer(0) if not g.is_valid else g
    )

    # 2. Drop micro-slivers that survived the area filter
    gdf = gdf[gdf.geometry.area > sliver_threshold_m2].copy()

    # 3. Smooth stair-step pixelation — preserve_topology prevents fragmentation
    gdf["geometry"] = gdf.geometry.simplify(
        tolerance=simplify_tolerance, preserve_topology=True
    )

    # 4. Final topology enforcement via GEOS make_valid
    gdf["geometry"] = gdf.geometry.make_valid()

    # 5. Assert all geometries are valid and non-empty before returning
    assert gdf.geometry.is_valid.all(), "make_valid() did not resolve all issues"
    assert not gdf.geometry.is_empty.any(), "Empty geometries after cleaning"

    return gdf.reset_index(drop=True)

Validate:

gdf_clean = clean_boundaries(gdf, simplify_tolerance=2.0)
assert gdf_clean.geometry.is_valid.all()
assert (gdf_clean.geometry.geom_type.isin(["Polygon", "MultiPolygon"])).all()

Step 4: Dissolve, validate, and export to GeoPackage

import pandas as pd


def validate_and_export(
    gdf: gpd.GeoDataFrame,
    output_path: str,
    dissolve_by: str = "class_val",
) -> None:
    """
    Dissolve adjacent same-class polygons, add metadata, and export to GeoPackage.
    """
    # Dissolve adjacent polygons sharing the same field class
    gdf = gdf.dissolve(by=dissolve_by, aggfunc="first").reset_index()

    # Recompute final area after dissolve (dissolve may merge geometries)
    gdf["final_area_ha"] = gdf.geometry.area / 10_000.0
    gdf["source"] = "geopandas_extraction"
    gdf["extracted_at"] = pd.Timestamp.now(tz="UTC").isoformat()

    # Assert reasonable area bounds (flag outliers > 500 ha as suspect)
    assert (gdf["final_area_ha"] < 500).all(), (
        "Boundary larger than 500 ha detected — check for mask inversion"
    )

    # GeoPackage: single-file, CRS-embedded, no field name truncation
    gdf.to_file(output_path, driver="GPKG", layer="field_boundaries")
    print(f"Exported {len(gdf)} field boundaries to {output_path}")
    print(gdf[["final_area_ha"]].describe())

GeoPackage (.gpkg) is the recommended format over Shapefile for field boundaries: it stores CRS metadata natively, supports column names longer than 10 characters (critical for agronomic attribute schemas), has no 2 GB size limit, and ships as a single file with no sidecar .prj or .dbf that can be accidentally lost. Before deploying the output to a machinery controller, visually verify alignment in QGIS against the source orthomosaic.

3. Key parameters and tuning

4. Handling edge cases and failure modes

Stair-step polygon edges: Occurs when the input mask resolution is coarser than 0.5 m/px (common with Sentinel-2 10 m imagery). Increase simplify_tolerance to 5–10 m or apply Gaussian smoothing to the mask array (scipy.ndimage.gaussian_filter) before morphological cleaning. For ingesting multispectral drone imagery at 5–10 cm resolution, a tolerance of 1.5–2.0 m is usually sufficient.

UTM zone boundary crossings: Farms straddling two UTM zones (e.g., a 40,000 ha operation crossing 6° longitude) produce coordinates in two different EPSG codes. Attempting to sjoin or overlay GeoDataFrames in mismatched UTM zones silently produces geometry errors. Always call gdf.to_crs(epsg=target) to force a single projection before any spatial operation. The how to convert WGS84 to UTM for farm mapping guide covers the full reprojection workflow.

Zero-area polygons after dissolve: If dissolve() merges non-adjacent polygons by the same class_val key, it can produce a MultiPolygon where some components are degenerate. Always call .explode(index_parts=False) after dissolve() if the downstream system (e.g., John Deere Operations Center) requires simple Polygons rather than MultiPolygons.

Mask inversion: A common error is processing an inverted mask where background is 1 and crop is 0. The shapes() generator traces every connected region including the entire raster border, producing one enormous polygon. Detect this with: assert gdf["final_area_ha"].max() < 500, "Possible mask inversion".

Shapefile encoding for legacy FMIS: Some farm management information systems (FMIS) still ingest Shapefiles with Latin-1 encoding for field names. GeoPandas defaults to UTF-8 when writing .shp. Pass encoding="latin-1" to to_file() for those integrations, and truncate attribute names to ≤10 characters manually before export.

5. Verification and output validation

Visual spot-check: Open the exported .gpkg in QGIS alongside the original orthomosaic or imagery. Enable transparency on the boundary layer and verify that edges align with visible crop rows, not offset by one or more tractor widths. A boundary offset of >3 m indicates a CRS mismatch between the mask and the export transform.

Histogram inspection: For multi-class masks, plot the distribution of final_area_ha per class:

import matplotlib.pyplot as plt

gdf_clean["final_area_ha"].hist(bins=30)
plt.xlabel("Field area (ha)")
plt.ylabel("Count")
plt.title("Field boundary area distribution")
plt.show()

Fields larger than 500 ha or smaller than 0.05 ha should be investigated individually — they are almost always mask artefacts rather than real agronomic parcels.

Known-value assertion on calibration target: If the farm has a reference field with known area from a certified surveyor (e.g., 42.7 ha from FSA records), assert within 2%:

ref_field_id = 5
ref_area_ha = 42.7
extracted_ha = float(gdf_clean.loc[gdf_clean["class_val"] == ref_field_id, "final_area_ha"].iloc[0])
assert abs(extracted_ha - ref_area_ha) / ref_area_ha < 0.02, (
    f"Extracted area {extracted_ha:.2f} ha differs from reference {ref_area_ha} ha by "
    f"{abs(extracted_ha - ref_area_ha) / ref_area_ha * 100:.1f}%"
)

Geometry validity check:

invalid = gdf_clean[~gdf_clean.geometry.is_valid]
if len(invalid) > 0:
    print(f"WARNING: {len(invalid)} invalid geometries in output")
    print(invalid[["class_val", "final_area_ha"]].to_string())

6. Integration with the broader pipeline

The GeoPackage produced by this workflow feeds directly into three downstream operations:

Variable rate prescription generation: The boundary polygon is the spatial envelope within which spatial interpolation for yield data applies kriging or IDW to yield monitor points. Boundaries with self-intersecting rings cause geopandas.clip() to fail or produce incorrect clipping extents, so topology validity is a hard prerequisite.

Shapefile export for machinery controllers: Before exporting a prescription to a John Deere GreenStar terminal or CNH AFS system, field boundaries must pass shapefile validation for farm equipment — specifically geometry type consistency checks that reject MultiPolygon layers. Call .explode(index_parts=False) to flatten MultiPolygons before the prescription overlay step.

Vegetation index masking: When calculating NDVI, NDRE, or SAVI from drone tiles, the field boundary polygon is used to clip the raster to the crop area before index computation, excluding roads, treelines, and irrigation infrastructure. This masking step is covered in band math and raster algebra in Python.

Scaling to multi-farm deployments: For organisations managing hundreds of farms, implement chunked processing using rasterio.windows to tile large orthomosaics, process each tile independently, and merge resulting GeoDataFrames using gpd.overlay(how="union") to handle cross-tile boundary segments. Mount cloud storage via rasterio’s virtual filesystem (/vsis3/bucket/path.tif) to eliminate local disk I/O bottlenecks in AWS or GCS environments.

FAQ

Why does my output have stair-step edges even after simplification? The mask resolution is likely above 0.5 m/px. Apply scipy.ndimage.gaussian_filter(mask.astype(float), sigma=1.5) and binarise at > 0.5 before morphological cleaning, then set simplify_tolerance=3.0. This pre-smooths the pixel boundary before vectorization and produces visibly cleaner curves.

Why does rasterio.features.shapes output pixel-space coordinates? The transform argument was omitted. Always pass transform=src.transform (the affine transform from rasterio.open()). Without it, shapes returns coordinates in column/row index space, which GeoPandas cannot interpret as geographic coordinates.

GeoPackage vs Shapefile for field boundaries? GeoPackage stores the CRS internally in the SQLite database and supports long column names. Shapefile splits metadata across 6+ sidecar files (.prj, .dbf, .shx, .cpg, .qpj, .xml) — any of which can be lost during transfer. Use GeoPackage for all internal storage; export to Shapefile only when a legacy FMIS or machinery terminal explicitly requires it.

Related

• Understanding CRS in Precision Agriculture — CRS selection, UTM zone assignment, and datum shifts for farm mapping
• Ingesting Multispectral Drone Imagery — loading and band-aligning MicaSense RedEdge-MX and DJI P4 Multispectral rasters before classification
• Orthomosaic Stitching Workflows — building the georeferenced mosaic from which segmentation masks are derived
• Shapefile Validation for Farm Equipment — geometry and schema checks required before uploading boundaries to machinery terminals
• Spatial Interpolation for Yield Data — using extracted field boundaries to clip and zone yield monitor point clouds