A crop classifier confuses two crops every year
On this page
Problem statement
The overall accuracy is 84.5% and a single pair of classes accounts for most of the errors. On the Flevoland parcels used on this page, 41 of the 50 red-onion fields are predicted as yellow onion; elsewhere it is winter wheat as winter barley, or a cover crop as grassland — and adding more trees, more features or a different classifier changes nothing.
That is the signature of a problem the model cannot solve: the two classes are genuinely similar in the features provided. Either the feature that separates them is missing, or the dates that would show it were cloudy, or the labels themselves are unreliable for one of the pair.
Quick answer
Read the confusion matrix, then look at the two classes' mean profiles:
import numpy as np, pandas as pd
from sklearn.metrics import confusion_matrix
print(f"overall accuracy {(pred == y).mean():.3f} on {len(y)} parcels")
labels = sorted(np.unique(y))
cm = pd.DataFrame(confusion_matrix(y, pred, labels=labels), index=labels, columns=labels)
off = cm.values.copy()
np.fill_diagonal(off, 0)
i, j = np.unravel_index(off.argmax(), off.shape)
print(f"largest confusion: {labels[i]} → {labels[j]} "
f"({off[i, j]} of {cm.values[i].sum()}, {off[i, j] / cm.values[i].sum():.1%})")
for a, b in [(labels[i], labels[j])]:
pa = X[y == a].mean(axis=0)
pb = X[y == b].mean(axis=0)
sep = np.abs(pa - pb) / np.sqrt((X[y == a].var(axis=0) + X[y == b].var(axis=0)) / 2)
print("best separating feature:", feature_names[int(np.nanargmax(sep))],
f"(effect size {np.nanmax(sep):.2f})")overall accuracy 0.845 on 890 parcels
largest confusion: red onion → yellow onion (41 of 50, 82.0%)
best separating feature: ndvi_07-02 (effect size 0.43)
The two onions peak at 0.43, on NDVI in early July: nothing in eleven dates of four indices tells them apart. An effect size below about 0.5 on every feature means the pair is not separable with what you have. Above 1.5 on one date says the class means differ, and that the information may be there for a feature to use — a different problem with a different fix. It does not promise the fix will work: Example 3 shows a pair at 1.58 that a targeted feature barely helps, because an effect size compares averages and ignores how many of the larger class overlap the smaller one.
One pair or general noise, mean profiles, critical dates
1. Confirm it is one pair, not general noise
A confusion matrix with one large off-diagonal cell is a different problem from one with errors spread everywhere. The first is a feature problem; the second is usually a labelling or a validation problem.
2. Plot the two classes' mean temporal profiles
If they overlap everywhere, no classifier can separate them. If they diverge on two dates in April, the model has the information and the features are burying it.
3. Check whether the separating date was cloudy
Two cereals differ most at heading and at senescence, which are narrow windows. If the nearest clear observation is three weeks away, the interpolated value at the critical date is an average of the two sides of the gap. Count the clear observations near each class's divergence.
4. Check the index for saturation over the critical window
Two cereals at heading are both at full canopy, where NDVI is flat. On a June 2025 scene over this same part of Flevoland, 50.7% of vegetation pixels were above NDVI 0.8, and across 379 dense arable parcels NDRE spread 0.231 between fields where NDVI spread 0.123.
5. Add a feature that targets the difference
A red-edge index at heading, a shortwave index at senescence, the slope of the falling limb, or the date each class reached a thermal threshold. One targeted feature usually beats a hundred generic ones.
6. Check the labels
Declared crop is the intended crop. A field declared as barley and sown with wheat is a permanently mislabelled training sample, and a pair of similar crops is exactly where declaration errors concentrate.
7. If it still cannot be separated, merge the classes
A "winter cereal" class that is 96% accurate is more useful than wheat and barley classes that are each 70% accurate. Merging is a legitimate result, and it should be reported as one.
A separability report for two crops
The data used on this page
Every 2025 parcel of nine crops, of at least a hectare, in the Dutch parcel register BRP Gewaspercelen within 5.55–5.75° E, 52.48–52.56° N in Flevoland — the box the crop phenology and NDVI saturation guides use. The crops include three pairs that could plausibly be confused: yellow and red onions, ware and seed potatoes, temporary and permanent grassland. The features are four indices — NDVI, NDRE, NDWI and a shortwave-infrared ratio — as each parcel's median over clear pixels (vegetation or bare soil in the scene classification) on eleven nearly cloud-free Sentinel-2 dates from March to September 2025, read at 20 m from Earth Search, with the parcel shrunk 20 m to keep edge pixels out. A date a parcel was clouded on is filled from its neighbours in time. The classifier is a 500-tree random forest scored with 5-fold cross-validation; each parcel is one sample.
import os
import geopandas as gpd
import numpy as np
import pandas as pd
import rasterio
import requests
from rasterio.enums import Resampling
from rasterio.features import rasterize
from rasterio.windows import bounds as window_bounds, from_bounds
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import StratifiedKFold, cross_val_predict
os.environ.update(AWS_NO_SIGN_REQUEST="YES", GDAL_DISABLE_READDIR_ON_OPEN="EMPTY_DIR")
CROPS = {"Tarwe, winter-": "winter wheat", "Aardappelen, consumptie": "ware potato",
"Aardappelen, poot NAK": "seed potato", "Uien, gele zaai-": "yellow onion",
"Uien, rode zaai-": "red onion", "Bieten, suiker-": "sugar beet", "Winterpeen, productie": "carrot",
"Grasland, tijdelijk": "temporary grass", "Grasland, blijvend": "permanent grass"}
SCENES = ["S2B_31UFU_20250307_0_L2A", "S2A_31UFU_20250403_0_L2A", "S2C_31UFU_20250501_0_L2A",
"S2A_31UFU_20250513_0_L2A", "S2B_31UFU_20250519_0_L2A", "S2A_31UFU_20250612_0_L2A",
"S2C_31UFU_20250630_0_L2A", "S2A_31UFU_20250702_0_L2A", "S2A_31UFU_20250814_1_L2A",
"S2C_31UFU_20250819_0_L2A", "S2B_31UFU_20250906_0_L2A"]
if not os.path.exists("parcel_series.csv"):
url, params, feats = ("https://api.pdok.nl/rvo/gewaspercelen/ogc/v1/collections/brpgewas/items",
{"f": "json", "bbox": "5.55,52.48,5.75,52.56", "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
brp = gpd.GeoDataFrame.from_features(feats, crs=4326)
brp["ha"] = brp.to_crs(28992).area / 1e4
parcels = brp[brp["gewas"].isin(CROPS) & (brp["ha"] >= 1)].reset_index(drop=True)
parcels["crop"], parcels["pid"] = parcels["gewas"].map(CROPS), np.arange(1, len(parcels) + 1)
items = {i["id"]: i for i in requests.post("https://earth-search.aws.element84.com/v1/search", timeout=120,
json={"collections": ["sentinel-2-l2a"], "ids": SCENES, "limit": 50}).json()["features"]}
utm = parcels.to_crs(32631)
box = utm.total_bounds + np.array([-40, -40, 40, 40])
rows = []
for scene in SCENES:
a = items[scene]["assets"]
with rasterio.open(a["swir16"]["href"]) as src: # the 20 m grid
win = from_bounds(*box, transform=src.transform).round_offsets().round_lengths()
shape, grid, extent = (int(win.height), int(win.width)), src.window_transform(win), window_bounds(win, src.transform)
swir = src.read(1, window=win).astype("float32")
band = {}
for k in ("red", "green", "nir", "rededge1"):
with rasterio.open(a[k]["href"]) as src:
band[k] = src.read(1, window=from_bounds(*extent, transform=src.transform), out_shape=shape,
resampling=Resampling.average).astype("float32")
with rasterio.open(a["scl"]["href"]) as src:
clear = np.isin(src.read(1, window=win), [4, 5])
ids = rasterize(zip(utm.geometry.buffer(-20), utm["pid"]), out_shape=shape, transform=grid, fill=0, dtype="int32")
ok = (ids > 0) & clear
red, green, nir, re1 = band["red"], band["green"], band["nir"], band["rededge1"]
with np.errstate(all="ignore"):
idx = {"ndvi": (nir - red) / (nir + red), "ndre": (nir - re1) / (nir + re1),
"ndwi": (green - nir) / (green + nir), "swir": (swir - nir) / (swir + nir)}
frame = pd.DataFrame({"pid": ids[ok], **{k: v[ok] for k, v in idx.items()}})
med = frame.groupby("pid").median()
med["date"] = items[scene]["properties"]["datetime"][:10]
rows.append(med.reset_index())
pd.concat(rows).to_csv("parcel_series.csv", index=False)
parcels.drop(columns="geometry").to_csv("parcels.csv", index=False)
series = pd.read_csv("parcel_series.csv")
wide = series.pivot(index="pid", columns="date", values=["ndvi", "ndre", "ndwi", "swir"])
feature_names = [f"{index}_{date[5:]}" for index, date in wide.columns]
feature_dates = [date for _, date in wide.columns]
wide.columns = feature_names
for index in ("ndvi", "ndre", "ndwi", "swir"): # fill clouded dates along time
cols = [c for c in feature_names if c.startswith(index)]
wide[cols] = wide[cols].interpolate(axis=1, limit_direction="both")
data = pd.read_csv("parcels.csv").set_index("pid").join(wide, how="inner")
X, y = data[feature_names].to_numpy(), data["crop"].to_numpy()
pred = cross_val_predict(RandomForestClassifier(500, n_jobs=4, random_state=0), X, y,
cv=StratifiedKFold(5, shuffle=True, random_state=0))
print(len(y), "parcels,", X.shape[1], "features;", data["crop"].value_counts().to_dict())890 parcels, 44 features; {'ware potato': 202, 'winter wheat': 181, 'yellow onion': 121, 'sugar beet': 108, 'temporary grass': 90, 'carrot': 67, 'red onion': 50, 'seed potato': 46, 'permanent grass': 25}
Example 1 — the separability report for a pair
import numpy as np, pandas as pd
def separability(X, y, a, b, feature_names, dates=None):
A, B = X[y == a], X[y == b]
mu_a, mu_b = A.mean(axis=0), B.mean(axis=0)
pooled = np.sqrt((A.var(axis=0) + B.var(axis=0)) / 2)
effect = np.abs(mu_a - mu_b) / np.where(pooled == 0, np.nan, pooled)
df = pd.DataFrame({"feature": feature_names,
f"mean_{a}": mu_a.round(3), f"mean_{b}": mu_b.round(3),
"effect_size": effect.round(3)})
if dates is not None:
df["date"] = dates
df = df.sort_values("effect_size", ascending=False)
print(df.head(8).to_string(index=False))
print(f"\nn({a}) = {len(A)}, n({b}) = {len(B)}; "
f"max effect size {np.nanmax(effect):.2f}")
if np.nanmax(effect) < 0.5:
print("! these classes are not separable with these features")
return dfonions = separability(X, y, "red onion", "yellow onion", feature_names, feature_dates)
print()
potatoes = separability(X, y, "seed potato", "ware potato", feature_names, feature_dates) feature mean_red onion mean_yellow onion effect_size date
ndvi_07-02 0.592 0.530 0.426 2025-07-02
swir_07-02 -0.318 -0.256 0.413 2025-07-02
ndre_07-02 0.430 0.387 0.370 2025-07-02
ndwi_07-02 -0.573 -0.534 0.357 2025-07-02
swir_06-30 -0.359 -0.303 0.339 2025-06-30
ndvi_06-30 0.627 0.577 0.319 2025-06-30
ndwi_06-30 -0.604 -0.572 0.273 2025-06-30
ndre_06-30 0.437 0.409 0.240 2025-06-30
n(red onion) = 50, n(yellow onion) = 121; max effect size 0.43
! these classes are not separable with these features
feature mean_seed potato mean_ware potato effect_size date
ndvi_08-19 0.205 0.564 1.579 2025-08-19
swir_08-19 0.137 -0.158 1.567 2025-08-19
ndwi_08-19 -0.297 -0.524 1.556 2025-08-19
swir_08-14 0.118 -0.198 1.496 2025-08-14
ndvi_08-14 0.233 0.608 1.492 2025-08-14
ndre_08-19 0.121 0.323 1.487 2025-08-19
ndwi_08-14 -0.312 -0.550 1.438 2025-08-14
ndre_08-14 0.127 0.367 1.365 2025-08-14
n(seed potato) = 46, n(ware potato) = 202; max effect size 1.58
The two pairs are opposite cases. The onions never separate: the same species, sown and lifted on the same calendar, in fields that look alike from space on every date. The potatoes separate strongly in August — seed potatoes have their haulm destroyed early so the tubers stay small, while ware potatoes are still green — so their confusion, if any, is a feature or a date problem, not a data problem.
Example 2 — is the critical date observed?
import numpy as np, pandas as pd
def observation_support(sep_df, observed_dates, window_days=7, top_n=3):
obs = pd.DatetimeIndex(observed_dates)
rows = []
for _, r in sep_df.head(top_n).iterrows():
d = pd.Timestamp(r["date"])
near = np.abs((obs - d).days) <= window_days
gaps = np.diff(np.sort(obs)).astype("timedelta64[D]").astype(int)
rows.append({"date": str(d.date()), "effect_size": r["effect_size"],
"clear_obs_within_7d": int(near.sum()),
"days_to_nearest": int(np.abs((obs - d).days).min()),
"largest_season_gap": int(gaps.max())})
df = pd.DataFrame(rows)
print(df.to_string(index=False))
if (df["clear_obs_within_7d"] == 0).any():
print("! the most separating date has no nearby clear observation")
return dfobserved = sorted(set(feature_dates))
support = observation_support(potatoes, observed) date effect_size clear_obs_within_7d days_to_nearest largest_season_gap
2025-08-19 1.579 2 0 43
2025-08-19 1.567 2 0 43
2025-08-19 1.556 2 0 43
A separating date with no clear observation within a week is an argument for adding a sensor — Landsat, or a radar series — rather than for tuning the model.
Example 3 — targeted features, and merging as a fallback
import numpy as np, pandas as pd
from sklearn.metrics import classification_report
def add_targeted_features(series, dates, window, index="ndre"):
"""Mean, slope and range of one index within a specific window."""
d = pd.DatetimeIndex(dates)
sel = (d >= pd.Timestamp(window[0])) & (d <= pd.Timestamp(window[1]))
sub = series[:, sel]
return {
f"{index}_mean_{window[0]}": np.nanmean(sub, axis=1),
f"{index}_slope_{window[0]}": np.nanmean(np.diff(sub, axis=1), axis=1),
f"{index}_range_{window[0]}": np.nanmax(sub, axis=1) - np.nanmin(sub, axis=1),
}
def merge_classes(y, pred, groups):
"""groups: {'winter cereal': ['winter wheat', 'winter barley'], ...}"""
remap = {old: new for new, olds in groups.items() for old in olds}
ym = np.array([remap.get(v, v) for v in y])
pm = np.array([remap.get(v, v) for v in pred])
print(classification_report(ym, pm, zero_division=0))
return ym, pmfrom sklearn.metrics import recall_score
for pair in (["red onion", "yellow onion"], ["seed potato", "ware potato"]):
r = recall_score(y, pred, labels=pair, average=None)
print(pair, "recall", r.round(2).tolist())
ym, pm = merge_classes(y, pred, {"onion": ["red onion", "yellow onion"]})['red onion', 'yellow onion'] recall [0.16, 0.85]
['seed potato', 'ware potato'] recall [0.43, 0.92]
precision recall f1-score support
carrot 0.91 0.90 0.90 67
onion 0.98 0.98 0.98 171
permanent grass 0.79 0.60 0.68 25
seed potato 0.59 0.43 0.50 46
sugar beet 1.00 0.96 0.98 108
temporary grass 0.86 0.87 0.86 90
ware potato 0.83 0.92 0.87 202
winter wheat 0.98 0.99 0.98 181
accuracy 0.91 890
macro avg 0.87 0.83 0.84 890
weighted avg 0.90 0.91 0.90 890
The potato pair is the other branch. They separate on 19 August at an effect size of 1.58, so a targeted feature — the NDVI mean, slope and range over July and August — should be the fix, and class weights should stop the 202 ware parcels outvoting the 46 seed parcels:
from sklearn.model_selection import StratifiedKFold, cross_val_predict
ndvi_cols = [i for i, n in enumerate(feature_names) if n.startswith("ndvi")]
extra = add_targeted_features(X[:, ndvi_cols], [feature_dates[i] for i in ndvi_cols],
("2025-07-02", "2025-08-19"), index="ndvi")
X_targeted = np.column_stack([X, *extra.values()])
for label, features, weight in (("as before", X, None), ("balanced weights", X, "balanced_subsample"),
("+ targeted features", X_targeted, None)):
p = cross_val_predict(RandomForestClassifier(500, n_jobs=4, random_state=0, class_weight=weight), features, y,
cv=StratifiedKFold(5, shuffle=True, random_state=0))
print(f"{label:22} seed potato recall {recall_score(y, p, labels=['seed potato'], average=None)[0]:.2f}")
aug19 = pd.Series(X[:, feature_names.index("ndvi_08-19")], index=y)
for crop in ("seed potato", "ware potato"):
v = aug19[crop]
print(f"{crop:12} NDVI < 0.4 on 19 August: {(v < 0.4).sum()} of {len(v)} ({(v < 0.4).mean():.0%})")as before seed potato recall 0.43
balanced weights seed potato recall 0.46
+ targeted features seed potato recall 0.46
seed potato NDVI < 0.4 on 19 August: 42 of 46 (91%)
ware potato NDVI < 0.4 on 19 August: 64 of 202 (32%)
Neither helps much, and the last two lines show why. Seed potatoes have their haulm destroyed early, and 91% of them are brown by 19 August — but so are 32% of the ware-potato fields, harvested early for the new-potato market. Sixty-four brown ware parcels against forty-two brown seed parcels: the feature is real, the averages differ enormously, and among the parcels it picks out, the larger class is still the majority. The fix here is an earlier date — before either kind is killed off — or information that is not optical at all, such as the certification register seed potatoes are grown under.
Reporting both the split and the merged result is the honest way to present it: the reader can see what was gained by merging and decide whether the distinction mattered to them.
Why more model capacity does not help
Why more model capacity does not help
A classifier can only separate classes that are separated in the feature space. If two crops' feature distributions overlap almost completely, the Bayes-optimal error rate is high, and a random forest, a gradient booster and a neural network will all converge on roughly the same accuracy. The flat response to model changes is itself the diagnosis.
Why two cereals are the hardest pair
Winter wheat and winter barley are sown within weeks of each other, reach full canopy at the same time and senesce within a fortnight of each other. The differences that exist — barley heads earlier and senesces earlier — are expressed in narrow windows, which optical satellites may or may not observe, and both crops are at full canopy for the period when observations are most reliable.
Why declared labels concentrate errors on similar pairs
A farmer who intended barley and sowed wheat rarely amends the declaration, and a declaration covering a parcel with two crops carries only one. Both errors fall disproportionately on similar crops, because those are the ones substituted for each other. A model trained on such labels learns the confusion.
Why merging is a result rather than a retreat
A classification is used for something: an area statistic, a rotation analysis, an input-demand estimate. If the user needs "winter cereal area" and the model gives it at 96%, the wheat-barley distinction was never required. Reporting the merged accuracy alongside the split one lets the user make that call rather than the modeller.
Off-diagonal cells, directions and imbalance
- One large off-diagonal cell is a feature problem; scattered errors are usually labels.
- Check both directions. Wheat → barley and barley → wheat can differ a lot.
- Class imbalance inflates one direction. Try balanced weights before concluding — for the potatoes they moved seed recall by three points.
- A large effect size is not separability. It compares means; count how many of the larger class sit in the smaller one's range.
- Radar is unaffected by cloud and separates cereals reasonably well.
- Sub-field mixtures cannot be resolved by a field-level classifier.
- Cover crops confuse everything in the shoulder seasons.
- Report per-class recall, not overall accuracy.
- Publish the merged result too when merging is the recommendation.
Internal links
- How to classify crop types from a satellite time series — the pipeline
- Crop phenology and growing seasons explained — where the classes diverge
- Vegetation indices explained: NDVI, EVI, NDRE and when each fails — the index over the critical window
- How to extract phenology metrics from an NDVI time series — targeted features
- How to handle class imbalance in a spatial classifier — the other reason one direction dominates
- How to evaluate a spatial model in Python — reading the metrics
- Field boundary data explained: where it comes from — where the labels come from
- Spatial leakage explained — if the accuracy is high and the field performance is not
FAQ
Why does my classifier always confuse the same two crops?
Because they are genuinely similar in the features provided. Check the effect size between their mean profiles: below about 0.5 everywhere, no model can separate them.
Will a better model fix it?
No. A flat response to model changes is the diagnosis, not a tuning problem.
What feature would separate two cereals?
A red-edge index at heading or the slope of the senescence limb — the narrow windows where they differ. One targeted feature beats a hundred generic ones.
What if the separating date was cloudy?
Then the information is not in the data. Add another sensor — Landsat, or a radar series — rather than tuning the model.
Could the labels be wrong?
Often. Declared crop is the intended crop, and substitution errors concentrate exactly on similar pairs.
Is merging the two classes acceptable?
Yes, and it is often the right answer. Here red onions were recalled 16% of the time and yellow onions 85%; merged, the onion class is recalled 98% of the time. Report both.