Time Series Analysis of Satellite Data
Extracting trends, seasonality, and phenology from multi-temporal imagery
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.
1. The Question
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
2. The Conceptual Model
Time Series Components
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
Phenological Metrics
Key dates in annual cycle:
- Start of Season (SOS): When NDVI crosses threshold (spring green-up)
- Peak of Season (POS): Maximum NDVI (summer)
- End of Season (EOS): When NDVI falls below threshold (senescence)
- Growing Season Length (GSL): EOS - SOS
Amplitude: Difference between peak and minimum NDVI
Integrated NDVI: Area under curve (proportional to productivity)
Harmonic Analysis
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.
One NDVI Record Carries Trend, Seasonality, And Events
The main analytical move is to separate the repeating seasonal wave from the slower long-term drift and the short-lived anomalies riding on top of both.
3. Building the Mathematical Model
Moving Average Smoothing
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)
Savitzky-Golay Filter
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)
Trend Extraction
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.
Harmonic Regression
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)
Phenology Extraction
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
4. Worked Example by Hand
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
Solution
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)
5. Computational Implementation
Below is an interactive time series analyzer.
Start of Season: --
Peak of Season: --
End of Season: --
Growing Season Length: -- days
Trend: -- NDVI/year
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!
6. Interpretation
Climate Change Signals
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
Drought Monitoring
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
Agricultural Phenology
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
7. What Could Go Wrong?
Cloud Contamination
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
Missing Data
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)
Sensor Drift
Long-term sensors (AVHRR, MODIS) degrade over time.
Calibration drift → false trend.
Solution:
- Cross-calibration between sensors
- Atmospheric correction updates
- Use surface reflectance products
Mixed Land Cover
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
8. Extension: TIMESAT
Software package for time series processing.
Functions:
- Gap-filling (interpolation)
- Smoothing (Savitzky-Golay, asymmetric Gaussian)
- Phenology extraction (multiple methods)
- Visualization
Three fitting methods:
- Asymmetric Gaussian: Separate rise/fall curves
- Double logistic: S-curves for green-up/senescence
- 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.
9. Math Refresher: Fourier Series
Periodic Function Decomposition
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
Application to NDVI
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)
Summary
- Time series analysis decomposes observations into trend, seasonality, and noise
- Phenological metrics: Start of season (SOS), peak (POS), end (EOS), growing season length (GSL)
- Harmonic analysis: Models seasonality using sine/cosine functions
- Trend extraction: Linear regression to detect long-term change
- Smoothing filters: Moving average or Savitzky-Golay to reduce noise
- Applications: Climate change detection, drought monitoring, crop calendars, forest health
- Challenges: Clouds, missing data, sensor drift, mixed pixels
- TIMESAT: Standard software for operational phenology extraction
- Global observations: Growing season lengthening, Arctic greening, NDVI anomalies indicate drought
- Critical for understanding vegetation dynamics and climate change impacts