Is terminal heat stressing the wheat during grain-fill this season?

Wheat spends its last few weeks filling its grain, pumping starch into the kernels before harvest. A burst of late-season heat in that window, called terminal heat stress, cuts the fill short and shrivels the harvest. This question watches two warning signs from space: how hot the fields are running, and whether the crop is browning off (drying down) earlier than it should.

Which data to use, and how

Use this dataset: MODIS MOD11A2 LST (short name MOD11A2). LST means land-surface temperature, how hot the actual ground and crop canopy are, measured from orbit. This product gives you one temperature map every 8 days at ~1 km resolution (each pixel is a 1 km square), going back to 2000. Start here because it’s a long, steady record with one number per pixel.

Then do four steps. The runnable code further down does all of this; this is just the idea.

1. Pick your field or region (a small latitude/longitude box around it) and find your grain-fill window: the few weeks when wheat is filling grain, usually Feb–March in the Indo-Gangetic plain, but set it from your own sowing date.
2. Map the heat. Pull MODIS LST over that window and see which fields run hottest. It gives a separate daytime and nighttime temperature.
3. Count the hot days. For a true heat threshold you need air temperature, not surface temperature, so bring in ERA5 2 m air temperature and count the days the daily high (Tmax) goes over a crop-specific cutoff. Get that cutoff from a local farming-weather (agromet) source. Don’t invent one.
4. Watch for early browning. Pull MODIS MOD13Q1 NDVI, a greenness index (high = lush green, falling = drying/dying), and check whether it turns down earlier this year than in a normal year. An early down-turn is the senescence (the crop’s natural end-of-life browning) arriving too soon, a possible heat-stress sign.

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

• How many hot days hit during the grain-fill window over your area (using ERA5 air temperature against a local threshold).
• Which fields run hottest: the canopy/surface heating pattern from MODIS LST (~1 km), or ECOSTRESS (~70 m) for a much finer field-scale view.
• Whether the crop is browning early: does NDVI turn down sooner than in a normal/reference year?
• This season vs a typical one, compared against the usual grain-fill window.
• Day vs night heat. MODIS gives both, and warm nights also stress grain-fill (the plant burns energy respiring instead of storing it).

What it can’t tell you

• An exact yield loss in kg/ha. These are stress indicators, not a yield model. They tell you something looks wrong, not how much grain you’ll lose. For numbers you need agronomic crop models plus ground data.
• One universal heat threshold. The “too hot” cutoff depends on the wheat variety and the growth stage, so there’s no single magic number to hard-code. Take it from local agromet guidance.
• Air temperature straight from the satellite. MODIS and ECOSTRESS measure the surface/canopy temperature, not the 2 m air temperature a weather station reports. For real threshold work use ERA5 or station data.
• Hour-by-hour heat spikes from MODIS alone. MODIS only flies over ~1–2 times a day, so it can miss a short afternoon spike. For sub-daily detail use ERA5 hourly data or a geostationary satellite (one parked over the same region).

Gotchas to watch for

• Surface temperature is not air temperature (LST ≠ air temp). MODIS/ECOSTRESS report the skin/canopy temperature, which on clear dry days can run several °C hotter than the 2 m air temperature. Hold LST up against an air-temp threshold and you’ll over-count hot days. Keep ERA5 or station air temperature for the actual threshold, and lean on LST for the spatial pattern of where it’s hottest.
• Clouds blank out the view. Optical LST and NDVI need a clear sky, so cloud (and the Feb–March haze/fog over the Indo-Gangetic plain) leaves gaps. Composite or gap-fill over a window rather than trusting a single day.
• The window must match the sowing date. Late-sown wheat fills grain later, so a fixed Feb–March window can miss the real vulnerable weeks entirely. Align the window to the local sowing calendar.
• The threshold isn’t one-size-fits-all. Terminal-heat cutoffs depend on cultivar and stage. Verify against IMD/agromet advisories rather than asserting a single number.
• 1 km pixels are mixed. A single MODIS 1 km pixel blends wheat with roads, villages, and other crops, smearing the reading. For smallholder fields, reach for HLS (30 m) or ECOSTRESS (~70 m) to get finer detail.

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, open the granules, map the grain-fill heat, count hot canopy days, compare this season’s NDVI down-turn to a normal year, and plot it. It needs a free Earthdata Login.

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

earthaccess.login(strategy="netrc")

aoi = (75, 28, 82, 30.5)      # Indo-Gangetic wheat belt as (W, S, E, N)
this_year, ref_year = 2026, 2024
grainfill = (f"{this_year}-02-10", f"{this_year}-03-31")  # SET FROM LOCAL SOWING DATE
lst_threshold_c = 33.0        # canopy-heat flag; SET FROM LOCAL AGROMET GUIDANCE

