K-Means Clustering for Yield Zone Delineation

TL;DR: Stack a kriged yield raster with optional EC and elevation bands into an (n_pixels × n_features) matrix, standardize the columns with StandardScaler, pick k between 3 and 5 with silhouette and elbow scores, fit scikit-learn KMeans, reshape the labels back onto the raster grid, sieve out speckle, and polygonize the result into controller-ready management zones.

Why K-Means for Management Zones

A kriged yield surface is continuous — every pixel carries a floating-point yield estimate — but a variable-rate controller cannot vary seed or fertilizer rate per pixel. It needs a small number of discrete zones, each mapped to one prescription rate. Management zone delineation is the classification step that converts the smooth interpolated surface produced by interpolating sparse yield monitor data with kriging into 3–5 spatially coherent management units.

K-Means is the workhorse here because it is fast, unsupervised, and multivariate: it can fuse yield with co-registered layers such as apparent soil electrical conductivity (EC), elevation, or a multi-year yield stability index into a single zonal map. The failure mode that ruins most first attempts is skipping feature standardization — yield in bu/ac has roughly 40× the numeric spread of EC in mS/m, so the algorithm silently clusters on yield alone. Get the scaling wrong and the EC layer might as well not be in the stack, which defeats the entire point of a multivariate zonation and can leave 30–40% of the field assigned to the wrong input rate.

The diagram below shows the full transform from stacked rasters to polygon zones.

K-Means management zone delineation pipeline Five-stage flow: stack yield, EC and elevation rasters; standardize features; choose k with silhouette; fit KMeans and reshape labels; sieve and polygonize zones. Stack layers yield · EC elevation Standardize StandardScaler Choose k silhouette · elbow Fit KMeans reshape labels Sieve + polygonize zone polygons 1. Assemble 2. Scale 3. Select 4. Cluster 5. Vectorize

This guide is part of Management Zone Classification Algorithms — see there for the full pipeline context, including Gaussian Mixture alternatives and multi-year stability zoning.

Prerequisites

These add scikit-learn on top of the raster stack used elsewhere in this section:

TEXT 
scikit-learn==1.5.0
rasterio==1.3.10
numpy==1.26.4
shapely==2.0.4

Install with:

BASH 
pip install scikit-learn==1.5.0 rasterio==1.3.10 numpy==1.26.4 shapely==2.0.4

Input requirements:

  • A kriged yield GeoTIFF (Band 1 = predicted yield), already in a projected metric CRS such as EPSG:32615 (UTM zone 15 N).
  • Optional co-variate rasters (soil EC, elevation) resampled to the exact same grid, transform, and shape as the yield raster. Misaligned grids produce meaningless per-pixel feature vectors.
  • A consistent nodata value on every input (this script treats numpy.nan as nodata).

Step-by-Step

Choosing k with a number rather than a preference Bars of silhouette score for four cluster counts on the same field: two zones score highest but are coarse, three and four remain defensible, and six drops sharply because the zones stop being distinct from one another. k = 2 silhouette 0.61 — coarse but stable k = 3 silhouette 0.58 — the usual choice k = 4 silhouette 0.49 — still defensible k = 6 silhouette 0.31 — zones stop being distinct Silhouette score, higher is better. The agronomic constraint usually binds first: every zone must be at least twice the implement width.

Step 1 — Stack the feature rasters

Read each aligned raster into a 2-D array and stack them along a new axis. The pixels where any feature is nodata must be dropped before clustering, so build a shared valid-data mask. Every layer must share the same transform and shape; assert it rather than trust it.

Step 2 — Standardize the feature columns

Reshape the stack to an (n_valid_pixels × n_features) matrix and fit StandardScaler. This centres each feature at mean 0 and unit variance so yield, EC, and elevation contribute equally to the Euclidean distance K-Means minimizes. Keep the fitted scaler — you never need to invert it here, but persisting it makes zone maps from different fields comparable.

Step 3 — Choose k with silhouette and elbow

