Extracting trends, seasonality, and phenology from multi-temporal imagery
2026-02-27
Landsat 5, launched in 1984 and operational until 2013, accumulated 29 years of continuous Earth observation — 2.5 million scenes, one every 16 days at 30-metre resolution. For any specific location, that record is a time series: a sequence of NDVI or reflectance values sampled through almost three decades of seasons, droughts, fires, floods, and land-use changes. This is a remarkable scientific resource. It means that for a forest stand in Saskatchewan, or a wetland on the Mackenzie Delta, or a cropping district in southern Alberta, we can reconstruct the detailed seasonal and interannual trajectory of vegetation greenness since the mid-1980s. Embedded in that trajectory are signals of browning during the 2002 drought, accelerated green-up as spring arrives earlier, and the distinctive dip-and-recovery signature of a bark beetle outbreak.
Time series analysis of satellite data extracts those signals from what is otherwise a noisy record of spectral reflectance. The key decomposition separates a time series into trend (long-term directional change), seasonality (annual cycle), and residual (anomalies and noise). The seasonal component — the dominant signal for most vegetated surfaces — can be modelled as a sum of sinusoids fitted by Fourier decomposition, enabling precise extraction of phenological metrics: the date of green-up, the peak NDVI, the date of senescence, and the integrated annual productivity. The trend component, once seasonality is removed, reveals whether long-term vegetation conditions are improving or declining. This model derives the decomposition framework, implements harmonic regression, and extracts phenological metrics from MODIS time series data.
How has the growing season changed over the past 20 years?
Time series analysis examines sequences of observations through time:
Satellite time series: - MODIS: 500m, 16-day composites, 2000-present - Landsat: 30m, 16-day revisit, 1984-present - Sentinel-2: 10m, 5-day revisit, 2015-present
Applications: - Phenology: Timing of green-up, peak, senescence - Productivity trends: Is vegetation increasing or decreasing? - Drought detection: NDVI anomalies below normal - Crop monitoring: Within-season growth patterns - Climate change: Lengthening/shortening growing seasons - Forest health: Detecting gradual decline vs. abrupt disturbance
The mathematical question: Given a noisy time series of NDVI values, how do we separate: 1. Trend: Long-term change 2. Seasonality: Repeating annual cycle 3. Noise: Random variation and clouds
Additive model:
y(t) = T(t) + S(t) + \varepsilon(t)
Where: - y(t) = observed value at time t - T(t) = trend component - S(t) = seasonal component - \varepsilon(t) = noise/residual
Example - NDVI time series: - Trend: Gradual greening (+0.001 NDVI/year) - Seasonal: Spring green-up, summer peak, autumn senescence - Noise: Clouds, atmospheric effects, sensor variations
Key dates in annual cycle:
Amplitude: Difference between peak and minimum NDVI
Integrated NDVI: Area under curve (proportional to productivity)
Model seasonal cycle with sine/cosine:
S(t) = A \sin(2\pi t / P + \phi)
Where: - A = amplitude - P = period (365 days) - \phi = phase (timing of peak)
Multiple harmonics capture complex seasonal patterns:
S(t) = \sum_{k=1}^{n} \left[a_k \cos(2\pi k t / P) + b_k \sin(2\pi k t / P)\right]
k=1: Annual cycle
k=2: Bi-annual (two peaks per year)
k=3: Tri-annual, etc.
Simple moving average (window size w):
\hat{y}_t = \frac{1}{w} \sum_{i=t-w/2}^{t+w/2} y_i
Reduces noise but blurs edges.
Typical: w = 3 (adjacent observations)
Fits polynomial to local window, evaluates at center.
Preserves peaks better than moving average.
Algorithm: For each point, fit polynomial of degree d to w neighbors, use fitted value.
Typical: w = 5, d = 2 (quadratic)
Linear trend:
T(t) = \beta_0 + \beta_1 t
Fit via least squares:
\beta_1 = \frac{\sum (t_i - \bar{t})(y_i - \bar{y})}{\sum (t_i - \bar{t})^2}
\beta_0 = \bar{y} - \beta_1 \bar{t}
Interpretation: - \beta_1 > 0: Increasing trend (greening) - \beta_1 < 0: Decreasing trend (browning) - \beta_1 = 0: No trend
Mann-Kendall test: Non-parametric trend significance test.
Fit annual cycle:
y(t) = \beta_0 + a_1\cos(\omega t) + b_1\sin(\omega t) + \varepsilon(t)
Where \omega = 2\pi / 365 (radians/day).
Amplitude:
A = \sqrt{a_1^2 + b_1^2}
Phase (peak timing):
\phi = \arctan(b_1 / a_1)
Day of year for peak:
\text{DOY}_{\text{peak}} = \phi \times \frac{365}{2\pi}
Add trend:
y(t) = \beta_0 + \beta_1 t + a_1\cos(\omega t) + b_1\sin(\omega t)
Threshold method:
Start of Season: First day when \text{NDVI} > T_{\text{green}}
End of Season: Last day when \text{NDVI} > T_{\text{green}}
Typical threshold: T_{\text{green}} = 0.3 or T_{\text{green}} = \text{NDVI}_{\min} + 0.2(\text{NDVI}_{\max} - \text{NDVI}_{\min})
Derivative method:
SOS: Maximum rate of increase
\frac{d\text{NDVI}}{dt} = \max
EOS: Maximum rate of decrease
\frac{d\text{NDVI}}{dt} = \min
Problem: Extract phenology from simplified annual NDVI time series.
Monthly NDVI values (2023):
| Month | DOY | NDVI |
|---|---|---|
| Jan | 15 | 0.25 |
| Feb | 45 | 0.28 |
| Mar | 75 | 0.35 |
| Apr | 105 | 0.50 |
| May | 135 | 0.68 |
| Jun | 165 | 0.75 |
| Jul | 195 | 0.72 |
| Aug | 225 | 0.65 |
| Sep | 255 | 0.50 |
| Oct | 285 | 0.35 |
| Nov | 315 | 0.28 |
| Dec | 345 | 0.25 |
Find: SOS, POS, EOS, GSL using threshold T = 0.4
Step 1: Find Start of Season
First month when NDVI > 0.4: - Jan: 0.25 < 0.4 - Feb: 0.28 < 0.4 - Mar: 0.35 < 0.4 - Apr: 0.50 > 0.4 ✓
SOS ≈ DOY 105 (mid-April)
Step 2: Find Peak of Season
Maximum NDVI: - Max = 0.75 in June
POS = DOY 165 (mid-June)
Step 3: Find End of Season
Last month when NDVI > 0.4: - Sep: 0.50 > 0.4 ✓ - Oct: 0.35 < 0.4
EOS ≈ DOY 255 (mid-September)
Step 4: Growing Season Length
\text{GSL} = \text{EOS} - \text{SOS} = 255 - 105 = 150 \text{ days}
Step 5: Integrated NDVI (productivity)
Trapezoidal integration (approximate area under curve):
\int \text{NDVI} \, dt \approx \sum_{i=1}^{n-1} \frac{\text{NDVI}_i + \text{NDVI}_{i+1}}{2} \times \Delta t
\approx \frac{0.25+0.28}{2}\times30 + \frac{0.28+0.35}{2}\times30 + \cdots
\approx 30 \times (0.265 + 0.315 + 0.425 + 0.590 + 0.715 + 0.735 + 0.685 + 0.575 + 0.425 + 0.315 + 0.265)
\approx 30 \times 5.31 = 159.3 \text{ NDVI-days}
Summary: - SOS: DOY 105 - POS: DOY 165
- EOS: DOY 255 - GSL: 150 days - Integrated NDVI: ~159 (productivity proxy)
Below is an interactive time series analyzer.
<label>
Location type:
<select id="location-type">
<option value="temperate" selected>Temperate Forest</option>
<option value="grassland">Grassland</option>
<option value="tropical">Tropical (Evergreen)</option>
<option value="cropland">Cropland (Double-crop)</option>
</select>
</label>
<label>
Show components:
<input type="checkbox" id="show-trend" checked> Trend
<input type="checkbox" id="show-seasonal" checked> Seasonal
<input type="checkbox" id="show-smoothed"> Smoothed
</label>
<label>
Phenology threshold:
<input type="range" id="pheno-threshold" min="0.2" max="0.6" step="0.05" value="0.4">
<span id="pheno-threshold-val">0.40</span>
</label>
<div class="timeseries-info">
<p><strong>Start of Season:</strong> <span id="sos-doy">--</span></p>
<p><strong>Peak of Season:</strong> <span id="pos-doy">--</span></p>
<p><strong>End of Season:</strong> <span id="eos-doy">--</span></p>
<p><strong>Growing Season Length:</strong> <span id="gsl-days">--</span> days</p>
<p><strong>Trend:</strong> <span id="trend-value">--</span> NDVI/year</p>
</div><canvas id="timeseries-canvas" width="700" height="400" style="border: 1px solid #ddd;"></canvas>Try this: - Temperate forest: Strong seasonal cycle, single peak - Grassland: Similar but smaller amplitude - Tropical evergreen: Minimal seasonality (always green) - Cropland: Double-crop pattern (two peaks per year!) - Show trend: Red dashed line (long-term change) - Show seasonal: Green line (trend + seasonality) - Show smoothed: Blue line (noise reduced) - Adjust threshold: Changes SOS/EOS detection - Phenology markers: Green (SOS), Purple (POS), Orange (EOS)
Key insight: Time series decomposition separates long-term trends from seasonal cycles, enabling climate change detection and phenology monitoring!
Global greening trend: - NDVI increased ~0.001-0.002/year (1982-2015) - Reasons: CO₂ fertilization, warmer temperatures, longer growing seasons - Arctic: Strongest greening (shrub expansion)
Growing season extension: - SOS earlier (1-2 weeks in 40 years) - EOS later (1 week) - Net: +10-20 days GSL in northern latitudes
NDVI anomaly:
A(t) = \text{NDVI}(t) - \text{NDVI}_{\text{climatology}}(t)
Where climatology = long-term average for that time of year.
Negative anomaly = drought stress
Example: 2012 US drought - NDVI anomaly < -0.15 across Corn Belt - Early warning of crop failure
Crop calendars from satellites: - SOS = planting date - POS = peak vegetative growth - EOS = harvest date
Applications: - Yield forecasting (early prediction) - Insurance claims (verify timing) - Double-cropping detection
Clouds reduce NDVI (appear as low values).
Spurious phenology: - Cloudy spring → delayed SOS (false) - Need cloud masking
Solutions: - Composite images (maximum NDVI in period) - Cloud screening (quality flags) - Gap-filling interpolation
Landsat: 16-day repeat but clouds may leave gaps.
Irregular time series → difficult to analyze.
Solutions: - Interpolation (linear, spline) - Harmonic fitting (smooth through gaps) - Data fusion (combine Landsat + Sentinel)
Long-term sensors (AVHRR, MODIS) degrade over time.
Calibration drift → false trend.
Solution: - Cross-calibration between sensors - Atmospheric correction updates - Use surface reflectance products
Pixel contains multiple types (e.g., 50% forest, 50% grass).
Phenology represents blend → hard to interpret.
Solution: - Use finer resolution - Classify first, then analyze per class - Sub-pixel unmixing
Software package for time series processing.
Functions: - Gap-filling (interpolation) - Smoothing (Savitzky-Golay, asymmetric Gaussian) - Phenology extraction (multiple methods) - Visualization
Three fitting methods: 1. Asymmetric Gaussian: Separate rise/fall curves 2. Double logistic: S-curves for green-up/senescence 3. Savitzky-Golay filter: Polynomial smoothing
Outputs: - Start/end of season - Length of season - Base/peak/amplitude - Small/large integrals (productivity metrics)
Widely used in operational phenology monitoring.
Any periodic function with period P can be represented as:
f(t) = a_0 + \sum_{k=1}^{\infty} \left[a_k\cos(k\omega t) + b_k\sin(k\omega t)\right]
Where \omega = 2\pi / P (fundamental frequency).
Coefficients:
a_0 = \frac{1}{P}\int_0^P f(t)\,dt
a_k = \frac{2}{P}\int_0^P f(t)\cos(k\omega t)\,dt
b_k = \frac{2}{P}\int_0^P f(t)\sin(k\omega t)\,dt
For annual cycle (P = 365 days):
First harmonic (k=1): Annual component
Second harmonic (k=2): Bi-annual (captures double-cropping)
Typically: 1-3 harmonics sufficient for most vegetation.
Amplitude of harmonic k:
A_k = \sqrt{a_k^2 + b_k^2}
Phase (timing of peak):
\phi_k = \arctan(b_k / a_k)