Cloud-Native Raster Workflows for Forestry in Python

The traditional shape of a satellite forestry workflow is: find the scenes, download 40 GB, unzip, reproject everything, then compute over a study area that occupies 2% of what you downloaded. Cloud-native geospatial inverts that: the imagery stays where it is, an index describes what exists, and your code reads only the bytes covering your stands. The concrete problem this page solves is a regional NDVI or burn-severity time series over a forest ownership, computed from public Sentinel-2 or Landsat archives without a bulk download and without a machine large enough to hold a year of imagery. It sits inside Ecological GIS Data Foundations in Python and feeds the analysis in Vegetation Index Calculation in Python.

Three technologies carry the whole pattern, and each solves one problem. The Cloud Optimized GeoTIFF (COG) makes a remote raster randomly readable, so a window costs a few HTTP range requests instead of a file transfer. SpatioTemporal Asset Catalog (STAC) makes the archive searchable by space, time and metadata without a proprietary API. Chunked arrays — xarray over dask — let you express an operation over a hundred scenes as if it were one array and have it execute in bounded memory.

Prerequisites

• Python ≥ 3.10 with rasterio ≥ 1.3, pystac-client ≥ 0.7, stackstac ≥ 0.5 or odc-stac, xarray ≥ 2023.1, and dask ≥ 2023.1.
• GDAL ≥ 3.6 built with curl support, so /vsicurl/ and the COG driver are available (gdalinfo --formats | grep COG).
• A study area as a vector layer in a known CRS — stand boundaries, an ownership, or a plot buffer.
• Network access to a STAC endpoint, and credentials if the collection is requester-pays.
• A decision about the analysis CRS and grid, made before any read, so every scene lands on one grid.
• GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR exported — without it every remote open issues a directory listing that can dominate the runtime.

If the CRS side of that list is not settled, resolve it first with Coordinate Reference Systems for Forestry; a lazy stack built on an ambiguous CRS reprojects wrongly at every chunk and the error is very hard to see.

What makes a GeoTIFF cloud-optimized

A COG is an ordinary GeoTIFF with two disciplines: the pixels are stored in internal tiles rather than strips, and the file carries an overview pyramid, with all of it laid out so the header comes first. A reader can then fetch the header, learn which tiles cover the window it wants, and issue byte-range requests for only those tiles.

The practical consequence for forestry is a change in cost model. Reading a 5 km × 5 km stand from a 110 km Sentinel-2 tile touches a few hundred internal tiles out of tens of thousands, so it costs megabytes rather than the whole scene. Reading a whole scene from a COG is no cheaper than reading a plain GeoTIFF — cloud-native pays when your area of interest is small relative to the archive, which is the normal case in forest management.

Stage 1 — Search the archive with STAC

STAC replaces per-provider APIs with one query grammar: a bounding box, a datetime range, a collection, and arbitrary property filters. Searching does not read pixels, so it is cheap enough to run interactively while you decide what you need. The details of building a robust query, paging results, and filtering on scene-level cloud cover are covered in Searching STAC Catalogs for Sentinel-2 with pystac-client.

What matters at this level is that a search returns items, each holding assets — one per band — whose hrefs point at COGs. Nothing is downloaded. The whole plan of the analysis is decided here: which scenes, which bands, which dates.

Stage 2 — Read only the window you need

With an asset URL in hand, a windowed read costs a few range requests. rasterio handles this transparently through GDAL’s virtual file system; the practical work is converting a study-area geometry into a pixel window in the raster’s own CRS. That conversion — and the surprisingly easy ways to get it wrong — is the subject of Reading Cloud Optimized GeoTIFFs with rasterio Windows.

The environment matters as much as the code. Three GDAL settings change remote read performance by an order of magnitude:

import os

import rasterio

GDAL_ENV = {
    # Without this, every open lists the containing directory — often the dominant cost.
    "GDAL_DISABLE_READDIR_ON_OPEN": "EMPTY_DIR",
    # Only these extensions are considered sidecars, so no speculative .aux.xml requests.
    "CPL_VSIL_CURL_ALLOWED_EXTENSIONS": ".tif,.TIF,.tiff",
    # A block cache large enough to hold the tiles a window touches.
    "GDAL_CACHEMAX": "512",
    "VSI_CACHE": "TRUE",
    "VSI_CACHE_SIZE": "134217728",
}


def open_remote(url: str):
    """Open a remote COG with the settings that make windowed reads fast."""
    with rasterio.Env(**GDAL_ENV):
        return rasterio.open(url)

Stage 3 — Stack scenes lazily

A time series is many items sharing a grid. Building it eagerly means holding every scene in memory; building it lazily means describing the stack and letting dask fetch chunks only as a computation demands them. stackstac and odc-stac both turn a STAC item collection into an xarray.DataArray with dimensions of time, band, y and x — the reprojection to a common grid is declared rather than performed.

The mechanics of chunk sizing, of why an ill-chosen chunk shape can make a lazy stack slower than a download, and of when to persist an intermediate, are in Chunked Raster Processing with xarray and dask.

The pattern for a forest NDVI time series is short enough to state here in full:

import stackstac
import pystac_client


def ndvi_series(bbox, start: str, end: str, epsg: int = 32633, res: int = 20):
    """Lazy NDVI time series over a bounding box; nothing is read until .compute()."""
    client = pystac_client.Client.open("https://earth-search.aws.element84.com/v1")
    search = client.search(
        collections=["sentinel-2-l2a"],
        bbox=bbox,
        datetime=f"{start}/{end}",
        query={"eo:cloud_cover": {"lt": 40}},
    )
    items = list(search.items())

    stack = stackstac.stack(
        items,
        assets=["red", "nir", "scl"],
        epsg=epsg,
        resolution=res,
        bounds_latlon=bbox,
        chunksize=1024,
    )
    red = stack.sel(band="red").astype("float32")
    nir = stack.sel(band="nir").astype("float32")
    scl = stack.sel(band="scl")

    # Keep only vegetation, bare soil, water and unclassified; drop cloud and shadow.
    clear = scl.isin([4, 5, 6, 7])
    ndvi = ((nir - red) / (nir + red)).where(clear)
    return ndvi.rename("ndvi")

