Serverless NDVI Tiling from Sentinel-2 Imagery

A serverless NDVI tiling pipeline reacts to a new Sentinel-2 scene by fanning out one function per tile: each worker window-reads the red (B04) and near-infrared (B08) bands directly from the source Cloud-Optimized GeoTIFF, computes NDVI = (NIR - red) / (NIR + red) with rasterio and numpy, and writes a catalogued NDVI COG tile. The governing rule is: never load a full 10980×10980 scene into one function — split it into fixed windows (for example 2048×2048) so each worker finishes in seconds inside a few hundred megabytes, well under the 15-minute AWS Lambda timeout and the 10 GB memory ceiling. This recipe sits within the broader discipline of event-driven geospatial processing patterns, specializing the fan-out-and-merge shape for single-scene raster math.

Why Serverless NDVI Tiling Matters for Geospatial Workloads

NDVI (Normalized Difference Vegetation Index) is the workhorse index for vegetation health, crop monitoring, and land-cover change. Computing it is arithmetically trivial — one subtraction and one division per pixel — but the data volume makes it an engineering problem. A single Sentinel-2 Level-2A tile covers 100×100 km at 10 m resolution, which is 10980×10980 pixels per band. Reading two full bands as uint16 and producing a float32 result touches roughly 240 MB just in raw arrays, and a naive whole-scene job flirts with both the memory ceiling and the timeout.

The serverless answer is to treat each scene as a grid of independent windows and process them in parallel. Because Sentinel-2 imagery is published as Cloud-Optimized GeoTIFFs, a worker can issue HTTP range requests to fetch only the bytes for its window rather than downloading the whole file — the same access pattern optimized in tuning HTTP range requests for COG reads on S3. This keeps each function’s /tmp footprint and memory small enough that the workload maps cleanly onto Lambda or Cloud Functions.

The two ingredients that determine whether the pipeline is fast or fragile are band I/O and memory sizing. Reading B04 and B08 efficiently depends on the COG’s internal tiling and overview structure, which is exactly the concern of chunked I/O for large satellite imagery and its Sentinel-2-specific companion optimizing chunked I/O for multi-band Sentinel-2 processing. Sizing the function so numpy has enough vCPU to run the division at speed is governed by the memory and CPU allocation model for raster workloads, where memory and CPU scale together.

The diagram below shows the fan-out from scene event to catalogued tiles:

Platform-by-Platform Limits

Where the per-tile compute runs is the main choice. These are the constraints that decide window size and parallelism.

On AWS, the natural big-scene escape hatch is a Step Functions map that fans out per-tile Lambdas and merges results — the same orchestration used in process a 10 GB GeoTIFF with Step Functions and Lambda. Azure Functions Consumption is the most constrained target: with 1.5 GB of memory and a 10-minute ceiling, keep windows to 1024×1024 or move to a container-based job. GCP’s higher memory and 60-minute ceiling make Cloud Functions 2nd gen forgiving for larger windows.

Step-by-Step Implementation

Step 1: Trigger on a New Scene and Enumerate the Grid

An S3 event (or Eventarc on GCS) fires when a new Level-2A scene lands. An orchestrator function opens the source COG’s profile, reads its width and height, and enumerates fixed windows without reading any pixels.

# enumerate_tiles.py — build the window grid from the COG header only
import rasterio
from rasterio.windows import Window

TILE = 2048  # window edge in pixels

def enumerate_windows(cog_uri: str) -> list[dict]:
    with rasterio.open(cog_uri) as src:  # reads header, not pixels
        w, h = src.width, src.height
    windows = []
    for row_off in range(0, h, TILE):
        for col_off in range(0, w, TILE):
            width = min(TILE, w - col_off)
            height = min(TILE, h - row_off)
            windows.append({
                "col_off": col_off, "row_off": row_off,
                "width": width, "height": height,
            })
    return windows  # each dict becomes one fan-out payload

Step 2: Fan Out One Worker per Tile

Publish each window as a message so the platform runs workers in parallel. On AWS this is an SQS send (or a Step Functions map); on GCP a Pub/Sub publish. The routing and concurrency controls are the same ones covered in SQS and Pub/Sub queue routing strategies.

# fanout.py — publish one message per window
import json, boto3

sqs = boto3.client("sqs")
QUEUE_URL = "https://sqs.us-east-1.amazonaws.com/123456789012/ndvi-tiles"

