Temporal Aggregation of Vegetation Indices: Production-Ready Python Workflows for Crop Monitoring

Q: Why does my seasonal mean contain NaN over pixels that clearly have data?

The most common cause is omitting skipna=True. By default, xarray's .mean() propagates NaN if any observation in the stack is masked, so a single cloudy date invalidates every pixel it touches. Pass skipna=True to all reduction methods and also check that your cloud mask aligns pixel-perfectly with the reflectance rasters — a one-pixel shift can corrupt valid data.

Q: My xarray stack runs out of memory even though I set chunks. What is wrong?

Calling .compute() or .values() anywhere before the final export materialises the entire dataset. Check that you are not accidentally triggering computation inside a loop or a print statement. Also verify that your chunks parameter creates tiles small enough — for float32 data, a 1024×1024 tile per time step is roughly 4 MB; multiply by the number of dates to estimate peak memory per worker.

Temporal aggregation of vegetation indices solves a core precision-agriculture problem: a single drone or satellite overpass captures one moment, but crop health decisions require a view of the whole growing season. This page shows you how to stack date-stamped NDVI, NDRE, and EVI rasters into robust seasonal statistics — growing-season means, peak-canopy maximums, and 90th-percentile stress thresholds — using xarray, rioxarray, and dask. The output is a set of analysis-ready GeoTIFF layers that feed directly into field zoning, yield forecasting, and anomaly detection.

This workflow sits in the Drone Imagery Processing & Vegetation Index Workflows pipeline, downstream of index calculation and upstream of threshold-based zone mapping.

Pipeline overview

GeoTIFF, float32 per-pixel NaN lazy / out-of-core skipna=True DEFLATE, EPSG:32754

Prerequisites

Requirement	Detail
Python	3.9 or later
`rasterio`	≥ 1.3
`xarray`	≥ 2023.1
`rioxarray`	≥ 0.15
`numpy`	≥ 1.24
`dask`	≥ 2023.1
Input rasters	Single-band float32 GeoTIFF, consistent CRS (state the EPSG explicitly — never assume), identical spatial resolution and pixel grid
Nodata convention	`NaN` preferred; if sensor outputs `-9999` or `65535`, convert during load
Sensor assumptions	Applies to any calibrated multispectral source — MicaSense RedEdge-MX, DJI P4 Multispectral, Parrot Sequoia, or Sentinel-2 exports — provided reflectance conversion is already complete
Upstream dependency	Raw digital numbers must already be converted to reflectance and index values produced via band math & raster algebra in Python before stacking

Install with:

BASH

pip install "rasterio>=1.3" "xarray>=2023.1" rioxarray "numpy>=1.24" "dask[array]>=2023.1"

1. Concept & algorithm

Single-date vegetation index snapshots are noisy. Phenological progression, solar angle variation, atmospheric haze, sensor temperature drift, and patchy cloud cover all shift pixel values between flights. Aggregating a time series into stable statistics — seasonal maximum canopy vigour, growing-season mean biomass, or 90th-percentile stress threshold — dampens this noise and reveals genuine crop-health signals.

Why these statistics matter agronomically:

Seasonal maximum corresponds to peak canopy closure, a robust proxy for above-ground biomass and the basis for yield-potential zone maps. MicaSense RedEdge-MX flights over wheat typically show NDRE seasonal maxima above 0.45 in high-yield zones and below 0.30 in stressed patches.
Growing-season mean smooths phenological asymmetry between early- and late-maturing paddocks, making it suitable for input into management-zone classification where you need comparable field-wide signals.
90th percentile (p90) captures peak stress episodes without being distorted by outlier cloudy dates. It is the recommended statistic for input to variable-rate nitrogen prescription models where over-application risk is high.

The algorithm is straightforward: mask invalid observations to NaN, stack along a time dimension, then apply a reduction function with skipna=True. The computational challenge is doing this without loading the entire time series into RAM — xarray with dask solves this through lazy, chunked evaluation.

2. Step-by-step implementation

Step 1 — Standardise CRS and pixel grid

