Delineated field boundaries merge fields or swallow roads

Problem statement

The segmentation ran, the polygons look plausible, and two things are wrong. Some adjacent fields with the same crop have become one polygon, and some polygons include the road, the ditch and the hedge between two fields.

Both are the same mechanism from opposite directions: the segmenter grouped pixels whose temporal profiles were similar enough, and "similar enough" is one parameter. Tighten it and fields split into tramline-sized fragments; loosen it and neighbouring wheat fields merge. There is no setting that does both, which is why the fix is in the features and the post-processing rather than in the threshold.

Quick answer

Diagnose by measuring what the segments contain, not by looking at them:

import numpy as np, geopandas as gpd

seg = segments.copy()
seg["area_ha"] = seg.area / 1e4
seg["compactness"] = 4 * np.pi * seg.area / (seg.length ** 2)

print(f"{len(seg):,} segments; median {seg.area_ha.median():.1f} ha, "
      f"max {seg.area_ha.max():.1f} ha")
print(f"over 3x the median area: {(seg.area_ha > 3 * seg.area_ha.median()).sum()} "
      "(candidate merges)")
print(f"compactness below 0.3: {(seg.compactness < 0.3).sum()} "
      "(candidate ribbons along roads)")

673 segments; median 8.6 ha, max 662.8 ha
over 3x the median area: 100 (candidate merges)
compactness below 0.3: 341 (candidate ribbons along roads)

The OpenStreetMap farmland layer for the same Flevoland box has 492 polygons with a median of 18.62 ha and a maximum of 234.4 ha; 16 are under half a hectare and 54 have a compactness below 0.3. The segmentation here has eleven times as many low-compactness shapes — partly real ribbons, partly the staircase edges of polygons traced from pixels. A delineation whose profile is very different from a reference is describing something other than fields.

Is the boundary in the data, is the index saturated, add dates

1. Check whether the boundary exists in the data at all

Two fields of the same crop, sown the same week, differ only slightly. Plot the temporal profiles of the two halves of a merged segment: if they are indistinguishable, no segmenter can separate them and the fix is more dates, a different index, or an external boundary source.

2. Check the index for saturation

A saturated index has no dynamic range, so real differences disappear. On a June 2025 scene over these fields, 50.7% of vegetation pixels were above NDVI 0.8, and across the 462 fields already at full canopy on 12 June, field-mean NDVI spanned 0.110 where NDRE spanned 0.263 (measured here). A red-edge index keeps its range.

3. Add the dates that separate fields

Emergence and senescence separate fields far better than peak canopy does, because sowing and harvest dates differ even when the crop does not. If the stack is peak-heavy, the merges are inevitable.

4. Over-segment deliberately, then merge

Set the segmentation scale small enough that fields split, then merge adjacent segments whose profiles are close. Merging is a rule you control; separating merged fields is not.

5. Cut the linear features out rather than tuning them away

Roads, hedges and ditches are in OpenStreetMap. Buffering them and erasing them from the segments removes road pixels from field polygons without touching the threshold — but check what it does to the rest. In Flevoland it removed 160 ha of road-edge fragments and cut the segments into nearly twice as many pieces, without reducing the share that spans two parcels (Example 2): here the parcels are separated by ditches OpenStreetMap mostly lacks, and the merges are same-crop fields, which no road layer separates.

6. Use an edge signal as well as a region signal

Segmenters that use only region similarity have no reason to place a boundary at a hedge. Adding a gradient layer — the temporal standard deviation, or the magnitude of the gradient of the first principal component — gives the segmenter an edge to snap to.

7. Filter the artefacts by shape

Ribbons have low compactness, slivers have small area, and segments that span two fields are unusually large. All three are one line each.

Testing a merged segment and erasing roads

The data used on this page

The Flevoland box from the delineation guide — 5.55–5.75°E, 52.48–52.56°N — with the same inputs: 19 cloud-screened Sentinel-2 NDRE dates from March to September 2025, reduced to four principal components and segmented with Felzenszwalb (scale=200, min_size=50); the 2025 BRP parcel register as the reference; and OpenStreetMap roads and waterways from the Geofabrik Flevoland extract.

import json
import os
import urllib.request

import geopandas as gpd
import numpy as np
import rasterio
import requests
from rasterio.features import rasterize, shapes
from rasterio.warp import transform_bounds
from rasterio.windows import from_bounds
from shapely.geometry import shape
from skimage.segmentation import felzenszwalb
from sklearn.decomposition import PCA

os.environ.update(AWS_NO_SIGN_REQUEST="YES", GDAL_DISABLE_READDIR_ON_OPEN="EMPTY_DIR")
BBOX = [5.55, 52.48, 5.75, 52.56]

