How to build management zones with clustering in Python

On this page

Problem statement

Clustering a field into zones is four lines of scikit-learn and a great many decisions: which layers, how to normalise them, how many clusters, how to smooth the result, and how to tell whether the zones mean anything.

The failure is always the same shape — zones that look convincing, follow the previous year's weather, and explain almost none of the variation in any other year. The fix is not a better algorithm; it is stable inputs, a held-out year and a test that produces a number.

This guide builds the zoning on a real field — the 22 ha corn field at Porto Felloni, near Ferrara, whose 2016–2018 yield maps are published on Zenodo (CC BY 4.0) — smooths it into something a machine can drive, and evaluates it on the year it was not fitted to.

Quick answer

import numpy as np
from sklearn.cluster import KMeans
from sklearn.preprocessing import StandardScaler

X = np.column_stack([mean_z.ravel(), stability.ravel()])
ok = np.isfinite(X).all(axis=1)

labels = np.full(X.shape[0], -1, dtype="int16")
km = KMeans(n_clusters=3, n_init=20, random_state=0).fit(StandardScaler().fit_transform(X[ok]))
labels[ok] = km.labels_

# order zones by productivity so the labels mean the same thing every run
order = np.argsort([mean_z.ravel()[labels == k].mean() for k in range(3)])
labels[ok] = np.vectorize({old: new for new, old in enumerate(order)}.get)(labels[ok])
zones = labels.reshape(mean_z.shape)
print("share of the field in each zone:", {k: f"{(zones == k).sum() / (zones >= 0).sum():.1%}" for k in range(3)})

share of the field in each zone: {0: '15.5%', 1: '30.2%', 2: '54.3%'}

Re-ordering the labels by productivity is not cosmetic. K-means labels are arbitrary, so without it zone 1 means something different every time the script runs, and a prescription written against last week's run applies the wrong rates.

Six steps; the two that are usually skipped are label ordering and validation.

Stable layers, per-year normalisation, k-means

1. Assemble layers that do not change between years

A multi-year yield or vegetation index mean, its temporal standard deviation, elevation and its derivatives where the field has relief, and soil electrical conductivity if a survey exists. Put every year on a common grid first. A layer earns its place by explaining something: this field is reclaimed land in the Po delta, where the 30 m Copernicus DEM spans −2.93 to −2.45 m, and adding it made the held-out result worse (Example 2).

2. Normalise each year within the field before combining

Absolute levels differ between seasons. A within-field z-score per year removes the season and keeps the pattern, which is the thing the zones are for.

3. Include stability as its own variable

The temporal standard deviation separates a consistently mediocre area from an alternately excellent and terrible one, which a mean cannot.

4. Choose the number of clusters, and justify it

Sweep two to six and look at the silhouette score, the variance explained in a held-out year, and whether the smallest zone is large enough to manage. They do not have to agree. On this field the held-out figure rose with every extra cluster, from 6.5% at two to 14.4% at six, while the smallest zone shrank from 24% of the field to under 1%. The held-out number stops the in-sample one from flattering you (54.7% at three clusters); the smallest zone and the machine are what cap k.

5. Order the labels and smooth the result

A modal filter removes speckle, a minimum area merges unmanageable fragments into their neighbours, and simplifying at roughly half the implement width removes detail the machine cannot follow. Each of the three has a tempting version that is wrong: a median filter is not a modal one, dropping small fragments leaves holes in the field, and simplifying each zone polygon on its own opens gaps and overlaps between them (Example 3).

6. Validate on a year that was not used

Fit the zones on some years, compute the variance explained in another. That number is the zoning's value, and it is usually much lower than the within-sample figure: here 7.2% of the 2018 yield variance against 54.7% of the 2016–17 pattern the zones were fitted to.

7. Publish the inputs with the zones

Which years, which layers, which parameters. A zone map without them cannot be reproduced or updated.

Measured on the Ferrara field: held out, more clusters still helped a little — the smallest zone, not the score, is what capped k.

Choosing k on a held-out year

The data used on this page

The three kriged yield maps in the Zenodo archive (2016, 2017, 2018; 2.2 m cells on one grid) and the Copernicus 30 m DEM resampled onto the same grid. The kriged grids carry the yield points' Gauss-Boaga eastings with the leading million dropped, which lands within about 37 m of UTM zone 32N (EPSG:32632) — close enough to sample a 30 m DEM, and the grid is treated as 32632 throughout. Zones are fitted on 2016 and 2017; 2018 is held out.