Temporal aggregation fails when rasters have mismatched projections, pixel offsets, or inconsistent bit depths. Before stacking, reproject every raster to the same projected CRS. Use an explicit EPSG code — never let rasterio infer the zone.

PYTHON

import rasterio
from rasterio.warp import calculate_default_transform, reproject, Resampling
from pathlib import Path

TARGET_CRS = "EPSG:32754"   # UTM zone 54 S — adjust to your field location
TARGET_RES = (2.0, 2.0)     # 2 m GSD matching MicaSense RedEdge-MX at 60 m AGL

def align_to_target(src_path: Path, dst_path: Path) -> None:
    """Reproject a VI raster to TARGET_CRS and TARGET_RES."""
    with rasterio.open(src_path) as src:
        transform, width, height = calculate_default_transform(
            src.crs, TARGET_CRS, src.width, src.height,
            *src.bounds, resolution=TARGET_RES
        )
        profile = src.profile.copy()
        profile.update(
            crs=TARGET_CRS,
            transform=transform,
            width=width,
            height=height,
            dtype="float32",
            nodata=float("nan"),
            compress="deflate",
        )
        with rasterio.open(dst_path, "w", **profile) as dst:
            reproject(
                source=rasterio.band(src, 1),
                destination=rasterio.band(dst, 1),
                src_transform=src.transform,
                src_crs=src.crs,
                dst_transform=transform,
                dst_crs=TARGET_CRS,
                resampling=Resampling.bilinear,
            )

aligned_dir = Path("data/vi_aligned")
aligned_dir.mkdir(parents=True, exist_ok=True)
for p in sorted(Path("data/vi_raw").glob("*.tif")):
    align_to_target(p, aligned_dir / p.name)
    print(f"Aligned: {p.name}")

Verify alignment before proceeding:

PYTHON

files = sorted(aligned_dir.glob("*.tif"))
with rasterio.open(files[0]) as ref:
    ref_shape, ref_crs, ref_transform = ref.shape, ref.crs, ref.transform
for f in files[1:]:
    with rasterio.open(f) as src:
        assert src.shape == ref_shape, f"Shape mismatch: {f.name}"
        assert src.crs == ref_crs, f"CRS mismatch: {f.name}"
        assert src.transform == ref_transform, f"Transform mismatch: {f.name}"
print("All rasters aligned ✓")

Step 2 — Apply cloud and quality masks

Optical imagery requires per-pixel quality filtering before any temporal statistic is computed. Apply your cloud and shadow masks — following the approach in cloud masking for agricultural imagery — and replace invalid pixels with float("nan"). This ensures skipna=True in the reduction step ignores them cleanly.

PYTHON

import numpy as np

def apply_qa_mask(vi_path: Path, mask_path: Path, out_path: Path) -> None:
    """
    Replace masked pixels with NaN.
    mask_path: single-band uint8 raster — 1 = valid, 0 = cloud/shadow/invalid.
    """
    with rasterio.open(vi_path) as vi_src, rasterio.open(mask_path) as mask_src:
        assert vi_src.shape == mask_src.shape, "VI / mask shape mismatch"
        profile = vi_src.profile.copy()
        profile.update(dtype="float32", nodata=float("nan"))
        vi = vi_src.read(1).astype("float32")
        mask = mask_src.read(1)           # 1 = valid, 0 = masked
        vi[mask == 0] = float("nan")
        with rasterio.open(out_path, "w", **profile) as dst:
            dst.write(vi, 1)

masked_dir = Path("data/vi_masked")
masked_dir.mkdir(parents=True, exist_ok=True)
for vi_path in sorted(aligned_dir.glob("*.tif")):
    mask_path = Path("data/masks") / vi_path.name
    apply_qa_mask(vi_path, mask_path, masked_dir / vi_path.name)

Always generate an observation-count layer alongside your aggregated output. Pixels with fewer than three valid observations across the growing season should be flagged or excluded from downstream modelling to avoid spurious agronomic signals.

Step 3 — Build a lazy xarray time-series stack

