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

What you can answer

• Cropping intensity per year — single / double / triple cropping from counted NDVI cycles
• Change in intensity over time — where cycles-per-year are rising across the record
• Spatial map of intensification — which sub-regions added a season
• Crop-calendar timing — when each green-up and senescence occurs through the year
• Coarse vs fine comparison — MODIS 250 m trend vs HLS 30 m for fragmented fields

What you can NOT answer with these alone

• Which crop is grown — NDVI cycles count seasons, not species; pair with crop-type maps or ground truth.
• Exact field boundaries at MODIS scale — 250 m pixels mix multiple smallholder plots; use HLS/Sentinel-2 30 m.
• Irrigated vs rain-fed cause without ancillary data — a second cycle implies irrigation/residual moisture but doesn’t prove the water source.
• Yield or production — cycle count is intensity, not output; combine with yield models and census data.

Code template

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

earthaccess.login(strategy="netrc")

# AOI: Punjab
aoi = (74, 29.5, 77, 32)
years = range(2010, 2026)

# Greenness / amplitude rules — tune and verify against local crop calendars
ndvi_threshold = 0.4    # 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 NDVI series (already time-ordered)."""
    smooth = savgol_filter(ndvi_1d, window_length=5, polyorder=2)
    peaks, _ = find_peaks(smooth, height=ndvi_threshold,
                          prominence=min_prominence)
    return len(peaks)

# 1. MODIS NDVI (16-day, 250 m) across the year range
granules = earthaccess.search_data(short_name="MOD13Q1",
                                   bounding_box=aoi,
                                   temporal=(f"{min(years)}-01-01",
                                             f"{max(years)}-12-31"))
# Open NDVI (scale ×0.0001), apply the VI QA/pixel-reliability mask,
# stack to a (time, y, x) DataArray, then per pixel per crop-year:
#   intensity[year] = cycles_per_year(ndvi_series_for_that_year)

# 2. Trend: regress intensity vs year per pixel → map of rising cropping
#    (positive slope = shifting toward double/triple cropping).

# 3. (Optional) HLS 30 m for smallholder fields where 250 m mixes plots
#    short_name="HLSS30" / "HLSL30"; same cycle-counting on a finer grid.

# 4. Mask to cropland (e.g. a land-cover product) so forest/grass green-up
#    is not miscounted as a crop cycle.

Expected output

• Per-year cropping-intensity map (1 / 2 / 3 crops) from counted NDVI cycles
• Intensity-trend map — slope of cycles-per-year over the record (intensifying vs stable)
• NDVI-vs-time curves for sample pixels showing one, two, or three green waves
• Summary: share of the AOI that gained a cropping season over the period

Caveats

• Mixed pixels at 250 m — a MODIS pixel can blend several fields with different calendars, smearing or double-counting cycles; use HLS 30 m for smallholder landscapes.
• Natural vegetation greens up too — forest, grassland, and weeds produce NDVI cycles; mask to cropland before counting.
• Cloud/monsoon gaps — kharif-season cloud cover thins the NDVI record; smooth and gap-fill, and use pixel-reliability QA.
• Threshold/amplitude rules are tunable, not universal — the greenness cutoff and prominence depend on crop and region; calibrate locally rather than asserting one value.
• Verify the shift against the census — cross-check any intensification trend with agricultural-census / district crop statistics before claiming a real change.

Cross-DAAC composition

LP DAAC only for the core workflow (MODIS MOD13, HLS) — single DAAC, single Earthdata auth. Optional cropland masks (e.g. national land-cover products) may add a second 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/