Resampling Irregular Satellite Time Series with xarray

Resample irregular, cloud-gapped satellite time series to a regular cadence with xarray: use .resample, interpolate_na with a gap limit, and rolling composites.

TL;DR Answer

Satellite revisit is irregular and clouds punch ragged holes in the record, so an optical time series arrives on an uneven, gap-riddled temporal axis that no model can consume directly. Bin it onto a fixed calendar cadence with da.resample(time="MS").median(), fill only the short revisit gaps with da.interpolate_na(dim="time", max_gap=...), and leave long gaps as NaN so you never fabricate a rainy-season signal that was never observed. Keep the y/x dimensions and the rioxarray CRS untouched throughout. This is one building block of the wider Temporal Aggregation for Time-Series Geodata workflow.

import xarray as xr

def resample_to_cadence(da, freq="MS", reducer="median", max_gap="45D"):
    """Resample a (time, y, x) DataArray to a regular cadence with a gap limit."""
    binned = getattr(da.resample(time=freq), reducer)()
    filled = binned.interpolate_na(dim="time", method="linear", max_gap=max_gap)
    return filled

Part of: Temporal Aggregation for Time-Series Geodata

Why Irregular Satellite Series Break Model Pipelines

A Sentinel-2 tile is nominally revisited every five days, but the usable record is far sparser. Cloud masking removes whole scenes, orbit geometry shifts acquisition times, and quality filters drop partial observations. What lands in your data cube is a (time, y, x) array whose time axis has irregular spacing — two clear views one week, then nothing for six weeks. Most machine learning models, and nearly every temporal feature extractor, assume a fixed sampling interval. Feed them a ragged axis and one of three things happens: the code errors on shape mismatches, silently pads with garbage, or learns the sampling pattern instead of the land surface signal.

The corruption is subtle because the raw values are all real. The problem is that observation density correlates with weather, and weather correlates with the thing you are trying to predict. A pixel with many clear observations in a dry month and almost none in a wet month encodes a seasonal bias that has nothing to do with vegetation or soil. When you compute a rolling mean or a seasonal amplitude over an uneven axis, dense periods dominate the statistic and sparse periods vanish. The model appears to fit, then fails on a region or year with a different cloud regime.

Resampling to a regular grid fixes this by decoupling the feature cadence from the acquisition cadence. Every output timestamp is a fixed calendar bin, and every bin is reduced the same way whether it held one observation or ten. The steps most affected downstream are:

  • Temporal feature extraction: Fourier terms, seasonal amplitudes, and slopes are only meaningful on an evenly spaced axis.
  • Multi-sensor fusion: Aligning Sentinel-2 with Landsat or radar requires both on the same regular grid first.
  • Sequence models: LSTMs and temporal transformers assume a constant step size between elements of the sequence.

This is the natural next step after aggregating daily satellite data to monthly features: once you know how to collapse a dense record, you need the same discipline for a sparse, uneven one.

Core Principles for Sound Temporal Resampling

  • Bin to a calendar frequency, not a rolling window. resample(time="MS") groups by month-start regardless of how many observations fall in each month. This is what converts an irregular axis into a regular DatetimeIndex.
  • Prefer median over mean for optical data. A single unmasked cloud edge or shadow inflates a monthly mean. The median is robust to those residual outliers that survive the cloud mask.
  • Never interpolate across long gaps. Pass max_gap to interpolate_na so only short revisit gaps are filled. A three-month gap has no honest linear path and must stay NaN.
  • Resample the time axis only. Leave y and x alone. Reducing along time keeps the spatial grid, CRS, and affine transform intact under the rioxarray .rio accessor.
  • Keep NaN meaningful. After filling, the remaining NaN values are a signal — they mark where the sensor genuinely saw nothing. Carry them forward as an explicit mask rather than zero-filling.
  • Validate before you write. Assert the output index is monotonic and regularly spaced, and that no NaN run exceeds the gap limit, before persisting features.

Production-Ready Code

The function below resamples a (time, y, x) DataArray to a fixed cadence, applies a robust reducer, fills only short gaps, and preserves the rioxarray spatial metadata. It validates inputs and re-stamps the CRS so the result is safe to write as a stack of GeoTIFFs.

import logging
import numpy as np
import xarray as xr
import rioxarray  # noqa: F401 — registers the .rio accessor

logger = logging.getLogger(__name__)