Sweep k from 2 to 7. For each, record the inertia (within-cluster sum of squares, the elbow signal) and the silhouette score (cluster separation, higher is better). Pick the k with the best silhouette, capped at 5 for agronomic tractability. On a subsampled pixel set the silhouette computation stays cheap; scoring the full field is unnecessary.

Step 4 — Fit K-Means and rebuild the raster

Fit KMeans with the chosen k, a fixed random_state, and n_init=10 so the result is reproducible and not stuck in a bad local optimum. Scatter the 1-D labels back into a full-size integer raster using the valid-data mask, filling nodata pixels with an out-of-range sentinel (-1).

Step 5 — Sieve speckle and polygonize

K-Means ignores spatial adjacency, so single noisy pixels get their own label. Remove clumps smaller than the equipment working width with rasterio.features.sieve, then vectorize the cleaned label raster with rasterio.features.shapes. The complete, runnable script:

PYTHON 
import numpy as np
import rasterio
from rasterio.features import sieve, shapes
from shapely.geometry import shape
from sklearn.preprocessing import StandardScaler
from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_score

# ── 1. Stack aligned feature rasters ───────────────────────────────────────
FEATURE_PATHS = {
    "yield": "yield_kriged.tif",   # Band 1 = predicted yield (bu/ac)
    "ec":    "soil_ec.tif",        # apparent EC (mS/m), same grid
    "elev":  "elevation.tif",      # metres, same grid
}

layers, profile = [], None
for name, path in FEATURE_PATHS.items():
    with rasterio.open(path) as src:
        arr = src.read(1, masked=True).filled(np.nan).astype(np.float64)
        if profile is None:
            profile, ref_transform, ref_shape = src.profile, src.transform, arr.shape
        else:
            assert arr.shape == ref_shape, f"{name}: shape {arr.shape} != {ref_shape}"
            assert src.transform == ref_transform, f"{name}: transform mismatch"
        layers.append(arr)

stack = np.stack(layers, axis=-1)                 # (rows, cols, n_features)
rows, cols, n_features = stack.shape
assert profile["crs"].is_projected, "Inputs must be in a projected metric CRS"

# Shared valid-data mask: a pixel is usable only if every feature is finite
valid = np.all(np.isfinite(stack), axis=-1)       # (rows, cols) bool
X = stack[valid]                                  # (n_valid, n_features)
assert X.shape[0] > 0, "No pixels valid across all feature layers"

# ── 2. Standardize so each layer contributes equal variance ────────────────
X_scaled = StandardScaler().fit_transform(X)

# ── 3. Choose k by silhouette (elbow via inertia for cross-check) ──────────
rng = np.random.default_rng(42)
sample_idx = rng.choice(X_scaled.shape[0], size=min(5000, X_scaled.shape[0]),
                        replace=False)
scores = {}
for k in range(2, 8):
    km = KMeans(n_clusters=k, n_init=10, random_state=42).fit(X_scaled)
    sil = silhouette_score(X_scaled[sample_idx], km.labels_[sample_idx])
    scores[k] = (sil, km.inertia_)
    print(f"k={k}  silhouette={sil:.3f}  inertia={km.inertia_:,.0f}")

best_k = min(max(scores, key=lambda k: scores[k][0]), 5)   # cap at 5 zones
print(f"Selected k = {best_k}")

# ── 4. Fit final model and scatter labels back onto the raster ─────────────
final = KMeans(n_clusters=best_k, n_init=10, random_state=42).fit(X_scaled)
# Order labels by mean yield so zone 0 = lowest-yielding, k-1 = highest
yield_col = list(FEATURE_PATHS).index("yield")
cluster_mean_yield = [X[final.labels_ == c, yield_col].mean() for c in range(best_k)]
rank = np.argsort(np.argsort(cluster_mean_yield))          # old label -> rank
ordered = rank[final.labels_].astype(np.int32)

zones = np.full((rows, cols), -1, dtype=np.int32)          # -1 = nodata
zones[valid] = ordered

# ── 5. Sieve speckle, then polygonize ──────────────────────────────────────
# Drop clumps smaller than ~0.25 ha at 10 m pixels (25 pixels)
zones_clean = sieve(zones, size=25, connectivity=8)