xarray with dask chunks lets you represent a multi-year archive without loading it into RAM. Parse acquisition dates from filenames (ISO-8601 format is strongly recommended: ndvi_2024-05-18.tif) and attach them as a time coordinate.

PYTHON

import xarray as xr
import rioxarray          # registers .rio accessor on DataArray/Dataset
import pandas as pd
from pathlib import Path

masked_dir = Path("data/vi_masked")

def date_from_filename(path: Path) -> pd.Timestamp:
    """Extract ISO-8601 date from filename stem, e.g. 'ndvi_2024-05-18'."""
    date_str = path.stem.split("_")[-1]   # '2024-05-18'
    return pd.Timestamp(date_str)

vi_paths = sorted(masked_dir.glob("*.tif"))
dates = [date_from_filename(p) for p in vi_paths]

# Open lazily — engine="rasterio" preserves CRS and affine metadata
arrays = [
    xr.open_dataarray(p, engine="rasterio", chunks={"x": 1024, "y": 1024})
    .squeeze("band", drop=True)         # drop the band dim (single-band rasters)
    .assign_coords(time=d)
    for p, d in zip(vi_paths, dates)
]

vi_stack = xr.concat(arrays, dim="time")
vi_stack = vi_stack.sortby("time")

# Sanity-check without triggering full compute
assert vi_stack.dims == ("time", "y", "x"), f"Unexpected dims: {vi_stack.dims}"
print(f"Stack shape: {vi_stack.sizes}  |  CRS: {vi_stack.rio.crs}")

The chunks={"x": 1024, "y": 1024} parameter delegates memory management to dask. At float32, each 1024×1024 tile is 4 MB — multiply by the number of dates to estimate peak memory per dask worker. Never chunk along the time dimension for reduction operations; chunking along x/y keeps each worker’s time series contiguous.

Step 4 — Apply aggregation logic

Slice to your growing season then apply reduction functions along time. Always pass skipna=True explicitly — the default varies by operation and numpy version.

PYTHON

# Clip to growing season — adjust for your hemisphere and crop calendar
SEASON_START = "2024-04-01"
SEASON_END   = "2024-09-30"

season = vi_stack.sel(time=slice(SEASON_START, SEASON_END))

# Count valid (non-NaN) observations per pixel
obs_count = season.count(dim="time")                                # int64

# Core seasonal statistics — all lazy until .compute() or export
season_max  = season.max(dim="time", skipna=True)                   # peak canopy
season_mean = season.mean(dim="time", skipna=True)                  # biomass proxy
season_p90  = season.quantile(0.90, dim="time", skipna=True)        # stress threshold

# Flag pixels with insufficient observations
MIN_OBS = 3
season_max  = season_max.where(obs_count >= MIN_OBS)
season_mean = season_mean.where(obs_count >= MIN_OBS)
season_p90  = season_p90.where(obs_count >= MIN_OBS)

# Validate value ranges — NDVI should sit in [-0.2, 1.0] for vegetated areas
assert float(season_max.min(skipna=True)) >= -0.5, "Unexpected low minimum — check masking"
assert float(season_max.max(skipna=True)) <= 1.05, "Unexpected high maximum — check calibration"
print("Value range check passed ✓")

For non-standard reductions — rolling phenological windows or Savitzky-Golay smoothing before reduction — use xarray’s .rolling() or .groupby("time.month") methods. Chain operations without calling .compute() until export.

Step 5 — Export to Cloud Optimised GeoTIFF

Write aggregated rasters with rioxarray’s .rio.to_raster() to guarantee CRS and affine metadata are preserved. Use DEFLATE compression (zlevel=6) and float32 to minimise storage without losing precision.

PYTHON

from pathlib import Path

output_dir = Path("output/seasonal_stats")
output_dir.mkdir(parents=True, exist_ok=True)

export_kwargs = {
    "compress": "deflate",
    "zlevel": 6,
    "tiled": True,
    "blockxsize": 512,
    "blockysize": 512,
    "dtype": "float32",
}