def fan_out(cog_red: str, cog_nir: str, windows: list[dict]) -> None:
    entries = []
    for i, win in enumerate(windows):
        entries.append({
            "Id": str(i),
            "MessageBody": json.dumps({"red": cog_red, "nir": cog_nir, "window": win}),
        })
        if len(entries) == 10:  # SQS batch limit
            sqs.send_message_batch(QueueUrl=QUEUE_URL, Entries=entries)
            entries = []
    if entries:
        sqs.send_message_batch(QueueUrl=QUEUE_URL, Entries=entries)

Step 3: Window-Read the B04 and B08 Bands

Each worker reads only its window from each single-band COG. Because B04 and B08 are both 10 m bands they share the same grid, so one Window maps identically to both.

# read_bands.py — windowed reads of red and NIR
import rasterio
from rasterio.windows import Window

def read_window(red_uri: str, nir_uri: str, win: dict):
    window = Window(win["col_off"], win["row_off"], win["width"], win["height"])
    # GDAL config for efficient COG range reads over S3/GCS is set via env vars:
    #   GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR
    #   CPL_VSIL_CURL_ALLOWED_EXTENSIONS=.tif
    #   VSI_CACHE=TRUE
    with rasterio.open(red_uri) as red_src:
        red = red_src.read(1, window=window)
        profile = red_src.profile
        red_nodata = red_src.nodata
    with rasterio.open(nir_uri) as nir_src:
        nir = nir_src.read(1, window=window)
        nir_nodata = nir_src.nodata
    return red, nir, profile, red_nodata, nir_nodata, window

Setting GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR and CPL_VSIL_CURL_ALLOWED_EXTENSIONS=.tif prevents GDAL from issuing wasteful directory listings against object storage on every open — a common cause of slow COG reads.

Step 4: Compute NDVI with rasterio and numpy

Cast to float32, suppress the divide-by-zero warning, and mask pixels where either band is nodata or the denominator is zero. The exact masking recipe is expanded in computing NDVI per tile with rasterio and numpy.

# compute_ndvi.py — float32 NDVI with divide-by-zero masking
import numpy as np

NDVI_NODATA = np.float32(-9999.0)

def compute_ndvi(red, nir, red_nodata, nir_nodata):
    red = red.astype("float32")
    nir = nir.astype("float32")
    denom = nir + red

    with np.errstate(divide="ignore", invalid="ignore"):
        ndvi = np.where(denom == 0, NDVI_NODATA, (nir - red) / denom)

    # Propagate source nodata into the NDVI mask
    if red_nodata is not None:
        ndvi = np.where(red == red_nodata, NDVI_NODATA, ndvi)
    if nir_nodata is not None:
        ndvi = np.where(nir == nir_nodata, NDVI_NODATA, ndvi)
    return ndvi.astype("float32")

Step 5: Write NDVI COG Tiles and Register a Catalog Entry

Write the window as a Cloud-Optimized GeoTIFF with the correct geotransform for the window offset, then register the tile in a catalog (a STAC item or a DynamoDB row) so downstream consumers can find it.

# write_tile.py — write an NDVI COG tile and catalog it
import rasterio
from rasterio.windows import transform as win_transform

def write_ndvi_tile(ndvi, src_profile, window, out_uri):
    profile = src_profile.copy()
    profile.update(
        dtype="float32",
        count=1,
        nodata=-9999.0,
        compress="LZW",
        predictor=3,          # floating-point predictor improves compression
        tiled=True,
        blockxsize=512,
        blockysize=512,
        driver="COG",
        width=window.width,
        height=window.height,
        transform=win_transform(window, src_profile["transform"]),
    )
    with rasterio.open(out_uri, "w", **profile) as dst:
        dst.write(ndvi, 1)

Measurement and Verification

Confirm each tile is a valid single-band float32 COG with sane NDVI statistics — values must fall in [-1, 1] for real pixels and equal the nodata sentinel elsewhere.

# verify_tile.py — check NDVI range and COG structure
import rasterio
import numpy as np

def verify(tile_uri: str):
    with rasterio.open(tile_uri) as ds:
        assert ds.count == 1 and ds.dtypes[0] == "float32"
        assert ds.nodata == -9999.0
        arr = ds.read(1, masked=True)
        valid = arr.compressed()
        return {
            "min": float(valid.min()),
            "max": float(valid.max()),
            "mean": round(float(valid.mean()), 3),
            "is_tiled": ds.profile.get("tiled", False),
            "blocksize": ds.block_shapes[0],
        }

Expected output for a healthy vegetated tile:

{
  "min": -0.21,
  "max": 0.89,
  "mean": 0.42,
  "is_tiled": true,
  "blocksize": [512, 512]
}

A max above 1.0 or below -1.0 means the bands were divided as integers (no float32 cast) or the wrong bands were read. Run rio cogeo validate tile.tif to confirm the output satisfies the COG spec — internal tiling and overviews present, IFD ordering correct.

