Harmonic regression on a seasonal raster

remote sensing

time series

terra

ecology tutorial

Fit a sine and cosine pair per pixel to a year of NDVI composites in R, map amplitude and peak day, and measure what the raw seasonal maximum gets wrong.

Author

Tidy Ecology

Published

2026-08-03

A window 1200 metres on a side over a lowland mosaic: rough grassland on free-draining gravels in the south, a block of oak and hornbeam in the north, and a ragged band of encroaching scrub between them. The satellite archive holds one year of 8-day NDVI composites over that window, 46 layers, on a 30 metre grid, so the cube is 40 pixels by 40 pixels with 46 observations at each. It is a deliberately small window: everything below runs on a laptop in under a minute.

The question is not what the NDVI was on any particular date. It is what the seasonal cycle looks like at each pixel: how far the vegetation swings between its winter floor and its summer ceiling, and when it peaks. Those two numbers, an amplitude and a phase, are what goes into the habitat model or the comparison between management compartments, and they have to be recovered from 46 noisy values per pixel, with clouds deleting some of them and the composites eight days apart to begin with. The obvious way to get the peak is to ask which layer is largest, which is one line of code and what most people write first. It is also the estimator this post is mainly about, because a fitted seasonal curve beats it by a factor of four in these data.

NDVI time series from a raster stack built a cube of this shape and took per-pixel summaries off it: means, maxima, an integral over the season. This post fits a model to the time axis instead, which is a different object, because a summary throws away the shape and a model parametrises it. Raster data in R with terra covers single-layer work and is not repeated. Filling cloud gaps in a satellite series is about reconstructing missing observations; the section on missing dates below asks only what an unfilled gap does to a fitted coefficient. Checking a remote sensing covariate takes the outputs of this post and asks whether they deserve to be in a model at all. Two older posts set the other boundaries. Phenological trends and temperature estimates timing from ground observation: somebody walked the transect and wrote down a date, and here the same quantity is read off a reflectance curve by a machine that never visited the site. And circular data and the von Mises distribution owns the arithmetic of angles, which matters because a phase is an angle and forgetting that is one of the standard mistakes in this workflow.

Six things get measured: how well amplitude, phase and mean are recovered when the truth is known; how the answer moves as harmonics are added; how far the raw seasonal maximum is from the truth compared with two fitted alternatives; what an unevenly missing set of dates does; which direction a season-shaped gap pushes the amplitude; and how large the per-pixel standard error on the amplitude is before anybody carries it into another model.

library(ggplot2)

te_pal <- list(forest = "#275139", green = "#2f8f63", sage = "#93a87f",
               clay = "#b5534e", gold = "#cda23f", line = "#dad9ca",
               ink = "#16241d", paper = "#f5f4ee")

theme_te <- function() {
  theme_minimal(base_size = 12) +
    theme(panel.grid.minor = element_blank(),
          panel.grid.major = element_line(colour = "#e7e6dc"),
          plot.background = element_rect(fill = "#f5f4ee", colour = NA),
          panel.background = element_rect(fill = "#f5f4ee", colour = NA),
          plot.title = element_text(face = "bold", colour = te_pal$ink),
          axis.title = element_text(colour = "#2c3a31"),
          axis.text = element_text(colour = "#2c3a31"))
}

A cube with a known answer inside it

Simulated data, because the point is to measure error and that needs a truth. Each pixel gets a mean NDVI, an amplitude, a peak day and a small second harmonic, all varying smoothly across the window through a woodland index running from open grass in the south to closed canopy in the north. The second harmonic is what makes the curve asymmetric: green-up in a deciduous canopy is fast and senescence slow, and a single cosine cannot express that.

suppressPackageStartupMessages(library(terra))

n_side <- 40
n_cell <- n_side^2
n_lyr <- 46
doy <- seq(1, 361, by = 8)
omega <- 2 * pi / 365

tmpl <- rast(nrows = n_side, ncols = n_side, xmin = 0, xmax = 1200,
             ymin = 0, ymax = 1200, crs = "EPSG:32634")
xy_cell <- crds(tmpl)
g_east <- xy_cell[, 1] / 1200; g_north <- xy_cell[, 2] / 1200
wood <- plogis(7 * (g_north - 0.5 + 0.28 * sin(3.1 * g_east)))
mean_true <- 0.34 + 0.16 * wood
amp_true <- 0.13 + 0.17 * wood + 0.02 * g_east
peak_par <- 150 + 36 * wood + 10 * g_north
amp2_true <- 0.008 + 0.022 * wood

day_fine <- 1:365
ang_fine <- omega * (matrix(day_fine, n_cell, 365, byrow = TRUE) - peak_par)
sig_fine <- mean_true + amp_true * cos(ang_fine) + amp2_true * sin(2 * ang_fine)
peak_true <- day_fine[max.col(sig_fine, ties.method = "first")]

ang_obs <- omega * (matrix(doy, n_cell, n_lyr, byrow = TRUE) - peak_par)
sig_obs <- mean_true + amp_true * cos(ang_obs) + amp2_true * sin(2 * ang_obs)

sd_noise <- 0.04
set.seed(20260803)
obs_mat <- sig_obs + matrix(rnorm(n_cell * n_lyr, 0, sd_noise), n_cell, n_lyr)

ndvi <- setValues(rast(tmpl, nlyr = n_lyr), obs_mat)
names(ndvi) <- sprintf("d%03d", doy)

print(c(pixels = ncell(ndvi), layers = nlyr(ndvi),
        composite_days = doy[2] - doy[1], first_doy = doy[1],
        last_doy = doy[n_lyr]))

        pixels         layers composite_days      first_doy       last_doy 
          1600             46              8              1            361

print(round(c(mean_lo = min(mean_true), mean_hi = max(mean_true),
              amp_lo = min(amp_true), amp_hi = max(amp_true),
              peak_lo = min(peak_true), peak_hi = max(peak_true),
              noise_sd = sd_noise), 4))

 mean_lo  mean_hi   amp_lo   amp_hi  peak_lo  peak_hi noise_sd 
  0.3455   0.4993   0.1361   0.3153 158.0000 206.0000   0.0400

The cube holds 1600 pixels and 46 layers, one composite every 8 days from day 1 to day 361. True amplitude runs from 0.136 in the driest grass to 0.315 in the closed canopy, and the true peak day, defined as the argmax of the noiseless annual curve, runs from day 158 to day 206. Per-observation noise has a standard deviation of 0.04, which is an ordinary figure for a composited vegetation index once residual cloud, aerosol and view-angle effects are in it.

One terra habit is worth stating first, because it fails silently. A SpatRaster plus a bare numeric vector recycles across layers, not across cells, so adding a 1600-element per-pixel vector to a 46-layer cube produces a cube with hundreds of layers and no error. The safe route is the one used above: build the values as a matrix of cells by layers, hand it to setValues(), check nlyr().

One pixel, one linear model

Write $y_t$ for the NDVI of a pixel at day of year $t$ and $\omega = 2\pi/365$. The one-harmonic model is

\[y_t = \mu + a\cos(\omega t) + b\sin(\omega t) + \varepsilon_t\]

which looks like a curve-fitting problem and is not one. The pair $\cos(\omega t)$ and $\sin(\omega t)$ are known numbers once the dates are known, so $\mu$, $a$ and $b$ enter linearly and the whole thing is lm(). No optimiser, no starting values, no convergence to worry about. That is why harmonic regression is the workhorse in this corner of remote sensing. Jakubauskas, Legates and Kastens (2002) used the amplitudes and phases of the first few harmonics of an AVHRR NDVI series as classification features, the direct ancestor of using them as ecological covariates; Zhu and Woodcock (2014) put the same form at the centre of continuous change detection, fitting it to a stable period and flagging a pixel when new observations fall outside the prediction interval, which only works because refitting everywhere is cheap; and Verbesselt, Hyndman, Newnham and Culvenor (2010) take the complementary route in BFAST, separating a seasonal component from a piecewise trend and testing for breaks in each.

The two coefficients are then rewritten as an amplitude and a phase, because that is what anybody wants to look at:

\[a\cos(\omega t) + b\sin(\omega t) = A\cos(\omega t - \phi), \qquad A = \sqrt{a^2 + b^2}, \qquad \phi = \operatorname{atan2}(b, a)\]

and the phase converts to a day of year as $\phi/\omega$, taken modulo 365. Adding more harmonics means adding $\cos(2\omega t)$, $\sin(2\omega t)$ and so on: still linear, still lm().

harm_design <- function(days, n_harm) {
  des <- matrix(1, length(days), 1 + 2 * n_harm)
  colnames(des) <- c("mean", as.vector(rbind(paste0("cos", seq_len(n_harm)),
                                             paste0("sin", seq_len(n_harm)))))
  for (j in seq_len(n_harm)) {
    des[, 2 * j] <- cos(2 * pi * j * days / 365)
    des[, 2 * j + 1] <- sin(2 * pi * j * days / 365)
  }
  des
}

pix_id <- cellFromXY(tmpl, cbind(315, 915))
pix_dat <- data.frame(ndvi = obs_mat[pix_id, ], cos1 = cos(omega * doy),
                      sin1 = sin(omega * doy))
fit_pix <- lm(ndvi ~ cos1 + sin1, data = pix_dat)
print(round(summary(fit_pix)$coefficients, 5))

            Estimate Std. Error   t value Pr(>|t|)
(Intercept)  0.49198    0.00618  79.61077        0
cos1        -0.28323    0.00870 -32.54016        0
sin1        -0.05196    0.00877  -5.92202        0

cf_pix <- coef(fit_pix)
amp_pix <- sqrt(cf_pix[2]^2 + cf_pix[3]^2)
phase_pix <- (atan2(cf_pix[3], cf_pix[2]) / omega) %% 365
lmfit_pix <- lm.fit(harm_design(doy, 1), obs_mat[pix_id, ])$coefficients

print(round(c(intercept = unname(cf_pix[1]), true_mean = mean_true[pix_id],
              amplitude = unname(amp_pix), true_amplitude = amp_true[pix_id],
              phase_day = unname(phase_pix), true_peak_day = peak_true[pix_id],
              woodland_index = wood[pix_id]), 4))

     intercept      true_mean      amplitude true_amplitude      phase_day 
        0.4920         0.4941         0.2880         0.2990       193.0409 
 true_peak_day woodland_index 
      203.0000         0.9631

print(round(c(max_gap_lm_vs_lmfit = max(abs(cf_pix - lmfit_pix))), 12))

max_gap_lm_vs_lmfit 
                  0

At this pixel the intercept comes out at 0.492 against a true mean of 0.4941, and the amplitude at 0.288 against a true 0.299. The phase converts to day 193 while the true peak of the pixel’s curve is day 203, a gap of 10 days that is not noise and will be taken apart later. lm() and lm.fit() on the same design return identical coefficients, which matters only because the per-pixel loop below uses the second one for speed.

The phase is an angle, and angles do not average

Once a phase has been converted to a day of year it looks like an ordinary number, and the moment somebody averages peak days over a management unit or differences them between years, the branch cut bites. Day 362 and day 3 are six days apart, not 359.

wrap_days <- c(356, 362, 3, 9, 14)
naive_mean <- mean(wrap_days)
wrap_ang <- omega * wrap_days
vec_mean <- (atan2(mean(sin(wrap_ang)), mean(cos(wrap_ang))) / omega) %% 365
resultant <- sqrt(mean(cos(wrap_ang))^2 + mean(sin(wrap_ang))^2)

print(round(c(arithmetic_mean_day = naive_mean, vector_mean_day = vec_mean,
              resultant_length = resultant, error_days = naive_mean - vec_mean), 3))

arithmetic_mean_day     vector_mean_day    resultant_length          error_days 
            148.800               2.802               0.990             145.998

Five pixels whose peaks straddle the new year have an arithmetic mean peak day of 148.8, which is in late May. The vector mean, taken by averaging the unit vectors and converting back with atan2(), is day 2.8, and it is the right answer: the arithmetic version is off by 146 days. The resultant length 0.99 measures the agreement, near one when the phases coincide and near zero when they cancel. Fisher (1993) is the standard treatment and the base R mechanics are in the circular post. Nothing in this window peaks near the new year, so the problem never surfaces here, but it is exactly what a semi-arid system with a winter growing season does, and it fails silently when it does.

Three panels side by side, each with day of year from zero to three hundred and sixty-five on the horizontal axis and NDVI from about two tenths to eight tenths on the vertical axis. In every panel forty-six round dots scatter closely around two nearly coincident curves, one pale sage and one dark green, that rise from a flat winter floor near 0.2 to a single summer hump and fall again. The left panel, labelled open grassland, has the shallowest hump, cresting near day one hundred and seventy at about 0.55. The middle panel, the scrub edge, crests near day one hundred and ninety-five at about 0.76. The right panel, closed canopy, has the tallest hump, cresting just past day two hundred at about 0.8, and its rise is steeper than its fall. — Figure 1: Three pixels from the window: open grassland, the scrub edge and closed canopy. Dots are the 46 composite observations, the pale line is the noiseless annual curve the observations were drawn from, and the dark line is the two-harmonic fit. The grassland pixel swings less and peaks earlier; the canopy pixel swings more, peaks later and is visibly asymmetric, rising faster than it falls.

