Is the growing season shifting — are plants greening up earlier or staying green longer?

Phenology is the calendar of plant life: the date leaves come out in spring (greenup), the date they turn and fall (senescence/dormancy), and how long the green stays in between. Satellites watch whole landscapes turn green and brown from orbit every year. So you can ask a simple question. Are plants waking up earlier, or staying green longer, than they used to?

Which data to use, and how

Use this dataset: MODIS MCD12Q2 (short name MCD12Q2), NASA’s ready-made phenology product. It has already found the greenup date for you. For every 500-metre patch of ground, once a year, it records the day the plants greened up, peaked, and went dormant. That’s why it’s the easiest start: the hard part is pre-computed. (Later you can add the raw greenness datasets below for more detail.)

Then do four steps (the runnable code further down does all of this; this is just the idea):

1. Pick your region. A small latitude/longitude box around it, and the range of years you want to compare.
2. Read the dates. For each year, pull two numbers out of MCD12Q2: the greenup day and the dormancy day. Growing-season length is just dormancy − greenup.
3. Average your region. Take the average greenup day across your whole box for each year, so you have one number per year instead of millions of pixels.
4. Fit a trend. Draw a straight line through “greenup day vs. year.” A downward slope means greenup is creeping earlier (days per decade). Do the same for dormancy and for season length.

That’s the whole method. Below, “How a scientist answers this” names the exact tools, and the Run it cell does it live on synthetic data so you can see each step work.

What you can find out

• Whether greenup is starting earlier over your chosen years (the trend in MCD12Q2’s Greenup day-of-year, either for the whole region or pixel-by-pixel).
• Whether plants stay green longer: later Senescence/Dormancy dates, and a longer growing season (Dormancy − Greenup).
• Whether peak greenness is shifting: the date of peak green (MidGreenup/Peak) and how much green there is (EVI_Amplitude/EVI_Area, EVI being a satellite “greenness” index).
• Whether the season is getting more productive, by adding the greenness datasets MOD15A2H (a leaf-area measure called LAI/FPAR) or MOD13Q1 (the NDVI/EVI greenness indices) and summing greenness across the season.
• Whether the shift tracks warmer springs. If you bring in a temperature record of your own, you can line the two up and see if earlier greenup matches warmer years.

What it can’t tell you

• Why the season shifted. The data shows that greenup moved, not what moved it: warming, more CO₂, or a farmer changing what they plant. Pinning down the cause needs extra climate or land-management data.
• How much carbon the plants took in. Greenness tells you there are leaves, not how much carbon they captured. For actual carbon uptake you need a different product (a “flux” product like MOD17 GPP, or ground towers).
• Which crop is which. Each pixel is 250–500 m wide and blends whole fields together, and MCD12Q2 reports one generic “land greening” date. It can’t tell corn from soybean.
• Timing finer than the satellite’s rhythm. The greenness datasets only refresh every 8 or 16 days, so a greenup date is never more precise than that window. Don’t read a 3-day shift as real.
• Whether two green-ups a year mean double-cropping. MCD12Q2 can report up to two green-up cycles per pixel, but it won’t tell you if the second one is a second planted crop or just noise.

Gotchas to watch for

• The dates aren’t day-of-year, so convert them first. MCD12Q2 stores greenup as “days since 1 January 1970,” not “day 105 of this year.” You have to convert. Also throw out the fill value 32767, which marks “no valid date here” (e.g. bare ground). Leave it in and it wrecks your average.
• There can be two cycles per pixel, so pick the right one. MCD12Q2 reports up to two green-up cycles a year (cycle 0 and cycle 1). For ordinary single-season areas use cycle 0. Only double-cropped or irrigated land fills in cycle 1.
• Very dense plants max out the greenness meter. The leaf-area measure (LAI) flattens out once a canopy is thick (around LAI 5–6), so a lush cornfield can look “equally green” across its peak even when it isn’t. Matters only if you compare peak-summer productivity in dense crops.
• The old satellite ages, so cross-check it. MODIS Terra’s sensor drifts late in the record, so a long-term trend can partly be the instrument aging, not the plants. The fix: confirm the trend against the newer VIIRS record (VNP15A2H) before believing a multi-decade signal.
• 23 years is short for climate claims. A two-decade record is brief for a climate trend. Report the trend with its confidence interval (the plausible range), not as a settled fact.

The real-data code