import os
import zipfile

import numpy as np
import rasterio
import requests
from rasterio.warp import Resampling, reproject

if not os.path.exists("FIELD_YIELD"):
    reply = requests.get("https://zenodo.org/records/14216165/files/FIELD_YIELD.zip?download=1", timeout=300)
    reply.raise_for_status()
    open("FIELD_YIELD.zip", "wb").write(reply.content)
    zipfile.ZipFile("FIELD_YIELD.zip").extractall(".")

CRS = "EPSG:32632"
years = {}
for year in (2016, 2017, 2018):
    with rasterio.open(f"FIELD_YIELD/yield{year}krig/{year}yieldkrig.tif") as src:
        years[year] = src.read(1, masked=True).filled(np.nan)
        transform = src.transform
inside = np.isfinite(years[2016])

elevation = np.full(inside.shape, np.nan, dtype="float32")
with rasterio.open("/vsicurl/https://copernicus-dem-30m.s3.amazonaws.com/"
                   "Copernicus_DSM_COG_10_N44_00_E012_00_DEM/Copernicus_DSM_COG_10_N44_00_E012_00_DEM.tif") as dem:
    reproject(rasterio.band(dem, 1), elevation, dst_transform=transform, dst_crs=CRS,
              resampling=Resampling.bilinear)
elevation[~inside] = np.nan

z = {y: (a - np.nanmean(a)) / np.nanstd(a) for y, a in years.items()}
fit = np.stack([z[2016], z[2017]])
mean_z, stability = np.nanmean(fit, axis=0), -np.nanstd(fit, axis=0)
holdout = z[2018]

print(f"{inside.sum():,} cells of {abs(transform.a):.2f} m = {inside.sum() * abs(transform.a * transform.e) / 1e4:.2f} ha")
for year, a in years.items():
    print(year, f"mean {np.nanmean(a):,.0f} kg/ha, sd {np.nanstd(a):,.0f}")
pairs = [(2016, 2017), (2016, 2018), (2017, 2018)]
print("year-to-year correlation:", {f"{a}-{b}": round(float(np.corrcoef(z[a][inside], z[b][inside])[0, 1]), 3) for a, b in pairs})
print(f"elevation {np.nanmin(elevation):.2f} to {np.nanmax(elevation):.2f} m, "
      f"correlation with 2018 yield {np.corrcoef(elevation[inside], years[2018][inside])[0, 1]:.3f}")

41,126 cells of 2.21 m = 20.01 ha
2016 mean 14,320 kg/ha, sd 1,483
2017 mean 14,496 kg/ha, sd 1,923
2018 mean 13,914 kg/ha, sd 2,071
year-to-year correlation: {'2016-2017': 0.396, '2016-2018': 0.321, '2017-2018': 0.327}
elevation -2.93 to -2.45 m, correlation with 2018 yield -0.022

Example 1 — stack, normalise and cluster

import numpy as np
from sklearn.cluster import KMeans
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import silhouette_score

def build_inputs(year_arrays: dict, extra: dict = None):
    z = {}
    for label, a in year_arrays.items():
        a = np.asarray(a, dtype="float32")
        z[label] = (a - np.nanmean(a)) / np.nanstd(a)
    zs = np.stack(list(z.values()))
    out = {"mean_z": np.nanmean(zs, axis=0), "stability": -np.nanstd(zs, axis=0)}
    for k, v in (extra or {}).items():
        v = np.asarray(v, dtype="float32")
        out[k] = (v - np.nanmean(v)) / np.nanstd(v)
    return out

def cluster_zones(inputs, use, k=3, seed=0):
    X = np.column_stack([inputs[c].ravel() for c in use])
    ok = np.isfinite(X).all(axis=1)
    Xs = StandardScaler().fit_transform(X[ok])
    km = KMeans(n_clusters=k, n_init=20, random_state=seed).fit(Xs)

    labels = np.full(X.shape[0], -1, dtype="int16")
    labels[ok] = km.labels_
    ref = inputs["mean_z"].ravel()
    order = np.argsort([np.nanmean(ref[labels == i]) for i in range(k)])
    remap = {old: new for new, old in enumerate(order)}
    labels[ok] = np.vectorize(remap.get)(labels[ok])

    sil = silhouette_score(Xs, km.labels_, sample_size=min(20000, ok.sum()),
                           random_state=seed)
    return labels.reshape(inputs["mean_z"].shape), {"k": k, "silhouette": round(float(sil), 3)}