Fitting every pixel with app()

terra::app() applies a function to the vector of layer values at each cell, which is exactly the shape a per-pixel time series has. The function returns a vector of coefficients, so the result is a multi-layer raster with one layer per coefficient. Two harmonics gives five: a mean and two cosine and sine pairs.

des2 <- harm_design(doy, 2)
harm_r <- app(ndvi, function(v) lm.fit(des2, v)$coefficients)
names(harm_r) <- colnames(des2)
print(harm_r)

class       : SpatRaster
size        : 40, 40, 5  (nrow, ncol, nlyr)
resolution  : 30, 30  (x, y)
extent      : 0, 1200, 0, 1200  (xmin, xmax, ymin, ymax)
coord. ref. : WGS 84 / UTM zone 34N (EPSG:32634)
source(s)   : memory
names       :     mean,      cos1,      sin1,      cos2,      sin2
min values  : 0.330629, -0.327367, -0.083869, -0.037725, -0.018976
max values  : 0.516825, -0.104717,   0.09443,  0.033952,  0.052607

cf2 <- values(harm_r)
cf_direct <- t(qr.coef(qr(des2), t(obs_mat)))
print(c(coef_layers = nlyr(harm_r),
        agrees_with_matrix_solve = max(abs(cf2 - cf_direct)) < 1e-12))

             coef_layers agrees_with_matrix_solve 
                       5                        1

The output carries 5 layers and reproduces a single matrix solve of the same design exactly, which is worth checking once and then relying on: the sweeps later in the post use the matrix form, because running app() inside a replicate loop would spend its time on raster overheads for no gain.

fine2 <- harm_design(day_fine, 2)
circ_gap <- function(a, b) ((a - b + 182.5) %% 365) - 182.5

amp_hat <- sqrt(cf2[, "cos1"]^2 + cf2[, "sin1"]^2)
phase_hat <- (atan2(cf2[, "sin1"], cf2[, "cos1"]) / omega) %% 365
peak_hat <- day_fine[max.col(cf2 %*% t(fine2), ties.method = "first")]

amp_r <- setValues(rast(tmpl), amp_hat)
peak_r <- setValues(rast(tmpl), peak_hat)
print(c(amp_layers = nlyr(amp_r), peak_layers = nlyr(peak_r)))

 amp_layers peak_layers 
          1           1

amp_err <- amp_hat - amp_true
peak_err <- circ_gap(peak_hat, peak_true)
mean_err <- cf2[, "mean"] - mean_true
amp_mae <- median(abs(amp_err))
amp_p90 <- unname(quantile(abs(amp_err), 0.9))
peak_mae <- median(abs(peak_err))
peak_p90 <- unname(quantile(abs(peak_err), 0.9))

print(round(c(mean_median_abs_err = median(abs(mean_err)),
              amp_median_abs_err = amp_mae, amp_p90_abs_err = amp_p90,
              amp_median_pct_err = 100 * median(abs(amp_err / amp_true)),
              amp_bias = mean(amp_err)), 5))

mean_median_abs_err  amp_median_abs_err     amp_p90_abs_err  amp_median_pct_err 
            0.00391             0.00555             0.01370             2.20723 
           amp_bias 
           -0.00033

print(round(c(peak_median_abs_err_days = peak_mae, peak_p90_days = peak_p90,
              peak_bias_days = mean(peak_err),
              peak_worst_days = max(abs(peak_err))), 3))

peak_median_abs_err_days            peak_p90_days           peak_bias_days 
                   3.000                    7.000                   -0.072 
         peak_worst_days 
                  17.000

Against a known truth the two-harmonic fit recovers the annual mean to a median absolute error of 0.0039 NDVI units and the amplitude to 0.0056, or 2.21 per cent of the amplitude being estimated, with a ninetieth percentile of 0.0137: a tail about two and a half times the median rather than an order of magnitude away. Peak day comes back with a median absolute error of 3 days, a ninetieth percentile of 7 days and a mean signed error of -0.07 days, the worst pixel being 17 days out.

Two square maps side by side on warm off-white paper, each twelve hundred metres across with easting and northing axes in metres. The left map shows the seasonal amplitude and grades smoothly from pale sage along the bottom edge, where the legend reads about 0.15, to dark forest green across the top at about 0.30. The pale zone reaches highest up the left and right edges and is squeezed lowest across the middle of the panel, so the transition bows downwards. The right map shows the peak day of the growing season on the same shape, gold at about one hundred and sixty days along the bottom and dark green at about two hundred days along the top. Both maps carry a fine grainy speckle everywhere, which is the pixel-to-pixel estimation noise. — Figure 2: The two coefficient maps that a downstream model would actually use, both estimated pixel by pixel from 46 noisy composites. Left, the amplitude of the first harmonic: the seasonal swing, largest in the closed canopy in the north of the window and smallest over the grassland in the south, with a broad transition where the scrub is. Right, the peak day taken as the argmax of the fitted two-harmonic curve, which follows the same pattern and spans about seven weeks across the window. The speckle in both is estimation noise, not vegetation.

How many harmonics

One harmonic is a symmetric hump. Two can lean. Three can carry a shoulder. The usual way to choose is a fit statistic, so start there and then ask whether the fit statistic is answering the question.

fit_harm <- function(n_harm, keep = seq_len(n_lyr), y_mat = obs_mat) {
  des <- harm_design(doy[keep], n_harm)
  yk <- y_mat[, keep, drop = FALSE]
  cf <- t(qr.coef(qr(des), t(yk)))
  rss <- rowSums((yk - cf %*% t(des))^2)
  list(cf = cf, des = des, rss = rss, npar = ncol(des),
       amp = sqrt(cf[, 2]^2 + cf[, 3]^2),
       peak = day_fine[max.col(cf %*% t(harm_design(day_fine, n_harm)),
                               ties.method = "first")],
       rse = sqrt(rss / (length(keep) - ncol(des))),
       r2 = 1 - rss / rowSums((yk - rowMeans(yk))^2))
}

k_grid <- 1:3
scan_out <- lapply(k_grid, fit_harm)
k_tab <- data.frame(
  harmonics = k_grid,
  parameters = sapply(scan_out, function(s) s$npar),
  resid_sd = sapply(scan_out, function(s) median(s$rse)),
  r_squared = sapply(scan_out, function(s) median(s$r2)),
  amp_mae = sapply(scan_out, function(s) median(abs(s$amp - amp_true))),
  amp_bias = sapply(scan_out, function(s) mean(s$amp - amp_true)),
  peak_mae = sapply(scan_out, function(s) mean(abs(circ_gap(s$peak, peak_true)))),
  peak_bias = sapply(scan_out, function(s) mean(circ_gap(s$peak, peak_true))))
print(round(k_tab, 5))

  harmonics parameters resid_sd r_squared amp_mae amp_bias peak_mae peak_bias
1         1          3  0.04305   0.95088 0.00555 -0.00031  9.54750  -9.54500
2         2          5  0.03963   0.96196 0.00555 -0.00033  3.31063  -0.07187
3         3          7  0.03959   0.96401 0.00556 -0.00034  5.60500  -0.14875

bic_mat <- sapply(scan_out, function(s)
  n_lyr * log(s$rss / n_lyr) + (s$npar + 1) * log(n_lyr))
aic_mat <- sapply(scan_out, function(s)
  n_lyr * log(s$rss / n_lyr) + 2 * (s$npar + 1))
bic_pick <- tabulate(max.col(-bic_mat, ties.method = "first"), 3)
aic_pick <- tabulate(max.col(-aic_mat, ties.method = "first"), 3)
print(rbind(harmonics = k_grid, bic_pixels = bic_pick, aic_pixels = aic_pick))

           [,1] [,2] [,3]
harmonics     1    2    3
bic_pixels  670  886   44
aic_pixels  320  988  292

The fit statistics behave the way fit statistics always do. Median residual standard deviation falls from 0.0430 at one harmonic to 0.0396 at two and 0.0396 at three, and the median R squared climbs from 0.9509 through 0.9620 to 0.9640: read that column alone and three harmonics looks harmless. The columns that matter say something else. Mean absolute error in the estimated peak day is 9.55 days at one harmonic, drops to 3.31 days at two, and rises again to 5.6 days at three: the third harmonic makes the estimate 69.3 per cent worse while making the residuals look better. What the third harmonic is fitting is the noise, and the argmax of a curve is a sharp functional of that curve, so a wiggle worth 0.00003 NDVI units of residual moves the crest by days.

The amplitude, by contrast, does not care. Median absolute error is 0.00555, 0.00555 and 0.00556 at one, two and three harmonics, and the bias is under 0.00034 throughout. That is not luck. On a complete year of evenly spaced dates the harmonic columns are mutually orthogonal, so leaving out the second harmonic cannot contaminate the first, and adding a third costs almost nothing in variance. Hold on to that sentence: the next section but one takes the even spacing away and the property goes with it.

Per-pixel information criteria do not settle it either. BIC picks one harmonic at 670 pixels, two at 886 and three at 44; AIC, being more permissive, picks three at 292 pixels, or 18.2 per cent of the window. Selecting per pixel would hand a third harmonic to hundreds of pixels that measurably do not want one, and produce a peak-day map whose estimator changes from place to place.

Two harmonics, fitted everywhere, is the recommendation these numbers support, and it is conditional on the shape of this signal: one growing season a year with a moderately asymmetric crest. A genuine double season, a spring flush and an autumn flush, needs at least the second harmonic to represent the two peaks and probably a third to place them, and a canopy that sits on a plateau for four months needs more terms than either. The way to find out is the comparison run above, against a truth or against withheld dates, not the R squared column.

A line chart on warm off-white paper with one, two and three harmonics on the horizontal axis and a vertical axis running from about thirty-five to one hundred per cent of the one-harmonic value. All three lines start together at one hundred on the left. A pale sage line labelled amplitude error stays flat along the top of the panel at one hundred across the whole width. A gold line labelled residual standard deviation drops to about ninety-two at two harmonics and stays there. A dark green line labelled peak day error plunges to about thirty-five at two harmonics and rebounds to about fifty-nine at three, making a deep V well below the other two. — Figure 3: Three summaries of the same three fits, each expressed as a percentage of its value at one harmonic so that they can be read on one scale. Residual standard deviation, the statistic that is easy to look at, improves by a few per cent and then stops. Amplitude error does not move at all, because on a complete regular year the harmonic terms do not interfere. Peak day error, the quantity the map is for, falls to a third at two harmonics and then gives most of that back at three.

Three ways to read the peak day, and the gap between them

Here is the comparison the post exists for. Three estimators of the same quantity, the day of year on which a pixel’s greenness crests:

the argmax of the raw series, that is the day of the composite with the largest NDVI;
the argmax of the fitted two-harmonic curve, evaluated on a daily grid;
the analytic phase of the first harmonic, $\phi/\omega$.

The first is the line most people write. The third is what you get if you read the phase coefficient off the model and call it the peak, which is common and is not the same as the second whenever the curve is asymmetric.

raw_peak <- doy[max.col(obs_mat, ties.method = "first")]
clean_peak <- doy[max.col(sig_obs, ties.method = "first")]

est_days <- list("raw seasonal maximum" = raw_peak,
                 "argmax of the fitted curve" = peak_hat,
                 "phase of the first harmonic" = phase_hat)
peak_tab <- do.call(rbind, lapply(names(est_days), function(nm) {
  e <- circ_gap(est_days[[nm]], peak_true)
  data.frame(estimator = nm, mean_abs_days = mean(abs(e)),
             median_abs_days = median(abs(e)), bias_days = mean(e),
             sd_days = sd(e), p90_abs_days = unname(quantile(abs(e), 0.9)),
             worse_8d = mean(abs(e) > 8))
}))
print(round(peak_tab[, -1], 4))

  mean_abs_days median_abs_days bias_days sd_days p90_abs_days worse_8d
1       13.6763         12.0000   -0.1188 17.0193      28.0000   0.6250
2        3.3106          3.0000   -0.0719  4.3187       7.0000   0.0562
3        9.5459          9.7902   -9.5435  2.3897      12.3119   0.7825

print(peak_tab$estimator)

[1] "raw seasonal maximum"        "argmax of the fitted curve" 
[3] "phase of the first harmonic"

raw_floor <- mean(abs(circ_gap(clean_peak, peak_true)))
gain_ratio <- peak_tab$mean_abs_days[1] / peak_tab$mean_abs_days[2]
print(round(c(noiseless_raw_argmax_error = raw_floor,
              quantisation_floor = (doy[2] - doy[1]) / 4,
              raw_over_fitted = gain_ratio, raw_sd_days = peak_tab$sd_days[1],
              fitted_sd_days = peak_tab$sd_days[2]), 4))

noiseless_raw_argmax_error         quantisation_floor 
                    2.1437                     2.0000 
           raw_over_fitted                raw_sd_days 
                    4.1310                    17.0193 
            fitted_sd_days 
                    4.3187