# 1. MODIS LST (MOD11A2, 8-day, 1 km) over the grain-fill window. The gridded
#    HDF-EOS opens with xarray; scale daytime LST (×0.02, in K) and go to °C.
lst_grans = earthaccess.search_data(short_name="MOD11A2", bounding_box=aoi, temporal=grainfill)
day_c = []
for fh in earthaccess.open(lst_grans):
    ds = xr.open_dataset(fh, engine="rasterio")          # HDF-EOS grid via GDAL
    lst_k = ds["LST_Day_1km"].where(ds["LST_Day_1km"] > 0) * 0.02
    day_c.append(lst_k.squeeze().values - 273.15)
lst_c = np.stack(day_c)                                   # (time, y, x) in °C

# 2. Map the spatial pattern of peak canopy heat, and count hot canopy days.
peak_lst = np.nanmax(lst_c, axis=0)                       # hottest 8-day value per pixel
hot_days = np.nansum(lst_c > lst_threshold_c, axis=0)     # 8-day periods over threshold
print(f"Peak grain-fill canopy LST: {np.nanmax(peak_lst):.1f} °C, "
      f"area-mean {np.nanmean(peak_lst):.1f} °C; "
      f"hottest pixel had {int(np.nanmax(hot_days))} of {lst_c.shape[0]} periods "
      f"above {lst_threshold_c:.0f} °C")

# 3. NDVI senescence timing (MOD13Q1, 250 m): build an area-mean greenness curve
#    for this season and a reference year; an earlier down-turn flags heat stress.
def ndvi_curve(year):
    grans = earthaccess.search_data(short_name="MOD13Q1", bounding_box=aoi,
                                    temporal=(f"{year}-01-01", f"{year}-04-30"))
    days, vals = [], []
    for fh in earthaccess.open(grans):
        ds = xr.open_dataset(fh, engine="rasterio")
        ndvi = (ds["_250m_16_days_NDVI"].where(ds["_250m_16_days_NDVI"] > -2000) * 1e-4)
        days.append(int(str(ds.attrs.get("RANGEBEGINNINGDATE", f"{year}-01-01"))[5:7]))
        vals.append(float(ndvi.mean()))
    return np.array(days), np.array(vals)
m_now, n_now = ndvi_curve(this_year)
m_ref, n_ref = ndvi_curve(ref_year)
peak_now, peak_ref = m_now[np.nanargmax(n_now)], m_ref[np.nanargmax(n_ref)]
verdict = "EARLY browning — possible terminal heat stress" if peak_now < peak_ref else "on schedule"
print(f"NDVI peak (greenness max) this year ~month {peak_now} vs ref ~month {peak_ref}: {verdict}")

# 4. Plot: heat map + the two NDVI curves side by side.
fig, (a, b) = plt.subplots(1, 2, figsize=(11, 4))
im = a.imshow(peak_lst, cmap="inferno"); a.set_title("Peak grain-fill LST (°C)")
fig.colorbar(im, ax=a, shrink=0.8)
b.plot(m_ref, n_ref, "o-", label=f"{ref_year} (reference)")
b.plot(m_now, n_now, "o-", label=f"{this_year}")
b.set_xlabel("month"); b.set_ylabel("area-mean NDVI"); b.set_title("Senescence timing"); b.legend()
plt.tight_layout(); plt.show()

# Cross-check before trusting it: LST is canopy SKIN temp, not 2 m air temp — count true
# heat days against ERA5 Tmax (or an IMD station) before calling a threshold breach (see gotchas).

Where the data comes from

Most of this is NASA data through one free login (Earthdata Login): the MODIS LST and NDVI products, ECOSTRESS, and HLS all come from NASA’s LP DAAC archive. The one outside piece is ERA5 air temperature, a non-NASA source (ECMWF/Copernicus, via the Climate Data Store). You need it because satellites can’t measure true 2 m air temperature directly, and that’s what heat thresholds are defined on.

Sources

• MODIS MOD11 LST: https://lpdaac.usgs.gov/products/mod11a2v061/
• MODIS MOD13 VI: https://lpdaac.usgs.gov/products/mod13q1v061/
• ECOSTRESS LST: https://lpdaac.usgs.gov/products/eco2lstev001/
• HLS: https://lpdaac.usgs.gov/products/hlsl30v002/
• ERA5: https://cds.climate.copernicus.eu/datasets/reanalysis-era5-single-levels
• IMD Agromet advisories (crop-stage thresholds): https://mausam.imd.gov.in/