# Trigger dask compute and write
season_max.rio.to_raster(output_dir / "ndvi_season_max.tif",  **export_kwargs)
season_mean.rio.to_raster(output_dir / "ndvi_season_mean.tif", **export_kwargs)
season_p90.rio.to_raster(output_dir / "ndvi_season_p90.tif",  **export_kwargs)
obs_count.astype("int16").rio.to_raster(output_dir / "ndvi_obs_count.tif",
                                         compress="deflate", dtype="int16")

print("Export complete.")

The tiled=True option with 512×512 blocks produces a Cloud Optimised GeoTIFF layout that supports efficient windowed reads by downstream tools and web tile servers.

3. Key parameters and tuning

Parameter	Type	Default	Agronomic Effect
`chunks` (`x`, `y`)	`int`	`1024`	Larger chunks reduce dask task-graph overhead but increase per-worker RAM. For DJI P4 Multispectral 3 cm imagery, 512 is safer.
`SEASON_START` / `SEASON_END`	`str` (ISO-8601)	Crop-specific	Determines which flights contribute. Extending past senescence inflates the mean with low post-harvest values.
`MIN_OBS`	`int`	`3`	Minimum cloud-free observations for a pixel to be included. Increase to 5 for high-cloud regions; decrease to 2 only for short-season crops with infrequent flights.
`quantile` `q`	`float`	`0.90`	p90 is the recommended stress-threshold input for nitrogen VRA. p50 (median) is more robust to outliers for zone classification.
`resampling`	`Resampling` enum	`bilinear`	Use `nearest` for categorical quality masks; `bilinear` for continuous reflectance rasters to avoid spectral ringing at edges.
`zlevel`	`int` (1–9)	`6`	Higher levels compress more but write slower. Level 6 is a good balance for archival COG storage.

4. Edge cases and failure modes

NaN propagation across the grid. If your seasonal mean or max is NaN everywhere despite the stack containing data, you omitted skipna=True. This is the most frequent bug. Explicitly pass it to every .mean(), .max(), and .quantile() call.

Misaligned CRS. Stacking rasters with different CRS silently produces nonsense spatial outputs because xarray.concat aligns on coordinate labels, not on geographic extent. Always run the alignment step (Step 1) before building the stack; assert CRS equality after opening.

Memory exhaustion during .quantile(). Quantile computation is more memory-intensive than mean or max because dask must sort values along the time dimension. If you hit OOM errors, reduce chunk size from 1024 to 512, or use method="nearest" instead of the default linear interpolation — this halves the intermediate array size.

Timestamp parsing failures. Mixed datetime formats (2024-05-18 vs 20240518 vs 2024_05_18) across a multi-year archive cause xr.concat to fail or, worse, silently mis-sort dates. Enforce ISO-8601 filenames at ingest or normalise with pd.to_datetime(date_str, infer_datetime_format=False).

UTM zone boundary crossings. Fields near UTM zone meridians (e.g., spanning zones 54 and 55) will have rasters projected into different EPSG codes by default. Always specify the same TARGET_CRS for every flight in the project; use the zone covering the majority of the field.

Cloud shadows mis-identified as crop stress. Shadow-masked NaN pixels cluster spatially and create artefacts that look like localised stress zones in the seasonal statistics. Inspect the obs_count raster before interpreting the mean — a low count in a ring or wedge pattern indicates a systematic shadow-masking issue rather than genuine crop variability.

5. Verification and output validation

After export, validate outputs without loading the full raster into memory:

PYTHON

import rasterio
import numpy as np

output_path = "output/seasonal_stats/ndvi_season_max.tif"

with rasterio.open(output_path) as src:
    # Metadata checks
    assert src.crs.to_epsg() == 32754, f"CRS mismatch: {src.crs}"
    assert src.dtypes[0] == "float32", f"Unexpected dtype: {src.dtypes[0]}"
    assert src.nodata is None or np.isnan(src.nodata), "nodata should be NaN"

    # Sample the centre of the raster for a sanity check
    cx, cy = src.width // 2, src.height // 2
    window = rasterio.windows.Window(cx - 64, cy - 64, 128, 128)
    sample = src.read(1, window=window)
    valid = sample[~np.isnan(sample)]

    assert len(valid) > 0, "Centre window is entirely NaN — check masking"
    assert valid.min() > -0.5, f"Min NDVI too low: {valid.min():.3f}"
    assert valid.max() < 1.05, f"Max NDVI too high: {valid.max():.3f}"
    print(f"Centre window: min={valid.min():.3f}  mean={valid.mean():.3f}  max={valid.max():.3f}")

