How to delineate field boundaries from imagery in Python

Problem statement

Where no agricultural parcel register exists, the boundaries have to come from the imagery — and a single date cannot provide them, because two adjacent fields growing the same crop at the same stage are the same colour.

What separates fields is time. Two fields differ because they were sown on different dates, with different varieties, and harvested in different weeks. A stack of dates turns an invisible boundary into a sharp one, and the delineation becomes a segmentation of a multi-band image whose bands are dates.

This guide builds that stack from real Sentinel-2 data, segments it, and checks the result against the artefacts that always appear.

Quick answer

Segment a temporal stack, not a single image:

import numpy as np
from skimage.segmentation import felzenszwalb

stack = np.dstack(ndre_by_date)                              # rows × cols × dates
stack = np.nan_to_num(stack, nan=float(np.nanmedian(stack)))

segments = felzenszwalb(stack, scale=200, sigma=1.0, min_size=250, channel_axis=2)
print(f"{segments.max() + 1:,} segments from {stack.shape[2]} dates")

scale controls how large the segments are and is the parameter to sweep. min_size is in pixels — at 20 m resolution, 250 pixels is 10 ha. That is too coarse even for Flevoland's large arable fields: on the run below it left 42% of segments spanning two registered parcels, against 31% at 50 pixels (2 ha).

6464 BRP parcels; {'Bouwland': 2659, 'Landschapselement': 2435, 'Grasland': 1321, 'Natuurterrein': 34, 'Braakland': 11, 'Overige': 4}

A clean temporal stack, a non-saturating index, segmentation

1. Build a clean temporal stack

Cloud-screened, gap-filled and on a common grid. Twenty dates across a season is plenty. The Flevoland series used below — Sentinel-2 scenes under 20% cloud from March to September 2025 — has 21 dates from 7 March to 6 September, and 19 of them are more than 60% clear over the area. The search returns 25 items for those 21 dates, because Earth Search lists four acquisitions twice; stack them unchecked and those dates count double.

2. Use an index that does not saturate

The stack's job is to separate fields, and NDVI is flat across the whole of mid-summer — on a June scene over these same fields, 50.7% of vegetation pixels were above 0.8. A red-edge index keeps its dynamic range and therefore keeps the boundaries.

3. Reduce the dates to a small number of components

Twenty correlated bands is more than any segmenter needs. A principal component transform reduces them to three or four that hold nearly all the variance, which is faster and less noisy.

4. Segment

Felzenszwalb, SLIC and watershed on a gradient all work. Felzenszwalb needs the least tuning and produces variable-sized segments, which suits fields. SLIC produces roughly equal-sized superpixels, which does not.

5. Merge segments that are the same field

Over-segmentation is easier to fix than under-segmentation. Merging adjacent segments whose temporal profiles are close recovers fields that were split by a tramline or a wet patch.

6. Filter the artefacts

Slivers, ribbons along roads and hedges, and segments that are mostly non-agricultural. Area and compactness catch most of them, but check both against a reference before choosing thresholds. In the Dutch parcel register (BRP) for this area, 55.4% of the 3,980 arable and grassland parcels are under half a hectare, yet together they hold only 2.6% of the farmed area; a 20 m delineation should not try to find them. Registered parcels of 2 ha or more have a median compactness of 0.66, and only 6% fall below 0.3.

Measure compactness on smoothed outlines. Polygons traced from pixels have staircase edges that inflate the perimeter: on the run below, 203 of 312 segments score below 0.3 as traced and 158 after a one-pixel simplification. Even smoothed, half the segments are less compact than 94% of registered parcels, because a segment that spans two fields is an irregular shape. A 0.3 threshold taken from the register would discard most of the output; low compactness here is a sign of merged fields, not of slivers.

7. Validate against something

Where any reference parcels exist — a register, a sample of hand-digitised fields, OSM — compare. The useful metrics are the fraction of reference parcels matched by one segment, the fraction split across several, and the fraction of segments that span two reference parcels. The first two move against the third as segments shrink, so report all of them. Example 4 computes them against the BRP.

Delineating fields in Flevoland

The data used on this page

The area is 5.55–5.75°E, 52.48–52.56°N in Flevoland: large arable fields with a public parcel register to check against. The register is the 2025 BRP (Basisregistratie Gewaspercelen) from PDOK's OGC API; the imagery is Sentinel-2 L2A from Element 84's Earth Search, read directly from the cloud-optimised GeoTIFFs.

import os
import geopandas as gpd
import requests

BBOX = [5.55, 52.48, 5.75, 52.56]