if not os.path.exists("stack.npz"):
    body = {"collections": ["sentinel-2-l2a"], "bbox": BBOX, "limit": 60,
            "datetime": "2025-03-01T00:00:00Z/2025-09-30T23:59:59Z",
            "query": {"eo:cloud_cover": {"lt": 20}}}
    request = urllib.request.Request("https://earth-search.aws.element84.com/v1/search",
                                     data=json.dumps(body).encode(),
                                     headers={"Content-Type": "application/json"})
    best = {}
    for f in json.loads(urllib.request.urlopen(request, timeout=120).read())["features"]:
        d = f["properties"]["datetime"][:10]
        if d not in best or f["properties"]["eo:cloud_cover"] < best[d]["properties"]["eo:cloud_cover"]:
            best[d] = f
    layers = []
    for d in sorted(best):
        a = best[d]["assets"]
        def read(key, step=1):
            with rasterio.open(a[key]["href"]) as src:
                win = from_bounds(*transform_bounds("EPSG:4326", src.crs, *BBOX), transform=src.transform)
                return src.read(1, window=win).astype("float32")[::step, ::step], src.window_transform(win)
        nir, _ = read("nir", 2)
        re1, transform = read("rededge1")
        scl, _ = read("scl")
        h, w = (min(x.shape[0] for x in (nir, re1, scl)), min(x.shape[1] for x in (nir, re1, scl)))
        nir, re1, clear = nir[:h, :w], re1[:h, :w], np.isin(scl[:h, :w], [4, 5, 6])
        if clear.mean() < 0.6:
            continue
        with np.errstate(invalid="ignore", divide="ignore"):
            layers.append(np.where(clear, (nir - re1) / (nir + re1), np.nan))
    np.savez_compressed("stack.npz", stack=np.dstack(layers), transform=np.array(transform)[:6])
saved = np.load("stack.npz")
stack, transform = saved["stack"], rasterio.Affine(*saved["transform"])

h, w, t = stack.shape
flat = stack.reshape(-1, t)
flat = np.where(np.isfinite(flat), flat, np.nanmedian(flat, axis=0))
comp = PCA(n_components=4).fit_transform(flat).reshape(h, w, -1)
comp = (comp - comp.min(axis=(0, 1))) / np.ptp(comp, axis=(0, 1))
labels = felzenszwalb(comp, scale=200, sigma=1.0, min_size=50, channel_axis=2)
segments = gpd.GeoDataFrame([{"seg": int(v), "geometry": shape(g)} for g, v in
                             shapes(labels.astype("int32"), transform=transform)],
                            crs=32631).dissolve("seg", as_index=False)

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").to_crs(32631)
parcels = brp[brp["category"].isin(["Bouwland", "Grasland"]) & (brp.area >= 2e4)].reset_index(drop=True)
parcel_id = rasterize(zip(parcels.geometry, range(1, len(parcels) + 1)), out_shape=(h, w),
                      transform=transform, fill=0, dtype="int32")

if not os.path.exists("flevoland.zip"):
    with requests.get("https://download.geofabrik.de/europe/netherlands/flevoland-260101-free.shp.zip",
                      stream=True, timeout=600) as reply:
        with open("flevoland.zip", "wb") as out:
            for chunk in reply.iter_content(1 << 20):
                out.write(chunk)
box = tuple(transform_bounds("EPSG:4326", "EPSG:4326", *BBOX))
roads = gpd.read_file("zip://flevoland.zip!gis_osm_roads_free_1.shp", bbox=box).to_crs(32631)
waterways = gpd.read_file("zip://flevoland.zip!gis_osm_waterways_free_1.shp", bbox=box).to_crs(32631)
print(f"{t} NDRE dates, {len(segments)} segments; {len(parcels):,} BRP parcels of 2 ha or more; "
      f"{len(roads):,} OSM roads, {len(waterways)} waterways")

19 NDRE dates, 673 segments; 1,361 BRP parcels of 2 ha or more; 3,837 OSM roads, 149 waterways

Example 1 — test whether a merged segment really is two fields

import numpy as np
from sklearn.cluster import KMeans

def split_test(stack, seg, seg_id, min_pixels=200):
    """Cluster the temporal profiles inside one segment into two."""
    mask = seg == seg_id
    if mask.sum() < min_pixels:
        return None
    X = stack[mask]
    ok = np.isfinite(X).all(axis=1)
    if ok.sum() < min_pixels:
        return None

    km = KMeans(n_clusters=2, n_init=10, random_state=0).fit(X[ok])
    a, b = X[ok][km.labels_ == 0], X[ok][km.labels_ == 1]
    sep = np.abs(a.mean(axis=0) - b.mean(axis=0))
    pooled = np.sqrt((a.var(axis=0) + b.var(axis=0)) / 2)
    effect = float(np.nanmax(sep / np.where(pooled == 0, np.nan, pooled)))

    return {"pixels": int(ok.sum()),
            "balance": round(float(min(len(a), len(b)) / ok.sum()), 3),
            "max_effect_size": round(effect, 2),
            "verdict": "two fields" if effect > 1.5 and min(len(a), len(b)) / ok.sum() > 0.2
                       else "one field"}

