Computing NDVI per Tile with Rasterio and NumPy

Computing NDVI for one tile is a windowed read of the red (B04) and NIR (B08) bands, a cast to float32, and ndvi = np.where(denom == 0, -9999, (nir - red) / denom) inside a np.errstate(divide='ignore', invalid='ignore') block. Cast before the arithmetic so uint16 bands cannot underflow or truncate, mask every pixel where the denominator is zero or either source band is nodata, and write the result as a float32 COG with a -9999 nodata tag. This is the single computational kernel that each per-tile worker in the pipeline runs.

Per-tile NDVI computation kernel Five-stage kernel: window-read red B04 and NIR B08 arrays, cast both to float32, compute NDVI with masked division guarding denominator equal to zero, propagate source nodata into the mask, and write the float32 NDVI COG tile with a -9999 nodata value. Read Windows B04 red B08 NIR uint16 Cast float32 no overflow no truncation Masked Division (nir-red)/(nir+red) guard denom=0 Nodata Mask propagate source nodata Write COG float32 nodata -9999 NDVI kernel — window in, float32 NDVI COG tile out

Context

The serverless NDVI tiling pipeline fans out one worker per window of a Sentinel-2 scene. This page zooms into the single function each worker calls: the NDVI kernel. It is short, but every line guards a specific failure. Skipping the float32 cast silently produces garbage; skipping the denominator guard raises RuntimeWarning and writes inf or nan over nodata; skipping source-nodata propagation paints spurious vegetation over water. Getting this kernel exactly right is what separates a tile catalog you can trust from one full of edge artifacts.

The windowed read that feeds this kernel is the same access pattern optimized in chunked I/O for large satellite imagery — you fetch only the bytes for the window, not the whole band. Sizing the function so numpy runs the division without swapping is governed by the memory and CPU allocation model for raster workloads, and if a scene ever outgrows a single Lambda the orchestration falls back to process a 10 GB GeoTIFF with Step Functions and Lambda.

Prerequisites

Before deploying the handler, confirm the following:

  • Runtime: Python 3.11 on AWS Lambda, rasterio==1.4.3 and numpy==2.1.3 from a runtime-matched layer built as in stripping unnecessary Python packages from AWS Lambda Layers
  • Memory: 1,769 MB or higher so a 2048×2048 window (three float32 arrays ≈ 48 MB plus GDAL buffers) fits comfortably with room for one full vCPU
  • IAM: s3:GetObject on the source Sentinel-2 bucket and s3:PutObject on the NDVI output bucket
  • Environment variables (GDAL COG read tuning and library data paths):
    code 
    GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR
CPL_VSIL_CURL_ALLOWED_EXTENSIONS=.tif
VSI_CACHE=TRUE
GDAL_DATA=/opt/share/gdal
PROJ_LIB=/opt/share/proj
  • Input contract: the message payload carries red and nir COG URIs (single-band B04 and B08) and a window dict of col_off, row_off, width, height

Implementation

The handler reads both windows, computes NDVI with fully guarded arithmetic, and writes a float32 COG tile with the correct windowed geotransform. Every masking step is commented to make the failure it prevents explicit.

python 
# handler.py — per-tile NDVI kernel for AWS Lambda
import os

import numpy as np
import rasterio
from rasterio.windows import Window
from rasterio.windows import transform as window_transform

NDVI_NODATA = np.float32(-9999.0)  # sentinel outside the valid [-1, 1] range


def compute_ndvi(red_u16: np.ndarray, nir_u16: np.ndarray,
                 red_nodata, nir_nodata) -> np.ndarray:
    """Return a float32 NDVI array with nodata pixels set to -9999."""
    # 1) Cast BEFORE arithmetic. uint16 subtraction underflows to huge values
    #    and uint16 division truncates every fractional NDVI to 0 or 1.
    red = red_u16.astype("float32")
    nir = nir_u16.astype("float32")

    denom = nir + red  # zero over nodata and deep shadow

    # 2) Scope the divide-by-zero suppression to exactly this division, then
    #    substitute the sentinel wherever the denominator is zero.
    with np.errstate(divide="ignore", invalid="ignore"):
        ndvi = np.where(denom == 0.0, NDVI_NODATA, (nir - red) / denom)

    # 3) Propagate source nodata: a pixel that was nodata in EITHER band
    #    cannot yield a meaningful NDVI, regardless of the denominator.
    if red_nodata is not None:
        ndvi = np.where(red == np.float32(red_nodata), NDVI_NODATA, ndvi)
    if nir_nodata is not None:
        ndvi = np.where(nir == np.float32(nir_nodata), NDVI_NODATA, ndvi)

    # 4) Clamp any residual out-of-range values from numerical noise.
    valid = ndvi != NDVI_NODATA
    np.clip(ndvi, -1.0, 1.0, out=ndvi, where=valid)
    return ndvi.astype("float32")