if not os.path.exists("brp.gpkg"):
    url = "https://api.pdok.nl/rvo/gewaspercelen/ogc/v1/collections/brpgewas/items"
    params, feats = {"f": "json", "bbox": ",".join(map(str, BBOX)), "limit": 1000}, []
    while url:
        page = requests.get(url, params=params, timeout=120).json()
        feats += page["features"]
        url, params = next((l["href"] for l in page["links"] if l["rel"] == "next"), None), None
    gpd.GeoDataFrame.from_features(feats, crs=4326).to_file("brp.gpkg")
brp = gpd.read_file("brp.gpkg")
print(len(brp), "BRP parcels;", brp["category"].value_counts().to_dict())

6464 BRP parcels; {'Bouwland': 2659, 'Landschapselement': 2435, 'Grasland': 1321, 'Natuurterrein': 34, 'Braakland': 11, 'Overige': 4}

Example 1 — build the stack from a STAC search

import json, urllib.request, os, numpy as np, rasterio
from rasterio.warp import transform_bounds
from rasterio.windows import from_bounds

os.environ.setdefault("GDAL_DISABLE_READDIR_ON_OPEN", "EMPTY_DIR")
os.environ.setdefault("AWS_NO_SIGN_REQUEST", "YES")
SCALE = 1e-4

def search(bbox, start, end, cloud=20, limit=60):
    body = {"collections": ["sentinel-2-l2a"], "bbox": bbox,
            "datetime": f"{start}T00:00:00Z/{end}T23:59:59Z", "limit": limit,
            "query": {"eo:cloud_cover": {"lt": cloud}}}
    req = urllib.request.Request("https://earth-search.aws.element84.com/v1/search",
        data=json.dumps(body).encode(), headers={"Content-Type": "application/json"})
    with urllib.request.urlopen(req, timeout=120) as r:
        return json.loads(r.read())["features"]

def window_read(href, bbox, step=1):
    with rasterio.open(href) as src:
        b = transform_bounds("EPSG:4326", src.crs, *bbox)
        win = from_bounds(*b, transform=src.transform)
        return (src.read(1, window=win).astype("float32")[::step, ::step],
                src.crs, src.window_transform(win))

def build_stack(bbox, start, end, min_clear=0.6):
    # Earth Search holds some acquisitions twice (reprocessed); keep one per date
    best = {}
    for f in search(bbox, start, end):
        d, cc = f["properties"]["datetime"][:10], f["properties"]["eo:cloud_cover"]
        if d not in best or cc < best[d]["properties"]["eo:cloud_cover"]:
            best[d] = f
    layers, dates, dropped = [], [], []
    for d in sorted(best):
        f = best[d]
        a = f["assets"]
        nir, crs, tr = window_read(a["nir"]["href"], bbox, 2)
        re1, _, tr20 = window_read(a["rededge1"]["href"], bbox, 1)
        scl, _, _ = window_read(a["scl"]["href"], bbox, 1)
        h = min(nir.shape[0], re1.shape[0], scl.shape[0])
        w = min(nir.shape[1], re1.shape[1], scl.shape[1])
        nir, re1, scl = nir[:h, :w] * SCALE, re1[:h, :w] * SCALE, scl[:h, :w]
        clear = np.isin(scl, [4, 5, 6])
        if clear.mean() < min_clear:
            dropped.append(d)
            continue
        with np.errstate(invalid="ignore", divide="ignore"):
            ndre = np.where(clear, (nir - re1) / (nir + re1), np.nan)
        layers.append(ndre)
        dates.append(d)
    print(f"{len(best)} dates, {len(dates)} kept, dropped below {min_clear:.0%} clear: {dropped}")
    return np.dstack(layers), dates, crs, tr20

stack, dates, crs, tr = build_stack(BBOX, "2025-03-01", "2025-09-30")
print(stack.shape, dates[0], "to", dates[-1])

21 dates, 19 kept, dropped below 60% clear: ['2025-04-26', '2025-05-09']
(470, 695, 19) 2025-03-07 to 2025-09-06

Dropping dates below 60% clear is a decision worth reporting: here it removed 2 of 21 dates, both in the spring gap between sowing and emergence. The duplicate check matters as much. Without it, the stack has 22 layers, 17 March, 6 April and 14 August each appear twice, and the PCA weights those dates double.

Example 2 — reduce, segment and vectorise

import numpy as np, geopandas as gpd
from sklearn.decomposition import PCA
from skimage.segmentation import felzenszwalb
from rasterio.features import shapes
from shapely.geometry import shape

def vectorise(labels, transform, crs):
    polys = [{"seg": int(v), "geometry": shape(g)}
             for g, v in shapes(labels.astype("int32"), transform=transform)]
    return gpd.GeoDataFrame(polys, crs=crs).dissolve("seg", as_index=False)