Run against the register, on every segment of at least 200 pixels that is mostly farmland:

import pandas as pd

rows = []
for sid in np.unique(labels):                     # labels: the segmentation raster from the setup
    inside = labels == sid
    ids = parcel_id[inside]
    ids = ids[ids > 0]
    if inside.sum() < 200 or len(ids) < 0.5 * inside.sum():
        continue
    counts = np.bincount(ids)
    first, second = np.argsort(counts)[::-1][:2]
    if counts[second] >= 0.2 * inside.sum():
        kind = ("two parcels, same crop" if parcels["gewas"][first - 1] == parcels["gewas"][second - 1]
                else "two parcels, different crops")
    elif counts[first] >= 0.9 * inside.sum():
        kind = "one parcel"
    else:
        continue
    result = split_test(stack, labels, sid)
    if result:
        rows.append({"segment": kind, **result})
tested = pd.DataFrame(rows)
print(tested.groupby("segment")["verdict"].value_counts().unstack(fill_value=0))
print(tested.groupby("segment")["max_effect_size"].min().round(2).to_string())

verdict                       one field  two fields
segment                                            
one parcel                           26           8
two parcels, different crops          4          29
two parcels, same crop               20          12
segment
one parcel                      1.96
two parcels, different crops    2.80
two parcels, same crop          2.52

The test finds most merges of two different crops and a minority of same-crop merges, and it flags about a quarter of single-parcel segments as two fields. The effect-size threshold of 1.5 does no work on a 19-date stack: taking the maximum over 19 dates, every segment exceeded it, so the verdict rested on the balance between the two clusters. Treat "two fields" as a candidate for review, and same-crop merges as mostly beyond what these data can separate.

Example 2 — erase linear features

import geopandas as gpd

def erase_linear(seg, roads, hedges=None, buffer_m=4, min_area_ha=0.5):
    lines = [roads] + ([hedges] if hedges is not None else [])
    mask = gpd.GeoDataFrame(
        geometry=[g.buffer(buffer_m) for gdf in lines for g in gdf.geometry],
        crs=seg.crs).dissolve()

    out = gpd.overlay(seg, mask, how="difference").explode(index_parts=False)
    out["area_ha"] = out.area / 1e4
    dropped = out[out.area_ha < min_area_ha]
    if len(dropped):
        print(f"erasing linear features split off {len(dropped)} fragments "
              f"({dropped.area_ha.sum():.2f} ha)")
    return out[out.area_ha >= min_area_ha]

Buffering by half the road width plus the pixel size is about right. Too small and the ribbon survives; too large and it eats the headland. What it does to the segmentation, measured against the register:

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
    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"  {one.mean():.1%} of parcels in one segment; {spans.size:,} farmed segments, "
          f"{spans.mean():.1%} spanning two parcels")


rural = roads[roads["fclass"].isin(["track", "unclassified", "tertiary", "secondary", "primary",
                                     "trunk", "motorway"])]
reference = parcels.assign(pid=np.arange(len(parcels)))
print(f"{len(rural):,} rural roads of {len(roads):,}")
for label, erased in (("as segmented", segments),
                      ("rural roads + drains erased", erase_linear(segments, rural, waterways))):
    erased = erased.reset_index(drop=True).assign(seg=lambda g: np.arange(len(g)))
    print(f"{label}: {len(erased):,} segments, median {(erased.area / 1e4).median():.1f} ha")
    validate(erased, reference)

1,000 rural roads of 3,837
erasing linear features split off 2741 fragments (160.10 ha)
as segmented: 673 segments, median 8.6 ha
  82.3% of parcels in one segment; 551 farmed segments, 29.4% spanning two parcels
rural roads + drains erased: 1,220 segments, median 5.7 ha
  82.3% of parcels in one segment; 766 farmed segments, 34.3% spanning two parcels

The erase did what it says — road-edge pixels are gone — and made the segmentation no better against the register: more pieces, a smaller median, and a slightly larger share of segments spanning two parcels. Flevoland's parcels are divided by ditches that OpenStreetMap mostly does not map (149 waterways in the whole box), and its merges are same-crop fields. Erase linear features when the ribbons are the problem, and merge the fragments afterwards; do not expect it to split fields.

