Filling Gaps in NDVI Time Series with Interpolation

To fill short cloud gaps without inventing a season, interpolate along the time axis with an explicit maximum gap:

import xarray as xr

regular = ndvi.resample(time="5D").mean()             # explicit, regular time axis
filled = regular.interpolate_na(dim="time", method="time", max_gap="15D")
was_filled = regular.isnull() & filled.notnull()      # the flag consumers need

Gap filling is the step that makes the series in Temporal Aggregation and Time-Series Analysis usable — and the step most likely to make it dishonest.

Why This Arises in Remote Sensing Workflows

Optical time series are full of holes. Sentinel-2 revisits every five days, but cloud removes half the observations in temperate regions and far more in the tropics, and the remaining acquisitions are irregularly spaced. Most downstream methods — phenology metrics, harmonic fits, change detection, anything that differences neighbouring dates — assume a regular series.

Filling those holes is therefore routine, and it is where remote sensing products most often stop being measurements. An interpolated point looks exactly like an observed one: same dtype, same range, same plot. Unless it is flagged, nothing downstream can tell the difference, and a phenology metric computed from a series that is 60 percent interpolated is reporting the interpolation method as much as the vegetation.

The honest approach is not to avoid filling but to bound it: interpolate only across gaps short enough that the underlying process really is smooth, leave longer gaps empty, and carry a flag so the distinction survives.

Where interpolation is information and where it is invention Across a ten-day gap in mid-season, a straight line closely follows the true curve. Across a two-month gap spanning green-up, the same straight line cuts the corner of the seasonal rise and reports a trajectory the vegetation never followed. Same method, two gap lengths 0.8 0.1 true seasonal curve (dashed) 10-day gap: line ≈ truth 60-day gap across green-up: the line misses the whole curve shape Jan Jun Dec max_gap is the parameter that separates these two cases — without it, both are filled identically. Dots are real observations; everything between them is a modelling choice.

Environment & Setup

Package Version Why
xarray ≥2023.6 resample, interpolate_na with max_gap
pandas ≥2.0 Datetime handling behind the time axis
numpy ≥1.23 Array operations and flags
scipy ≥1.11 Optional: Savitzky-Golay smoothing after filling 
pip install "xarray>=2023.6" "pandas>=2.0" "scipy>=1.11"

Complete Working Example

This function regularises the series, fills only short gaps, optionally fills longer ones from a climatology, and returns a flag array alongside the result.

Which fill method for which gap Short gaps in a smooth signal take time-weighted interpolation. Gaps spanning a phenological transition need a shorter limit or no fill at all. Long gaps take a climatological substitute, and systematically empty bins should be left empty and documented. Gap length decides the method, not convenience gap method flag value ≤ 15 days, mid-season time-weighted interpolation 1 — interpolated ≤ 15 days across green-up shorten the limit, or leave 1 or 3, decided per site 15–60 days day-of-year climatology 2 — climatology every year, same months leave empty, document why 3 — missing, deliberately The flag column is what lets a consumer reconstruct which rows were observations. 
import numpy as np
import xarray as xr


def fill_series(
    ndvi: xr.DataArray,
    *,
    step: str = "5D",
    max_gap: str = "15D",
    climatology: xr.DataArray | None = None,
    clim_max_gap: str = "60D",
) -> tuple[xr.DataArray, xr.DataArray]:
    """Return (filled series, fill flag) where flag is 0 observed, 1 interpolated, 2 climatology."""
    regular = ndvi.resample(time=step).mean(skipna=True)

    observed = regular.notnull()

    # Time-weighted interpolation respects irregular spacing within the regular axis
    interpolated = regular.interpolate_na(dim="time", method="time", max_gap=max_gap)
    flag = xr.where(observed, 0, xr.where(interpolated.notnull(), 1, np.nan))

    filled = interpolated
    if climatology is not None:
        # Only where a longer gap remains: substitute the day-of-year climatology
        doy = filled["time"].dt.dayofyear
        clim_aligned = climatology.sel(dayofyear=doy, method="nearest")
        still_missing = filled.isnull()
        long_enough = still_missing & _gap_length(regular) <= xr.DataArray(
            np.timedelta64(int(clim_max_gap.rstrip("D")), "D"))
        filled = xr.where(still_missing, clim_aligned, filled)
        flag = xr.where(still_missing & filled.notnull(), 2, flag)

    flag = flag.fillna(3)                     # 3 = still missing, deliberately
    filled.attrs.update(fill_max_gap=max_gap, fill_step=step,
                        flag_meaning="0 observed, 1 interpolated, 2 climatology, 3 missing")
    return filled.astype("float32"), flag.astype("uint8")