out_profile = profile.copy()
out_profile.update(dtype="int32", count=1, nodata=-1)
with rasterio.open("management_zones.tif", "w", **out_profile) as dst:
    dst.write(zones_clean, 1)

polys = []
for geom, value in shapes(zones_clean, mask=(zones_clean >= 0),
                          transform=ref_transform):
    if value >= 0:
        polys.append({"zone": int(value), "geometry": shape(geom)})
print(f"Polygonized into {len(polys)} zone parts across {best_k} zones")

# ── Inline verification ────────────────────────────────────────────────────
labels, counts = np.unique(zones_clean[zones_clean >= 0], return_counts=True)
print(dict(zip(labels.tolist(), counts.tolist())))
assert len(labels) == best_k, f"Expected {best_k} zones, got {len(labels)}"
assert counts.min() > 25, "A zone is smaller than the sieve threshold — re-check sieve"
assert set(labels.tolist()) == set(range(best_k)), "Zone labels are not contiguous 0..k-1"
print("Zone delineation verified.")

The verification block asserts three things that catch the common failures: the number of distinct labels equals best_k (no zone was fully sieved away), no surviving zone is below the speckle threshold, and the labels form a contiguous 0..k-1 range after the yield-ordering step so downstream prescription tables can index them directly.

Gotchas & Edge Cases

  • Unscaled features silently collapse the zonation onto the highest-range layer. If you skip StandardScaler, yield (0–250 bu/ac) dominates EC (0–60 mS/m) and the clustering ignores the soil signal entirely. Always scale, and sanity-check by confirming each K-Means group differs in EC as well as yield.
  • Salt-and-pepper zones are a spatial artefact, not a clustering error. K-Means has no concept of adjacency, so noisy pixels scatter. The sieve(size=25, connectivity=8) call merges sub-quarter-hectare clumps into their neighbours; raise size to match your applicator’s working width if zones still look grainy.
  • k too high produces zones no controller can act on. A silhouette that keeps rising past k=6 usually reflects noise structure, not agronomy. Cap k at 5 and prefer the coarser map — the min(..., 5) guard enforces this.
  • Misaligned co-variate grids poison every feature vector. Resample EC and elevation to the exact yield-raster transform first; the per-layer transform and shape asserts fail loudly rather than clustering garbage.
Why standardisation is not optional before clustering Two panels comparing clustering on raw and standardised layers. With raw units, yield in kilograms per hectare spans thousands while conductivity spans tens, so the distance metric is dominated by one layer; standardising lets every layer contribute comparably. Unstandardised inputs Yield in kg/ha, EC in mS/m, elevation in m Yield spans thousands; EC spans tens. Distance is dominated by the largest range. The clustering is a yield map with extra steps. Standardised inputs Each layer centred and scaled Every layer contributes comparably. PCA can then drop what is redundant. Zones reflect all the evidence.

Frequently Asked Questions

Why must I scale features before running K-Means on yield and EC layers?

K-Means minimizes squared Euclidean distance, so a feature with a large numeric range dominates the clustering. Yield in bushels per acre spans 0 to 250 while soil EC spans 0 to 60 mS/m, so without standardization the zones would track yield alone and ignore EC. Fit StandardScaler on the stacked feature columns so each layer contributes equal variance.

How many management zones should I create?

Three to five zones is the practical range for most fields because variable-rate controllers and input budgets cannot act on finer distinctions. Use the silhouette score to pick k objectively, but cap it at five even when a higher k scores marginally better, since agronomically meaningless zones increase application complexity without a yield return.

Why do my zones look like scattered speckle instead of contiguous blocks?

K-Means clusters in feature space and ignores spatial adjacency, so noisy pixels produce a salt-and-pepper pattern that no controller can follow. Apply a majority filter or a connected-component sieve to remove clumps smaller than the equipment working width before polygonizing, which merges isolated pixels into their surrounding zone.

Parent section Management Zone Classification for Precision Agriculture