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.

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.

Array memory for a single 2,048 by 2,048 pixel windowBar chart of array memory for one 2,048 by 2,048 pixel window: 8 mebibytes for a band read as uint16, 16 mebibytes for the same window cast to float32, 48 mebibytes for the kernel's three float32 arrays together, and 96 mebibytes if the same kernel is run in float64.What one 2,048 px window actually costs in RAMB04 window as read, uint168 MiBSame window cast to float3216 MiBKernel working set, float3248 MiBSame kernel in float6496 MiB0MiB for one 2,048 x 2,048 window2,048 x 2,048 x 4 bytes is 16 MiB per float32 array, and the kernel holds three of them — red, NIR and the result. float64 doubles every rowto resolve a range that only ever spans -1 to 1.
The whole kernel working set is 48 MiB — roughly a thirty-fifth of the 1,769 MB the handler is sized at. That memory tier is bought for the vCPU that comes with it, not because the arrays need the room.

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.

The guarded order of operations inside the NDVI kernelFour numbered steps: cast both bands to float32 before any arithmetic, guard the zero denominator inside an errstate block, propagate the source nodata of either band into the result, and clamp the surviving values into the valid minus one to one range.Four steps, each guarding one specific failure1Cast both bands to float32uint16 subtraction underflows to a huge positive value and uint16 division truncates every fraction to0 or 1.astype('float32')2Guard the denominatornir + red is zero over nodata and deep shadow, so np.where substitutes the sentinel for the wholebranchnp.errstate(...)3Propagate source nodataa pixel that was nodata in either band cannot yield a meaningful index whatever the denominator says-9999.04Clamp the survivorsnumerical noise can leave a valid pixel a hair outside the range NDVI is defined onnp.clip(-1, 1)
Reversing steps two and three is harmless only while both write the same sentinel. Reversing steps one and two is never harmless: the arithmetic has already been done in uint16 by the time the guard runs.
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.

Kernel behaviour on a vegetation pixel, a zero denominator and a nodata pixelThree panels, one per branch of the kernel. A vegetation pixel with red 1000 and NIR 3000 yields 0.5. A pixel where both bands are zero has a zero denominator and yields the minus 9999 sentinel. A pixel whose value equals the band's declared nodata is masked to the same sentinel by the nodata propagation step.The three pixels the unit test has to exerciseVegetation pixelred = 1000, NIR = 3000denominator = 4000(3000 - 1000) / 4000 = 0.5asserted to within 1e-6Zero denominatorred = 0, NIR = 0denominator = 0np.where picks the sentinel branchresult -9999.0Source nodatathe band declares nodata = 0masked whatever the denominator isthe mask is applied after thedivisionresult -9999.0All three assertions must hold together: an exact 0.5 proves the float32 division is real-valued, and both sentinels provethe guards fire.
A kernel that passes only the first panel looks correct on farmland and writes inf across every water body and scene edge — which is exactly the class of bug that survives eyeballing a preview.
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