Crop phenology and growing seasons explained
On this page
Problem statement
A field is not one object. In March it is bare soil, in June a closed canopy, in September a stubble — and any analysis that samples it once is an analysis of a moment, not of a crop.
Phenology is the study of that curve: when the crop emerges, how fast it develops, when it peaks, when it senesces. Nearly everything useful in agricultural remote sensing is a property of the curve rather than of any single image, and the two things that make it hard are that the curve is sampled irregularly by cloud and that its calendar shifts by weeks between years and between fields.
On a real Sentinel-2 series over an 11.4 ha spring-barley field in Flevoland in 2025, the field-median NDVI sat at 0.163 on 7 March, rose to 0.959 on 19 May, and was back to 0.163 by 20 July, after harvest. Only 39 of the 47 dates with any view of the field were at least 60% clear. By late September the stubble had greened again, to 0.723: a second rise that a single-season extractor reads as part of the barley.
Quick answer
Build the curve first, extract metrics second:
import pandas as pd
usable = season[season.clear_frac >= 0.6]
crop = usable[usable.date.between("2025-03-01", "2025-09-15")] # before the stubble regrows
s = (crop.set_index("date")["ndvi_med"]
.reindex(pd.date_range(crop.date.min(), crop.date.max(), freq="D"))
.interpolate("time")
.rolling(15, center=True, min_periods=5).mean())
peak_value = s.max()
peak_date = s.idxmax()
amplitude = peak_value - s.min()
sos = s[(s >= s.min() + 0.2 * amplitude) & (s.index < peak_date)].index.min()
eos = s[(s >= s.min() + 0.2 * amplitude) & (s.index > peak_date)].index.max()
print(f"start {sos.date()}, peak {peak_date.date()} at {peak_value:.3f}, end {eos.date()}")start 2025-03-28, peak 2025-05-18 at 0.950, end 2025-07-07
Interpolating to a daily grid and smoothing is not cosmetic. Even in this well-observed season the usable dates leave a 24-day hole across the peak, and any metric computed from the observations alone is a metric about when the sky was clear. The window stops at 15 September on purpose; the reason is in step 7.
The metrics, a moving calendar, irregular sampling
1. Know the metrics and what each is for
- Start of season — emergence or green-up, usually a threshold on the rising limb.
- Peak value and peak date — maximum greenness and when it occurred.
- End of season — senescence or harvest on the falling limb.
- Length of season — end minus start.
- Amplitude — peak minus the seasonal minimum, a proxy for how much biomass was grown.
- Integral — the area under the curve above the baseline, the best single proxy for accumulated production.
2. Expect the calendar to move
The same crop in the same field emerges weeks apart between years. Comparing fields or years by calendar date rather than by phenological stage is the commonest error in this area, and thermal time is the usual fix — Growing degree days explained.
3. Expect irregular sampling
Optical satellites see the ground only when it is clear. In the reference series, 39 usable dates across seven months still left a 24-day gap from 19 May to 12 June, across the peak, and an 18-day gap from 19 August to 6 September. A cloudier year, or a sensor with a longer revisit, leaves far bigger holes.
4. Interpolate deliberately, and say how
Linear interpolation on a daily grid plus a smoothing window is the simple approach. Savitzky–Golay, double-logistic fitting and harmonic regression are the standard alternatives; each imposes a shape, and the choice changes the metrics.
5. Choose thresholds and state them
A start of season at 20% of amplitude and one at 50% are different dates: 28 March against 6 April on the reference field, and 7 July against 27 June at the end. Neither is right; the undocumented one is wrong.
6. Work per field, not per pixel, where you can
Field-median series are far less noisy than pixel series, and the field is the management unit anyway. Pixel-level phenology matters for within-field variation and for fields with more than one crop.
7. Beware the double peak
A cover crop before a main crop, a cut grassland, a catch crop or regrowth after harvest — all produce two peaks in one season, and a single-peak extractor reports a nonsense season length. On the reference field the stubble greened to 0.723 by 28 September. Run on all of March to September, the same extractor put the end of season on 30 September and the length at 186 days; cut at 15 September, it gave 7 July and 101 days.
The seasonal curve of a spring-barley field
The data used on this page
One real field: the largest spring-barley parcel (Gerst, zomer-) declared in 2025 in the Dutch parcel register BRP Gewaspercelen within 5.55–5.75° E, 52.48–52.56° N, east of Dronten in Flevoland. Its series comes from every Sentinel-2 L2A scene in Earth Search over the parcel from 1 March to 30 September 2025. For each date the setup records the share of the parcel that the scene classification calls vegetation or bare soil (clear_frac), and the median NDVI (B8, B4) and NDRE (B8, B5) over those pixels. Where two scenes share a date, the clearer one is kept. Dates with no clear pixel are dropped. Reading about a hundred scenes takes several minutes the first time; after that season.csv is reused.
import os
import geopandas as gpd
import numpy as np
import pandas as pd
import rasterio
import requests
from rasterio.features import geometry_mask
from rasterio.windows import from_bounds
os.environ.update(AWS_NO_SIGN_REQUEST="YES", GDAL_DISABLE_READDIR_ON_OPEN="EMPTY_DIR",
CPL_VSIL_CURL_ALLOWED_EXTENSIONS=".tif")
if not os.path.exists("field.gpkg"):
url = "https://api.pdok.nl/rvo/gewaspercelen/ogc/v1/collections/brpgewas/items"
feats, params = [], {"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
brp.to_file("brp.gpkg")
brp[brp.gewas.str.startswith("Gerst, zomer")].nlargest(1, "ha").to_file("field.gpkg")
field = gpd.read_file("field.gpkg")
def window(href, bounds, shape=None):
with rasterio.open(href) as src:
win = from_bounds(*bounds, transform=src.transform)
extra = {"out_shape": shape} if shape else {}
return src.read(1, window=win, boundless=True, **extra).astype("float32"), src.window_transform(win)
if not os.path.exists("season.csv"):
body = {"collections": ["sentinel-2-l2a"], "intersects": field.geometry.iloc[0].__geo_interface__,
"datetime": "2025-03-01T00:00:00Z/2025-09-30T23:59:59Z", "limit": 100}
items = []
while body:
page = requests.post("https://earth-search.aws.element84.com/v1/search", json=body, timeout=120).json()
items += page["features"]
body = next((l["body"] for l in page["links"] if l["rel"] == "next"), None)
rows = []
for day in sorted({i["properties"]["datetime"][:10] for i in items}):
best = None
for item in (i for i in items if i["properties"]["datetime"][:10] == day):
a = item["assets"]
utm = field.to_crs(item["properties"]["proj:epsg"])
b = utm.total_bounds + np.array([-20, -20, 20, 20])
red, transform = window(a["red"]["href"], b)
inside = ~geometry_mask(utm.geometry, red.shape, transform)
scl, _ = window(a["scl"]["href"], b, red.shape)
clear = inside & np.isin(scl, [4, 5]) # vegetation, bare soil
if best is None or clear.sum() > best[1].sum():
best = (item, clear, inside, b, red)
item, clear, inside, b, red = best
if not clear.any():
continue
# DN / 10,000 is reflectance: no offset to remove, whatever the item's
# boa_offset_applied flag says (water reads 1.5–2.6% on 7 and 10 March either way)
nir, _ = window(item["assets"]["nir"]["href"], b, red.shape)
re1, _ = window(item["assets"]["rededge1"]["href"], b, red.shape)
r, n, e = red[clear], nir[clear], re1[clear]
rows.append({"date": day, "scene": item["id"],
"clear_frac": round(float(clear.sum() / inside.sum()), 3),
"ndvi_med": round(float(np.median((n - r) / (n + r))), 3),
"ndre_med": round(float(np.median((n - e) / (n + e))), 3)})
pd.DataFrame(rows).to_csv("season.csv", index=False)
season = pd.read_csv("season.csv", parse_dates=["date"])
print(f"{field.gewas.iloc[0]} parcel, {field.ha.iloc[0]:.1f} ha; "
f"{len(season)} dates with a clear view of it, {season.date.min().date()} to {season.date.max().date()}")Gerst, zomer- parcel, 11.4 ha; 47 dates with a clear view of it, 2025-03-07 to 2025-09-30
Example 1 — the real seasonal curve
usable = season[season.clear_frac >= 0.6].reset_index(drop=True)
print(f"{len(season)} dates with any clear view of the field, {len(usable)} at least 60% clear")
gaps = pd.DataFrame({"from": usable.date.shift(), "to": usable.date,
"days": usable.date.diff().dt.days}).dropna().nlargest(4, "days")
for g in gaps.itertuples():
print(f" gap {g._1:%d %b} → {g.to:%d %b}: {g.days:.0f} days")
half_months = (usable.set_index("date")
.resample("SMS")[["ndvi_med", "ndre_med"]]
.agg(["count", "median"]).dropna())
half_months.columns = ["n", "ndvi", "n_", "ndre"]
print(half_months[["n", "ndvi", "ndre"]].round(3).to_string())
print()
print(season.set_index("date").loc[["2025-03-07", "2025-05-19", "2025-07-20", "2025-09-28"],
["clear_frac", "ndvi_med", "ndre_med"]].to_string())47 dates with any clear view of the field, 39 at least 60% clear
gap 19 May → 12 Jun: 24 days
gap 19 Aug → 06 Sep: 18 days
gap 13 Sep → 28 Sep: 15 days
gap 02 Jul → 15 Jul: 13 days
n ndvi ndre
date
2025-03-01 2 0.164 0.111
2025-03-15 2 0.190 0.126
2025-04-01 7 0.580 0.383
2025-04-15 2 0.934 0.726
2025-05-01 5 0.943 0.758
2025-05-15 2 0.954 0.803
2025-06-01 2 0.911 0.678
2025-06-15 3 0.865 0.620
2025-07-01 1 0.355 0.199
2025-07-15 4 0.192 0.119
2025-08-01 3 0.198 0.110
2025-08-15 1 0.226 0.133
2025-09-01 3 0.242 0.142
2025-09-15 2 0.599 0.308
clear_frac ndvi_med ndre_med
date
2025-03-07 1.0 0.163 0.103
2025-05-19 1.0 0.959 0.816
2025-07-20 1.0 0.163 0.091
2025-09-28 1.0 0.723 0.331
Most half-months had two to five usable dates, but the first half of July — harvest — and the second half of August had one each. The plateau from mid-April to mid-June is NDVI saturating, not a canopy standing still: NDRE keeps moving underneath it (Example 3).
Example 2 — extract the metrics with the thresholds explicit
import numpy as np, pandas as pd
def phenology(dates, values, threshold=0.2, smooth_days=15):
s = (pd.Series(np.asarray(values, dtype=float), index=pd.DatetimeIndex(dates))
.reindex(pd.date_range(min(dates), max(dates), freq="D"))
.interpolate("time")
.rolling(smooth_days, center=True, min_periods=smooth_days // 3).mean()
.dropna())
base, peak = float(s.min()), float(s.max())
amp = peak - base
peak_date = s.idxmax()
level = base + threshold * amp
rising = s[(s >= level) & (s.index <= peak_date)]
falling = s[(s >= level) & (s.index >= peak_date)]
sos = rising.index.min() if len(rising) else None
eos = falling.index.max() if len(falling) else None
return {
"baseline": round(base, 3), "peak": round(peak, 3),
"amplitude": round(amp, 3), "peak_date": str(peak_date.date()),
"sos": str(sos.date()) if sos is not None else None,
"eos": str(eos.date()) if eos is not None else None,
"length_days": (eos - sos).days if sos is not None and eos is not None else None,
"integral": round(float((s - base).clip(lower=0).sum()), 1),
"threshold": threshold, "smooth_days": smooth_days,
}
whole = phenology(usable.date, usable.ndvi_med)
crop = usable[usable.date.between("2025-03-01", "2025-09-15")]
print("all of March–September:", {k: whole[k] for k in ("sos", "peak_date", "eos", "length_days")})
print("barley season only: ", phenology(crop.date, crop.ndvi_med))all of March–September: {'sos': '2025-03-28', 'peak_date': '2025-05-18', 'eos': '2025-09-30', 'length_days': 186}
barley season only: {'baseline': 0.169, 'peak': 0.95, 'amplitude': 0.781, 'peak_date': '2025-05-18', 'sos': '2025-03-28', 'eos': '2025-07-07', 'length_days': 101, 'integral': 67.1, 'threshold': 0.2, 'smooth_days': 15}
Returning the threshold and the smoothing window inside the result is the point: a phenology metric without them cannot be compared with anybody else's.
Example 3 — show how much the threshold matters
for t in (0.1, 0.2, 0.3, 0.5):
m = phenology(crop.date, crop.ndvi_med, threshold=t)
print(f"threshold {t:.1f}: sos {m['sos']}, eos {m['eos']}, "
f"length {m['length_days']} days")
plateau = usable[usable.ndvi_med > 0.9]
print(f"\nNDVI above 0.9 on {len(plateau)} dates, {plateau.date.min():%d %b} to {plateau.date.max():%d %b}: "
f"NDVI {plateau.ndvi_med.min():.3f}–{plateau.ndvi_med.max():.3f}, "
f"NDRE {plateau.ndre_med.min():.3f}–{plateau.ndre_med.max():.3f}")threshold 0.1: sos 2025-03-23, eos 2025-07-14, length 113 days
threshold 0.2: sos 2025-03-28, eos 2025-07-07, length 101 days
threshold 0.3: sos 2025-03-30, eos 2025-07-03, length 95 days
threshold 0.5: sos 2025-04-06, eos 2025-06-27, length 82 days
NDVI above 0.9 on 10 dates, 26 Apr to 13 Jun: NDVI 0.923–0.959, NDRE 0.659–0.816
A 0.1 and a 0.5 threshold differ by 14 days at the start and 17 at the end, and the season length runs from 82 to 113 days. That is larger than most of the differences anyone is trying to detect. On the plateau NDVI moved 0.036 across ten dates while NDRE moved 0.157.
The curve, not the image, is the unit
Why the curve, not the image, is the unit
A single image separates bare soil from vegetation and very little else — two crops at the same growth stage look alike. What distinguishes them is when they green up, how long they stay green and how fast they senesce, which is why crop classification from a time series works and classification from one date does not.
Why cloud makes this a modelling problem
Optical phenology depends on observations that arrive when the sky allows. A 24-day gap over the peak means the peak date is interpolated rather than observed, and the confidence in that metric should reflect it. Reporting the number of clear observations within each window is the minimum honest accompaniment to any phenology metric.
Why thermal time beats calendar time
Development rate depends on temperature, so a cold spring delays everything by a fortnight. Comparing fields or years on the calendar mixes weather with management; comparing them on accumulated degree days removes most of the weather and leaves the management.
Why saturation distorts the curve's shape
From 26 April to 13 June the reference series stayed between 0.923 and 0.959 — flat. The canopy was not flat; NDVI was. Over the same ten dates NDRE ranged from 0.659 to 0.816. That flattening makes the peak date arbitrary, flattens the apparent peak, and makes the integral insensitive to exactly the period when most biomass is being produced. A red-edge index gives a curve with the peak still in it.
Double cropping, cut grassland and perennials
- Double cropping gives two peaks. Detect them rather than fitting one.
- Cut grassland has several peaks a year and no single season.
- Perennials never reach a low baseline. The amplitude metric fails on orchards.
- Winter crops span the new year. Define the season, not the calendar year.
- Shadow and haze depress an index and look like early senescence.
- Snow gives a negative index and must be masked, not interpolated across.
- Smoothing shifts the peak by roughly half the window on an asymmetric curve.
- Report the observation count in every window the metric depends on.
Internal links
- How to extract phenology metrics from an NDVI time series — the implementation in full
- Vegetation indices explained: NDVI, EVI, NDRE and when each fails — why the peak is flat
- Growing degree days explained — the thermal clock
- How to classify crop types from a satellite time series — the main consumer of the curve
- How to build an NDVI time series for a polygon in Python — assembling the series
- How to resample a time series in xarray — regular grids from irregular observations
- How to build a cloud-free composite in Python — filling the gaps differently
- NDVI stops responding in a dense canopy — the flat peak
FAQ
What is crop phenology?
The timing of a crop's development through a season — emergence, green-up, peak, senescence and harvest — and the metrics derived from a vegetation index curve that tracks it.
Why do I need a time series rather than one image?
Because crops at the same growth stage look alike on one date. What separates them is the shape and timing of the curve.
How do I handle cloudy gaps?
Interpolate onto a daily grid and smooth, and report how many clear observations fall in each window. A real, fairly clear season still had a 24-day gap across the peak.
What threshold defines the start of season?
Whatever you state. On a real barley field a 20% and a 50% threshold differed by 9 days at the start and 10 at the end.
Should I work per pixel or per field?
Per field for management decisions — the series is far less noisy. Per pixel for within-field variation.
Why is my season curve flat at the top?
NDVI saturation. From 26 April to 13 June a real barley field sat between 0.923 and 0.959 while its NDRE ranged from 0.659 to 0.816.