The raw seasonal maximum is off by 13.68 days on average, with a ninetieth percentile of 28 days and 62.5 per cent of pixels more than a composite interval from the truth. The fitted curve’s argmax averages 3.31 days, a factor of 4.13 better, from exactly the same 46 numbers.

Two separate things hurt the raw estimator, and the numbers separate them. Run the same argmax on the noiseless signal and the error is still 2.14 days, because an 8-day product can only return one of 46 dates and the average distance from a uniformly placed true peak to the nearest of them is 2 days. That is the compositing floor, and no care with the raw series gets under it. Noise then takes the raw estimator from 2.14 days to 13.68, because the top of a seasonal curve is nearly flat: over the three composites either side of the crest the signal changes by less than the observation noise, so which is largest is close to a coin toss. The fitted curve has neither problem. It is evaluated on a daily grid, so no compositing floor applies, and all 46 observations contribute to where the crest is put rather than just the winner.

The third estimator is the interesting failure, and it is the most precise of the three. Its standard deviation is 2.39 days against 4.32 for the fitted argmax, because it uses two coefficients rather than four and is a smooth function of them instead of a maximum. It is also biased by -9.54 days and puts 78.2 per cent of pixels more than a composite interval from the truth. The bias is less a defect than a statement that it estimates a different quantity: the phase of the first harmonic is the crest of the symmetric part of the curve, and when the second harmonic leans the curve towards a slow senescence the crest moves later while the phase stays put. Here the lean is worth 9.5 days, a fifth of the whole grassland-to-woodland difference the map is meant to show.

The rule out of this section: fit the model, then find the maximum of the fitted curve numerically. Do not read the peak off the phase coefficient unless the curve is known to be symmetric, and do not take the argmax of the raw composites at all. Harmonic fitting is in any case not the usual route to a phenological date: Zhang and colleagues (2003) fit piecewise logistic functions and take the extremes of the rate of change of curvature as the transitions, which suits a fast green-up better than a low order Fourier series does.

Three overlapping density curves on warm off-white paper, horizontal axis the error in the estimated peak day running from about minus fifty to plus fifty days, vertical axis density. A vertical dashed line stands at zero. A broad low pale sage curve labelled raw seasonal maximum spreads across the whole panel with a rounded top near zero. A tall narrow dark green curve labelled argmax of the fitted curve rises steeply at zero. An even taller and narrower gold curve labelled phase of the first harmonic peaks about ten days left of zero, with almost no overlap with the dashed line. — Figure 4: Error in the estimated peak day for all 1600 pixels, in days, for the three estimators. Zero is the truth. The raw seasonal maximum is unbiased and very wide, because the crest of the curve is flat compared with the observation noise. The argmax of the fitted curve is unbiased and narrow. The phase of the first harmonic is narrower still and sits about ten days early, because it measures the crest of the symmetric part of the curve rather than the crest of the curve.

When the dates go missing

Everything above assumed 46 evenly spaced composites. Real archives do not deliver that. Cloud, snow, sensor outages and quality flags remove observations, and they do not remove them at random across the year: in a temperate summer it is convective cloud, in a monsoon system it is a solid three months of nothing. Roerink, Menenti and Verhoef (2000) built HANTS around this problem, wrapping the same regression in an iterative rejection loop that discards observations lying below the fitted curve, on the grounds that cloud almost always depresses a vegetation index, and then predicting the missing dates from the fit. That reconstruction is the gap-filling post’s subject; the question here is narrower. Three scenarios, all dropping the same number of dates so that sample size is held fixed and only the pattern changes: nine at random, nine consecutive composites lost across the spring green-up, and nine consecutive composites lost across the senescence, each run over 50 noise replicates.

spring_out <- which(doy >= 97 & doy <= 161)
autumn_out <- which(doy >= 193 & doy <= 257)
keep_full <- seq_len(n_lyr)
keep_spring <- setdiff(keep_full, spring_out)
keep_autumn <- setdiff(keep_full, autumn_out)
n_drop <- length(spring_out)
print(c(dropped_per_scenario = n_drop, retained = n_lyr - n_drop,
        spring_gap_from = doy[min(spring_out)], spring_gap_to = doy[max(spring_out)],
        autumn_gap_from = doy[min(autumn_out)], autumn_gap_to = doy[max(autumn_out)]))

dropped_per_scenario             retained      spring_gap_from 
                   9                   37                   97 
       spring_gap_to      autumn_gap_from        autumn_gap_to 
                 161                  193                  257

gap_once <- function(y_mat, keep, n_harm) {
  f <- fit_harm(n_harm, keep, y_mat)
  c(amp_rel_bias = 100 * mean(f$amp / amp_true - 1),
    amp_mae = mean(abs(f$amp - amp_true)),
    peak_mae = mean(abs(circ_gap(f$peak, peak_true))), condition = kappa(f$des))
}

n_rep_gap <- 50
set.seed(4471)
gap_acc <- list()
for (b in seq_len(n_rep_gap)) {
  y_b <- sig_obs + matrix(rnorm(n_cell * n_lyr, 0, sd_noise), n_cell, n_lyr)
  keep_scatter <- sort(sample.int(n_lyr, n_lyr - n_drop))
  keeps <- list(complete = keep_full, scattered = keep_scatter,
                spring = keep_spring, autumn = keep_autumn)
  for (nm in names(keeps)) for (kk in 1:2) {
    gap_acc[[length(gap_acc) + 1]] <- data.frame(
      scenario = nm, harmonics = kk, t(gap_once(y_b, keeps[[nm]], kk)))
  }
}
gap_all <- do.call(rbind, gap_acc)
gap_tab <- aggregate(cbind(amp_rel_bias, amp_mae, peak_mae, condition) ~
                       scenario + harmonics, gap_all, mean)
gap_tab <- gap_tab[order(gap_tab$harmonics, match(gap_tab$scenario,
                    c("complete", "scattered", "spring", "autumn"))), ]
print(round(gap_tab[, -1], 4))

  harmonics amp_rel_bias amp_mae peak_mae condition
2         1       0.0470  0.0066   9.6285    1.3108
3         1      -0.0071  0.0076   9.6258    1.3792
4         1       3.2672  0.0105  11.6313    1.7843
1         1      -3.3844  0.0113  11.6809    1.7792
6         2       0.0463  0.0066   3.3146    1.3574
7         2       0.0580  0.0075   3.7491    1.5081
8         2       0.0861  0.0087   5.1927    3.5916
5         2       0.1151  0.0090   5.0670    3.4672

print(gap_tab$scenario)

[1] "complete"  "scattered" "spring"    "autumn"    "complete"  "scattered"
[7] "spring"    "autumn"

gval <- function(sc, kk, col) gap_tab[[col]][gap_tab$scenario == sc &
                                               gap_tab$harmonics == kk]

Dropping the dates at random costs precision and nothing else. With two harmonics the mean absolute amplitude error goes from 0.00663 on the complete year to 0.00750 with nine scattered dates missing, and peak day error from 3.31 to 3.75 days, while the relative bias in the amplitude is 0.058 per cent, which is nothing. A correctly specified least squares fit does not care which subset of dates it gets, only how many and how spread out.

A season-shaped gap is a different animal, and what it does depends on the number of harmonics. With two, the correct model here, the amplitude bias is still 0.086 per cent for the spring gap and 0.115 per cent for the autumn one. What is paid instead is conditioning and precision: the condition number of the design rises from 1.36 to 3.59, peak day error from 3.31 to 5.19 days, and amplitude error by a factor of 1.31.

With one harmonic the same gap produces a bias. On the complete year the omitted second harmonic was orthogonal to the retained columns and cost nothing: 0.047 per cent. Take out the spring and the amplitude comes back 3.27 per cent too high; take out the senescence instead and it comes back -3.38 per cent too low. Same number of observations, same noise, opposite sign of error, and the direction is set purely by which side of the crest is missing.

sweep_len <- 9
sweep_start <- seq(1, n_lyr, by = 2)
gap_sweep <- do.call(rbind, lapply(sweep_start, function(st) {
  drop_idx <- ((st - 1 + seq_len(sweep_len) - 1) %% n_lyr) + 1
  keep_i <- setdiff(seq_len(n_lyr), drop_idx)
  mid_day <- ((doy[st] + (sweep_len - 1) * (doy[2] - doy[1]) / 2 - 1) %% 365) + 1
  do.call(rbind, lapply(1:2, function(kk) {
    f <- fit_harm(kk, keep_i, sig_obs)
    data.frame(mid_day = mid_day, harmonics = kk,
               amp_rel_bias = 100 * mean(f$amp / amp_true - 1),
               frac_over = mean(f$amp > amp_true))
  }))
}))
gap_sweep <- gap_sweep[order(gap_sweep$harmonics, gap_sweep$mid_day), ]
s1 <- gap_sweep[gap_sweep$harmonics == 1, ]
worst_hi <- s1[which.max(s1$amp_rel_bias), ]
worst_lo <- s1[which.min(s1$amp_rel_bias), ]
s2 <- gap_sweep[gap_sweep$harmonics == 2, ]
print(round(s1[, c("mid_day", "amp_rel_bias", "frac_over")], 4))

   mid_day amp_rel_bias frac_over
43       4       0.5518    0.5406
45      20       2.5044    1.0000
1       33       3.4227    1.0000
3       49       3.2428    1.0000
5       65       2.2642    1.0000
7       81       1.3671    1.0000
9       97       1.3013    1.0000
11     113       2.1440    1.0000
13     129       3.2139    1.0000
15     145       3.5604    1.0000
17     161       2.6277    0.8219
19     177       0.6098    0.5969
21     193      -1.6647    0.2087
23     209      -3.2104    0.0000
25     225      -3.4726    0.0000
27     241      -2.6704    0.0000
29     257      -1.6208    0.0000
31     273      -1.1601    0.0000
33     289      -1.6059    0.0000
35     305      -2.5881    0.0000
37     321      -3.2906    0.0000
39     337      -3.0361    0.0175
41     353      -1.6142    0.2894

print(round(c(worst_over_pct = worst_hi$amp_rel_bias,
              worst_over_at_day = worst_hi$mid_day,
              worst_under_pct = worst_lo$amp_rel_bias,
              worst_under_at_day = worst_lo$mid_day,
              largest_two_harmonic_bias_pct = max(abs(s2$amp_rel_bias))), 4))

               worst_over_pct             worst_over_at_day 
                       3.5604                      145.0000 
              worst_under_pct            worst_under_at_day 
                      -3.4726                      225.0000 
largest_two_harmonic_bias_pct 
                       0.0000

Sweeping the block of 9 missing composites right around the year, on the noiseless signal so that only the geometry is left, gives the whole picture. With one harmonic the relative amplitude bias runs from +3.56 per cent when the block is centred on day 145 to -3.47 per cent when it is centred on day 225, and at both extremes 100 per cent of pixels are pushed the same way. With two harmonics the largest bias anywhere in the sweep is 4.1e-14 per cent, which is zero to machine precision.

The mechanism generalises past this example. Least squares does not bias the coefficients of a correctly specified model whatever the design is, so every bias here is the projection of the part of the signal the model does not carry onto the columns it does. On a complete evenly spaced year that projection is exactly zero by orthogonality of the harmonic basis, which is why a model one term short can look harmless. Remove a contiguous block of dates, the columns stop being orthogonal over the retained set, and the missing term leaks into the retained ones with a sign that follows its own sign over the retained window. Missing the green-up removes the part of the year where the asymmetry pulls the curve down, so the fitted swing is too large; missing the senescence removes where it pushes the curve up, so the swing is too small.

That is the finding of this section, and it is unpleasant. The number of harmonics stops being a variance question and becomes a bias question the moment the missing dates cluster in a season, which is the normal case. Choosing the model order on a complete year, or on the fit statistics, says nothing about how it will behave on the pixels where the cloud sat.

A line chart on warm off-white paper. The horizontal axis is the middle day of a missing block of nine composites, running from a few days after the new year to about day three hundred and fifty. The vertical axis is the relative amplitude error in per cent, from about minus four to plus four, with a horizontal reference line at zero. A gold line of round dots labelled one harmonic climbs from about plus half a per cent at the left edge to a hump of about plus three and a half near day thirty-five, dips to about plus one and a third near day ninety, rises to its maximum of about plus three and a half at day one hundred and forty-five, then falls steeply, crossing zero a little before day one hundred and ninety, reaching its minimum of about minus three and a half near day two hundred and twenty-five, easing back to about minus one near day two hundred and seventy, dipping again to about minus three and a third near day three hundred and twenty, and ending near minus one and a half. A dark green line labelled two harmonics is flat and indistinguishable from the zero line across the whole panel. — Figure 5: Relative error in the estimated seasonal amplitude when a block of nine consecutive composites is missing, against the middle of the missing block, computed on the noiseless signal so that only the geometry of the design is in it. With one harmonic the error swings from over three per cent too high when the green-up is missing to over three per cent too low when the senescence is missing. With two harmonics, the correct model here, the line sits on zero throughout.