inputs = build_inputs({2016: years[2016], 2017: years[2017]}, {"elevation": elevation})
labels, meta = cluster_zones(inputs, ["mean_z", "stability"], k=3)
print(meta, "cells per zone:", np.bincount(labels[labels >= 0]).tolist())

{'k': 3, 'silhouette': 0.386} cells per zone: [6363, 12437, 22326]

The re-ordering is there because of what the next block shows. The same data, the same k, five seeds — and how many cells get the same raw k-means label as seed 0:

from sklearn.cluster import KMeans
from sklearn.preprocessing import StandardScaler

X = np.column_stack([inputs[c].ravel() for c in ["mean_z", "stability"]])
ok = np.isfinite(X).all(axis=1)
Xs = StandardScaler().fit_transform(X[ok])
raw = [KMeans(n_clusters=3, n_init=20, random_state=seed).fit(Xs).labels_ for seed in range(5)]
print("raw labels, share equal to seed 0:", [round(float((r == raw[0]).mean()), 3) for r in raw])

ordered = [cluster_zones(inputs, ["mean_z", "stability"], k=3, seed=seed)[0] for seed in range(5)]
print("ordered labels, share equal to seed 0:", [round(float((o == ordered[0])[inside].mean()), 3) for o in ordered])

raw labels, share equal to seed 0: [1.0, 0.0, 0.543, 0.0, 0.302]
ordered labels, share equal to seed 0: [1.0, 1.0, 0.999, 1.0, 0.999]

The fits are the same clustering; only the integers differ, and with them which zone a prescription calls "1". Ordered by mean productivity, the five runs agree on all but 0.1% of cells.

Measured on the Ferrara field, three zones fitted on 2016–17.

Example 2 — choose k on a held-out year

import numpy as np, pandas as pd

def variance_explained(labels, target):
    lab, y = labels.ravel(), np.asarray(target).ravel()
    ok = (lab >= 0) & np.isfinite(y)
    lab, y = lab[ok], y[ok]
    ss_total = ((y - y.mean()) ** 2).sum()
    ss_within = sum(((y[lab == k] - y[lab == k].mean()) ** 2).sum()
                    for k in np.unique(lab))
    return float(1 - ss_within / ss_total)

def choose_k(fit_years, holdout, extra, use, ks=range(2, 7)):
    rows = []
    inputs = build_inputs(fit_years, extra)
    for k in ks:
        labels, meta = cluster_zones(inputs, use, k=k)
        sizes = np.bincount(labels[labels >= 0])
        rows.append({
            "k": k, "silhouette": meta["silhouette"],
            "in_sample": round(variance_explained(labels, inputs["mean_z"]), 3),
            "held_out": round(variance_explained(labels, holdout), 3),
            "smallest_zone_share": round(float(sizes.min() / sizes.sum()), 3),
        })
    df = pd.DataFrame(rows)
    print(df.to_string(index=False))
    return df

for use in (["mean_z", "stability"], ["mean_z", "stability", "elevation"]):
    print(" + ".join(use))
    choose_k({2016: years[2016], 2017: years[2017]}, holdout, {"elevation": elevation}, use)

mean_z + stability
 k  silhouette  in_sample  held_out  smallest_zone_share
 2       0.459      0.342     0.065                0.243
 3       0.386      0.547     0.072                0.155
 4       0.388      0.619     0.078                0.077
 5       0.379      0.706     0.104                0.063
 6       0.365      0.793     0.144                0.009
mean_z + stability + elevation
 k  silhouette  in_sample  held_out  smallest_zone_share
 2       0.336      0.342     0.065                0.241
 3       0.277      0.284     0.051                0.172
 4       0.264      0.473     0.073                0.131
 5       0.262      0.527     0.077                0.070
 6       0.255      0.611     0.078                0.063

Read the two variance columns against each other. In-sample rises steeply with k — 34.2% to 79.3% — because that is what adding clusters does by construction. Held out, on 2018, the zones explain 6.5% to 14.4%: real, small, and still rising at six, so on this field it does not pick k for you. What does is the last column: at six clusters the smallest zone is 0.9% of the field, about 0.2 ha scattered in pieces no machine can follow. Adding elevation lowered the held-out figure at every k above two, because on a field with 48 cm of relief the DEM contributes 30 m noise and nothing else.

