XArray GroupBy for Temporal Raster Aggregation: Monthly Means at Scale

A decade of Landsat 8 imagery over eastern Australia: roughly 12,000 scenes, each a 7.5GB GeoTIFF stack at 30m resolution. Total: approximately 90TB. Loading that into a NumPy array is not an option. Even loading a single band for a single year (1.2TB) would exhaust most workstation RAM.

XArray solves this with lazy evaluation backed by Dask. You can define aggregation operations over the entire dataset — monthly means, annual maxima, decadal trends — and have them execute in chunks that fit in memory, in parallel across CPU cores, without writing intermediate files.

# The Problem: Temporal Raster Stacks

A “temporal raster stack” is a 3D array with dimensions (time, y, x) — or 4D with (time, band, y, x). The time dimension represents acquisition dates; each (y, x) position is a pixel. Common use cases:

NDVI time series: compute monthly mean vegetation index from daily Landsat/Sentinel scenes
Sea surface temperature: monthly climatologies from decades of MODIS data
Rainfall accumulation: aggregate daily CHIRPS rainfall to monthly or seasonal totals
Land cover change: identify pixels that changed class between years

These operations follow a common pattern: group by time period → apply aggregation → return reduced array.

# Setting Up a Temporal Stack with XArray and Dask

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

# ── Simulate a 3-year daily NDVI stack: (1095, 2000, 2000) ────────────────────
# In practice, load with stackstac or xr.open_mfdataset from GeoTIFF/Zarr
time_index = pd.date_range("2020-01-01", periods=365*3, freq="D")

# Lazy Dask-backed array — chunks of 30 days × 500px × 500px
ndvi = xr.DataArray(
    data=np.random.default_rng(0).uniform(-0.1, 0.9, (len(time_index), 2000, 2000)).astype(np.float32),
    dims=["time", "y", "x"],
    coords={
        "time": time_index,
        "y": np.linspace(-27.0, -30.0, 2000),
        "x": np.linspace(152.5, 155.5, 2000),
    },
    name="ndvi"
)

# Chunk it for lazy computation (in production, load directly as chunked)
ndvi_dask = ndvi.chunk({"time": 30, "y": 500, "x": 500})
print(ndvi_dask)
# <xarray.DataArray 'ndvi' (time: 1095, y: 2000, x: 2000)>
# dask.array<...> chunks=(30, 500, 500)
# Size: ~8.8 GB on disk, loaded lazily

# GroupBy: Computing Monthly Means

# ── Monthly mean NDVI — entirely lazy until .compute() ────────────────────────
monthly_mean = ndvi_dask.groupby("time.month").mean(dim="time")
# Result shape: (12, 2000, 2000) — one layer per calendar month
# This is still lazy — no computation has run

# Trigger computation
import time
t0 = time.perf_counter()
result = monthly_mean.compute()
elapsed = time.perf_counter() - t0
print(f"Monthly means (3yr × 2000×2000 daily): {elapsed:.1f}s")
# → ~8.2s on 8-core machine with Dask threading scheduler
print(result.shape)   # (12, 2000, 2000)

The .groupby("time.month") accessor groups by the month component of the time coordinate (1–12), then .mean(dim="time") reduces across all times within each group. The result is a 12-layer array of climatological monthly means.

# Resample: Calendar-Consistent Aggregation

resample is a time-aware GroupBy that uses Pandas time offsets, producing one output per calendar period:

# ── Monthly mean using resample (more calendar-accurate than groupby.month) ───
monthly_resampled = ndvi_dask.resample(time="1MS").mean()   # "1MS" = month start
# Result: one layer per month, labelled with the first day of each month
# e.g. 2020-01-01, 2020-02-01, ..., 2022-12-01 → 36 layers

# ── Quarterly maximum ─────────────────────────────────────────────────────────
quarterly_max = ndvi_dask.resample(time="QS").max()
# 12 layers for 3 years (4 quarters × 3 years)

# ── Annual mean ───────────────────────────────────────────────────────────────
annual_mean = ndvi_dask.resample(time="AS").mean()
# 3 layers

Common time frequency aliases: "D" (day), "W" (week), "MS" (month start), "ME" (month end), "QS" (quarter start), "AS" (year start).

# Loading Real Data: xr.open_mfdataset

For real workflows, data comes from multiple GeoTIFF files — one per scene or per date. open_mfdataset opens them all lazily as a single concatenated array:

import xarray as xr
import glob