The Run it cell above runs this method on synthetic data with no login. Below is the whole recipe against the real archive: search MCD12Q2, open the HDF-EOS phenology layers, convert the dates, region-average each year, fit the trend in days-per-decade, and plot it. It needs a free Earthdata Login.

import earthaccess, numpy as np, xarray as xr
import datetime as dt
import matplotlib.pyplot as plt

earthaccess.login(strategy="netrc")

aoi = (-96.0, 39.0, -87.0, 44.0)        # your region as (W, S, E, N) — here, the US Corn Belt
years = range(2001, 2022)               # one MCD12Q2 granule set per year

# 1. One greenup/dormancy/length number per year. MCD12Q2 stores Greenup & Dormancy
#    as "days since 1970-01-01"; convert to day-of-year, drop the 32767 fill, use
#    cycle 0 (band 0) for single-season land, and area-average over the box.
def yearly_phenology(y):
    grans = earthaccess.search_data(
        short_name="MCD12Q2", bounding_box=aoi,
        temporal=(f"{y}-01-01", f"{y}-12-31"))
    gd, dd = [], []
    for fh in earthaccess.open(grans):
        ds = xr.open_dataset(fh, engine="rasterio", group="MODIS_Grid_2D")
        green = ds["Greenup"].isel(band=0).values.astype("float64")
        dorm = ds["Dormancy"].isel(band=0).values.astype("float64")
        green[green == 32767] = np.nan       # 32767 = "no valid date" (bare ground)
        dorm[dorm == 32767] = np.nan
        # days-since-1970 -> day-of-year within year y
        jan1 = (dt.date(y, 1, 1) - dt.date(1970, 1, 1)).days
        gd.append(np.nanmean(green - jan1))
        dd.append(np.nanmean(dorm - jan1))
    return np.nanmean(gd), np.nanmean(dd)

greenup, dormancy = zip(*[yearly_phenology(y) for y in years])
greenup, dormancy = np.array(greenup), np.array(dormancy)
length = dormancy - greenup             # season length = dormancy - greenup

# 2. Fit a straight line of region-mean day vs year; slope*10 = days per decade.
yrs = np.array(list(years))
def trend(series):
    ok = ~np.isnan(series)
    slope = np.polyfit(yrs[ok], series[ok], 1)[0]
    return slope * 10.0                  # days per decade

g_dec, d_dec, l_dec = trend(greenup), trend(dormancy), trend(length)

# 3. One-line verdict.
verb = "earlier" if g_dec < 0 else "later"
print(f"Greenup is moving {verb}: {g_dec:+.1f} days/decade | "
      f"dormancy {d_dec:+.1f} d/dec | season length {l_dec:+.1f} d/dec")

# 4. Plot all three series with their fitted lines.
fig, ax = plt.subplots()
for series, lab in [(greenup, "greenup"), (dormancy, "dormancy"), (length, "length")]:
    ok = ~np.isnan(series)
    ax.plot(yrs[ok], series[ok], "o-", label=lab)
    ax.plot(yrs[ok], np.poly1d(np.polyfit(yrs[ok], series[ok], 1))(yrs[ok]), "--", lw=1)
ax.set_xlabel("year"); ax.set_ylabel("day of year"); ax.legend()
plt.tight_layout(); plt.show()

# Before you trust it: re-run against the VIIRS record (VNP15A2H) — MODIS Terra's sensor
# drifts late in the record, so confirm the trend survives the instrument handoff (see gotchas).

Where the data comes from

This is the simple case: everything lives in one place under one free login. All four products (MCD12Q2, MOD15A2H, MOD13Q1, and the VIIRS continuation VNP15A2H) come from NASA’s LP DAAC (the Land Processes archive) via Earthdata Login. The only real work is stitching time together (MODIS hands off to VIIRS around 2012) and opening the slightly fiddly HDF-EOS file format. No juggling multiple logins. If you later add a climate record from another archive, see recipes/r01-three-daac-composition.mdx for the multi-archive pattern.

Sources + further reading

• MCD12Q2 (Land Cover Dynamics / phenology) user guide: https://lpdaac.usgs.gov/products/mcd12q2v061/
• MOD15A2H (LAI/FPAR) user guide: https://lpdaac.usgs.gov/products/mod15a2hv061/
• MOD13Q1 (Vegetation Indices) user guide: https://lpdaac.usgs.gov/products/mod13q1v061/
• VNP15A2H (VIIRS LAI/FPAR) user guide: https://lpdaac.usgs.gov/products/vnp15a2hv002/
• earthaccess docs: https://earthaccess.readthedocs.io/