The coefficients are estimates, and they have standard errors

The amplitude map is a covariate now. It will go into a species distribution model, a comparison of grazing compartments, a regression of bird abundance on seasonal productivity. What travels into that model is a single number per pixel, and that number came out of a regression with 46 observations and 5 parameters, so it has a standard error the map does not show. The amplitude is a nonlinear function of two coefficients, so the delta method supplies it. With $V$ the covariance matrix of the fitted coefficients,

\[\operatorname{Var}(A) \approx \frac{a^2 V_{aa} + 2ab\,V_{ab} + b^2 V_{bb}}{a^2 + b^2}\]

and the same construction with the derivatives of $\operatorname{atan2}$ gives the standard error of the phase.

resid_full <- obs_mat - cf2 %*% t(des2)
df_resid <- n_lyr - ncol(des2)
s2_pix <- rowSums(resid_full^2) / df_resid
xtx_inv <- solve(crossprod(des2))
v_aa <- xtx_inv["cos1", "cos1"]; v_bb <- xtx_inv["sin1", "sin1"]
v_ab <- xtx_inv["cos1", "sin1"]
a_hat <- cf2[, "cos1"]; b_hat <- cf2[, "sin1"]
var_amp <- s2_pix * (a_hat^2 * v_aa + 2 * a_hat * b_hat * v_ab +
                       b_hat^2 * v_bb) / amp_hat^2
se_amp <- sqrt(var_amp)
var_phase <- s2_pix * (b_hat^2 * v_aa - 2 * a_hat * b_hat * v_ab +
                         a_hat^2 * v_bb) / amp_hat^4
se_phase_days <- sqrt(var_phase) / omega

cover90 <- mean(abs(amp_err) < qt(0.95, df_resid) * se_amp)
print(round(c(median_se_amp = median(se_amp), min_se_amp = min(se_amp),
              max_se_amp = max(se_amp), actual_sd_of_amp_error = sd(amp_err),
              coverage_of_90pct_interval = cover90), 5))

             median_se_amp                 min_se_amp 
                   0.00823                    0.00501 
                max_se_amp     actual_sd_of_amp_error 
                   0.01200                    0.00829 
coverage_of_90pct_interval 
                   0.90625

print(round(c(median_rel_se_pct = 100 * median(se_amp / amp_hat),
              min_rel_se_pct = 100 * min(se_amp / amp_hat),
              max_rel_se_pct = 100 * max(se_amp / amp_hat),
              median_se_phase_days = median(se_phase_days)), 4))

   median_rel_se_pct       min_rel_se_pct       max_rel_se_pct 
              3.1205               1.5897               7.8691 
median_se_phase_days 
              1.8271

The median per-pixel standard error on the amplitude is 0.00823, and the actual standard deviation of the amplitude error across the window is 0.00829, so the reported uncertainty is honest: a nominal 90 per cent interval covers the truth at 90.6 per cent of pixels. As a fraction of the amplitude being estimated it runs from 1.59 per cent in the closed canopy to 7.87 per cent over the thinnest grass, because the standard error is nearly constant across the window while the amplitude is not. The companion standard error on the phase has a median of 1.83 days.

That heterogeneity is the reason to keep the standard error rather than admire it. A covariate measured with error attenuates the coefficient it carries, by the reliability ratio: the variance of the true covariate over the variance of the true covariate plus the measurement error variance.

grass_only <- wood < 0.15
rel_all <- var(amp_true) / (var(amp_true) + mean(var_amp))
rel_grass <- var(amp_true[grass_only]) /
  (var(amp_true[grass_only]) + mean(var_amp[grass_only]))

print(round(c(pixels_all = n_cell, pixels_grass = sum(grass_only),
              sd_true_amp_all = sd(amp_true),
              sd_true_amp_grass = sd(amp_true[grass_only])), 5))

       pixels_all      pixels_grass   sd_true_amp_all sd_true_amp_grass 
        1.600e+03         1.230e+02         5.324e-02         9.800e-03

print(round(c(reliability_all = rel_all,
              attenuation_pct_all = 100 * (1 - rel_all),
              reliability_grass = rel_grass,
              attenuation_pct_grass = 100 * (1 - rel_grass)), 4))

      reliability_all   attenuation_pct_all     reliability_grass 
               0.9763                2.3715                0.5763 
attenuation_pct_grass 
              42.3727

Over the whole window the reliability ratio is 0.9763, so a slope fitted on the amplitude map would be attenuated by about 2.37 per cent. That sounds like permission to ignore the problem, and it is not, because the ratio is a property of the sample and not of the method. Restrict the analysis to the 123 grassland pixels, where the true amplitude has a standard deviation of 0.00980 instead of 0.05324, and the reliability falls to 0.5763: an attenuation of 42.4 per cent, enough to turn a real effect into a marginal one. The measurement error did not change. The contrast asked of it did.

Carrying the point estimate forward on its own treats a noisy quantity as known. The cheap partial fix is to carry se_amp alongside amp_hat, so the downstream model can at least weight the pixels, or report the reliability ratio and the attenuation. What that will not fix is a bias like the one in the previous section, which no standard error describes.

A scatter plot on warm off-white paper with true amplitude from about 0.15 to 0.31 on the horizontal axis and estimated amplitude on the same range vertically. Eighty short vertical error bars with points at their centres lie along a dashed diagonal line of equality, scattered narrowly either side of it. The bars are all about the same length. Their colour runs from rust red at the left-hand, low amplitude end to dark green at the right-hand end, keyed to a legend giving the standard error as a percentage of the amplitude and running from about two at the dark end to about seven at the rust end. — Figure 6: Estimated seasonal amplitude against the truth for a systematic sample of 80 pixels, with plus and minus one standard error from the delta method, coloured by the standard error as a percentage of the amplitude. The dashed line is equality. The intervals are of nearly constant length across the range, so the same absolute uncertainty is a small fraction of a woodland amplitude and a large fraction of a grassland one.

The honest limit

The truth here was generated as a sum of two harmonics and then fitted with harmonics, which is the friendliest possible case and inflates every recovery number above. The amplitude came back to 2.21 per cent and the peak day to 3 days because the model could represent the signal exactly. On real reflectance it never can, so the errors reported here are lower bounds. The one place the simulation was not friendly is the section on missing dates, where a model one term short of the truth was fitted on purpose, and that is exactly where the only real bias in the post appeared.

The harmonic form itself is a strong assumption. A single cosine is symmetric about its crest, which no growing season is. Two harmonics can lean but they cannot produce a fast green-up followed by a long slow senescence and a flat winter, which is the actual shape of a deciduous canopy, and pushing the order up to get it runs into the result of the harmonics section: more terms fit the noise and move the peak. A genuine double season is the worse case, because a low order fit to two peaks does not fail loudly. It returns a single smooth hump somewhere between them, with a plausible amplitude and a peak day corresponding to nothing that happened on the ground, and nothing in the residuals of a 40 by 40 window at this noise level would make that obvious. The general warning has been measured on real data: Atkinson, Jeganathan, Dash and Atzberger (2012) compared four smoothing models, harmonic analysis among them, and found the estimated dates depend on which was used, while White and colleagues (2009) put a number on it for North America, where ten methods on the same record gave mean start-of-spring dates spanning around three weeks, with the disagreement varying by land cover. A peak day is a property of a fitted curve, and the curve is a choice.

The gap analysis used contiguous blocks of missing dates, a caricature of what cloud does. Real masks are patchy, differ from pixel to pixel within a scene, and lose observations at a rate correlated with the surface state, so the missingness is not ignorable the way the random-drop scenario assumed: the direction of the bias should carry over, the size should not be read across. And nothing here has atmospheric correction, sensor drift, changing view geometry or georegistration error in it, the cube having arrived calibrated and aligned. Each of those can put a seasonal or spatial pattern into the coefficients that has nothing to do with vegetation, and a sub-pixel misregistration in particular puts a false amplitude signal along every land cover boundary, which is where the interesting ecology usually is.

Where to go next

The three coefficient rasters produced here, mean, amplitude and peak day, are covariates, and the checking post is the next step before they go anywhere: whether they vary enough to be worth modelling, and whether their errors are spatially correlated, are questions nothing above answers. If the missing dates in your own cube are the problem rather than a sideshow, the gap-filling post treats reconstruction properly. And if what you want is a phenological date to set against ground records, the phenology posts are where that comparison happens, with the warning from here attached: the satellite date and the observer date are two different measurements of a partly shared thing.

References

Jakubauskas ME, Legates DR, Kastens JH 2002 Computers and Electronics in Agriculture 37(1-3):127-139 (10.1016/S0168-1699(02)00116-3)

Roerink GJ, Menenti M, Verhoef W 2000 International Journal of Remote Sensing 21(9):1911-1917 (10.1080/014311600209814)

Zhang X, Friedl MA, Schaaf CB, Strahler AH, Hodges JCF, Gao F, Reed BC, Huete A 2003 Remote Sensing of Environment 84(3):471-475 (10.1016/S0034-4257(02)00135-9)

White MA, de Beurs KM, Didan K, et al 2009 Global Change Biology 15(10):2335-2359 (10.1111/j.1365-2486.2009.01910.x)

Verbesselt J, Hyndman R, Newnham G, Culvenor D 2010 Remote Sensing of Environment 114(1):106-115 (10.1016/j.rse.2009.08.014)

Atkinson PM, Jeganathan C, Dash J, Atzberger C 2012 Remote Sensing of Environment 123:400-417 (10.1016/j.rse.2012.04.001)

Zhu Z, Woodcock CE 2014 Remote Sensing of Environment 144:152-171 (10.1016/j.rse.2014.01.011)

Fisher NI 1993 Statistical Analysis of Circular Data. Cambridge University Press. ISBN 978-0-521-56890-6.

