Filling Gaps in Satellite Time Series Before Aggregation

Fill cloud gaps in satellite time series with xarray: interpolation, harmonic fitting and Savitzky-Golay smoothing, plus a validity mask so imputed values never masquerade as observations.

Cloud removes 30–70% of optical acquisitions in most of the world, so a “monthly mean NDVI” is really a mean over whichever days happened to be clear. Before aggregating, fill the gaps deliberately — linear interpolation for short ones, a harmonic fit for long ones — and carry a mask that records what was filled. The mask matters as much as the fill: without it, a pixel reconstructed from two clear scenes is indistinguishable from one observed thirty times.

This page covers gap filling. For the aggregation that follows, see aggregating daily satellite data to monthly features and the wider workflow in Temporal Aggregation for Time Series Geodata.

Across a long gap, a straight line removes the season A vegetation index time series rises to a summer peak and falls again. A cloud gap covers the peak. Linear interpolation draws a straight chord between the last clear observation before the gap and the first after it, cutting the peak off entirely. A harmonic fit follows the seasonal curve through the gap and recovers a peak close to the true one. The gap always falls where it hurts — peak growing season is the cloudiest NDVI MarMayJulSepNov cloud gap linear: peak lost harmonic: peak recovered Peak NDVI and its date are usually the two strongest features a crop model has.

Why This Fails in Geospatial ML Pipelines

Cloud is not random. It clusters in the growing season, in mountains, and in the tropics — exactly where and when the signal you want is strongest. So a naive mean over available observations is biased in a spatially structured way: dry lowland pixels are averaged over many clear days, humid upland pixels over a handful of unusual clear ones. A model then learns “cloudiness” as if it were a land-cover property, and any spatial cross-validation that separates regions will expose it as an unexplained regional performance gap.

The second failure is chord-cutting. Linear interpolation across a six-week gap draws a straight line between the last pre-gap and the first post-gap observation. If that gap covers the seasonal peak — which it usually does, because peak growth coincides with convective cloud — the reconstructed series has no peak at all. Features derived from it, such as maximum NDVI or the day of maximum, are then systematically wrong in the direction that most reduces class separability.

The third is provenance loss. Once a filled array is written to disk as a plain float raster, nothing distinguishes an observation from an estimate. Downstream aggregation weights them identically, the model treats them identically, and no drift monitor can tell that a region’s predictions degraded because its imagery got cloudier rather than because the world changed.

Core Principles

  • Mask before filling. Cloud, shadow, snow and saturated pixels are not data.
  • Choose the filler by gap length. Interpolate short gaps; fit a seasonal model across long ones.
  • Never extrapolate past the series. Leading and trailing gaps stay NaN.
  • Refuse to fill hopeless pixels. Below a minimum valid count, emit NaN and let the model handle missing.
  • Emit a companion mask and a valid count. Propagate both through every aggregation.
  • Fit the seasonal model per pixel, not per scene. Phenology varies with elevation and crop.

Production-Ready Code

from __future__ import annotations

import logging
import numpy as np
import xarray as xr

logger = logging.getLogger(__name__)


def mask_invalid(da: xr.DataArray, scl: xr.DataArray) -> xr.DataArray:
    """Set cloud, shadow, snow and saturated pixels to NaN using the scene classification."""
    bad = scl.isin([0, 1, 3, 8, 9, 10, 11])          # nodata, saturated, shadow, cloud, snow
    out = da.where(~bad)
    logger.info("masked %.1f%% of observations", 100 * float(bad.mean()))
    return out


def harmonic_design(doy: np.ndarray, n_harmonics: int = 2) -> np.ndarray:
    """Design matrix for a seasonal harmonic model: intercept plus sin/cos pairs."""
    t = 2 * np.pi * doy / 365.25
    cols = [np.ones_like(t)]
    for k in range(1, n_harmonics + 1):
        cols += [np.sin(k * t), np.cos(k * t)]
    return np.column_stack(cols)