def resample_time_series(
    da: xr.DataArray,
    freq: str = "MS",
    reducer: str = "median",
    max_gap: str = "45D",
    interp_method: str = "linear",
) -> xr.DataArray:
    """Resample an irregular satellite series to a regular cadence with a gap limit.

    Parameters
    ----------
    da : xr.DataArray
        Array with dims (time, y, x). ``time`` may be irregularly spaced.
    freq : str, default "MS"
        Pandas offset alias for the target cadence (e.g. "MS" month-start, "10D").
    reducer : str, default "median"
        Reduction applied within each bin. One of "median", "mean", "max".
    max_gap : str, default "45D"
        Longest gap (as a timedelta string) that interpolation is allowed to bridge.
        Gaps longer than this stay NaN.
    interp_method : str, default "linear"
        Interpolation method passed to interpolate_na.

    Returns
    -------
    xr.DataArray
        Array on a regular DatetimeIndex with the same y/x grid and CRS.
    """
    if "time" not in da.dims:
        raise ValueError("DataArray must have a 'time' dimension.")
    if not {"y", "x"}.issubset(da.dims):
        raise ValueError("DataArray must have spatial 'y' and 'x' dimensions.")
    if reducer not in {"median", "mean", "max"}:
        raise ValueError(f"Unsupported reducer: {reducer!r}")
    if da.sizes["time"] < 2:
        raise ValueError("Need at least two time steps to resample.")

    # Preserve spatial metadata before any operation that might drop it.
    crs = da.rio.crs

    # 1. Bin the irregular axis onto a fixed calendar cadence.
    resampler = da.resample(time=freq)
    binned = getattr(resampler, reducer)()
    logger.info(
        "Resampled %d irregular steps to %d regular bins at freq=%s",
        da.sizes["time"], binned.sizes["time"], freq,
    )

    # 2. Fill only short gaps; long gaps remain NaN by design.
    filled = binned.interpolate_na(
        dim="time", method=interp_method, max_gap=np.timedelta64(_to_days(max_gap), "D")
    )

    # 3. Re-stamp the CRS in case a reduction cleared the .rio metadata.
    if crs is not None:
        filled = filled.rio.write_crs(crs)

    filled.attrs.update(da.attrs)
    return filled


def _to_days(gap: str) -> int:
    """Convert a simple '<n>D' timedelta string to an integer number of days."""
    gap = gap.strip().upper()
    if not gap.endswith("D"):
        raise ValueError("max_gap must be given in days, e.g. '45D'.")
    return int(gap[:-1])

interpolate_na accepts max_gap as a numpy.timedelta64, so the helper keeps the public API a friendly string while passing the exact type xarray expects. If you prefer, pass a pandas offset directly — the key point is that the gap limit is explicit and auditable rather than an implicit “fill everything”.

Step-by-Step Walkthrough: An Irregular NDVI Cube to Monthly Composites

1. Load the cube and inspect the time axis

Start from a cloud-masked NDVI stack. If you have not masked clouds yet, do that first — resampling over cloud-contaminated pixels poisons every bin. See masking clouds before computing NDVI for the scene-classification approach, and how to calculate NDVI and EVI with rasterio for building the index itself.

import xarray as xr
import numpy as np

da = xr.open_dataarray("ndvi_masked_2025.nc")  # dims: (time, y, x)

# Inspect the spacing between consecutive acquisitions.
gaps = np.diff(da["time"].values).astype("timedelta64[D]")
print("acquisitions:", da.sizes["time"])
print("median gap:", np.median(gaps), "max gap:", gaps.max())
acquisitions: 41
median gap: 5 days  max gap: 38 days

Forty-one clear observations over a year, spaced anywhere from 5 to 38 days apart — a classic irregular optical record.

2. Resample to month-start with a median reducer

monthly = da.resample(time="MS").median()
print(monthly["time"].values[:3])
['2025-01-01T00:00:00.000000000'
 '2025-02-01T00:00:00.000000000'
 '2025-03-01T00:00:00.000000000']

The output has exactly twelve steps on a clean month-start index. A month with no clear acquisition — say a fully clouded July — becomes an all-NaN slice rather than being skipped, which keeps the axis regular.

3. Fill short gaps, keep long gaps as NaN

filled = monthly.interpolate_na(
    dim="time", method="linear", max_gap=np.timedelta64(45, "D")
)

With monthly bins, a single missing month (about 31 days) is bridged, but two consecutive missing months (about 62 days) exceed the 45-day limit and stay NaN. That boundary is a design choice: widen max_gap if your model tolerates more inference, tighten it if you need every feature grounded in a real observation.