# List of GeoTIFFs, one per acquisition date
tif_files = sorted(glob.glob("landsat_ndvi/LC08_*_NDVI.tif"))

# Open all as a single dataset — lazy, no data loaded yet
ds = xr.open_mfdataset(
    tif_files,
    engine="rasterio",            # requires rioxarray
    concat_dim="time",
    combine="nested",
    parallel=True,                # open file headers in parallel
    preprocess=lambda ds: ds.assign_coords(
        time=pd.Timestamp(ds.encoding['source'].split("_")[3])
    )
)

ndvi = ds["band_1"].rename({"band_1": "ndvi"})
ndvi_chunked = ndvi.chunk({"time": 30, "y": 512, "x": 512})

With parallel=True, file headers are opened concurrently. Actual data only loads when you call .compute().

# Masking Cloudy Pixels Before Aggregation

Real satellite data requires quality masking before temporal aggregation. Landsat provides a QA_PIXEL band:

import xarray as xr
import numpy as np

# Load QA_PIXEL alongside NDVI
qa = ds["QA_PIXEL"]

# Bit 3: cloud shadow. Bit 4: snow. Bit 5: cloud
cloud_mask = (qa & 0b00111000).astype(bool)   # True = cloudy/shadowed/snow

# Mask NDVI: set cloudy pixels to NaN
ndvi_masked = ndvi_chunked.where(~cloud_mask)

# Now monthly mean ignores NaN (masked) values
monthly_mean_masked = ndvi_masked.resample(time="1MS").mean(skipna=True)

skipna=True (the default) means the aggregation computes the mean over only the unmasked observations — equivalent to cloud-free compositing.

# Dask Scheduler Choices

The default Dask scheduler for XArray is the threaded scheduler, which works well for I/O-bound and NumPy operations:

# Default: threaded scheduler (good for most raster work)
result = monthly_mean.compute()

# For CPU-intensive operations on large arrays, try the distributed scheduler:
from dask.distributed import Client

client = Client(n_workers=4, threads_per_worker=2, memory_limit="8GB")
result = monthly_mean.compute()
client.close()

The distributed scheduler adds overhead for small computations but can be significantly faster for large ones by enabling more precise memory management and task graph visualisation.

# Percentile and Custom Aggregations

Not all aggregations are built-in. Custom reductions work with .reduce():

# 90th percentile NDVI per pixel per month (vegetation peak estimate)
def percentile_90(x, axis):
    return np.nanpercentile(x, 90, axis=axis)

monthly_p90 = ndvi_masked.resample(time="1MS").reduce(percentile_90)

For complex aggregations, apply_ufunc provides full control while preserving lazy evaluation:

import numpy as np
import xarray as xr

def pixel_trend_slope(ndvi_stack):
    """Compute linear trend slope per pixel using least squares."""
    # ndvi_stack: (time,) float array
    t = np.arange(len(ndvi_stack), dtype=float)
    valid = ~np.isnan(ndvi_stack)
    if valid.sum() < 2:
        return np.nan
    slope, _ = np.polyfit(t[valid], ndvi_stack[valid], 1)
    return slope

# Apply pixel-wise, preserving spatial dims
trend_map = xr.apply_ufunc(
    pixel_trend_slope,
    annual_mean,
    input_core_dims=[["time"]],
    vectorize=True,
    dask="parallelized",
    output_dtypes=[float]
)
result_trend = trend_map.compute()

# Writing Results

# Write to Zarr (best for continued XArray work)
monthly_mean_masked.to_zarr("monthly_ndvi.zarr", mode="w")

# Write to multi-band GeoTIFF via rioxarray
import rioxarray
monthly_mean_masked.rio.write_crs("EPSG:32755").to_raster("monthly_ndvi.tif")

# Write individual months as separate TIFFs
for i, month_data in enumerate(monthly_mean_masked):
    month_label = f"{month_data.time.dt.strftime('%Y-%m').item()}"
    month_data.rio.to_raster(f"monthly_ndvi_{month_label}.tif")

# Memory Management Rule of Thumb

Chunk size should be roughly 100–500 MB in memory after decompression. For float32 data:

A (30, 512, 512) chunk = 30 × 512 × 512 × 4 bytes ≈ 31 MB ✓
A (1, 2048, 2048) chunk = 1 × 2048 × 2048 × 4 bytes ≈ 16 MB ✓
A (1095, 2000, 2000) chunk (unchunked time dim) ≈ 8.8 GB ✗ — OOM

The time dimension should always be chunked when doing temporal aggregations, unless your entire time series fits in memory.