How to classify crop types from a satellite time series
On this page
Problem statement
Crop classification from one image is land-cover classification with extra steps: at peak canopy, wheat, barley and grass all look like dense green vegetation. What distinguishes them is the shape of the season — when they green up, how long they stay green, how fast they senesce — so the features have to come from a time series.
Two things then decide the accuracy, and neither is the classifier. The first is whether the features capture the season rather than the weather, which is why irregular cloud-driven sampling has to be resampled onto a common calendar. The second is whether the validation respects space, because fields are large and a random pixel split puts pixels from the same field in both training and test, which inflates the accuracy by a wide margin.
Quick answer
Interpolate onto a fixed calendar, extract features per field, and split by field:
import numpy as np, pandas as pd
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import GroupKFold, cross_val_score
# one row per field, columns = index value on each of a fixed set of dates
X = features.drop(columns=["field_id", "crop"]).values
y = features["crop"].values
groups = features["field_id"].values
clf = RandomForestClassifier(n_estimators=500, min_samples_leaf=2, random_state=0)
scores = cross_val_score(clf, X, y, groups=groups, cv=GroupKFold(5), n_jobs=-1)
print(f"field-grouped accuracy: {scores.mean():.3f} ± {scores.std():.3f}")GroupKFold on the field identifier is the difference between a defensible number and a meaningless one. On the Flevoland fields used below, a plain KFold on pixels reports 0.997 against 0.960 with each field kept whole (Example 3), and all of that difference is leakage.
What you get — from The data used on this page
57 dates, (515, 774) pixels, 867 fields: {'ware potato': 204, 'winter wheat': 183, 'onion': 173, 'grass': 129, 'sugar beet': 109, 'carrot': 69}
Step-by-step solution
1. Build a cloud-screened time series per field
Field medians rather than pixel values: they are far less noisy, the field is the management unit, and the resulting table is small enough to model comfortably.
2. Resample onto a fixed calendar
Observations arrive when the sky allows — over the Flevoland fields below, 57 Sentinel-2 passes from March to October gave each field between 34 and 50 usable dates, and every field a gap of at least three weeks somewhere in the season. Interpolating onto, say, ten-day steps from 1 March to 31 October gives every field the same feature vector, which a classifier requires.
3. Use several indices, not just one
NDVI for cover, a red-edge index for biomass at closure, and optionally a shortwave index for senescence and residue. NDVI alone is flat for weeks at the peak — on a June scene over these same fields, 50.7% of vegetation pixels were above 0.8 — which removes exactly the period when crops diverge most in biomass.
4. Add phenology metrics as features
Start of season, peak date, peak value, length and integral summarise the curve in a handful of numbers that transfer between years better than raw dates do.
5. Add thermal time if you have it
Expressing the phenology metrics in accumulated degree days rather than calendar days removes most of the between-year weather variation, which is what makes a model trained in one year work in the next.
6. Split by field, and preferably by region
GroupKFold on the field is the minimum. If the model is to be applied to a different area, a spatial block split is the honest test — Spatial cross-validation explained.
7. Report the confusion matrix, not the accuracy
Overall accuracy hides the pair that matters. Two cereals confused with each other is a different problem from a cereal confused with grassland, and only the matrix shows which.
Code examples
The data used on this page
The Flevoland box of the crop phenology guide — 5.55–5.75° E, 52.48–52.56° N — and every 2025 parcel of at least a hectare of six crops in the Dutch register BRP Gewaspercelen: winter wheat, ware potato, onion (yellow and red together), sugar beet, carrot and grass (temporary and permanent together). The imagery is every Sentinel-2 L2A pass from Earth Search between March and October 2025 with less than 60% cloud — one per date, the clearest — read at 20 m, with NDVI and NDRE set to NaN wherever the scene classification is not vegetation or bare soil. That leaves genuinely irregular series: clouded dates differ from field to field. Parcels are shrunk by 20 m before they are rasterised, so edge pixels do not mix crops. Building the stack reads about sixty scenes and takes a few minutes the first time.
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
os.environ.update(AWS_NO_SIGN_REQUEST="YES", GDAL_DISABLE_READDIR_ON_OPEN="EMPTY_DIR")
CROPS = {"Tarwe, winter-": "winter wheat", "Aardappelen, consumptie": "ware potato", "Uien, gele zaai-": "onion",
"Uien, rode zaai-": "onion", "Bieten, suiker-": "sugar beet", "Winterpeen, productie": "carrot",
"Grasland, tijdelijk": "grass", "Grasland, blijvend": "grass"}
if not os.path.exists("stack.npz"):
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
fields = brp[brp["gewas"].isin(CROPS) & (brp["ha"] >= 1)].reset_index(drop=True).to_crs(32631)
fields["crop"], fields["field_id"] = fields["gewas"].map(CROPS), np.arange(1, len(fields) + 1)
fields.to_file("fields.gpkg")
found = requests.post("https://earth-search.aws.element84.com/v1/search", timeout=120, json={
"collections": ["sentinel-2-l2a"], "bbox": [5.55, 52.48, 5.75, 52.56], "limit": 300,
"datetime": "2025-03-01T00:00:00Z/2025-10-31T23:59:59Z", "query": {"eo:cloud_cover": {"lt": 60}}}).json()["features"]
best = {}
for item in found:
day = item["properties"]["datetime"][:10]
if day not in best or item["properties"]["eo:cloud_cover"] < best[day]["properties"]["eo:cloud_cover"]:
best[day] = item
box = fields.total_bounds + np.array([-40, -40, 40, 40])
ndvi, ndre = [], []
for day in sorted(best):
assets = best[day]["assets"]
with rasterio.open(assets["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)
def read(key):
with rasterio.open(assets[key]["href"]) as src:
return src.read(1, window=from_bounds(*extent, transform=src.transform), out_shape=shape,
resampling=Resampling.average).astype("float32")
red, nir, re1 = read("red"), read("nir"), read("rededge1")
with rasterio.open(assets["scl"]["href"]) as src:
clear = np.isin(src.read(1, window=win), [4, 5])
with np.errstate(all="ignore"):
ndvi.append(np.where(clear, (nir - red) / (nir + red), np.nan).astype("float32"))
ndre.append(np.where(clear, (nir - re1) / (nir + re1), np.nan).astype("float32"))
ids = rasterize(zip(fields.geometry.buffer(-20), fields["field_id"]), out_shape=shape, transform=grid,
fill=0, dtype="int32")
np.savez_compressed("stack.npz", ndvi=np.stack(ndvi, 2), ndre=np.stack(ndre, 2), ids=ids,
dates=np.array(sorted(best)), transform=np.array(grid)[:6])
stack = np.load("stack.npz")
ndvi_stack, ndre_stack, ids = stack["ndvi"], stack["ndre"], stack["ids"]
dates = pd.to_datetime(stack["dates"])
transform = rasterio.Affine(*stack["transform"])
fields = gpd.read_file("fields.gpkg")
print(len(dates), "dates,", ndvi_stack.shape[:2], "pixels,", len(fields), "fields:",
fields["crop"].value_counts().to_dict())57 dates, (515, 774) pixels, 867 fields: {'ware potato': 204, 'winter wheat': 183, 'onion': 173, 'grass': 129, 'sugar beet': 109, 'carrot': 69}
Example 1 — per-field series onto a fixed calendar
import numpy as np, pandas as pd, geopandas as gpd, rasterio
from rasterio.features import rasterize
def field_series(index_stack, dates, transform, shape, fields, min_pixels=25):
"""index_stack: rows × cols × dates. Returns one row per field per date."""
ids = rasterize(((g, i + 1) for i, g in enumerate(fields.geometry)),
out_shape=shape, transform=transform, fill=0, dtype="int32")
rows = []
for i in np.unique(ids):
if i == 0:
continue
sel = ids == i
if sel.sum() < min_pixels:
continue
for d, layer in zip(dates, np.moveaxis(index_stack, 2, 0)):
v = layer[sel]
v = v[np.isfinite(v)]
if v.size < min_pixels // 2:
continue
rows.append({"field_id": int(i), "date": pd.Timestamp(d),
"value": float(np.median(v)), "n": int(v.size)})
return pd.DataFrame(rows)
def to_calendar(series, start, end, freq="10D"):
grid = pd.date_range(start, end, freq=freq)
out = {}
for fid, g in series.groupby("field_id"):
s = (g.set_index("date")["value"]
.reindex(pd.DatetimeIndex(g["date"]).union(grid)).sort_index()
.interpolate("time", limit_direction="both").reindex(grid))
out[fid] = s.values
return pd.DataFrame(out, index=grid).Tshrunk = fields.set_geometry(fields.geometry.buffer(-20))
ndre_series = field_series(ndre_stack, dates, transform, ids.shape, shrunk)
per_field = ndre_series.groupby("field_id")["date"].agg(["count", lambda d: d.sort_values().diff().dt.days.max()])
per_field.columns = ["clear_dates", "longest_gap_days"]
print(per_field.describe().loc[["min", "50%", "max"]].round(0).to_string())
ndre_cal = to_calendar(ndre_series, "2025-03-01", "2025-10-31")
ndvi_cal = to_calendar(field_series(ndvi_stack, dates, transform, ids.shape, shrunk), "2025-03-01", "2025-10-31")
print(ndre_cal.shape, "fields × calendar steps; NaN left:", int(ndre_cal.isna().sum().sum())) clear_dates longest_gap_days
min 34.0 22.0
50% 42.0 24.0
max 50.0 30.0
(776, 25) fields × calendar steps; NaN left: 0
Building the union index before reindexing is what makes the interpolation use the real observation dates rather than only the grid. The first version had two faults. It called g["date"].union(grid), which fails — a pandas Series has no union; the dates have to be a DatetimeIndex first. And it interpolated without limit_direction="both": the calendar starts on 1 March, six days before the first pass, time interpolation does not extrapolate, and every one of the 776 fields got a NaN in its first column — which the dropna() before modelling then turned into zero fields.
Example 2 — features from the curve
import numpy as np, pandas as pd
def curve_features(row, dates, prefix="ndre"):
v = np.asarray(row, dtype="float32")
ok = np.isfinite(v)
if ok.sum() < 4:
return {}
base, peak = float(np.nanmin(v)), float(np.nanmax(v))
amp = peak - base
peak_i = int(np.nanargmax(v))
level = base + 0.2 * amp
rising = np.flatnonzero(ok & (v >= level) & (np.arange(len(v)) <= peak_i))
falling = np.flatnonzero(ok & (v >= level) & (np.arange(len(v)) >= peak_i))
return {
f"{prefix}_base": base, f"{prefix}_peak": peak, f"{prefix}_amp": amp,
f"{prefix}_peak_doy": dates[peak_i].dayofyear,
f"{prefix}_sos_doy": dates[rising[0]].dayofyear if len(rising) else np.nan,
f"{prefix}_eos_doy": dates[falling[-1]].dayofyear if len(falling) else np.nan,
f"{prefix}_integral": float(np.nansum(np.clip(v - base, 0, None))),
f"{prefix}_slope_up": float(np.nanmax(np.diff(v[:peak_i + 1]))) if peak_i > 0 else np.nan,
f"{prefix}_slope_down": float(np.nanmin(np.diff(v[peak_i:]))) if peak_i < len(v) - 1 else np.nan,
}grid = ndre_cal.columns
curves = pd.DataFrame([curve_features(ndre_cal.loc[f], grid, "ndre") | curve_features(ndvi_cal.loc[f], grid, "ndvi")
for f in ndre_cal.index], index=ndre_cal.index)
crop_of = fields.set_index(fields.index + 1)["crop"] # field_series numbers fields from 1
summary = curves.join(crop_of.rename("crop")).groupby("crop")[["ndre_peak_doy", "ndre_slope_down"]].median()
print(summary.round(3).to_string()) ndre_peak_doy ndre_slope_down
crop
carrot 250.0 -0.119
grass 110.0 -0.266
onion 200.0 -0.167
sugar beet 200.0 -0.118
ware potato 170.0 -0.160
winter wheat 140.0 -0.272
The slope features earn their place, though not always for the reason expected. The textbook contrast — a cereal senesces sharply and a grass does not — fails here: grass in Flevoland is mown, and its steepest fall is as sharp as winter wheat's (−0.266 against −0.272 per ten-day step). What separates them is when: grass peaks around day 110, wheat around day 140. Shape features are useful precisely because they combine rate and timing; any one of them alone can mislead.
Example 3 — evaluate honestly
import numpy as np, pandas as pd
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import GroupKFold, cross_val_predict
from sklearn.metrics import confusion_matrix, classification_report
def evaluate(X, y, groups, n_splits=5):
clf = RandomForestClassifier(n_estimators=500, min_samples_leaf=2,
class_weight="balanced_subsample", random_state=0)
pred = cross_val_predict(clf, X, y, groups=groups,
cv=GroupKFold(n_splits), n_jobs=-1)
labels = sorted(np.unique(y))
cm = pd.DataFrame(confusion_matrix(y, pred, labels=labels),
index=labels, columns=labels)
print(cm.to_string())
print()
print(classification_report(y, pred, labels=labels, zero_division=0))
off = cm.values.copy()
np.fill_diagonal(off, 0)
i, j = np.unravel_index(off.argmax(), off.shape)
print(f"largest confusion: {labels[i]} predicted as {labels[j]} "
f"({off[i, j]} of {cm.values[i].sum()})")
return cmfeatures = pd.concat([ndre_cal.add_prefix("ndre_"), ndvi_cal.add_prefix("ndvi_"), curves], axis=1).dropna()
X = features.to_numpy()
y = crop_of.loc[features.index].to_numpy()
cm = evaluate(X, y, groups=features.index.to_numpy()) carrot grass onion sugar beet ware potato winter wheat
carrot 52 0 1 0 3 0
grass 0 83 0 1 0 1
onion 1 0 89 0 4 0
sugar beet 1 1 0 98 3 0
ware potato 0 0 4 3 181 0
winter wheat 0 0 0 0 0 168
precision recall f1-score support
carrot 0.96 0.93 0.95 56
grass 0.99 0.98 0.98 85
onion 0.95 0.95 0.95 94
sugar beet 0.96 0.95 0.96 103
ware potato 0.95 0.96 0.96 188
winter wheat 0.99 1.00 1.00 168
accuracy 0.97 694
macro avg 0.97 0.96 0.96 694
weighted avg 0.97 0.97 0.97 694
largest confusion: onion predicted as ware potato (4 of 94)
That is one row per field, so every fold already keeps fields whole. The leakage shows up when the rows are pixels — the setup most tutorials use:
from sklearn.model_selection import GroupKFold, KFold, cross_val_score
clear = [k for k in range(len(dates)) if np.isfinite(ndvi_stack[..., k][ids > 0]).mean() > 0.97]
pixel_features = np.concatenate([ndvi_stack[..., clear], ndre_stack[..., clear]], axis=2)
rng = np.random.default_rng(0)
rows, labels, groups = [], [], []
for fid in np.unique(ids[ids > 0]):
where = np.argwhere(ids == fid)
if len(where) < 40:
continue
pick = where[rng.choice(len(where), 40, replace=False)]
values = pixel_features[pick[:, 0], pick[:, 1]]
values = values[np.isfinite(values).all(axis=1)]
rows.append(values); labels += [crop_of[fid]] * len(values); groups += [fid] * len(values)
Xp, yp, gp = np.vstack(rows), np.array(labels), np.array(groups)
clf = RandomForestClassifier(n_estimators=300, min_samples_leaf=2, n_jobs=4, random_state=0)
print(f"{len(Xp):,} pixels from {len(np.unique(gp))} fields, {len(clear)} near-clear dates")
print("random pixel folds:", cross_val_score(clf, Xp, yp, cv=KFold(5, shuffle=True, random_state=0)).mean().round(3))
print("fields kept whole: ", cross_val_score(clf, Xp, yp, groups=gp, cv=GroupKFold(5)).mean().round(3))27,398 pixels from 717 fields, 19 near-clear dates
random pixel folds: 0.997
fields kept whole: 0.96
Printing the largest off-diagonal cell by name is a small thing that focuses the next iteration far better than an accuracy figure does.
Explanation
Why the calendar has to be fixed
A classifier needs the same feature in the same column for every sample. Raw observations arrive on different dates for different fields — different orbits, different cloud — so a raw stack has a column that means "the fourth clear date", which is not a variable. Interpolating onto a fixed grid makes column four mean "1 April" for everybody.
Why field-grouped validation matters so much
A 25-hectare field at 20 m resolution is about 600 pixels, all of the same crop, all highly correlated. A random pixel split puts most fields in both sets, so the model is tested on pixels it has effectively seen. The reported accuracy then measures memorisation: 0.997 on the Flevoland pixels, against 0.960 once each field sits in one fold only. Nearly all of the remaining 4% of errors were hidden by the random split.
Why thermal time transfers between years
A model trained on calendar features learns the date wheat peaks in the training year. In a cold year wheat peaks later, and the model fails. A model trained on thermal features learns the accumulated temperature at which wheat peaks, which is a property of the crop rather than of the calendar. Expressing the phenology metrics in degree days is the cheapest way to make a model portable across seasons.
Why the confusion matrix names the fix
Two winter cereals confused with each other need a feature that separates them — often a single early-spring or late-senescence date, or a red-edge index at booting. A cereal confused with grassland needs a better mask or a longer series. The overall accuracy cannot tell you which problem you have, and the matrix can.
Edge cases or notes
- Class imbalance is severe. A minor crop with thirty fields needs weighting or resampling.
- The training labels come from a declaration, which is the intended crop, not always the grown one.
- Mixed fields exist. A field-median feature for a field with two crops is neither.
- Cover crops confuse everything in the shoulder seasons.
- Grassland is a continuum, not a class.
- Different years need different models unless the features are thermal.
- Report per-class recall. Overall accuracy is dominated by the commonest crop.
- Record the dates and the calendar grid with the model.
Internal links
- Crop phenology and growing seasons explained — the curve the features summarise
- Vegetation indices explained: NDVI, EVI, NDRE and when each fails — which index to stack
- A crop classifier confuses two crops every year — reading the confusion matrix
- Growing degree days explained — making features transfer between years
- Spatial leakage explained — why a random split inflates the score
- How to do spatial cross-validation in Python — the block split
- How to handle class imbalance in a spatial classifier — minor crops
- How to build an NDVI time series for a polygon in Python — assembling the series
FAQ
Can I classify crops from a single image?
No. At peak canopy most crops look alike. What separates them is the shape of the season, which needs a time series.
How do I handle irregular observation dates?
Interpolate each field's series onto a fixed calendar — ten-day steps through the season — so every field has the same feature vector.
Why is my accuracy so high in testing and poor in practice?
Almost certainly a random pixel split. Fields are large and internally correlated; use GroupKFold on the field identifier.
Which classifier should I use?
A random forest or gradient boosting is fine. The classifier is rarely the limiting factor; the features and the validation are.
How do I make a model work in a different year?
Express the phenology features in accumulated degree days rather than calendar days.
What should I report?
The confusion matrix and per-class recall. Overall accuracy is dominated by the commonest crop and hides the pair that is actually failing.