print("Output validation passed ✓")

Run gdalinfo output/seasonal_stats/ndvi_season_max.tif to confirm the DEFLATE compression, 512×512 tiling, and embedded CRS. Cross-check the seasonal mean over known calibration targets (e.g., a spectrally stable panel or a harvested headland) — if the mean deviates more than 0.05 NDVI units from your hand-measured reference, revisit the radiometric calibration step.

Inspect the obs_count.tif histogram: for a 6-month growing season with fortnightly flights you expect 6–12 observations per pixel in cloud-free regions. Pixels consistently below 3 indicate a systematic cloud or shadow masking problem, not a data gap.

6. Integration with the broader pipeline

The outputs of this workflow — seasonal_max.tif, seasonal_mean.tif, seasonal_p90.tif, and obs_count.tif — are the direct inputs to threshold-based zone classification. Feed seasonal_mean.tif into the threshold mapping for crop health step to classify fields into low, medium, and high-vigour zones. The obs_count.tif layer should be used to mask out data-sparse pixels before zone boundaries are finalised, otherwise shapefile outputs will contain unreliable zones.

Upstream dependency: Raw VI rasters must already have reflectance conversion and band math applied — see band math & raster algebra in Python — and cloud contamination removed as described in cloud masking for agricultural imagery before passing to this aggregation step.

Storage format recommendation: Export seasonal statistics as Cloud Optimised GeoTIFF for direct consumption by web-served field dashboards and by Python tools using windowed reads. For multi-year archival, Zarr with chunked DEFLATE compression enables efficient partial reads when querying a single season from a ten-year archive without reading the whole dataset.

CRS note: All outputs must share the same projected CRS (state the EPSG code in every rio.write_crs() call). Downstream tools that join VI statistics with vector field boundaries — using the approach in field boundary extraction with GeoPandas — require exact CRS alignment for spatial joins to resolve correctly. Validate with understanding CRS in precision agriculture if you encounter join failures.

Frequently asked questions

Why does my seasonal mean contain NaN over pixels that clearly have data? The most common cause is omitting skipna=True. By default, xarray’s .mean() propagates NaN if any observation in the stack is masked, so a single cloudy date invalidates every pixel it touches. Pass skipna=True to all reduction methods and also check that your cloud mask aligns pixel-perfectly with the reflectance rasters — a one-pixel shift can corrupt valid data.

My xarray stack runs out of memory even though I set chunks. What is wrong? Calling .compute() or .values() anywhere before the final export materialises the entire dataset. Check that you are not accidentally triggering computation inside a loop or a print statement. Also verify that your chunks parameter creates tiles small enough — for float32 data, a 1024×1024 tile per time step is roughly 4 MB; multiply by the number of dates to estimate peak memory per worker.

How many valid observations per pixel do I need for a reliable seasonal statistic? As a practical minimum, three to five cloud-free observations spread across the growing season provide a defensible seasonal mean for most broadacre crops. Fewer than three observations should be flagged in your output obs_count raster and excluded from downstream zone classification or yield modelling.

Drone Imagery Processing & Vegetation Index Workflows — parent overview covering the full UAV-to-prescription pipeline
Band Math & Raster Algebra in Python — upstream step that produces the per-date VI rasters this workflow aggregates
Cloud Masking for Agricultural Imagery — produces the QA masks applied in Step 2
Threshold Mapping for Crop Health — consumes the seasonal statistics produced here to generate zone polygons
Field Boundary Extraction with GeoPandas — for clipping aggregated rasters to paddock geometries before zone export

Temporal Aggregation of Vegetation Indices: Production-Ready Python Workflows for Crop Monitoring

Related on this site