def delineate(stack, transform, crs, n_components=4, scale=200, min_size=250):
    h, w, t = stack.shape
    flat = stack.reshape(-1, t)
    med = np.nanmedian(flat, axis=0)
    flat = np.where(np.isfinite(flat), flat, med)

    pca = PCA(n_components=min(n_components, t)).fit(flat)
    comp = pca.transform(flat).reshape(h, w, -1)
    print("explained variance:", np.round(pca.explained_variance_ratio_, 3))

    comp = (comp - comp.min(axis=(0, 1))) / np.ptp(comp, axis=(0, 1))
    seg = felzenszwalb(comp, scale=scale, sigma=1.0, min_size=min_size,
                       channel_axis=2)

    g = vectorise(seg, transform, crs)
    g["area_ha"] = g.area / 1e4
    smooth = g.simplify(transform.a)            # staircase edges inflate the perimeter
    g["compactness"] = 4 * np.pi * smooth.area / (smooth.length ** 2)
    print(f"{len(g):,} segments; median {g.area_ha.median():.1f} ha, "
          f"{(g.area_ha < 0.5).sum()} under 0.5 ha, "
          f"{(g.compactness < 0.3).sum()} with compactness below 0.3")
    return g, seg

segs, seg = delineate(stack, tr, crs)
segs_fine, seg_fine = delineate(stack, tr, crs, min_size=50)

explained variance: [0.549 0.179 0.143 0.055]
312 segments; median 22.6 ha, 0 under 0.5 ha, 158 with compactness below 0.3
explained variance: [0.549 0.179 0.143 0.055]
673 segments; median 8.6 ha, 0 under 0.5 ha, 280 with compactness below 0.3

Four components hold 93% of the variance of 19 dates. The first run's median segment is 22.6 ha where the median registered parcel of 2 ha or more is 5.8 ha, which is the first sign that min_size=250 is merging fields.

Example 3 — merge over-segmented neighbours

import numpy as np, geopandas as gpd
from scipy import ndimage

def merge_similar(seg, stack, threshold=0.05):
    """Merge adjacent segments whose mean temporal profiles are close."""
    n = seg.max() + 1
    profiles = np.stack([
        np.nanmean(stack[seg == i], axis=0) if (seg == i).any() else np.zeros(stack.shape[2])
        for i in range(n)])

    parent = np.arange(n)
    def find(a):
        while parent[a] != a:
            parent[a] = parent[parent[a]]
            a = parent[a]
        return a

    # adjacency from a one-pixel dilation of each label boundary
    right = np.c_[seg[:, :-1].ravel(), seg[:, 1:].ravel()]
    down = np.c_[seg[:-1, :].ravel(), seg[1:, :].ravel()]
    pairs = np.unique(np.vstack([right, down]), axis=0)
    pairs = pairs[pairs[:, 0] != pairs[:, 1]]

    for a, b in pairs:
        d = np.nanmean(np.abs(profiles[a] - profiles[b]))
        if d < threshold:
            ra, rb = find(a), find(b)
            if ra != rb:
                parent[ra] = rb

    merged = np.vectorize(find)(seg)
    print(f"{n:,} segments → {len(np.unique(merged)):,} after merging at {threshold}")
    return merged

merged = merge_similar(seg_fine, stack)
segs_merged = vectorise(merged, tr, crs)

673 segments → 597 after merging at 0.05

Merging on the mean absolute difference between temporal profiles is crude and effective. The threshold is the parameter to sweep, and the right value depends on the index's dynamic range — which is another reason to use one that does not saturate.

Example 4 — validate against the parcel register

import numpy as np, geopandas as gpd

def validate(segs, parcels, min_ha=2, match=0.8, span=0.2):
    """Share of parcels inside one segment, and of segments spanning two parcels."""
    p = parcels[parcels.area >= min_ha * 1e4]
    ov = gpd.overlay(p[["pid", "geometry"]], segs[["seg", "geometry"]],
                     how="intersection", keep_geom_type=True)
    ov["a"] = ov.area
    one = ov.groupby("pid")["a"].max() / p.set_index("pid").area >= match
    seg_area = segs.set_index("seg").area
    by_seg = ov.groupby("seg")["a"]
    farmed = by_seg.sum() / seg_area >= 0.5          # ignore segments over towns and water
    second = by_seg.apply(lambda a: np.sort(a.values)[-2] if len(a) > 1 else 0) / seg_area
    spans = second[farmed[farmed].index] >= span
    print(f"{len(p):,} parcels ≥ {min_ha} ha: {one.mean():.1%} in one segment, "
          f"{1 - one.mean():.1%} split; {spans.size:,} farmed segments, "
          f"{spans.mean():.1%} spanning two parcels")

