Computing an NDVI Difference Between Two Dates

To difference NDVI between two dates, align the grids, mask both scenes, compute the index on each, and subtract:

import rioxarray
import xarray as xr
from rasterio.enums import Resampling

before = rioxarray.open_rasterio("ndvi_2022-07-12.tif", masked=True).squeeze(drop=True)
after = rioxarray.open_rasterio("ndvi_2023-07-15.tif", masked=True).squeeze(drop=True)
after = after.rio.reproject_match(before, resampling=Resampling.bilinear)

d_ndvi = after - before          # negative = loss, positive = gain
d_ndvi.rio.to_raster("dndvi_2022_2023.tif", driver="COG", compress="DEFLATE")

This is the minimal correct version of the workflow in Change Detection and Differencing Workflows.

Why This Arises in Remote Sensing Workflows

An NDVI difference is the workhorse of vegetation monitoring: it detects harvest, clearance, dieback, irrigation change, and the greening that follows a wet season. It is popular because it is interpretable — everyone understands “the vegetation index went down” — and because the index cancels much of the illumination difference that would otherwise dominate a raw band difference.

What it does not cancel is the difference in what was observed. Two dates rarely have the same clear-sky footprint, and a pixel that was cloudy on one date has no defined change value at all. Treating those pixels as zero change understates loss; treating them as valid data lets cloud edges masquerade as change. The only defensible option is to carry an explicit validity mask through the whole calculation, which is most of what the code below does.

The assessable area is an intersection, not a scene The before date is clear except for a cloud in the north-west; the after date is clear except for a band of cloud across the south. Change can only be computed where both were clear, which here is about sixty percent of the tile. The remaining area is unassessed, not unchanged. before: clear-sky mask after: clear-sky mask assessable = both cloud cloud 61% of the tile white = unassessed, not "no change" Report the assessable fraction with the product: a map covering 61% of a tile means something different from one covering 98%.

Environment & Setup

Package Version Why
rioxarray ≥0.15 Masked reads, reproject_match, COG writing
xarray ≥2023.6 Labelled arithmetic
rasterio ≥1.3.0 Resampling enums and I/O
numpy ≥1.23 Statistics for the verification step 
pip install "rioxarray>=0.15" "xarray>=2023.6" "rasterio>=1.3.0"

Complete Working Example

This function takes band paths and quality layers for both dates and returns a signed difference plus the metadata that makes it interpretable.

Statistics that qualify a difference layer The assessable fraction says how much of the tile could be compared at all. The median over stable ground exposes a radiometric offset between the dates. Sigma sets the threshold scale. Together they decide whether the difference layer is worth thresholding. Read these three before thresholding anything statistic healthy value what a bad value means assessable fraction > 0.6 of the tile one date was mostly cloudy median over stable ground within ±0.02 a baseline or atmospheric offset sigma over stable ground 0.02 – 0.05 above this, registration is suspect valid pixels in both dates similar count per date one date contributes almost nothing A difference layer whose median is far from zero maps the offset, not the change. 
import numpy as np
import rioxarray  # noqa: F401
import xarray as xr
from rasterio.enums import Resampling

CLOUDY_SCL = (0, 1, 3, 8, 9, 10, 11)     # no-data, defective, shadow, cloud ×3, snow


def ndvi_from(red_path: str, nir_path: str, scl_path: str, reference: xr.DataArray | None = None):
    """Masked NDVI for one date, optionally matched to a reference grid."""
    red = rioxarray.open_rasterio(red_path, masked=True).squeeze(drop=True).astype("float32")
    nir = rioxarray.open_rasterio(nir_path, masked=True).squeeze(drop=True).astype("float32")
    scl = rioxarray.open_rasterio(scl_path).squeeze(drop=True)

    if reference is not None:
        red = red.rio.reproject_match(reference, resampling=Resampling.bilinear)
        nir = nir.rio.reproject_match(reference, resampling=Resampling.bilinear)
    # The class layer is 20 m and categorical: nearest, always
    scl = scl.rio.reproject_match(red, resampling=Resampling.nearest)

    clear = ~scl.isin(list(CLOUDY_SCL))
    red, nir = red.where(clear), nir.where(clear)

    denom = nir + red
    ndvi = xr.where(denom != 0, (nir - red) / denom, np.nan).astype("float32")
    return ndvi, clear


def ndvi_difference(before: dict, after: dict) -> tuple[xr.DataArray, dict]:
    """Signed NDVI difference (after − before) over pixels clear on both dates."""
    ndvi_b, clear_b = ndvi_from(before["red"], before["nir"], before["scl"])
    ndvi_a, clear_a = ndvi_from(after["red"], after["nir"], after["scl"], reference=ndvi_b)
    clear_a = clear_a.rio.reproject_match(ndvi_b, resampling=Resampling.nearest)

    assessable = clear_b & clear_a
    diff = (ndvi_a - ndvi_b).where(assessable).astype("float32")

    stats = {
        "assessable_fraction": float(assessable.mean()),
        "median_difference": float(diff.median(skipna=True)),
        "sigma": float(diff.std(skipna=True)),
        "date_before": before["date"],
        "date_after": after["date"],
        "convention": "after_minus_before: negative = vegetation loss",
    }
    diff.attrs.update(stats)
    return diff, stats