Failure Modes and Debugging

ValueError: Integer overflow or NDVI pinned to 0 and 1: The bands were subtracted and divided as uint16 before casting. Integer division truncates every non-integer NDVI to 0. Always .astype("float32") both bands before the arithmetic.

Out-of-memory kill on the worker: The window is too large for the configured memory, or the code read all bands instead of just B04 and B08. Drop the window to 2048×2048, read exactly two bands, and confirm memory with the memory and CPU allocation model for raster workloads.

COG reads are slow / thousands of GET requests: GDAL is listing the bucket directory or fetching overviews it does not need. Set GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR and tune the read window to the COG’s internal block size, as in tuning HTTP range requests for COG reads on S3.

Nodata bleeding into NDVI as valid values: Source nodata (often 0 in Sentinel-2 L2A) was not masked, so pixels over water or scene edges compute a spurious NDVI. Propagate both bands’ nodata into the output mask before writing.

Timeout on a whole-scene job: A single Lambda tried to process the entire 10980×10980 scene and hit the 15-minute ceiling. Fan out per-tile as in Step 2, or move the orchestration to a Step Functions map for larger scenes.

Cost and Scaling Considerations

Per-tile fan-out turns a serial, timeout-prone job into an embarrassingly parallel one: 30 windows of a scene finish in the wall-clock time of a single window because they run concurrently, bounded only by account concurrency. Cost is dominated by function duration, and since vCPU scales with memory, a larger memory tier that finishes the numpy division in half the time often costs the same or less overall while returning results faster — run the tradeoff math in the memory and CPU allocation model for raster workloads.

I/O, not compute, is usually the real cost driver: each worker issues range GETs against object storage, and unoptimized reads can multiply request counts tenfold. Keep the window aligned to the COG’s internal tiling and cache reads within the invocation. For continent-scale mosaics that exceed the Lambda ceiling, run the identical rasterio+numpy code inside a Cloud Run job or a Step Functions map, both of which lift the 15-minute limit while preserving the per-tile parallelism. The output side benefits from the same tiling discipline as the input: write NDVI as an LZW-compressed, internally tiled float32 COG so downstream consumers can range-read it just as efficiently as you read Sentinel-2.

Frequently Asked Questions

Which Sentinel-2 bands do I read to compute NDVI?

NDVI uses the red and near-infrared bands. On Sentinel-2 these are B04 (red, 665 nm) and B08 (NIR, 842 nm), both at 10 m resolution, which is why they align pixel-for-pixel without resampling. Read only these two bands from the COG using windowed range requests — never download the full multi-band scene, which wastes memory and /tmp space.

Why fan out one Lambda per tile instead of processing the whole scene?

A full Sentinel-2 tile is 10980×10980 pixels per band. Holding both bands plus the float32 NDVI result in memory approaches 1.5 GB, and a single Lambda would risk both the 15-minute timeout and the memory ceiling. Splitting the scene into fixed windows lets many Lambdas run in parallel, each finishing in seconds within a few hundred MB — faster and cheaper than one large invocation.

How do I avoid divide-by-zero when NIR and red are both zero?

The NDVI denominator (nir + red) is zero over nodata and deep-shadow pixels. Compute NDVI with numpy using np.errstate to suppress the warning and np.where to substitute a sentinel value where the denominator is zero, then mask those pixels as nodata. Cast the bands to float32 before the division so integer overflow and truncation cannot occur.

Should the NDVI compute run in Lambda, Cloud Run, or Cloud Functions?

For per-tile windows that finish in seconds, AWS Lambda or GCP Cloud Functions 2nd gen are ideal — event-driven, massively parallel, billed per millisecond. For very large scenes or whole-continent mosaics that exceed the 15-minute Lambda ceiling, run the same rasterio+numpy code on Cloud Run or in a Step Functions map, which removes the hard timeout and lets a container hold more state.

Related

• Computing NDVI per Tile with Rasterio and NumPy — the exact windowed read and float32 masking wrapped for a Lambda handler
• Chunked I/O for Large Satellite Imagery — read only the bytes each tile needs from a COG
• Optimizing Chunked I/O for Multi-Band Sentinel-2 Processing — band-aware read patterns for Sentinel-2 scenes
• Memory and CPU Allocation for Raster Workloads — size each tile worker so numpy has the vCPU it needs
• Process a 10 GB GeoTIFF with Step Functions and Lambda — the orchestration to reach for when a scene outgrows a single Lambda

Back to Serverless Geospatial Pipeline Recipes