Example 3 — profile the result against a reference

import numpy as np, geopandas as gpd, pandas as pd

def compare_profile(seg, reference, crs=28992):
    rows = []
    for name, g in (("delineated", seg), ("reference", reference)):
        g = g.to_crs(crs)
        a = g.area / 1e4
        c = 4 * np.pi * g.area / (g.length ** 2)
        rows.append({"layer": name, "n": len(g),
                     "median_ha": round(float(a.median()), 2),
                     "max_ha": round(float(a.max()), 1),
                     "under_0.5_ha": int((a < 0.5).sum()),
                     "compact_below_0.3": int((c < 0.3).sum()),
                     "total_ha": round(float(a.sum()))})
    df = pd.DataFrame(rows)
    print(df.to_string(index=False))
    return df


compare_profile(segments, parcels)          # registered arable and grass parcels of 2 ha or more

     layer    n  median_ha  max_ha  under_0.5_ha  compact_below_0.3  total_ha
delineated  673       8.60   662.7             0                341     13064
 reference 1361       5.83    33.2             0                 79      9294

Comparing the shape profile with a reference layer of the same area catches both failure modes at once: merges push the median and maximum area up, ribbons push the low-compactness count up. Here the largest segment is twenty times the largest parcel, and low-compactness shapes are four times as common. The reference is the register's parcels of 2 ha or more, because a 20 m segmentation cannot see the rest: the full register is dominated by field margins and slivers, with a median parcel of 0.3 ha.

One parameter behind both failures

Why the same parameter causes both failures

A region-merging segmenter groups adjacent pixels whose values differ by less than a threshold. The distance between two fields of different crops is large; between two fields of the same crop it is small; between a field and the tramline running through it, smaller still. Those three distances are not separated by any single threshold, so one setting always merges something that should be split or splits something that should be merged.

Why more dates help and a better algorithm does not

Two fields of the same crop differ in when they were sown and when they were harvested, which is expressed at emergence and senescence and almost nowhere else. A stack that samples those periods contains the boundary; one that samples only the peak does not, and no segmenter can find information that is not there.

Why linear features should be erased rather than segmented around

A hedge or a road is three to eight metres wide, which at 10–20 m resolution is a partial pixel with a mixed spectral signature that resembles neither field. Asking a segmenter to place a boundary on such a pixel is asking for trouble; erasing a known vector layer is deterministic, reversible and uses information the imagery does not contain.

Why over-segmentation is the safe direction

A merged segment has lost information irretrievably: the mean profile of two fields is not either field, and every statistic computed on it is a weighted average of two crops. A split field retains all of it, and merging adjacent segments by profile similarity is a rule with a parameter you can sweep and validate.

Tramlines, two-crop fields and grassland

• Tramlines split fields at fine scales; merging recovers them.
• A field with two crops genuinely has an internal boundary.
• Grassland has little seasonal signal and delineates poorly whatever you do.
• Irrigation circles are easy; small irregular fields are not.
• min_size is in pixels. Convert to hectares before setting it.
• Erasing roads splits a field in two where a track crosses it; merge afterwards.
• Compare the shape profile with a reference rather than judging by eye.
• Record the dates used. A delineation is a statement about one season.

Internal links

• How to delineate field boundaries from imagery in Python — the pipeline this fixes
• Field boundary data explained: where it comes from — whether to delineate at all
• Vegetation indices explained: NDVI, EVI, NDRE and when each fails — why saturation causes merges
• NDVI stops responding in a dense canopy — the flat peak
• How to simplify geometry in GeoPandas — cleaning the polygons
• How to fix gaps and overlaps in a polygon coverage — after erasing
• How to download OpenStreetMap data in Python — getting the roads and hedges
• How to convert raster to vector in Python — the vectorisation step

FAQ

Why did two fields merge into one polygon?

Because their temporal profiles were too similar for the threshold — usually the same crop sown at the same time, or a stack dominated by saturated peak-season dates.

Why does my polygon include the road?

Because a road pixel at 10–20 m resolution is a mixed pixel that resembles neither field. Erase a buffered road layer rather than trying to tune it away.

Can I just change the segmentation threshold?

No. The same parameter controls both failures in opposite directions. Over-segment and merge instead.

How do I know whether a merged segment is really two fields?

Cluster the temporal profiles inside it into two and measure the separation. An effect size above about 1.5 on at least one date with both halves substantial means two fields.

Which dates matter most?

Emergence and senescence. Sowing and harvest dates differ between fields even when the crop does not; peak canopy does not separate them.

How do I check the result?

Compare the area and compactness profile with a reference layer for the same area. Merges raise the median area; ribbons raise the low-compactness count.