def _gap_length(series: xr.DataArray) -> xr.DataArray:
    """Length of the NaN run each timestep belongs to, as a timedelta."""
    valid = series.notnull()
    # Forward and backward fill of the last valid timestamp gives the enclosing gap
    times = series["time"]
    last_valid = xr.where(valid, times, np.datetime64("NaT")).ffill("time")
    next_valid = xr.where(valid, times, np.datetime64("NaT")).bfill("time")
    return (next_valid - last_valid)


if __name__ == "__main__":
    filled, flag = fill_series(ndvi_cube, step="5D", max_gap="15D")

    total = int(flag.size)
    print("observed:", int((flag == 0).sum()) / total)
    print("interpolated:", int((flag == 1).sum()) / total)
    print("still missing:", int((flag == 3).sum()) / total)

    filled.rio.to_raster("ndvi_filled.tif", driver="COG", compress="DEFLATE")
    flag.rio.to_raster("ndvi_fill_flag.tif", driver="COG", compress="DEFLATE", dtype="uint8")

Writing the flag as its own raster is the part that is easy to skip and expensive to add later. It costs one uint8 layer per date and it is the only thing that lets a consumer — or a reviewer — distinguish an observation from a reconstruction.

Variant Patterns

1. Climatological fill for long gaps

Where a gap spans a whole season, a straight line is indefensible but a long-term average for that day of year is often reasonable, because it converges on typical behaviour rather than drifting.

# Build a day-of-year climatology from several years of the same pixel stack
climatology = (history.groupby("time.dayofyear").median("time")
                      .rolling(dayofyear=5, center=True, min_periods=1).mean())

filled, flag = fill_series(ndvi_cube, max_gap="15D", climatology=climatology)

The climatology should be smoothed, or day-to-day noise in the historical sample becomes structure in the fill. It should also be flagged separately from interpolation, because its error behaves differently: interpolation errs toward the local trajectory, climatology toward the average year.

2. Choosing the interpolation method

Four interpolation methods on the same gap Linear interpolation ignores the actual spacing of observations and is wrong when they are irregular. Time-weighted interpolation accounts for the spacing and is the safe default. Nearest produces steps. Spline follows curvature but can overshoot beyond the observed range, which for a bounded index is a defect. Same two anchor points, four fills method behaviour use when method="time" weights by real spacing the default for irregular series method="linear" weights by index position only on an already-regular axis method="nearest" step function categorical series, never NDVI method="cubic" / spline follows curvature, can overshoot clip to [−1, 1] afterwards An index bounded by ±1 must never leave that range: clip after any higher-order interpolation.

3. Filling in space rather than time

For a single date with small cloud holes, spatial interpolation from surrounding pixels can be better than temporal interpolation, because neighbouring land is often more similar than the same pixel two weeks earlier. It is only defensible for small holes in homogeneous cover, and it should be flagged distinctly, since its failure mode — smearing a field boundary — is spatial rather than temporal.

Reporting a Filled Series

Three fields make a filled series interpretable, and all three come free from the code above.

The fill fraction per pixel, which is the share of timesteps that were reconstructed. A phenology metric computed where that fraction exceeds roughly half is describing the fill method, and consumers need to be able to exclude those pixels.

The maximum gap actually filled, as distinct from the configured limit. A series with a 15-day limit whose longest filled gap was 5 days is in much better shape than one that used the full allowance everywhere.

The method per filled point, which is what the flag layer encodes. Interpolated and climatological values have different uncertainty, and collapsing them into one “filled” category discards that.

Publishing these alongside the series is the difference between a product a reviewer can assess and one they have to take on trust — the same argument for carrying observation counts made in Creating Monthly NDVI Composites with xarray resample.

Common Errors

interpolate_na fills the entire year

max_gap was not set, so every gap is filled regardless of length. Always set it, and set it from the process being measured rather than from convenience.

Filled values fall outside ±1

A cubic or spline method overshot. Clip after interpolating, and prefer method="time" for bounded indices.

The series is filled but downstream results got worse

Interpolated points are being used by an algorithm that assumes observations — a threshold-crossing detector, for example. Feed such algorithms the gapped series and the flag, not the filled one.

Frequently Asked Questions

Q: How long a gap can I interpolate across? As long as the underlying process is smooth over that interval. For canopy greenness in mid-season, two to three weeks is defensible; across a green-up or senescence transition, even ten days can invent a trajectory the plant never followed.

Q: Should interpolated values go into the analytical product? Preferably not. Keep the gapped series as the analytical layer and publish the filled version as a clearly labelled derived product, because algorithms that threshold values will otherwise trigger on points with no observation behind them.

Q: Is Savitzky-Golay smoothing a gap-filling method? No. It smooths an already-complete series and will happily smooth across NaNs by ignoring them, which quietly changes the shape of the curve. Fill first, with a method and a limit you can defend, then smooth if the analysis needs it.