def fill_pixel(values: np.ndarray, doy: np.ndarray, max_linear_gap: int = 2,
               min_valid: int = 6, n_harmonics: int = 2) -> tuple[np.ndarray, np.ndarray]:
    """Fill one pixel's series. Returns (filled values, filled-flag mask).

    Short runs of missing steps are interpolated; longer ones are predicted from a
    harmonic fit over the observed steps. Leading and trailing gaps are left NaN.
    """
    valid = np.isfinite(values)
    filled = values.astype("float64").copy()
    was_filled = np.zeros_like(valid)

    if valid.sum() < min_valid:
        return np.full_like(filled, np.nan), was_filled

    # 1. Linear interpolation for interior gaps, applied everywhere first.
    idx = np.arange(len(values))
    interp = np.interp(idx, idx[valid], values[valid])

    # 2. Identify run lengths of missing steps.
    gaps = ~valid
    runs, start = [], None
    for i, g in enumerate(gaps):
        if g and start is None:
            start = i
        elif not g and start is not None:
            runs.append((start, i))
            start = None
    if start is not None:
        runs.append((start, len(gaps)))

    # 3. Fit harmonics once; use them for the long gaps only.
    design = harmonic_design(doy[valid], n_harmonics)
    coeffs, *_ = np.linalg.lstsq(design, values[valid], rcond=None)
    seasonal = harmonic_design(doy, n_harmonics) @ coeffs

    for a, b in runs:
        if a == 0 or b == len(values):
            continue                                   # never extrapolate
        source = interp if (b - a) <= max_linear_gap else seasonal
        filled[a:b] = source[a:b]
        was_filled[a:b] = True

    filled[~valid & ~was_filled] = np.nan
    return filled, was_filled


def fill_stack(da: xr.DataArray, max_linear_gap: int = 2,
               min_valid: int = 6) -> tuple[xr.DataArray, xr.DataArray]:
    """Apply fill_pixel across a (time, y, x) DataArray, returning values and a mask."""
    doy = da["time"].dt.dayofyear.values.astype("float64")

    def _apply(series):
        out, flag = fill_pixel(series, doy, max_linear_gap, min_valid)
        return np.concatenate([out, flag.astype("float64")])

    stacked = xr.apply_ufunc(
        _apply, da,
        input_core_dims=[["time"]], output_core_dims=[["time2"]],
        vectorize=True, dask="parallelized",
        dask_gufunc_kwargs={"output_sizes": {"time2": 2 * da.sizes["time"]}},
        output_dtypes=["float64"],
    )
    n = da.sizes["time"]
    values = stacked.isel(time2=slice(0, n)).rename({"time2": "time"}).assign_coords(time=da.time)
    flags = stacked.isel(time2=slice(n, None)).rename({"time2": "time"}).assign_coords(time=da.time)
    return values, flags.astype("bool")

Step-by-Step Walkthrough

Step 1 — Load and mask

import rioxarray
import xarray as xr

ndvi = xr.open_dataarray("ndvi_2026.nc", chunks={"time": -1, "y": 512, "x": 512})
scl = xr.open_dataarray("scl_2026.nc", chunks={"time": -1, "y": 512, "x": 512})
clean = mask_invalid(ndvi, scl)

valid_count = clean.notnull().sum("time")
print(f"median valid observations per pixel: {float(valid_count.median()):.0f} "
      f"of {clean.sizes['time']}")

Step 2 — Look at the gap structure before choosing a strategy

import numpy as np

sample = clean.isel(y=slice(0, 200), x=slice(0, 200)).compute()
per_pixel_gaps = (~np.isfinite(sample)).sum("time")
print(f"pixels with zero gaps: {float((per_pixel_gaps == 0).mean()):.1%}")
print(f"pixels below the min_valid floor: {float((valid_count < 6).mean()):.1%}")

If a large share of pixels sit below the floor, gap filling is not the right tool — the answer is a longer compositing window or a radar source that sees through cloud.

Step 3 — Fill, and keep the mask

filled, filled_flag = fill_stack(clean, max_linear_gap=2, min_valid=6)
filled = filled.compute()
filled_flag = filled_flag.compute()

print(f"{float(filled_flag.mean()):.1%} of all time steps were reconstructed")

Step 4 — Aggregate, carrying provenance

monthly = filled.resample(time="1MS").mean()
monthly_obs = (~filled_flag & filled.notnull()).resample(time="1MS").sum()
monthly_filled_frac = filled_flag.resample(time="1MS").mean()