Example 3 — smooth into drivable polygons

import numpy as np, geopandas as gpd
from scipy import ndimage
from rasterio.features import shapes
from shapely.geometry import shape

def zones_to_polygons(labels, transform, crs, boom_width_m=24,
                      min_area_ha=1.0, modal_cells=5):
    smooth = ndimage.median_filter(labels.astype("int16"), size=modal_cells)
    polys = [{"zone": int(v), "geometry": shape(g)}
             for g, v in shapes(smooth, mask=smooth >= 0, transform=transform)]
    g = (gpd.GeoDataFrame(polys, crs=crs)
         .dissolve("zone", as_index=False)
         .explode(index_parts=False))
    g["area_ha"] = g.area / 1e4

    small = g[g.area_ha < min_area_ha]
    if len(small):
        print(f"dropping {len(small)} fragments below {min_area_ha} ha "
              f"({small.area_ha.sum():.2f} ha total)")
    g = g[g.area_ha >= min_area_ha]
    g["geometry"] = g.geometry.simplify(boom_width_m / 2).buffer(0)
    out = g.dissolve("zone", as_index=False)
    out["area_ha"] = out.area / 1e4
    print(out[["zone", "area_ha"]].round(2).to_string(index=False))
    return out

import shapely

polygons = zones_to_polygons(labels, transform, CRS)
field_ha = inside.sum() * abs(transform.a * transform.e) / 1e4
print(f"field {field_ha:.2f} ha, zones cover {polygons.union_all().area / 1e4:.2f} ha; "
      f"zones present: {sorted(polygons['zone'])}")

dropping 253 fragments below 1.0 ha (5.27 ha total)
 zone  area_ha
    1     5.54
    2     9.79
field 20.01 ha, zones cover 15.27 ha; zones present: [1, 2]

That is not a zoning of the field. The low zone on this field exists only as scattered patches, none of them a hectare, so dropping fragments below the minimum deletes zone 0 altogether and leaves 4.74 ha of the field uncovered by any zone, and so by any rate. Three changes fix it: a true modal filter (a median of zone numbers is not the most common zone, and at the field edge it mixes in the -1 outside), rasterio.features.sieve to merge small patches into their neighbours instead of dropping them, and shapely.coverage_simplify, which simplifies the shared boundaries once rather than each polygon on its own:

from rasterio.features import sieve


def modal(labels, size=5):
    def most_common(window):
        window = window[window >= 0].astype(int)
        return np.bincount(window).argmax() if len(window) else -1
    return ndimage.generic_filter(labels.astype("int16"), most_common, size=size)


def zones_to_coverage(labels, transform, crs, boom_width_m=24, min_area_ha=1.0, modal_cells=5):
    field = labels >= 0
    smooth = modal(labels, modal_cells)
    cells = int(round(min_area_ha * 1e4 / abs(transform.a * transform.e)))
    merged = sieve(smooth + 1, size=cells, mask=field).astype("int16") - 1
    merged[~field] = -1
    g = gpd.GeoDataFrame([{"zone": int(v), "geometry": shape(p)}
                          for p, v in shapes(merged, mask=field, transform=transform)], crs=crs)
    g = g.dissolve("zone", as_index=False)
    g["geometry"] = shapely.coverage_simplify(g.geometry.values, boom_width_m / 2)
    g["area_ha"] = (g.area / 1e4).round(2)
    return g


def overlap_m2(g):
    parts = list(g.geometry)
    return sum(a.intersection(b).area for i, a in enumerate(parts) for b in parts[i + 1:])


smooth = ndimage.median_filter(labels.astype("int16"), size=5)
print("median filter: cells outside the field given a zone:", int((smooth[~inside] >= 0).sum()),
      "| field cells set to -1:", int((smooth[inside] < 0).sum()))
for min_ha in (0.25, 1.0):
    g = zones_to_coverage(labels, transform, CRS, min_area_ha=min_ha)
    each = g.assign(geometry=g.geometry.simplify(12))
    print(f"min {min_ha} ha: {dict(zip(g['zone'], g['area_ha']))}, parts {[len(getattr(p, 'geoms', [p])) for p in g.geometry]}, "
          f"covers {g.union_all().area / 1e4:.2f} ha, overlap {overlap_m2(g):.0f} m² "
          f"| simplified one by one: overlap {overlap_m2(each):,.0f} m², covers {each.union_all().area / 1e4:.2f} ha")