def handler(event, context):
    red_uri = event["red"]
    nir_uri = event["nir"]
    w = event["window"]
    out_uri = event["out"]

    window = Window(w["col_off"], w["row_off"], w["width"], w["height"])

    # Windowed reads: only the bytes for this tile leave object storage.
    with rasterio.open(red_uri) as red_src:
        red = red_src.read(1, window=window)
        red_nodata = red_src.nodata
        base_transform = red_src.transform
        crs = red_src.crs
    with rasterio.open(nir_uri) as nir_src:
        nir = nir_src.read(1, window=window)
        nir_nodata = nir_src.nodata

    ndvi = compute_ndvi(red, nir, red_nodata, nir_nodata)

    profile = {
        "driver": "COG",
        "dtype": "float32",
        "count": 1,
        "width": w["width"],
        "height": w["height"],
        "crs": crs,
        "transform": window_transform(window, base_transform),
        "nodata": float(NDVI_NODATA),
        "compress": "LZW",
        "predictor": 3,          # floating-point predictor for better ratios
        "blocksize": 512,
    }
    with rasterio.open(out_uri, "w", **profile) as dst:
        dst.write(ndvi, 1)

    valid = ndvi[ndvi != NDVI_NODATA]
    return {
        "out": out_uri,
        "valid_pixels": int(valid.size),
        "ndvi_min": round(float(valid.min()), 4) if valid.size else None,
        "ndvi_max": round(float(valid.max()), 4) if valid.size else None,
    }

The order of operations matters: cast first, guard the denominator second, propagate source nodata third. Reversing the last two steps would let a divide-by-zero result leak through before the nodata mask overwrites it — harmless here because both write the same sentinel, but fragile if you later change the fill value.

Verification

Run the kernel against a synthetic pair that exercises every branch: normal vegetation, a zero-denominator pixel, and a source-nodata pixel.

python 
# verify_kernel.py — unit-test the NDVI kernel branches
import numpy as np
from handler import compute_ndvi, NDVI_NODATA

red = np.array([[1000, 0, 0]], dtype="uint16")   # veg, zero-denom, nodata
nir = np.array([[3000, 0, 0]], dtype="uint16")

out = compute_ndvi(red, nir, red_nodata=0, nir_nodata=0)
print(out)
# Vegetation pixel: (3000-1000)/(3000+1000) = 0.5
# Zero-denominator and nodata pixels: -9999.0
assert abs(out[0, 0] - 0.5) < 1e-6
assert out[0, 1] == NDVI_NODATA
assert out[0, 2] == NDVI_NODATA
print("NDVI kernel OK")

Expected output:

code 
[[ 5.0000e-01 -9.9990e+03 -9.9990e+03]]
NDVI kernel OK

The first value is exactly 0.5, confirming the float32 division is real-valued rather than integer-truncated. Both the zero-denominator pixel and the nodata pixel resolve to the -9999 sentinel.

Gotchas and Edge Cases

  • np.where still evaluates the divide-by-zero branch. np.where(denom == 0, sentinel, a / b) computes a / b for the whole array before selecting, so the division still emits a warning over zero denominators — that is exactly why the np.errstate block wraps it. The warning is suppressed, not the computation; the inf/nan results are simply discarded by the where.
  • Do not store NaN as the COG nodata. Some downstream integer or tile-server tooling mishandles NaN nodata. A concrete -9999.0 sentinel outside the [-1, 1] range is portable and unambiguous. Never use 0, which is a legitimate NDVI value over bare soil and water margins.
  • float32 is enough; float64 doubles memory for no benefit. NDVI ranges over [-1, 1] and float32 resolves it to roughly seven significant digits — far finer than the sensor’s radiometric precision. Using float64 doubles array memory and can push a large window over the function’s memory ceiling.
  • Reprojected or resampled bands break pixel alignment. B04 and B08 are both native 10 m, so their windows align exactly. If you substitute a 20 m band (for example B05 red-edge) you must resample it to the 10 m grid first, or the subtraction operates on misregistered pixels.

Frequently Asked Questions

Why cast the bands to float32 before computing NDVI?

Sentinel-2 bands are unsigned 16-bit integers. Subtracting them can underflow below zero and wrap to a huge positive number, and integer division truncates every fractional NDVI to 0 or 1. Casting both arrays to float32 first makes the subtraction signed and the division real-valued, and float32 is precise enough for NDVI’s [-1, 1] range while using half the memory of float64.

How do I suppress the numpy divide-by-zero warning without hiding real bugs?

Wrap only the division in a with np.errstate(divide='ignore', invalid='ignore') block, and immediately replace the invalid results using np.where on the denominator being zero. This scopes the suppression to the one operation that legitimately divides by zero over nodata pixels, so genuine numerical errors elsewhere in the handler still raise.

What nodata value should the NDVI output use?

Use a sentinel outside NDVI’s valid [-1, 1] range, conventionally -9999.0, and set it as the nodata tag on the output raster. Do not use NaN as the stored nodata for a COG that downstream integer tools may read, and never use 0, which is a valid NDVI value over bare soil and water edges.

Back to Serverless NDVI Tiling from Sentinel-2 Imagery