features = xr.Dataset({
    "ndvi_mean": monthly,
    "ndvi_n_obs": monthly_obs,             # how many REAL observations backed it
    "ndvi_filled_frac": monthly_filled_frac,
})
features.to_netcdf("monthly_features.nc")

Giving the model ndvi_n_obs alongside ndvi_mean lets it learn to discount reconstructed months — and gives you a covariate that a drift monitor can watch, in the sense described in model drift detection for geospatial inference.

Two identical means, two very different levels of evidence Two pixels report the same monthly mean vegetation index. The first is backed by twenty-eight clear observations, the second by two clear observations and a reconstruction of the rest. Without an observation count column the model cannot distinguish them. Carry the evidence, not only the estimate pixelndvi_meanndvi_n_obsfilled_frac lowland A0.62 280.07 upland B0.62 20.93 Drop the last two columns and these rows are identical to the model. Keep them and the model can learn to trust one more than the other — and a monitor can alert when the cloud regime shifts.

Verification

import numpy as np


def test_short_gap_is_interpolated_linearly():
    doy = np.arange(1, 366, 10, dtype="float64")
    values = 0.5 + 0.3 * np.sin(2 * np.pi * doy / 365.25)
    obs = values.copy()
    obs[10] = np.nan                        # one missing step
    filled, flag = fill_pixel(obs, doy, max_linear_gap=2)
    assert flag[10] and abs(filled[10] - values[10]) < 0.02


def test_long_gap_uses_the_seasonal_shape():
    """A harmonic fill must recover a peak that linear interpolation would flatten."""
    doy = np.arange(1, 366, 10, dtype="float64")
    values = 0.3 + 0.4 * np.sin(2 * np.pi * (doy - 80) / 365.25)
    obs = values.copy()
    peak = np.arange(14, 24)
    obs[peak] = np.nan
    filled, flag = fill_pixel(obs, doy, max_linear_gap=2)
    linear = np.interp(np.arange(len(obs)), np.flatnonzero(np.isfinite(obs)),
                       obs[np.isfinite(obs)])
    assert filled[peak].max() > linear[peak].max(), "harmonic fill should keep the peak"


def test_leading_and_trailing_gaps_stay_nan():
    doy = np.arange(1, 366, 10, dtype="float64")
    obs = np.full(len(doy), 0.5)
    obs[:3] = np.nan
    obs[-3:] = np.nan
    filled, flag = fill_pixel(obs, doy)
    assert np.isnan(filled[:3]).all() and np.isnan(filled[-3:]).all()
    assert not flag[:3].any()


def test_hopeless_pixel_returns_all_nan():
    doy = np.arange(1, 366, 10, dtype="float64")
    obs = np.full(len(doy), np.nan)
    obs[:3] = 0.4
    filled, flag = fill_pixel(obs, doy, min_valid=6)
    assert np.isnan(filled).all() and not flag.any()

For a data-level check, hold out a set of clear observations, fill as if they were cloudy, and compare. The residual distribution tells you the real uncertainty of your reconstruction — a number worth reporting alongside the features, and one that is usually much larger than teams expect for long gaps.

Measure your fill error rather than assuming it Reconstruction error is plotted against gap length in acquisitions. Linear interpolation error rises steeply once gaps exceed two or three acquisitions. The harmonic fit's error rises much more slowly, so the two curves cross at a small gap length, which is where the switch should be made. Hold-out error against gap length gap length (acquisitions) → RMSE linear harmonic switch here The crossing point is a property of your phenology, not a universal constant.

FAQ

Should I fill gaps before or after aggregating to monthly features?

Before, provided the month has enough real observations to justify it. Aggregating first hides how much of each month was cloud. Either way, carry a valid-observation count so the model can discount reconstructed months.

Is linear interpolation good enough for a vegetation index?

For gaps of one or two acquisitions, yes. Across weeks it cuts the corner off the phenological curve and under-estimates the seasonal peak — precisely the feature most crop models rely on. Switch to a harmonic fit beyond the crossing point you measure.

How do I stop imputed values from being treated as observations?

Emit a companion mask, propagate it through aggregation, and expose a per-feature observation count. Without it, a pixel reconstructed from two clear scenes looks exactly as reliable as one observed thirty times.

Part of: Temporal Aggregation for Time Series Geodata Part of: Spatial Feature Engineering for Machine Learning