parcels = brp[brp["category"].isin(["Bouwland", "Grasland"])].to_crs(crs).reset_index(drop=True)
parcels["pid"] = np.arange(len(parcels))
for name, s in [("min_size 250", segs), ("min_size 50", segs_fine), ("50, merged", segs_merged)]:
    print(name, end=": ")
    validate(s, parcels)

min_size 250: 1,361 parcels ≥ 2 ha: 91.0% in one segment, 9.0% split; 273 farmed segments, 41.8% spanning two parcels
min_size 50: 1,361 parcels ≥ 2 ha: 82.3% in one segment, 17.7% split; 551 farmed segments, 29.4% spanning two parcels
50, merged: 1,361 parcels ≥ 2 ha: 83.2% in one segment, 16.8% split; 496 farmed segments, 30.6% spanning two parcels

A parcel counts as matched when one segment covers 80% of it, and a segment spans two parcels when the second-largest parcel inside it fills a fifth of it. Both thresholds are choices; report them with the numbers.

Why a single date cannot work

Why a single date cannot work

Field boundaries are management boundaries, and management is only visible when it differs. Two neighbouring fields of winter wheat sown a week apart are indistinguishable at peak canopy and clearly different at emergence and at senescence. A single image samples one moment and therefore finds only the boundaries that happen to be expressed then.

Why the index choice is not cosmetic here

A saturated index has no dynamic range, so adjacent fields that differ by real biomass have the same value and no boundary appears. On 12 June 2025, 462 of the 867 registered fields over 1 ha of wheat, potato, onion, sugar beet, carrot and grass in this area already had a mean NDVI above 0.8. Across those 462 fields, field-mean NDVI spanned 0.800–0.911, a range of 0.110, while NDRE spanned 0.515–0.778, a range of 0.263 — more than twice as much room for a boundary to show.

Why over-segmentation is the safer failure

The Flevoland run shows both failures. Dropping min_size from 250 to 50 pixels cut the segments spanning two parcels from 41.8% to 29.4%, and raised the parcels split across segments from 9.0% to 17.7%. Of the 104 segments still spanning two parcels at min_size=250 after merging, 46 join two parcels declared with the same crop, which is the case a temporal stack separates only if their dates differ. Merging by profile at 0.05 changed little here: it rejoined 76 segments but moved the split share by under a point.

A field split into three segments can be merged with a rule. Three fields merged into one segment cannot be separated without going back to the imagery, and any statistic computed on it is a weighted average of three different crops. Tune towards smaller segments and merge.

Why to validate on match fractions rather than on an overall accuracy

"Accuracy" for a delineation is ill-defined, because the reference itself is a choice of parcel type. The useful numbers are how many reference parcels are covered by exactly one segment, how many are split, and how many segments straddle two parcels. Those three name the failure modes directly.

Tramlines, grassland and irrigation circles

• Tramlines split a field in the segmentation; merging fixes it.
• Grassland has no strong seasonal signal and delineates poorly.
• Irrigation circles are easy, hedges are hard.
• Mixed cropping within a field produces genuine internal boundaries.
• The stack must be on one grid. Reproject before stacking, not after.
• Gap-fill before PCA. NaNs propagate through the decomposition.
• min_size is in pixels. Convert to hectares before choosing it.
• Record the dates used. A delineation is a statement about a season.

Internal links

• Field boundary data explained: where it comes from — whether you need to delineate at all
• Delineated field boundaries merge fields or swallow roads — the failure modes
• Vegetation indices explained: NDVI, EVI, NDRE and when each fails — why the index matters here
• How to classify crop types from a satellite time series — the next step
• How to classify land cover in Python — the pixel-based alternative
• How to convert raster to vector in Python — vectorising the segments
• How to build a STAC-based data cube with xarray — assembling the stack
• How to mask clouds in Sentinel-2 imagery in Python — screening the dates

FAQ

Can I delineate fields from a single image?

No. Adjacent fields with the same crop at the same stage are the same colour. Boundaries appear across a season, not within one date.

Which index should I stack?

One that does not saturate — a red-edge index. NDVI is flat across most of mid-summer, which is a large part of the season contributing nothing.

How many dates do I need?

Fifteen to twenty across a season is ample, provided they span emergence and senescence as well as the peak. The Flevoland run used 19.

Which segmentation algorithm?

Felzenszwalb needs the least tuning and produces variable-sized segments, which suits fields. SLIC's equal-sized superpixels do not.

Should I aim for over- or under-segmentation?

Over. Merging adjacent segments with similar temporal profiles is a rule; separating merged fields is not.

How do I know if it worked?

Compare with any reference parcels you have, and report three numbers: parcels matched by one segment, parcels split across several, and segments spanning two parcels.