Note what has not happened: no scene was read, no reprojection ran, and the result is a description of an array. A subsequent .median(dim="time").compute() fetches only the chunks that median needs.

Stage 4 — Keep provenance with the result

A cloud-native workflow removes the one artefact that used to prove what you analysed: the downloaded file. Reproducibility now depends on recording the query and the item identifiers, because the archive can and does change — scenes are reprocessed, baselines shift, and a collection’s definition of a band can be revised.

import json
import hashlib
from datetime import datetime, timezone


def write_provenance(items, params: dict, out_path: str) -> str:
    """Record exactly which STAC items produced a result, and how."""
    ids = sorted(item.id for item in items)
    record = {
        "generated_utc": datetime.now(timezone.utc).isoformat(timespec="seconds"),
        "n_items": len(ids),
        "item_ids": ids,
        "item_id_digest": hashlib.sha256("\n".join(ids).encode()).hexdigest(),
        "parameters": params,
    }
    with open(out_path, "w", encoding="utf-8") as fh:
        json.dump(record, fh, indent=2)
    return record["item_id_digest"]

Store that JSON beside every derived raster. When a result cannot be reproduced a year later, the digest tells you immediately whether the input set changed or your code did — a distinction that is otherwise nearly impossible to recover.

The Python library ecosystem

python -m pip install "rasterio>=1.3" "pystac-client>=0.7" "stackstac>=0.5" \
    "xarray>=2023.1" "dask[array]>=2023.1" "rioxarray>=0.15" geopandas

Production principles

Declare the grid once. Every scene should be resampled to one CRS, resolution and origin at stack construction, not per operation. A stack whose grid is decided implicitly by the first item will silently shift when that item changes.

Pin the collection, not just the dates. sentinel-2-l2a from two providers can differ in processing baseline and scaling. Record the catalogue URL and collection id in the provenance file.

Cache deliberately. Repeated interactive work over the same stands should persist an intermediate — a clipped, gridded, cloud-masked stack written once as a local COG or Zarr store — rather than re-fetching. The cloud-native argument is against bulk downloads, not against caching what you have already decided you need.

Fail on the CRS, never coerce it. If an item’s CRS is not what the stack declared, that is a data problem worth an exception. Silently reprojecting a stray scene produces a time series with one misaligned date, which looks exactly like real change.

Budget requests, not just bytes. A stack of 200 dates × 3 bands is 600 remote datasets. Opening them serially can take longer than the reads. Use a thread pool for opens and keep GDAL_DISABLE_READDIR_ON_OPEN set.

When Cloud-Native Is the Wrong Choice

The pattern has a cost, and it is worth naming the cases where paying it is a mistake.

The first is a study area that covers most of a scene. Windowed reads win because they transfer a fraction of the archive; when the fraction approaches one, they transfer the same bytes with more requests and more latency. A national analysis at 20 m over a country that fits inside two Sentinel-2 tiles is faster to download once.

The second is repeated interactive work over the same small area. The first read of a remote scene pays for the header and the tile index, and that cost recurs in every new process. Persisting a clipped, gridded, cloud-masked stack locally after the first pass turns an afternoon of exploration from a sequence of network round trips into local array work. Cloud-native is an argument against bulk downloads, not against caching what you have already decided you need.

The third is anything that must be defensible over years. Archives are not immutable: items are reprocessed, occasionally withdrawn, and a collection’s band definitions can be revised. A published result that reruns against a live catalogue may not reproduce, and without recorded item identifiers there is no way to tell whether the code or the data moved. For work of that kind, persist the exact stack used and treat the cloud-native path as how you found the data rather than how you store it.

The fourth is a compute-bound analysis. The lazy-stack pattern is optimised for I/O-bound work, where latency dominates and a threaded scheduler is the right tool. A per-pixel model fit over a large region is bound by CPU instead, and the same code will be slower than a straightforwardly parallel loop over materialised chunks, because the scheduler is solving the wrong problem.

Provenance in Practice

Because the downloaded file no longer exists as evidence, a cloud-native workflow has to manufacture its own. The minimum record is small: the catalogue URL, the collection identifier, the full query, the sorted list of item identifiers with a digest over them, the declared CRS, resolution and grid origin, and the code version that produced the result. Written as a JSON sidecar next to every derived raster, that record answers the question that otherwise cannot be answered a year later — whether a result that will not reproduce changed because the analysis changed or because the archive did.

It is also worth recording what was excluded. A cloud mask drops observations, a cloud-cover filter drops scenes, and a date window drops seasons; a composite built from three clear observations and one built from thirty are rendered identically and mean very different things. Shipping the per-pixel clear-observation count alongside the composite is the cheapest uncertainty layer available, and it is the one most often left out.

Frequently Asked Questions

Related

• Reading Cloud Optimized GeoTIFFs with rasterio Windows — the windowed read in detail
• Searching STAC Catalogs for Sentinel-2 with pystac-client — building queries that return the scenes you meant
• Chunked Raster Processing with xarray and dask — chunk sizing and bounded-memory execution
• Vegetation Index Calculation in Python — the analysis these stacks usually feed
• Coordinate Reference Systems for Forestry — settle the grid before building a stack

Up: Ecological GIS Data Foundations in Python