median filter: cells outside the field given a zone: 6 | field cells set to -1: 5
min 0.25 ha: {0: 0.74, 1: 6.95, 2: 12.33}, parts [2, 2, 2], covers 20.02 ha, overlap 0 m² | simplified one by one: overlap 711 m², covers 19.32 ha
min 1.0 ha: {1: 7.23, 2: 12.78}, parts [1, 1], covers 20.01 ha, overlap 0 m² | simplified one by one: overlap 486 m², covers 19.50 ha

Now the zones tile the field with no gaps and no overlaps. The price of a 1 ha minimum is still visible, though: it merges the low zone into its neighbours entirely, and the field becomes a two-zone field. A 0.25 ha minimum keeps a 0.74 ha low zone in two pieces. Which is right is an agronomic decision, but it should be made on purpose, not by a filter. Simplifying at half the boom width is the right order of magnitude: detail finer than the implement cannot be applied, so keeping it only enlarges the file and gives the controller more boundaries to cross.

Measured in Example 3. The minimum area decides whether the low zone survives at all.

Why k-means, and ordered labels

Why k-means and not something cleverer

The inputs are a handful of continuous, roughly elliptical variables and the output has to be a small number of compact groups. K-means does that, it is fast, it is deterministic with a fixed seed and enough restarts, and every agronomist has seen it before. Fuzzy c-means is the traditional alternative and gives a membership rather than a label, which is useful if you want a continuous rate surface rather than discrete zones.

Why label ordering matters more than it sounds

Cluster labels are arbitrary integers. Two runs with different seeds, different input order or a different scikit-learn version can swap zone 0 and zone 2. A prescription file written against one run and applied after another then puts the high rate where the low rate belongs. Sorting by a meaningful quantity — mean productivity — makes the label a property of the field rather than of the run.

Why in-sample variance explained is not evidence

Adding clusters always reduces within-cluster variance, so in-sample explanatory power rises with k by construction, all the way to one cluster per pixel. Only a year that was not used in the fitting tests whether the pattern is persistent, and that is the number worth quoting.

Why the smallest zone's size is a constraint

A zone of half a hectare scattered in twenty fragments cannot be managed: the machine cannot switch rates fast enough, the agronomy is not different at that scale, and the sampling needed to characterise it costs more than the response is worth. Setting a minimum area before clustering — or enforcing it afterwards — is part of the specification, not a tidy-up.

A prescription written against one run and applied after another puts the high rate in the wrong place.

One grid, scaling and n_init

Put every layer on one grid before stacking, with the same mask.
Scale before clustering. K-means is sensitive to units.
Set n_init high. Ten is the default; twenty is cheap insurance.
Headlands often form their own zone. Decide whether that is useful.
Fuzzy membership suits a continuous rate surface better than hard labels.
Check the zone means differ agronomically, not just statistically.
Re-fit as years accumulate, and record which years each version used.
Keep the fitted model if you will zone other fields the same way.

Internal links

Management zones explained — what to build zones from and why
How to write a variable-rate prescription map — turning zones into rates
Yield monitor data explained: what the numbers really are — cleaning the input
Spatial clustering explained — the algorithms
How to lay out a soil sampling grid in Python — sampling by zone
How to analyse an on-farm strip trial in Python — finding the response
How to calculate slope and aspect in Python — terrain inputs
DBSCAN returns one cluster or all noise — when a density method is the wrong tool

FAQ

How do I build management zones in Python?

Stack normalised multi-year layers, cluster with k-means, order the labels by productivity, smooth with a modal filter and a minimum area, and validate on a year you held out.

How many zones should I make?

Sweep two to six and pick on held-out variance explained and the smallest zone's size. Two or three usually wins.

Why do my zone numbers change between runs?

Because k-means labels are arbitrary. Sort them by mean productivity so zone 0 is always the poorest.

Which algorithm should I use?

K-means for discrete zones, fuzzy c-means if you want a continuous membership that becomes a smooth rate surface.

How do I know the zones are worth using?

Compute the variance they explain in a year that was not used to fit them. In-sample variance explained always rises with the number of clusters and proves nothing.

How much should I simplify the polygons?

To about half the implement width. Finer detail cannot be applied and only makes the file harder for the controller.