4. Optionally smooth with a rolling composite

When you want a less jumpy signal for slope or amplitude features, chain a centred rolling median over the regular axis. Because the axis is now regular, the window is unambiguous.

smoothed = filled.rolling(time=3, center=True, min_periods=2).median()

5. Write the regular stack back out

smoothed.rio.write_crs(da.rio.crs, inplace=True)
smoothed.rio.to_raster("ndvi_monthly_2025.tif")

The CRS and affine transform ride along because we never touched y or x. The result is a twelve-band GeoTIFF on a regular monthly grid, ready for feature extraction.

Visualising the Resampling Flow

The diagram below traces the transformation: an irregular, gap-riddled input axis on the left, calendar binning in the middle, and a regular output cadence on the right where short gaps are interpolated and a long gap is deliberately left empty.

Resampling an irregular satellite time series to a regular cadence Irregularly spaced acquisitions on the left are grouped into fixed monthly bins in the middle, producing a regular output axis on the right where a one-month gap is interpolated and a two-month gap stays empty. Irregular input resample(time="MS") Regular output long cloud gap Jan: median(3) Feb: median(3) Mar: NaN Apr: NaN May: median(2) interp NaN kept Short gap (Mar) is interpolated under max_gap; long gap (Apr) exceeds the limit and stays NaN.

Verification

Never trust that a resample “worked” — assert the two properties a downstream model relies on: a strictly regular DatetimeIndex, and no NaN run longer than the gap limit. The checks below fail loudly so a bad cube can never reach feature extraction.

import numpy as np
import pandas as pd
import xarray as xr


def assert_regular_cadence(da: xr.DataArray, freq: str = "MS") -> None:
    """Assert the time axis is a regular DatetimeIndex at the expected frequency."""
    idx = pd.DatetimeIndex(da["time"].values)
    expected = pd.date_range(idx[0], idx[-1], freq=freq)
    if not idx.equals(expected):
        raise AssertionError(
            f"Time axis is not regular at freq={freq}: "
            f"{len(idx)} steps vs {len(expected)} expected."
        )


def assert_gap_limit(da: xr.DataArray, max_gap_steps: int) -> None:
    """Assert no all-NaN run along time exceeds max_gap_steps consecutive bins."""
    all_nan = da.isnull().all(dim=("y", "x")).values  # True where a whole slice is NaN
    run = 0
    for is_nan in all_nan:
        run = run + 1 if is_nan else 0
        if run > max_gap_steps:
            raise AssertionError(
                f"Found a NaN run of {run} bins, exceeding the "
                f"allowed {max_gap_steps}."
            )


assert_regular_cadence(filled, freq="MS")
assert_gap_limit(filled, max_gap_steps=1)  # monthly bins, 45D gap ≈ 1 bin
print("Resampled cube passes cadence and gap-limit checks.")

The assert_regular_cadence check catches a silently dropped or duplicated bin, and assert_gap_limit confirms that interpolation respected the ceiling you set — long gaps really did stay empty. Wire both into your pipeline tests so a change in cloud cover or source cadence surfaces as a failed assertion, not a quietly biased feature. From here, the regular cube feeds straight into the seasonal and rolling statistics covered across Temporal Aggregation for Time-Series Geodata.

FAQ

Why not just interpolate every gap to fill the whole series?

Interpolating across arbitrarily long gaps invents data the sensor never observed. A three-month cloud gap over a monsoon forest has no reliable linear path between the last clear observation and the next. Passing max_gap to interpolate_na caps how far interpolation reaches, so short revisit gaps are filled but long ones stay NaN and can be masked or dropped downstream. This keeps fabricated values out of your training features.

Does xarray resample handle irregularly spaced timestamps?

Yes. xarray’s .resample groups observations into fixed calendar bins (for example "MS" for month-start) regardless of how unevenly the source timestamps are spaced. Two acquisitions in one week and none the next both fall into the correct monthly bin. The reducer you chain — mean, median, or max — collapses each bin to a single value on a regular DatetimeIndex, which is exactly the grid a model expects.

How do I keep the spatial CRS and transform intact during resampling?

Resample only along the time dimension and never rename or drop the y and x coordinates. rioxarray stores the CRS and affine transform in the .rio accessor, which survives time-axis reductions because those operations leave the spatial dims untouched. After resampling, call rio.write_crs to re-stamp the CRS if any operation cleared it, then verify with da.rio.crs before writing a GeoTIFF.

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