--- title: "Harmonic regression on a seasonal raster" description: "Fit a sine and cosine pair per pixel to a year of NDVI composites in R, map amplitude and peak day, and measure what the raw seasonal maximum gets wrong." date: "2026-08-03 10:00" categories: [R, remote sensing, time series, terra, ecology tutorial] image: thumbnail.png image-alt: "Three overlapping density curves on warm off-white paper showing the error in the estimated peak day of the growing season, in days, for sixteen hundred pixels. A vertical dashed line marks zero error. A broad low pale sage curve labelled raw seasonal maximum spreads from about minus thirty to plus forty days with a flat top near zero. A tall dark green curve labelled argmax of the fitted curve rises to a sharp point at zero. The tallest and narrowest curve, in gold, is labelled phase of the first harmonic and peaks about ten days to the left of zero, barely touching the dashed line." --- A window 1200 metres on a side over a lowland mosaic: rough grassland on free-draining gravels in the south, a block of oak and hornbeam in the north, and a ragged band of encroaching scrub between them. The satellite archive holds one year of 8-day NDVI composites over that window, 46 layers, on a 30 metre grid, so the cube is 40 pixels by 40 pixels with 46 observations at each. It is a deliberately small window: everything below runs on a laptop in under a minute. The question is not what the NDVI was on any particular date. It is what the seasonal cycle looks like at each pixel: how far the vegetation swings between its winter floor and its summer ceiling, and when it peaks. Those two numbers, an amplitude and a phase, are what goes into the habitat model or the comparison between management compartments, and they have to be recovered from 46 noisy values per pixel, with clouds deleting some of them and the composites eight days apart to begin with. The obvious way to get the peak is to ask which layer is largest, which is one line of code and what most people write first. It is also the estimator this post is mainly about, because a fitted seasonal curve beats it by a factor of four in these data. [NDVI time series from a raster stack](../ndvi-time-series-from-a-raster-stack/) built a cube of this shape and took per-pixel summaries off it: means, maxima, an integral over the season. This post fits a model to the time axis instead, which is a different object, because a summary throws away the shape and a model parametrises it. [Raster data in R with terra](../terra-raster-basics/) covers single-layer work and is not repeated. [Filling cloud gaps in a satellite series](../filling-cloud-gaps-in-a-satellite-series/) is about reconstructing missing observations; the section on missing dates below asks only what an unfilled gap does to a fitted coefficient. [Checking a remote sensing covariate](../checking-a-remote-sensing-covariate/) takes the outputs of this post and asks whether they deserve to be in a model at all. Two older posts set the other boundaries. [Phenological trends and temperature](../phenological-trends-and-temperature/) estimates timing from ground observation: somebody walked the transect and wrote down a date, and here the same quantity is read off a reflectance curve by a machine that never visited the site. And [circular data and the von Mises distribution](../circular-data-and-von-mises/) owns the arithmetic of angles, which matters because a phase is an angle and forgetting that is one of the standard mistakes in this workflow. Six things get measured: how well amplitude, phase and mean are recovered when the truth is known; how the answer moves as harmonics are added; how far the raw seasonal maximum is from the truth compared with two fitted alternatives; what an unevenly missing set of dates does; which direction a season-shaped gap pushes the amplitude; and how large the per-pixel standard error on the amplitude is before anybody carries it into another model. ```{r setup} #| message: false #| warning: false library(ggplot2) te_pal <- list(forest = "#275139", green = "#2f8f63", sage = "#93a87f", clay = "#b5534e", gold = "#cda23f", line = "#dad9ca", ink = "#16241d", paper = "#f5f4ee") theme_te <- function() { theme_minimal(base_size = 12) + theme(panel.grid.minor = element_blank(), panel.grid.major = element_line(colour = "#e7e6dc"), plot.background = element_rect(fill = "#f5f4ee", colour = NA), panel.background = element_rect(fill = "#f5f4ee", colour = NA), plot.title = element_text(face = "bold", colour = te_pal$ink), axis.title = element_text(colour = "#2c3a31"), axis.text = element_text(colour = "#2c3a31")) } ``` ## A cube with a known answer inside it Simulated data, because the point is to measure error and that needs a truth. Each pixel gets a mean NDVI, an amplitude, a peak day and a small second harmonic, all varying smoothly across the window through a woodland index running from open grass in the south to closed canopy in the north. The second harmonic is what makes the curve asymmetric: green-up in a deciduous canopy is fast and senescence slow, and a single cosine cannot express that. ```{r cube} suppressPackageStartupMessages(library(terra)) n_side <- 40 n_cell <- n_side^2 n_lyr <- 46 doy <- seq(1, 361, by = 8) omega <- 2 * pi / 365 tmpl <- rast(nrows = n_side, ncols = n_side, xmin = 0, xmax = 1200, ymin = 0, ymax = 1200, crs = "EPSG:32634") xy_cell <- crds(tmpl) g_east <- xy_cell[, 1] / 1200; g_north <- xy_cell[, 2] / 1200 wood <- plogis(7 * (g_north - 0.5 + 0.28 * sin(3.1 * g_east))) mean_true <- 0.34 + 0.16 * wood amp_true <- 0.13 + 0.17 * wood + 0.02 * g_east peak_par <- 150 + 36 * wood + 10 * g_north amp2_true <- 0.008 + 0.022 * wood day_fine <- 1:365 ang_fine <- omega * (matrix(day_fine, n_cell, 365, byrow = TRUE) - peak_par) sig_fine <- mean_true + amp_true * cos(ang_fine) + amp2_true * sin(2 * ang_fine) peak_true <- day_fine[max.col(sig_fine, ties.method = "first")] ang_obs <- omega * (matrix(doy, n_cell, n_lyr, byrow = TRUE) - peak_par) sig_obs <- mean_true + amp_true * cos(ang_obs) + amp2_true * sin(2 * ang_obs) sd_noise <- 0.04 set.seed(20260803) obs_mat <- sig_obs + matrix(rnorm(n_cell * n_lyr, 0, sd_noise), n_cell, n_lyr) ndvi <- setValues(rast(tmpl, nlyr = n_lyr), obs_mat) names(ndvi) <- sprintf("d%03d", doy) print(c(pixels = ncell(ndvi), layers = nlyr(ndvi), composite_days = doy[2] - doy[1], first_doy = doy[1], last_doy = doy[n_lyr])) print(round(c(mean_lo = min(mean_true), mean_hi = max(mean_true), amp_lo = min(amp_true), amp_hi = max(amp_true), peak_lo = min(peak_true), peak_hi = max(peak_true), noise_sd = sd_noise), 4)) ``` The cube holds `r ncell(ndvi)` pixels and `r nlyr(ndvi)` layers, one composite every `r doy[2] - doy[1]` days from day `r doy[1]` to day `r doy[n_lyr]`. True amplitude runs from `r round(min(amp_true), 3)` in the driest grass to `r round(max(amp_true), 3)` in the closed canopy, and the true peak day, defined as the argmax of the noiseless annual curve, runs from day `r min(peak_true)` to day `r max(peak_true)`. Per-observation noise has a standard deviation of `r sd_noise`, which is an ordinary figure for a composited vegetation index once residual cloud, aerosol and view-angle effects are in it. One `terra` habit is worth stating first, because it fails silently. A SpatRaster plus a bare numeric vector recycles across layers, not across cells, so adding a 1600-element per-pixel vector to a 46-layer cube produces a cube with hundreds of layers and no error. The safe route is the one used above: build the values as a matrix of cells by layers, hand it to `setValues()`, check `nlyr()`. ## One pixel, one linear model Write $y_t$ for the NDVI of a pixel at day of year $t$ and $\omega = 2\pi/365$. The one-harmonic model is $$y_t = \mu + a\cos(\omega t) + b\sin(\omega t) + \varepsilon_t$$ which looks like a curve-fitting problem and is not one. The pair $\cos(\omega t)$ and $\sin(\omega t)$ are known numbers once the dates are known, so $\mu$, $a$ and $b$ enter linearly and the whole thing is `lm()`. No optimiser, no starting values, no convergence to worry about. That is why harmonic regression is the workhorse in this corner of remote sensing. Jakubauskas, Legates and Kastens (2002) used the amplitudes and phases of the first few harmonics of an AVHRR NDVI series as classification features, the direct ancestor of using them as ecological covariates; Zhu and Woodcock (2014) put the same form at the centre of continuous change detection, fitting it to a stable period and flagging a pixel when new observations fall outside the prediction interval, which only works because refitting everywhere is cheap; and Verbesselt, Hyndman, Newnham and Culvenor (2010) take the complementary route in BFAST, separating a seasonal component from a piecewise trend and testing for breaks in each. The two coefficients are then rewritten as an amplitude and a phase, because that is what anybody wants to look at: $$a\cos(\omega t) + b\sin(\omega t) = A\cos(\omega t - \phi), \qquad A = \sqrt{a^2 + b^2}, \qquad \phi = \operatorname{atan2}(b, a)$$ and the phase converts to a day of year as $\phi/\omega$, taken modulo 365. Adding more harmonics means adding $\cos(2\omega t)$, $\sin(2\omega t)$ and so on: still linear, still `lm()`. ```{r one-pixel} harm_design <- function(days, n_harm) { des <- matrix(1, length(days), 1 + 2 * n_harm) colnames(des) <- c("mean", as.vector(rbind(paste0("cos", seq_len(n_harm)), paste0("sin", seq_len(n_harm))))) for (j in seq_len(n_harm)) { des[, 2 * j] <- cos(2 * pi * j * days / 365) des[, 2 * j + 1] <- sin(2 * pi * j * days / 365) } des } pix_id <- cellFromXY(tmpl, cbind(315, 915)) pix_dat <- data.frame(ndvi = obs_mat[pix_id, ], cos1 = cos(omega * doy), sin1 = sin(omega * doy)) fit_pix <- lm(ndvi ~ cos1 + sin1, data = pix_dat) print(round(summary(fit_pix)$coefficients, 5)) cf_pix <- coef(fit_pix) amp_pix <- sqrt(cf_pix[2]^2 + cf_pix[3]^2) phase_pix <- (atan2(cf_pix[3], cf_pix[2]) / omega) %% 365 lmfit_pix <- lm.fit(harm_design(doy, 1), obs_mat[pix_id, ])$coefficients print(round(c(intercept = unname(cf_pix[1]), true_mean = mean_true[pix_id], amplitude = unname(amp_pix), true_amplitude = amp_true[pix_id], phase_day = unname(phase_pix), true_peak_day = peak_true[pix_id], woodland_index = wood[pix_id]), 4)) print(round(c(max_gap_lm_vs_lmfit = max(abs(cf_pix - lmfit_pix))), 12)) ``` At this pixel the intercept comes out at `r round(coef(fit_pix)[1], 4)` against a true mean of `r round(mean_true[pix_id], 4)`, and the amplitude at `r round(amp_pix, 4)` against a true `r round(amp_true[pix_id], 4)`. The phase converts to day `r round(phase_pix, 1)` while the true peak of the pixel's curve is day `r peak_true[pix_id]`, a gap of `r round(peak_true[pix_id] - phase_pix, 1)` days that is not noise and will be taken apart later. `lm()` and `lm.fit()` on the same design return identical coefficients, which matters only because the per-pixel loop below uses the second one for speed. ### The phase is an angle, and angles do not average Once a phase has been converted to a day of year it looks like an ordinary number, and the moment somebody averages peak days over a management unit or differences them between years, the branch cut bites. Day 362 and day 3 are six days apart, not 359. ```{r branch-cut} wrap_days <- c(356, 362, 3, 9, 14) naive_mean <- mean(wrap_days) wrap_ang <- omega * wrap_days vec_mean <- (atan2(mean(sin(wrap_ang)), mean(cos(wrap_ang))) / omega) %% 365 resultant <- sqrt(mean(cos(wrap_ang))^2 + mean(sin(wrap_ang))^2) print(round(c(arithmetic_mean_day = naive_mean, vector_mean_day = vec_mean, resultant_length = resultant, error_days = naive_mean - vec_mean), 3)) ``` Five pixels whose peaks straddle the new year have an arithmetic mean peak day of `r round(naive_mean, 1)`, which is in late May. The vector mean, taken by averaging the unit vectors and converting back with `atan2()`, is day `r round(vec_mean, 1)`, and it is the right answer: the arithmetic version is off by `r round(naive_mean - vec_mean, 1)` days. The resultant length `r round(resultant, 3)` measures the agreement, near one when the phases coincide and near zero when they cancel. Fisher (1993) is the standard treatment and the base R mechanics are in the circular post. Nothing in this window peaks near the new year, so the problem never surfaces here, but it is exactly what a semi-arid system with a winter growing season does, and it fails silently when it does. ```{r fig-pixels} #| echo: false #| fig-cap: "Three pixels from the window: open grassland, the scrub edge and closed canopy. Dots are the 46 composite observations, the pale line is the noiseless annual curve the observations were drawn from, and the dark line is the two-harmonic fit. The grassland pixel swings less and peaks earlier; the canopy pixel swings more, peaks later and is visibly asymmetric, rising faster than it falls." #| fig-alt: "Three panels side by side, each with day of year from zero to three hundred and sixty-five on the horizontal axis and NDVI from about two tenths to eight tenths on the vertical axis. In every panel forty-six round dots scatter closely around two nearly coincident curves, one pale sage and one dark green, that rise from a flat winter floor near 0.2 to a single summer hump and fall again. The left panel, labelled open grassland, has the shallowest hump, cresting near day one hundred and seventy at about 0.55. The middle panel, the scrub edge, crests near day one hundred and ninety-five at about 0.76. The right panel, closed canopy, has the tallest hump, cresting just past day two hundred at about 0.8, and its rise is steeper than its fall." #| fig-width: 7.8 #| fig-height: 3.6 show_ids <- c(cellFromXY(tmpl, cbind(315, 135)), cellFromXY(tmpl, cbind(615, 555)), cellFromXY(tmpl, cbind(315, 915))) show_lab <- c("open grassland", "scrub edge", "closed canopy") des2_show <- harm_design(doy, 2) fine2_show <- harm_design(day_fine, 2) pt_df <- do.call(rbind, Map(function(id, lb) data.frame(day = doy, ndvi = obs_mat[id, ], panel = lb), show_ids, show_lab)) ln_df <- do.call(rbind, Map(function(id, lb) { cf <- lm.fit(des2_show, obs_mat[id, ])$coefficients rbind(data.frame(day = day_fine, ndvi = sig_fine[id, ], kind = "the truth", panel = lb), data.frame(day = day_fine, ndvi = as.vector(fine2_show %*% cf), kind = "two-harmonic fit", panel = lb)) }, show_ids, show_lab)) pt_df$panel <- factor(pt_df$panel, levels = show_lab) ln_df$panel <- factor(ln_df$panel, levels = show_lab) ggplot(pt_df, aes(day, ndvi)) + geom_point(size = 1.2, colour = te_pal$ink, alpha = 0.75) + geom_line(data = ln_df, aes(colour = kind, linewidth = kind)) + scale_colour_manual(values = c("the truth" = te_pal$sage, "two-harmonic fit" = te_pal$forest), name = NULL) + scale_linewidth_manual(values = c("the truth" = 1.6, "two-harmonic fit" = 0.7), name = NULL) + scale_x_continuous(breaks = c(1, 120, 240, 360)) + facet_wrap(~panel) + labs(x = "day of year", y = "NDVI", title = "One seasonal curve per pixel") + theme_te() + theme(legend.position = "bottom", plot.margin = margin(8, 14, 4, 8)) ``` ## Fitting every pixel with app() `terra::app()` applies a function to the vector of layer values at each cell, which is exactly the shape a per-pixel time series has. The function returns a vector of coefficients, so the result is a multi-layer raster with one layer per coefficient. Two harmonics gives five: a mean and two cosine and sine pairs. ```{r app-fit} des2 <- harm_design(doy, 2) harm_r <- app(ndvi, function(v) lm.fit(des2, v)$coefficients) names(harm_r) <- colnames(des2) print(harm_r) cf2 <- values(harm_r) cf_direct <- t(qr.coef(qr(des2), t(obs_mat))) print(c(coef_layers = nlyr(harm_r), agrees_with_matrix_solve = max(abs(cf2 - cf_direct)) < 1e-12)) ``` The output carries `r nlyr(harm_r)` layers and reproduces a single matrix solve of the same design exactly, which is worth checking once and then relying on: the sweeps later in the post use the matrix form, because running `app()` inside a replicate loop would spend its time on raster overheads for no gain. ```{r recovery} fine2 <- harm_design(day_fine, 2) circ_gap <- function(a, b) ((a - b + 182.5) %% 365) - 182.5 amp_hat <- sqrt(cf2[, "cos1"]^2 + cf2[, "sin1"]^2) phase_hat <- (atan2(cf2[, "sin1"], cf2[, "cos1"]) / omega) %% 365 peak_hat <- day_fine[max.col(cf2 %*% t(fine2), ties.method = "first")] amp_r <- setValues(rast(tmpl), amp_hat) peak_r <- setValues(rast(tmpl), peak_hat) print(c(amp_layers = nlyr(amp_r), peak_layers = nlyr(peak_r))) amp_err <- amp_hat - amp_true peak_err <- circ_gap(peak_hat, peak_true) mean_err <- cf2[, "mean"] - mean_true amp_mae <- median(abs(amp_err)) amp_p90 <- unname(quantile(abs(amp_err), 0.9)) peak_mae <- median(abs(peak_err)) peak_p90 <- unname(quantile(abs(peak_err), 0.9)) print(round(c(mean_median_abs_err = median(abs(mean_err)), amp_median_abs_err = amp_mae, amp_p90_abs_err = amp_p90, amp_median_pct_err = 100 * median(abs(amp_err / amp_true)), amp_bias = mean(amp_err)), 5)) print(round(c(peak_median_abs_err_days = peak_mae, peak_p90_days = peak_p90, peak_bias_days = mean(peak_err), peak_worst_days = max(abs(peak_err))), 3)) ``` Against a known truth the two-harmonic fit recovers the annual mean to a median absolute error of `r sprintf("%.4f", median(abs(mean_err)))` NDVI units and the amplitude to `r sprintf("%.4f", amp_mae)`, or `r round(100 * median(abs(amp_err / amp_true)), 2)` per cent of the amplitude being estimated, with a ninetieth percentile of `r sprintf("%.4f", amp_p90)`: a tail about two and a half times the median rather than an order of magnitude away. Peak day comes back with a median absolute error of `r peak_mae` days, a ninetieth percentile of `r peak_p90` days and a mean signed error of `r round(mean(peak_err), 2)` days, the worst pixel being `r round(max(abs(peak_err)), 0)` days out. ```{r fig-maps} #| echo: false #| fig-cap: "The two coefficient maps that a downstream model would actually use, both estimated pixel by pixel from 46 noisy composites. Left, the amplitude of the first harmonic: the seasonal swing, largest in the closed canopy in the north of the window and smallest over the grassland in the south, with a broad transition where the scrub is. Right, the peak day taken as the argmax of the fitted two-harmonic curve, which follows the same pattern and spans about seven weeks across the window. The speckle in both is estimation noise, not vegetation." #| fig-alt: "Two square maps side by side on warm off-white paper, each twelve hundred metres across with easting and northing axes in metres. The left map shows the seasonal amplitude and grades smoothly from pale sage along the bottom edge, where the legend reads about 0.15, to dark forest green across the top at about 0.30. The pale zone reaches highest up the left and right edges and is squeezed lowest across the middle of the panel, so the transition bows downwards. The right map shows the peak day of the growing season on the same shape, gold at about one hundred and sixty days along the bottom and dark green at about two hundred days along the top. Both maps carry a fine grainy speckle everywhere, which is the pixel-to-pixel estimation noise." #| fig-width: 8.0 #| fig-height: 3.9 to_df <- function(r) { d <- as.data.frame(r, xy = TRUE) names(d)[3] <- "z" d } map_panel <- function(dat, lo, hi, ttl, lname) { ggplot(dat, aes(x, y, fill = z)) + geom_raster() + coord_equal(expand = FALSE) + scale_fill_gradient(low = lo, high = hi, name = lname) + labs(x = "easting (m)", y = "northing (m)", title = ttl) + theme_te() + theme(plot.title = element_text(size = 10), axis.text = element_text(size = 7), axis.title = element_text(size = 8), legend.key.width = grid::unit(0.34, "cm"), legend.text = element_text(size = 7.5), legend.title = element_text(size = 8), plot.margin = margin(4, 6, 2, 4)) } grid::grid.newpage() grid::grid.rect(gp = grid::gpar(fill = te_pal$paper, col = NA)) grid::pushViewport(grid::viewport(layout = grid::grid.layout(1, 2))) print(map_panel(to_df(amp_r), te_pal$sage, te_pal$forest, "seasonal amplitude", "NDVI"), vp = grid::viewport(layout.pos.row = 1, layout.pos.col = 1)) print(map_panel(to_df(peak_r), te_pal$gold, te_pal$forest, "peak day of the fitted curve", "day"), vp = grid::viewport(layout.pos.row = 1, layout.pos.col = 2)) ``` ## How many harmonics One harmonic is a symmetric hump. Two can lean. Three can carry a shoulder. The usual way to choose is a fit statistic, so start there and then ask whether the fit statistic is answering the question. ```{r harmonic-scan} fit_harm <- function(n_harm, keep = seq_len(n_lyr), y_mat = obs_mat) { des <- harm_design(doy[keep], n_harm) yk <- y_mat[, keep, drop = FALSE] cf <- t(qr.coef(qr(des), t(yk))) rss <- rowSums((yk - cf %*% t(des))^2) list(cf = cf, des = des, rss = rss, npar = ncol(des), amp = sqrt(cf[, 2]^2 + cf[, 3]^2), peak = day_fine[max.col(cf %*% t(harm_design(day_fine, n_harm)), ties.method = "first")], rse = sqrt(rss / (length(keep) - ncol(des))), r2 = 1 - rss / rowSums((yk - rowMeans(yk))^2)) } k_grid <- 1:3 scan_out <- lapply(k_grid, fit_harm) k_tab <- data.frame( harmonics = k_grid, parameters = sapply(scan_out, function(s) s$npar), resid_sd = sapply(scan_out, function(s) median(s$rse)), r_squared = sapply(scan_out, function(s) median(s$r2)), amp_mae = sapply(scan_out, function(s) median(abs(s$amp - amp_true))), amp_bias = sapply(scan_out, function(s) mean(s$amp - amp_true)), peak_mae = sapply(scan_out, function(s) mean(abs(circ_gap(s$peak, peak_true)))), peak_bias = sapply(scan_out, function(s) mean(circ_gap(s$peak, peak_true)))) print(round(k_tab, 5)) bic_mat <- sapply(scan_out, function(s) n_lyr * log(s$rss / n_lyr) + (s$npar + 1) * log(n_lyr)) aic_mat <- sapply(scan_out, function(s) n_lyr * log(s$rss / n_lyr) + 2 * (s$npar + 1)) bic_pick <- tabulate(max.col(-bic_mat, ties.method = "first"), 3) aic_pick <- tabulate(max.col(-aic_mat, ties.method = "first"), 3) print(rbind(harmonics = k_grid, bic_pixels = bic_pick, aic_pixels = aic_pick)) ``` The fit statistics behave the way fit statistics always do. Median residual standard deviation falls from `r sprintf("%.4f", k_tab$resid_sd[1])` at one harmonic to `r sprintf("%.4f", k_tab$resid_sd[2])` at two and `r sprintf("%.4f", k_tab$resid_sd[3])` at three, and the median R squared climbs from `r sprintf("%.4f", k_tab$r_squared[1])` through `r sprintf("%.4f", k_tab$r_squared[2])` to `r sprintf("%.4f", k_tab$r_squared[3])`: read that column alone and three harmonics looks harmless. The columns that matter say something else. Mean absolute error in the estimated peak day is `r round(k_tab$peak_mae[1], 2)` days at one harmonic, drops to `r round(k_tab$peak_mae[2], 2)` days at two, and rises again to `r round(k_tab$peak_mae[3], 2)` days at three: the third harmonic makes the estimate `r round(100 * (k_tab$peak_mae[3] / k_tab$peak_mae[2] - 1), 1)` per cent worse while making the residuals look better. What the third harmonic is fitting is the noise, and the argmax of a curve is a sharp functional of that curve, so a wiggle worth `r sprintf("%.5f", k_tab$resid_sd[2] - k_tab$resid_sd[3])` NDVI units of residual moves the crest by days. The amplitude, by contrast, does not care. Median absolute error is `r sprintf("%.5f", k_tab$amp_mae[1])`, `r sprintf("%.5f", k_tab$amp_mae[2])` and `r sprintf("%.5f", k_tab$amp_mae[3])` at one, two and three harmonics, and the bias is under `r sprintf("%.5f", max(abs(k_tab$amp_bias)))` throughout. That is not luck. On a complete year of evenly spaced dates the harmonic columns are mutually orthogonal, so leaving out the second harmonic cannot contaminate the first, and adding a third costs almost nothing in variance. Hold on to that sentence: the next section but one takes the even spacing away and the property goes with it. Per-pixel information criteria do not settle it either. BIC picks one harmonic at `r bic_pick[1]` pixels, two at `r bic_pick[2]` and three at `r bic_pick[3]`; AIC, being more permissive, picks three at `r aic_pick[3]` pixels, or `r round(100 * aic_pick[3] / n_cell, 1)` per cent of the window. Selecting per pixel would hand a third harmonic to hundreds of pixels that measurably do not want one, and produce a peak-day map whose estimator changes from place to place. Two harmonics, fitted everywhere, is the recommendation these numbers support, and it is conditional on the shape of this signal: one growing season a year with a moderately asymmetric crest. A genuine double season, a spring flush and an autumn flush, needs at least the second harmonic to represent the two peaks and probably a third to place them, and a canopy that sits on a plateau for four months needs more terms than either. The way to find out is the comparison run above, against a truth or against withheld dates, not the R squared column. ```{r fig-harmonics} #| echo: false #| fig-cap: "Three summaries of the same three fits, each expressed as a percentage of its value at one harmonic so that they can be read on one scale. Residual standard deviation, the statistic that is easy to look at, improves by a few per cent and then stops. Amplitude error does not move at all, because on a complete regular year the harmonic terms do not interfere. Peak day error, the quantity the map is for, falls to a third at two harmonics and then gives most of that back at three." #| fig-alt: "A line chart on warm off-white paper with one, two and three harmonics on the horizontal axis and a vertical axis running from about thirty-five to one hundred per cent of the one-harmonic value. All three lines start together at one hundred on the left. A pale sage line labelled amplitude error stays flat along the top of the panel at one hundred across the whole width. A gold line labelled residual standard deviation drops to about ninety-two at two harmonics and stays there. A dark green line labelled peak day error plunges to about thirty-five at two harmonics and rebounds to about fifty-nine at three, making a deep V well below the other two." #| fig-width: 7.4 #| fig-height: 4.2 metric_lab <- c("residual standard deviation", "amplitude error", "peak day error") harm_long <- do.call(rbind, Map(function(v, m) data.frame(harmonics = k_grid, rel = 100 * v / v[1], metric = m), list(k_tab$resid_sd, k_tab$amp_mae, k_tab$peak_mae), metric_lab)) harm_long$metric <- factor(harm_long$metric, levels = metric_lab) ggplot(harm_long, aes(harmonics, rel, colour = metric)) + geom_line(linewidth = 0.9) + geom_point(size = 2.6) + scale_colour_manual(values = c(te_pal$gold, te_pal$sage, te_pal$forest), name = NULL) + scale_x_continuous(breaks = k_grid) + labs(x = "harmonics in the model", y = "per cent of the one-harmonic value", title = "The statistic that improves is not the one you need") + theme_te() + theme(legend.position = "bottom", plot.margin = margin(8, 14, 4, 8)) ``` ## Three ways to read the peak day, and the gap between them Here is the comparison the post exists for. Three estimators of the same quantity, the day of year on which a pixel's greenness crests: 1. the argmax of the raw series, that is the day of the composite with the largest NDVI; 2. the argmax of the fitted two-harmonic curve, evaluated on a daily grid; 3. the analytic phase of the first harmonic, $\phi/\omega$. The first is the line most people write. The third is what you get if you read the phase coefficient off the model and call it the peak, which is common and is not the same as the second whenever the curve is asymmetric. ```{r peak-estimators} raw_peak <- doy[max.col(obs_mat, ties.method = "first")] clean_peak <- doy[max.col(sig_obs, ties.method = "first")] est_days <- list("raw seasonal maximum" = raw_peak, "argmax of the fitted curve" = peak_hat, "phase of the first harmonic" = phase_hat) peak_tab <- do.call(rbind, lapply(names(est_days), function(nm) { e <- circ_gap(est_days[[nm]], peak_true) data.frame(estimator = nm, mean_abs_days = mean(abs(e)), median_abs_days = median(abs(e)), bias_days = mean(e), sd_days = sd(e), p90_abs_days = unname(quantile(abs(e), 0.9)), worse_8d = mean(abs(e) > 8)) })) print(round(peak_tab[, -1], 4)) print(peak_tab$estimator) raw_floor <- mean(abs(circ_gap(clean_peak, peak_true))) gain_ratio <- peak_tab$mean_abs_days[1] / peak_tab$mean_abs_days[2] print(round(c(noiseless_raw_argmax_error = raw_floor, quantisation_floor = (doy[2] - doy[1]) / 4, raw_over_fitted = gain_ratio, raw_sd_days = peak_tab$sd_days[1], fitted_sd_days = peak_tab$sd_days[2]), 4)) ``` The raw seasonal maximum is off by `r round(peak_tab$mean_abs_days[1], 2)` days on average, with a ninetieth percentile of `r round(peak_tab$p90_abs_days[1], 0)` days and `r round(100 * peak_tab$worse_8d[1], 1)` per cent of pixels more than a composite interval from the truth. The fitted curve's argmax averages `r round(peak_tab$mean_abs_days[2], 2)` days, a factor of `r round(gain_ratio, 2)` better, from exactly the same 46 numbers. Two separate things hurt the raw estimator, and the numbers separate them. Run the same argmax on the noiseless signal and the error is still `r round(raw_floor, 2)` days, because an 8-day product can only return one of 46 dates and the average distance from a uniformly placed true peak to the nearest of them is `r (doy[2] - doy[1]) / 4` days. That is the compositing floor, and no care with the raw series gets under it. Noise then takes the raw estimator from `r round(raw_floor, 2)` days to `r round(peak_tab$mean_abs_days[1], 2)`, because the top of a seasonal curve is nearly flat: over the three composites either side of the crest the signal changes by less than the observation noise, so which is largest is close to a coin toss. The fitted curve has neither problem. It is evaluated on a daily grid, so no compositing floor applies, and all 46 observations contribute to where the crest is put rather than just the winner. The third estimator is the interesting failure, and it is the most precise of the three. Its standard deviation is `r round(peak_tab$sd_days[3], 2)` days against `r round(peak_tab$sd_days[2], 2)` for the fitted argmax, because it uses two coefficients rather than four and is a smooth function of them instead of a maximum. It is also biased by `r round(peak_tab$bias_days[3], 2)` days and puts `r round(100 * peak_tab$worse_8d[3], 1)` per cent of pixels more than a composite interval from the truth. The bias is less a defect than a statement that it estimates a different quantity: the phase of the first harmonic is the crest of the symmetric part of the curve, and when the second harmonic leans the curve towards a slow senescence the crest moves later while the phase stays put. Here the lean is worth `r round(abs(peak_tab$bias_days[3]), 1)` days, a fifth of the whole grassland-to-woodland difference the map is meant to show. The rule out of this section: fit the model, then find the maximum of the fitted curve numerically. Do not read the peak off the phase coefficient unless the curve is known to be symmetric, and do not take the argmax of the raw composites at all. Harmonic fitting is in any case not the usual route to a phenological date: Zhang and colleagues (2003) fit piecewise logistic functions and take the extremes of the rate of change of curvature as the transitions, which suits a fast green-up better than a low order Fourier series does. ```{r fig-peak} #| echo: false #| fig-cap: "Error in the estimated peak day for all 1600 pixels, in days, for the three estimators. Zero is the truth. The raw seasonal maximum is unbiased and very wide, because the crest of the curve is flat compared with the observation noise. The argmax of the fitted curve is unbiased and narrow. The phase of the first harmonic is narrower still and sits about ten days early, because it measures the crest of the symmetric part of the curve rather than the crest of the curve." #| fig-alt: "Three overlapping density curves on warm off-white paper, horizontal axis the error in the estimated peak day running from about minus fifty to plus fifty days, vertical axis density. A vertical dashed line stands at zero. A broad low pale sage curve labelled raw seasonal maximum spreads across the whole panel with a rounded top near zero. A tall narrow dark green curve labelled argmax of the fitted curve rises steeply at zero. An even taller and narrower gold curve labelled phase of the first harmonic peaks about ten days left of zero, with almost no overlap with the dashed line." #| fig-width: 7.5 #| fig-height: 4.3 err_df <- do.call(rbind, lapply(names(est_days), function(nm) data.frame(err = circ_gap(est_days[[nm]], peak_true), estimator = nm))) err_df$estimator <- factor(err_df$estimator, levels = names(est_days)) ggplot(err_df, aes(err, colour = estimator, fill = estimator)) + geom_density(alpha = 0.22, linewidth = 0.8, adjust = 1.2) + geom_vline(xintercept = 0, linetype = "dashed", colour = te_pal$ink, linewidth = 0.5) + scale_colour_manual(values = c(te_pal$sage, te_pal$forest, te_pal$gold), name = NULL) + scale_fill_manual(values = c(te_pal$sage, te_pal$green, te_pal$gold), name = NULL) + labs(x = "error in estimated peak day (days)", y = "density", title = "Fitting the curve beats reading the maximum") + theme_te() + theme(legend.position = "bottom", plot.margin = margin(8, 14, 4, 8)) ``` ## When the dates go missing Everything above assumed 46 evenly spaced composites. Real archives do not deliver that. Cloud, snow, sensor outages and quality flags remove observations, and they do not remove them at random across the year: in a temperate summer it is convective cloud, in a monsoon system it is a solid three months of nothing. Roerink, Menenti and Verhoef (2000) built HANTS around this problem, wrapping the same regression in an iterative rejection loop that discards observations lying below the fitted curve, on the grounds that cloud almost always depresses a vegetation index, and then predicting the missing dates from the fit. That reconstruction is the gap-filling post's subject; the question here is narrower. Three scenarios, all dropping the same number of dates so that sample size is held fixed and only the pattern changes: nine at random, nine consecutive composites lost across the spring green-up, and nine consecutive composites lost across the senescence, each run over 50 noise replicates. ```{r gap-scenarios} spring_out <- which(doy >= 97 & doy <= 161) autumn_out <- which(doy >= 193 & doy <= 257) keep_full <- seq_len(n_lyr) keep_spring <- setdiff(keep_full, spring_out) keep_autumn <- setdiff(keep_full, autumn_out) n_drop <- length(spring_out) print(c(dropped_per_scenario = n_drop, retained = n_lyr - n_drop, spring_gap_from = doy[min(spring_out)], spring_gap_to = doy[max(spring_out)], autumn_gap_from = doy[min(autumn_out)], autumn_gap_to = doy[max(autumn_out)])) gap_once <- function(y_mat, keep, n_harm) { f <- fit_harm(n_harm, keep, y_mat) c(amp_rel_bias = 100 * mean(f$amp / amp_true - 1), amp_mae = mean(abs(f$amp - amp_true)), peak_mae = mean(abs(circ_gap(f$peak, peak_true))), condition = kappa(f$des)) } n_rep_gap <- 50 set.seed(4471) gap_acc <- list() for (b in seq_len(n_rep_gap)) { y_b <- sig_obs + matrix(rnorm(n_cell * n_lyr, 0, sd_noise), n_cell, n_lyr) keep_scatter <- sort(sample.int(n_lyr, n_lyr - n_drop)) keeps <- list(complete = keep_full, scattered = keep_scatter, spring = keep_spring, autumn = keep_autumn) for (nm in names(keeps)) for (kk in 1:2) { gap_acc[[length(gap_acc) + 1]] <- data.frame( scenario = nm, harmonics = kk, t(gap_once(y_b, keeps[[nm]], kk))) } } gap_all <- do.call(rbind, gap_acc) gap_tab <- aggregate(cbind(amp_rel_bias, amp_mae, peak_mae, condition) ~ scenario + harmonics, gap_all, mean) gap_tab <- gap_tab[order(gap_tab$harmonics, match(gap_tab$scenario, c("complete", "scattered", "spring", "autumn"))), ] print(round(gap_tab[, -1], 4)) print(gap_tab$scenario) gval <- function(sc, kk, col) gap_tab[[col]][gap_tab$scenario == sc & gap_tab$harmonics == kk] ``` Dropping the dates at random costs precision and nothing else. With two harmonics the mean absolute amplitude error goes from `r sprintf("%.5f", gval("complete", 2, "amp_mae"))` on the complete year to `r sprintf("%.5f", gval("scattered", 2, "amp_mae"))` with nine scattered dates missing, and peak day error from `r round(gval("complete", 2, "peak_mae"), 2)` to `r round(gval("scattered", 2, "peak_mae"), 2)` days, while the relative bias in the amplitude is `r sprintf("%.3f", gval("scattered", 2, "amp_rel_bias"))` per cent, which is nothing. A correctly specified least squares fit does not care which subset of dates it gets, only how many and how spread out. A season-shaped gap is a different animal, and what it does depends on the number of harmonics. With two, the correct model here, the amplitude bias is still `r sprintf("%.3f", gval("spring", 2, "amp_rel_bias"))` per cent for the spring gap and `r sprintf("%.3f", gval("autumn", 2, "amp_rel_bias"))` per cent for the autumn one. What is paid instead is conditioning and precision: the condition number of the design rises from `r round(gval("complete", 2, "condition"), 2)` to `r round(gval("spring", 2, "condition"), 2)`, peak day error from `r round(gval("complete", 2, "peak_mae"), 2)` to `r round(gval("spring", 2, "peak_mae"), 2)` days, and amplitude error by a factor of `r round(gval("spring", 2, "amp_mae") / gval("complete", 2, "amp_mae"), 2)`. With one harmonic the same gap produces a bias. On the complete year the omitted second harmonic was orthogonal to the retained columns and cost nothing: `r sprintf("%.3f", gval("complete", 1, "amp_rel_bias"))` per cent. Take out the spring and the amplitude comes back `r sprintf("%.2f", gval("spring", 1, "amp_rel_bias"))` per cent too high; take out the senescence instead and it comes back `r sprintf("%.2f", gval("autumn", 1, "amp_rel_bias"))` per cent too low. Same number of observations, same noise, opposite sign of error, and the direction is set purely by which side of the crest is missing. ```{r gap-sweep} sweep_len <- 9 sweep_start <- seq(1, n_lyr, by = 2) gap_sweep <- do.call(rbind, lapply(sweep_start, function(st) { drop_idx <- ((st - 1 + seq_len(sweep_len) - 1) %% n_lyr) + 1 keep_i <- setdiff(seq_len(n_lyr), drop_idx) mid_day <- ((doy[st] + (sweep_len - 1) * (doy[2] - doy[1]) / 2 - 1) %% 365) + 1 do.call(rbind, lapply(1:2, function(kk) { f <- fit_harm(kk, keep_i, sig_obs) data.frame(mid_day = mid_day, harmonics = kk, amp_rel_bias = 100 * mean(f$amp / amp_true - 1), frac_over = mean(f$amp > amp_true)) })) })) gap_sweep <- gap_sweep[order(gap_sweep$harmonics, gap_sweep$mid_day), ] s1 <- gap_sweep[gap_sweep$harmonics == 1, ] worst_hi <- s1[which.max(s1$amp_rel_bias), ] worst_lo <- s1[which.min(s1$amp_rel_bias), ] s2 <- gap_sweep[gap_sweep$harmonics == 2, ] print(round(s1[, c("mid_day", "amp_rel_bias", "frac_over")], 4)) print(round(c(worst_over_pct = worst_hi$amp_rel_bias, worst_over_at_day = worst_hi$mid_day, worst_under_pct = worst_lo$amp_rel_bias, worst_under_at_day = worst_lo$mid_day, largest_two_harmonic_bias_pct = max(abs(s2$amp_rel_bias))), 4)) ``` Sweeping the block of `r sweep_len` missing composites right around the year, on the noiseless signal so that only the geometry is left, gives the whole picture. With one harmonic the relative amplitude bias runs from `r sprintf("%+.2f", worst_hi$amp_rel_bias)` per cent when the block is centred on day `r round(worst_hi$mid_day, 0)` to `r sprintf("%+.2f", worst_lo$amp_rel_bias)` per cent when it is centred on day `r round(worst_lo$mid_day, 0)`, and at both extremes `r round(100 * max(worst_hi$frac_over, 1 - worst_lo$frac_over), 0)` per cent of pixels are pushed the same way. With two harmonics the largest bias anywhere in the sweep is `r formatC(max(abs(s2$amp_rel_bias)), format = "e", digits = 1)` per cent, which is zero to machine precision. The mechanism generalises past this example. Least squares does not bias the coefficients of a correctly specified model whatever the design is, so every bias here is the projection of the part of the signal the model does not carry onto the columns it does. On a complete evenly spaced year that projection is exactly zero by orthogonality of the harmonic basis, which is why a model one term short can look harmless. Remove a contiguous block of dates, the columns stop being orthogonal over the retained set, and the missing term leaks into the retained ones with a sign that follows its own sign over the retained window. Missing the green-up removes the part of the year where the asymmetry pulls the curve down, so the fitted swing is too large; missing the senescence removes where it pushes the curve up, so the swing is too small. That is the finding of this section, and it is unpleasant. The number of harmonics stops being a variance question and becomes a bias question the moment the missing dates cluster in a season, which is the normal case. Choosing the model order on a complete year, or on the fit statistics, says nothing about how it will behave on the pixels where the cloud sat. ```{r fig-gaps} #| echo: false #| fig-cap: "Relative error in the estimated seasonal amplitude when a block of nine consecutive composites is missing, against the middle of the missing block, computed on the noiseless signal so that only the geometry of the design is in it. With one harmonic the error swings from over three per cent too high when the green-up is missing to over three per cent too low when the senescence is missing. With two harmonics, the correct model here, the line sits on zero throughout." #| fig-alt: "A line chart on warm off-white paper. The horizontal axis is the middle day of a missing block of nine composites, running from a few days after the new year to about day three hundred and fifty. The vertical axis is the relative amplitude error in per cent, from about minus four to plus four, with a horizontal reference line at zero. A gold line of round dots labelled one harmonic climbs from about plus half a per cent at the left edge to a hump of about plus three and a half near day thirty-five, dips to about plus one and a third near day ninety, rises to its maximum of about plus three and a half at day one hundred and forty-five, then falls steeply, crossing zero a little before day one hundred and ninety, reaching its minimum of about minus three and a half near day two hundred and twenty-five, easing back to about minus one near day two hundred and seventy, dipping again to about minus three and a third near day three hundred and twenty, and ending near minus one and a half. A dark green line labelled two harmonics is flat and indistinguishable from the zero line across the whole panel." #| fig-width: 7.5 #| fig-height: 4.2 gap_sweep$model <- factor(ifelse(gap_sweep$harmonics == 1, "one harmonic", "two harmonics"), levels = c("one harmonic", "two harmonics")) ggplot(gap_sweep, aes(mid_day, amp_rel_bias, colour = model)) + geom_hline(yintercept = 0, colour = te_pal$line, linewidth = 0.6) + geom_line(linewidth = 0.9) + geom_point(size = 1.9) + scale_colour_manual(values = c("one harmonic" = te_pal$gold, "two harmonics" = te_pal$forest), name = NULL) + labs(x = "middle day of the missing block", y = "amplitude error (per cent)", title = "Which season is missing sets the sign of the error") + theme_te() + theme(legend.position = "bottom", plot.margin = margin(8, 14, 4, 8)) ``` ## The coefficients are estimates, and they have standard errors The amplitude map is a covariate now. It will go into a species distribution model, a comparison of grazing compartments, a regression of bird abundance on seasonal productivity. What travels into that model is a single number per pixel, and that number came out of a regression with 46 observations and `r scan_out[[2]]$npar` parameters, so it has a standard error the map does not show. The amplitude is a nonlinear function of two coefficients, so the delta method supplies it. With $V$ the covariance matrix of the fitted coefficients, $$\operatorname{Var}(A) \approx \frac{a^2 V_{aa} + 2ab\,V_{ab} + b^2 V_{bb}}{a^2 + b^2}$$ and the same construction with the derivatives of $\operatorname{atan2}$ gives the standard error of the phase. ```{r amp-se} resid_full <- obs_mat - cf2 %*% t(des2) df_resid <- n_lyr - ncol(des2) s2_pix <- rowSums(resid_full^2) / df_resid xtx_inv <- solve(crossprod(des2)) v_aa <- xtx_inv["cos1", "cos1"]; v_bb <- xtx_inv["sin1", "sin1"] v_ab <- xtx_inv["cos1", "sin1"] a_hat <- cf2[, "cos1"]; b_hat <- cf2[, "sin1"] var_amp <- s2_pix * (a_hat^2 * v_aa + 2 * a_hat * b_hat * v_ab + b_hat^2 * v_bb) / amp_hat^2 se_amp <- sqrt(var_amp) var_phase <- s2_pix * (b_hat^2 * v_aa - 2 * a_hat * b_hat * v_ab + a_hat^2 * v_bb) / amp_hat^4 se_phase_days <- sqrt(var_phase) / omega cover90 <- mean(abs(amp_err) < qt(0.95, df_resid) * se_amp) print(round(c(median_se_amp = median(se_amp), min_se_amp = min(se_amp), max_se_amp = max(se_amp), actual_sd_of_amp_error = sd(amp_err), coverage_of_90pct_interval = cover90), 5)) print(round(c(median_rel_se_pct = 100 * median(se_amp / amp_hat), min_rel_se_pct = 100 * min(se_amp / amp_hat), max_rel_se_pct = 100 * max(se_amp / amp_hat), median_se_phase_days = median(se_phase_days)), 4)) ``` The median per-pixel standard error on the amplitude is `r sprintf("%.5f", median(se_amp))`, and the actual standard deviation of the amplitude error across the window is `r sprintf("%.5f", sd(amp_err))`, so the reported uncertainty is honest: a nominal 90 per cent interval covers the truth at `r round(100 * cover90, 1)` per cent of pixels. As a fraction of the amplitude being estimated it runs from `r round(100 * min(se_amp / amp_hat), 2)` per cent in the closed canopy to `r round(100 * max(se_amp / amp_hat), 2)` per cent over the thinnest grass, because the standard error is nearly constant across the window while the amplitude is not. The companion standard error on the phase has a median of `r round(median(se_phase_days), 2)` days. That heterogeneity is the reason to keep the standard error rather than admire it. A covariate measured with error attenuates the coefficient it carries, by the reliability ratio: the variance of the true covariate over the variance of the true covariate plus the measurement error variance. ```{r reliability} grass_only <- wood < 0.15 rel_all <- var(amp_true) / (var(amp_true) + mean(var_amp)) rel_grass <- var(amp_true[grass_only]) / (var(amp_true[grass_only]) + mean(var_amp[grass_only])) print(round(c(pixels_all = n_cell, pixels_grass = sum(grass_only), sd_true_amp_all = sd(amp_true), sd_true_amp_grass = sd(amp_true[grass_only])), 5)) print(round(c(reliability_all = rel_all, attenuation_pct_all = 100 * (1 - rel_all), reliability_grass = rel_grass, attenuation_pct_grass = 100 * (1 - rel_grass)), 4)) ``` Over the whole window the reliability ratio is `r round(rel_all, 4)`, so a slope fitted on the amplitude map would be attenuated by about `r round(100 * (1 - rel_all), 2)` per cent. That sounds like permission to ignore the problem, and it is not, because the ratio is a property of the sample and not of the method. Restrict the analysis to the `r sum(grass_only)` grassland pixels, where the true amplitude has a standard deviation of `r sprintf("%.5f", sd(amp_true[grass_only]))` instead of `r sprintf("%.5f", sd(amp_true))`, and the reliability falls to `r round(rel_grass, 4)`: an attenuation of `r round(100 * (1 - rel_grass), 1)` per cent, enough to turn a real effect into a marginal one. The measurement error did not change. The contrast asked of it did. Carrying the point estimate forward on its own treats a noisy quantity as known. The cheap partial fix is to carry `se_amp` alongside `amp_hat`, so the downstream model can at least weight the pixels, or report the reliability ratio and the attenuation. What that will not fix is a bias like the one in the previous section, which no standard error describes. ```{r fig-se} #| echo: false #| fig-cap: "Estimated seasonal amplitude against the truth for a systematic sample of 80 pixels, with plus and minus one standard error from the delta method, coloured by the standard error as a percentage of the amplitude. The dashed line is equality. The intervals are of nearly constant length across the range, so the same absolute uncertainty is a small fraction of a woodland amplitude and a large fraction of a grassland one." #| fig-alt: "A scatter plot on warm off-white paper with true amplitude from about 0.15 to 0.31 on the horizontal axis and estimated amplitude on the same range vertically. Eighty short vertical error bars with points at their centres lie along a dashed diagonal line of equality, scattered narrowly either side of it. The bars are all about the same length. Their colour runs from rust red at the left-hand, low amplitude end to dark green at the right-hand end, keyed to a legend giving the standard error as a percentage of the amplitude and running from about two at the dark end to about seven at the rust end." #| fig-width: 7.4 #| fig-height: 4.4 pick <- seq(1, n_cell, length.out = 80) se_df <- data.frame(truth = amp_true[pick], est = amp_hat[pick], se = se_amp[pick], rel = 100 * se_amp[pick] / amp_hat[pick]) ggplot(se_df, aes(truth, est, colour = rel)) + geom_abline(slope = 1, intercept = 0, linetype = "dashed", colour = te_pal$ink, linewidth = 0.5) + geom_linerange(aes(ymin = est - se, ymax = est + se), linewidth = 0.7) + geom_point(size = 1.8) + scale_colour_gradient(low = te_pal$forest, high = te_pal$clay, name = "standard error\n(per cent of amplitude)") + labs(x = "true seasonal amplitude", y = "estimated seasonal amplitude", title = "The same error bar means different things across the map") + theme_te() + theme(legend.position = "right", legend.title = element_text(size = 8), legend.text = element_text(size = 8), plot.margin = margin(8, 14, 4, 8)) ``` ## The honest limit The truth here was generated as a sum of two harmonics and then fitted with harmonics, which is the friendliest possible case and inflates every recovery number above. The amplitude came back to `r round(100 * median(abs(amp_err / amp_true)), 2)` per cent and the peak day to `r peak_mae` days because the model could represent the signal exactly. On real reflectance it never can, so the errors reported here are lower bounds. The one place the simulation was not friendly is the section on missing dates, where a model one term short of the truth was fitted on purpose, and that is exactly where the only real bias in the post appeared. The harmonic form itself is a strong assumption. A single cosine is symmetric about its crest, which no growing season is. Two harmonics can lean but they cannot produce a fast green-up followed by a long slow senescence and a flat winter, which is the actual shape of a deciduous canopy, and pushing the order up to get it runs into the result of the harmonics section: more terms fit the noise and move the peak. A genuine double season is the worse case, because a low order fit to two peaks does not fail loudly. It returns a single smooth hump somewhere between them, with a plausible amplitude and a peak day corresponding to nothing that happened on the ground, and nothing in the residuals of a 40 by 40 window at this noise level would make that obvious. The general warning has been measured on real data: Atkinson, Jeganathan, Dash and Atzberger (2012) compared four smoothing models, harmonic analysis among them, and found the estimated dates depend on which was used, while White and colleagues (2009) put a number on it for North America, where ten methods on the same record gave mean start-of-spring dates spanning around three weeks, with the disagreement varying by land cover. A peak day is a property of a fitted curve, and the curve is a choice. The gap analysis used contiguous blocks of missing dates, a caricature of what cloud does. Real masks are patchy, differ from pixel to pixel within a scene, and lose observations at a rate correlated with the surface state, so the missingness is not ignorable the way the random-drop scenario assumed: the direction of the bias should carry over, the size should not be read across. And nothing here has atmospheric correction, sensor drift, changing view geometry or georegistration error in it, the cube having arrived calibrated and aligned. Each of those can put a seasonal or spatial pattern into the coefficients that has nothing to do with vegetation, and a sub-pixel misregistration in particular puts a false amplitude signal along every land cover boundary, which is where the interesting ecology usually is. ## Where to go next The three coefficient rasters produced here, mean, amplitude and peak day, are covariates, and the checking post is the next step before they go anywhere: whether they vary enough to be worth modelling, and whether their errors are spatially correlated, are questions nothing above answers. If the missing dates in your own cube are the problem rather than a sideshow, the gap-filling post treats reconstruction properly. And if what you want is a phenological date to set against ground records, the phenology posts are where that comparison happens, with the warning from here attached: the satellite date and the observer date are two different measurements of a partly shared thing. ## References Jakubauskas ME, Legates DR, Kastens JH 2002 Computers and Electronics in Agriculture 37(1-3):127-139 (10.1016/S0168-1699(02)00116-3) Roerink GJ, Menenti M, Verhoef W 2000 International Journal of Remote Sensing 21(9):1911-1917 (10.1080/014311600209814) Zhang X, Friedl MA, Schaaf CB, Strahler AH, Hodges JCF, Gao F, Reed BC, Huete A 2003 Remote Sensing of Environment 84(3):471-475 (10.1016/S0034-4257(02)00135-9) White MA, de Beurs KM, Didan K, et al 2009 Global Change Biology 15(10):2335-2359 (10.1111/j.1365-2486.2009.01910.x) Verbesselt J, Hyndman R, Newnham G, Culvenor D 2010 Remote Sensing of Environment 114(1):106-115 (10.1016/j.rse.2009.08.014) Atkinson PM, Jeganathan C, Dash J, Atzberger C 2012 Remote Sensing of Environment 123:400-417 (10.1016/j.rse.2012.04.001) Zhu Z, Woodcock CE 2014 Remote Sensing of Environment 144:152-171 (10.1016/j.rse.2014.01.011) Fisher NI 1993 Statistical Analysis of Circular Data. Cambridge University Press. ISBN 978-0-521-56890-6. ## Related tutorials - [NDVI time series from a raster stack](../ndvi-time-series-from-a-raster-stack/) - [Filling cloud gaps in a satellite series](../filling-cloud-gaps-in-a-satellite-series/) - [Checking a remote sensing covariate](../checking-a-remote-sensing-covariate/) - [Phenological trends and temperature](../phenological-trends-and-temperature/)