Is my region shifting from one crop a year to two or three?

When a field is planted, it turns green as the crop grows, then browns off at harvest. That’s one full “green wave” per crop. Satellites measure greenness from orbit as a number called NDVI. Count how many green waves a place goes through in a year and you know whether it grows one crop, two, or three. When that count climbs over the years, farmers and economists call it intensification.

Which data to use, and how

Use this dataset: MODIS MOD13Q1 (short name MOD13Q1), a greenness (NDVI/EVI) snapshot taken every 16 days at 250 m going back to 2000. The long record and ready-made greenness numbers make it the easiest place to start. If your fields are small and crowded together, swap in the sharper HLS (HLSS30 / HLSL30) at 30 m later, but start coarse.

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 years you care about.
2. Get the green wave. For each pixel, line up its NDVI readings through the year into a time series, and smooth out the jitter from clouds and noise.
3. Count the waves per year. Each green-up-then-brown-off bump above a greenness cutoff counts as roughly one crop. One bump is single-cropping, two is double, three is triple.
4. Watch the count over time. Track cycles-per-year across the record and map where it’s rising. Up means adding crop seasons (intensifying); flat means no change.

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 crops per year your region grows: single, double, or triple cropping, from counted NDVI cycles.
• Whether intensity is rising over time, by finding where cycles-per-year are climbing across the record.
• A map of intensification showing which sub-regions added a growing season.
• The crop calendar: when each green-up and harvest happens through the year.
• Coarse vs fine: the MODIS 250 m trend compared with sharper HLS 30 m for fragmented fields.

What it can’t tell you

• Which crop is actually grown. NDVI cycles count seasons, not species. Wheat, rice, and sugarcane all just look like “green.” To name the crop, pair this with crop-type maps or ground truth.
• Exact field boundaries (at MODIS scale). A 250 m pixel is bigger than many smallholder plots, so it blends several fields together. For sharp boundaries drop to HLS/Sentinel-2 at 30 m.
• Whether the water came from irrigation or rain. A second green wave suggests irrigation or leftover soil moisture, but it doesn’t prove the source. You need extra water/irrigation data for that.
• Yield or total production. Counting crops tells you how often a field is farmed (intensity), not how much it produced. For output, combine with yield models and census data.

Gotchas to watch for

• One pixel can cover several fields (mixed pixels). At 250 m, a single MODIS pixel may blend fields on different planting calendars, smearing the green waves together or double-counting them, so your crop count comes out blurry. For fragmented smallholder landscapes, drop to HLS at 30 m.
• Forests and grass green up too. Natural vegetation turns green and brown on its own schedule, producing NDVI cycles that aren’t crops and inflating your count. Mask to cropland with a land-cover map before you start counting.
• Clouds punch holes in the record. During the monsoon (the kharif growing season) heavy cloud hides the ground and thins the data, and a real green wave can vanish into the gap. Smooth and gap-fill the series, and lean on the built-in quality flag (pixel-reliability QA) that marks dodgy readings.
• The cutoffs are knobs, not laws. The greenness threshold and the minimum bump size (prominence) you use to call something “a crop” depend on the crop and region, so numbers borrowed from elsewhere will mis-count. Calibrate them against a local crop calendar rather than trusting one fixed value.
• Check the result against the census. Before claiming a region really shifted from one crop to two, cross-check your trend with the agricultural census / district crop statistics. If they disagree, trust the ground data and re-examine the code.

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 MOD13Q1, open the NDVI grids, count green-wave cycles per pixel per year, trend cycles-per-year over the record, and plot it. It needs a free Earthdata Login.

import earthaccess, numpy as np, xarray as xr
from scipy.signal import savgol_filter, find_peaks
import matplotlib.pyplot as plt

earthaccess.login(strategy="netrc")

aoi = (74, 29.5, 77, 32)     # your region as (W, S, E, N) — here, Punjab
years = range(2010, 2025)

# Greenness / amplitude rules — tune and verify against a local crop calendar.
ndvi_threshold = 0.4    # a peak must exceed this to count as a crop cycle
min_prominence = 0.15   # min NDVI rise/fall so noise isn't counted as a crop

def cycles_per_year(ndvi_1d):
    """Count crop cycles in one year's (time-ordered) NDVI series."""
    s = savgol_filter(ndvi_1d, window_length=5, polyorder=2)
    peaks, _ = find_peaks(s, height=ndvi_threshold, prominence=min_prominence)
    return len(peaks)

# 1. Search MOD13Q1 (16-day NDVI, 250 m) across the whole year range, then open
#    each granule's NDVI grid (scale 1e-4) and screen with pixel-reliability QA.
def ndvi_stack(y):
    grans = earthaccess.search_data(
        short_name="MOD13Q1", bounding_box=aoi,
        temporal=(f"{y}-01-01", f"{y}-12-31"))
    frames = []
    for fh in earthaccess.open(grans):
        ds = xr.open_dataset(fh, engine="rasterio")  # HDF-EOS subdatasets
        ndvi = ds["_250m_16_days_NDVI"].squeeze() * 1e-4
        qa = ds["_250m_16_days_pixel_reliability"].squeeze()
        frames.append(ndvi.where(qa == 0))           # 0 = good data
    return xr.concat(frames, dim="time").sortby("time")

# 2. Count cycles per pixel per year -> a (year, y, x) intensity cube.
intensity = []
for y in years:
    cube = ndvi_stack(y).interpolate_na("time").bfill("time").ffill("time")
    n = xr.apply_ufunc(cycles_per_year, cube, input_core_dims=[["time"]],
                       vectorize=True, dask="allowed")
    intensity.append(n)
intensity = xr.concat(intensity, dim=xr.DataArray(list(years), dims="year"))

# 3. Trend: regress cycles-per-year vs year per pixel; positive slope = shifting
#    toward double/triple cropping. Region mean gives the headline verdict.
yr = np.array(list(years), float)
mean_int = intensity.mean(dim=("y", "x"))
slope = np.polyfit(yr, mean_int.values, 1)[0]
verdict = "intensifying" if slope > 0.02 else "stable"
print(f"Punjab cropping intensity: {mean_int.values[0]:.2f} -> "
      f"{mean_int.values[-1]:.2f} crops/yr, slope {slope:+.3f}/yr ({verdict})")

# 4. Plot the regional cropping-intensity trend.
plt.plot(yr, mean_int.values, "o-")
plt.ylabel("mean crops / year"); plt.xlabel("year")
plt.title("Cropping intensity trend"); plt.tight_layout(); plt.show()

# Before you trust it: mask to cropland (forest/grass green up too) and cross-check
# the trend against the district agricultural census (see the gotchas above).

Where the data comes from

The whole core workflow comes from NASA through one free login (Earthdata Login). The MODIS greenness (MOD13) and the sharper HLS data both live at the LP DAAC archive, so there’s a single sign-in to manage. The only extra piece is an optional cropland mask (a national land-cover map) to screen out forest and grass, which may come from a separate source.

Sources

• MODIS MOD13 VI: https://lpdaac.usgs.gov/products/mod13q1v061/
• HLS: https://lpdaac.usgs.gov/products/hlss30v002/
• MODIS VI User Guide: https://lpdaac.usgs.gov/documents/621/MOD13_User_Guide_V61.pdf
• Cropping-intensity background (NDVI cycle counting): https://lpdaac.usgs.gov/