if __name__ == "__main__":
    diff, stats = ndvi_difference(
        {"red": "B04_2022.tif", "nir": "B08_2022.tif", "scl": "SCL_2022.tif", "date": "2022-07-12"},
        {"red": "B04_2023.tif", "nir": "B08_2023.tif", "scl": "SCL_2023.tif", "date": "2023-07-15"},
    )
    print(stats)
    diff.rio.write_nodata(np.nan, inplace=True)
    diff.rio.to_raster("dndvi.tif", driver="COG", compress="DEFLATE", blocksize=512)

The assessable_fraction and median_difference in the returned statistics are not decoration. A median far from zero means the two dates are not radiometrically comparable — different baselines, different atmosphere, or genuinely different seasons — and thresholding that difference will produce a map of the offset rather than a map of change.

Variant Patterns

1. Composite pairs instead of single scenes

Where cloud makes single dates unreliable, difference two composites built from short windows around each anniversary.

import xarray as xr

before_stack = xr.concat([ndvi_from(**p)[0] for p in before_scenes], dim="time")
after_stack = xr.concat([ndvi_from(**p, reference=before_stack[0])[0] for p in after_scenes], dim="time")

before_c = before_stack.median("time", skipna=True)
after_c = after_stack.median("time", skipna=True)
count_b = before_stack.count("time")           # carry the observation count
count_a = after_stack.count("time")

diff = (after_c - before_c).where((count_b >= 2) & (count_a >= 2))

Requiring at least two contributing observations per side is what stops a composite pixel resting on a single cloudy-day reading — the reliability argument made in Creating Monthly NDVI Composites with xarray resample.

2. Relative rather than absolute difference

For low-biomass areas, an absolute NDVI drop of 0.1 means something very different than it does in dense canopy. A relative difference normalises for the starting level.

rel = (after_c - before_c) / before_c.where(before_c > 0.15)     # guard sparse vegetation

The guard matters: dividing by a near-zero baseline produces enormous relative changes over bare ground and deserts, which then dominate the map.

3. Difference chunked, for large areas

How many arrays are alive at once A naive two-date difference holds red, near-infrared and the class layer for both dates plus two index arrays and the result — nine arrays at once. Chunking keeps the same nine arrays but at chunk size, so peak memory is set by the chunk rather than by the scene. Arrays alive during a two-date NDVI difference whole-scene 9 × 480 MB ≈ 4.3 GB peak — one worker, one tile chunked (1024²) 9 × 4 MB ≈ 36 MB peak per thread — scales to a cluster Same nine arrays either way; the only variable is how much of each is resident. Open with chunks={"x": 1024, "y": 1024} and the arithmetic above is unchanged. 
before = rioxarray.open_rasterio("B04_2022.tif", masked=True, chunks={"x": 1024, "y": 1024})

Passing chunks is the only change needed; the arithmetic stays identical and the computation becomes lazy, as covered in Tuning Dask Chunk Sizes for Raster Cubes.

Sanity Checks Before Thresholding

Three cheap checks catch the failures that a threshold would otherwise dress up as findings.

The median over the whole assessable area should be near zero for an anniversary pair. A median of −0.08 means a systematic difference — often a processing baseline change — and thresholding it will flag a large fraction of the tile as loss.

Self-differencing should return zero. Run the same pipeline with the same date as both inputs; anything other than an all-zero result reveals a bug in the alignment or masking path, and it takes seconds to run.

Stable water should not change. Water has near-constant, strongly negative NDVI, so a lake that shows change is a registration, masking or scaling problem rather than a hydrological one — the same use of a known target described in Spectral Index Calculation Pipelines.

Common Errors

The difference is dominated by a constant offset

Different processing baselines between the dates, so the reflectance scaling differs. Check the scale and offset on both, and see Handling nodata and Scale Factors in Band Math.

ValueError: cannot align objects with join='exact'

The two dates are on different grids and reproject_match was skipped. Align the second date to the first before any arithmetic.

The output has no CRS

Arithmetic dropped the accessor state. Re-attach with diff.rio.write_crs(before.rio.crs, inplace=True) before writing.

Frequently Asked Questions

Q: Should I difference NDVI or compute NDVI of the differenced bands? Difference the indices. NDVI of differenced bands is not a meaningful quantity, because the normalisation no longer refers to a single acquisition’s reflectance.

Q: Why is my difference layer mostly NaN? The valid footprint is the intersection of both dates’ clear-sky masks. Two 40-percent-cloudy scenes can leave well under half the tile assessable, which is expected and should be reported rather than hidden.

Q: What sign convention should the output use? After minus before, so negative means vegetation loss and positive means gain. Whatever you choose, write it into the file’s metadata — sign confusion is the most common misreading of a change product.