Raster-Vector Overlay Techniques for Forestry and Ecological Workflows in Python

Forest inventory and ecological monitoring constantly require pairing a continuous raster surface — a LiDAR canopy height model, a Sentinel-2 reflectance tile, a modelled climate grid — with discrete vector boundaries such as stands, plots, riparian buffers, or conservation easements. The recurring technical challenge is computing a defensible per-polygon summary of the raster (mean canopy height per stand, NDVI percentiles per management unit, biomass density per ownership parcel) without introducing the silent geometric and statistical errors that creep in when grids and polygons are misaligned. This page builds a reproducible, auditable Python overlay pipeline that a conservation agency or research team can run unattended across thousands of stands. It is one component of the broader Ecological GIS Data Foundations in Python workflow, which orchestrates the upstream and downstream steps this overlay depends on.

A worked scenario anchors the rest of the page: you have a 1 m LiDAR-derived canopy height raster covering a watershed and a shapefile of ~4,000 forest stand polygons, and you need mean and 95th-percentile canopy height for every stand, exported with the original attribute table intact.

Prerequisites

Confirm each of these before running the pipeline. Mismatches here are the dominant cause of wrong-but-plausible output later.

• Python 3.10+ with rasterio ≥ 1.3, geopandas ≥ 0.14, rasterstats ≥ 0.19, shapely ≥ 2.0, and numpy ≥ 1.24 installed in a clean environment.
• GDAL/PROJ available through the wheels (the rasterio and geopandas binary wheels bundle compatible versions — avoid mixing a system GDAL with pip-installed bindings).
• The raster has a defined, valid CRS and an explicit nodata value (not None) recorded in its metadata.
• Vector boundaries carry a defined CRS; if it differs from the raster you will reproject the vector to the raster’s CRS, never the reverse for a single overlay.
• Datum realisations are reconciled — see Coordinate Reference Systems for Forestry — so that NAD27/NAD83/ITRF shifts do not offset polygons from the grid by metres.
• You know whether area-weighted statistics are required; if so, the raster should be in an equal-area projection with square pixels.

Concept Background: What “Overlay” Means for a Grid and a Polygon

A raster is a regular matrix of cells, each anchored to ground coordinates by an affine transform; a vector polygon is continuous Euclidean geometry. Zonal extraction reconciles the two by deciding which cells belong to each polygon and how partial cells are weighted.

For a polygon $Z$ overlaid on raster cells indexed by $i$, the zonal mean is the coverage-weighted average of cell values $v_i$:

$${\bar{v}}_Z=\frac{{\sum }_{i\in Z}{w}_i{v}_i}{{\sum }_{i\in Z}{w}_i},w_i=\frac{a_i}{A_{\text{cell}}}$$

where $a_i$ is the area of cell $i$ that falls inside the polygon and $A_{\text{cell}}$ is the full cell area, so $w_i\in [0,1]$. Two simplifications dominate practice:

• Centroid sampling (all_touched=False): a cell counts fully ($w_i=1$) only if its centroid is inside the polygon, otherwise $w_i=0$. This is the standard, unbiased choice for forestry plots larger than a few cells.
• All-touched (all_touched=True): every cell intersected by the polygon counts fully. This inflates small or thin polygons (riparian strips, narrow buffers) by pulling in edge cells and should be reserved for cases where you would rather over-include than miss a feature.

The three rules differ only in how they treat the boundary cells — the cells the polygon edge cuts through — and that single decision drives the bias between methods:

True fractional weighting ($0<w_i<1$) requires a coverage-fraction engine such as exactextract; for most stand-level inventory the centroid rule is accurate and far cheaper. The bias from centroid sampling shrinks as polygon area grows relative to cell size, which is why grid resolution and minimum mapping unit should be chosen together.

Step-by-Step Python Pipeline

The pipeline runs in four ordered stages: harmonise CRS, repair topology, mask and extract, then export. Each step below is runnable in isolation and composes into the final driver.

Step 1 — Harmonise the coordinate framework

Treat the raster CRS as authoritative for a single overlay and reproject the vector to match it. Reprojecting the raster instead would resample every cell and blur the values you are about to summarise.

import geopandas as gpd
import rasterio

def align_to_raster(vector_path: str, raster_path: str) -> gpd.GeoDataFrame:
    """Load vector boundaries and reproject them onto the raster's CRS."""
    gdf = gpd.read_file(vector_path)
    if gdf.crs is None:
        raise ValueError(f"{vector_path} has no CRS; assign the correct EPSG before overlay.")
    with rasterio.open(raster_path) as src:
        raster_crs = src.crs
    if raster_crs is None:
        raise ValueError(f"{raster_path} has no CRS; the raster cannot be georeferenced.")
    return gdf.to_crs(raster_crs)

Step 2 — Repair vector topology

Stand polygons exported from desktop GIS routinely carry self-intersections, slivers, and unclosed rings. These trigger silent failures during sampling — null statistics or misattributed cells — so repair them before they reach the raster. The same discipline underpins Merging Shapefiles and Rasters in Python, where mismatched grid topology is the recurring failure.

import geopandas as gpd
from shapely.validation import make_valid

def clean_vector(gdf: gpd.GeoDataFrame) -> gpd.GeoDataFrame:
    """Fix invalid geometries and drop empty or degenerate features."""
    gdf = gdf.copy()
    invalid = ~gdf.geometry.is_valid
    if invalid.any():
        gdf.loc[invalid, "geometry"] = gdf.loc[invalid, "geometry"].apply(make_valid)
    # make_valid can emit GeometryCollections; keep only (multi)polygon parts.
    gdf["geometry"] = gdf.geometry.apply(
        lambda g: g if g.geom_type in ("Polygon", "MultiPolygon") else g.buffer(0)
    )
    gdf = gdf[~gdf.geometry.is_empty].copy()
    return gdf

Step 3 — Mask, extract, and aggregate

With clean, aligned inputs the core overlay proceeds through masking (clip the grid to the vector extent to cut memory and edge artefacts) and zonal extraction (per-polygon statistics). rasterstats.zonal_stats reads the source raster directly and intersects each geometry natively.

import geopandas as gpd
import rasterio
from rasterio.mask import mask
from rasterstats import zonal_stats

def run_forestry_overlay(raster_path: str, gdf: gpd.GeoDataFrame) -> gpd.GeoDataFrame:
    """Per-stand canopy statistics from a raster, with attributes preserved."""
    with rasterio.open(raster_path) as src:
        raster_crs = src.crs
        src_nodata = src.nodata          # capture before the dataset closes
        # Clip once to confirm overlap and bound memory; out_image unused downstream.
        mask(src, gdf.geometry, crop=True)

    if src_nodata is None:
        raise ValueError(
            "Raster nodata is None; rasterstats would treat fill as a real value. "
            "Set an explicit nodata before overlay."
        )

    stats = zonal_stats(
        vectors=gdf,
        raster=raster_path,
        stats=["mean", "std", "min", "max", "count"],
        percentiles=[95],                # adds a "percentile_95" field per zone
        nodata=src_nodata,
        all_touched=False,               # centroid sampling — standard for stand plots
        geojson_out=False,
    )

    out = gdf.copy()
    for key in ("mean", "std", "min", "max", "count", "percentile_95"):
        out[f"chm_{key}"] = [s.get(key) for s in stats]
    return gpd.GeoDataFrame(out, geometry=gdf.geometry, crs=raster_crs)

The percentiles=[95] argument adds a percentile_95 field to each result, giving the 95th-percentile canopy height the scenario asked for alongside the mean. Because attributes are joined back onto the original GeoDataFrame, every stand identifier and ownership field travels with its new statistics.

Step 4 — Compose and export

def overlay_to_file(raster_path: str, vector_path: str, output_path: str) -> gpd.GeoDataFrame:
    gdf = clean_vector(align_to_raster(vector_path, raster_path))
    result = run_forestry_overlay(raster_path, gdf)
    # GeoPackage preserves geometry + dtypes; CSV drops geometry but is portable.
    result.to_file(output_path, driver="GPKG")
    print(f"Overlay complete: {len(result)} stands written to {output_path}.")
    return result

Configuring rasterstats for Forestry Output

zonal_stats is controlled by a small set of parameters that materially change the result. Pin them deliberately rather than relying on defaults.

zonal_stats(
    vectors=gdf,                       # GeoDataFrame or path; CRS must match raster
    raster="chm.tif",                  # source raster, read window-by-window
    stats=["mean", "std", "count"],    # named reducers over each zone
    percentiles=[90, 95],              # appended as percentile_90 / percentile_95
    nodata=-9999,                      # MUST match the raster fill value
    all_touched=False,                 # centroid rule; True over-includes edge cells
    categorical=False,                 # True for class rasters (e.g. land-cover counts)
    boundless=True,                    # tolerate zones extending past the raster edge
)

For categorical inputs such as a forest-type or fire-fuel class raster, set categorical=True and zonal_stats returns a per-zone count of each class code instead of summary statistics — useful when an overlay feeds a fuel or habitat-class tally rather than a continuous metric.

Validation and Verification

Never trust an overlay you have not checked for spatial overlap and nodata handling. Run this guard before the expensive extraction; it fails loudly if the inputs cannot produce valid statistics.

import geopandas as gpd
import rasterio
from shapely.geometry import box

def check_overlay_alignment(raster_path: str, vector_path: str) -> bool:
    """Return True only if vector and raster extents overlap in a shared CRS."""
    with rasterio.open(raster_path) as src:
        raster_box = box(*src.bounds)
        raster_crs = src.crs
        if src.nodata is None:
            raise ValueError("Raster nodata is undefined; set it before overlay.")
        if abs(src.res[0] - src.res[1]) > 1e-9:
            print(f"Warning: non-square pixels {src.res}; area-weighting needs care.")
    gdf = gpd.read_file(vector_path).to_crs(raster_crs)
    if not raster_box.intersects(box(*gdf.total_bounds)):
        raise ValueError(
            "No spatial overlap after CRS alignment.\n"
            f"Raster bounds: {tuple(round(b, 2) for b in raster_box.bounds)}\n"
            f"Vector bounds: {tuple(round(b, 2) for b in gdf.total_bounds)}"
        )
    return True

Three post-run sanity checks confirm the numbers are real:

1. Coverage — every polygon should have chm_count > 0. A zero count means the stand fell on nodata or off-grid; investigate before reporting.
2. Plausible range — chm_max must not exceed the raster’s physical ceiling (a 1 m CHM should not report 400 m heights; if it does, nodata leaked into the statistics).
3. Spot-check against a manual zonal mean — clip one stand in a desktop GIS and compare; agreement within rounding confirms the centroid rule and nodata mask behaved as expected.

Failure Modes and Gotchas

• CRS mismatch slips through — if the vector lacks a CRS, to_crs cannot reproject it and the geometries land at the wrong ground location, producing all-nodata zones. Validate CRS presence at ingestion, not after.
• nodata left as None — rasterstats then counts the fill value (often -9999 or 0) as real data, dragging means toward the fill. Always set an explicit nodata on the raster.
• all_touched=True on thin polygons — riparian buffers and stand edges absorb neighbouring cells and report inflated means; reserve all-touched for cases where omission is worse than over-inclusion.
• Non-square pixels in area-weighted work — rectangular cells break the simple area weighting and silently bias hectare-normalised metrics; resample to square pixels first or use a coverage-fraction engine.
• NaN propagation — a single NaN in a band with nodata=None poisons mean/std for the whole zone; mask or set nodata before reducing.
• Memory blow-up on big rasters — calling mask(..., crop=True) on a continental tile materialises the whole window in RAM; window-read or tile instead (below).

Performance and Scale Notes

The watershed scenario (4,000 stands, one tile) runs comfortably in memory, but regional inventories spanning hundreds of tiles need a different strategy.

• Let zonal_stats window-read. Passing the raster path (not an in-memory array) lets rasterstats read only the window each polygon needs, keeping memory flat regardless of total raster size. Avoid pre-masking the whole raster into RAM.
• Partition by tile, parallelise by zone group. Split the vector by raster tile with a spatial join, then process tile groups across processes with concurrent.futures.ProcessPoolExecutor. Zonal statistics are embarrassingly parallel because polygons are independent.
• Use a spatial index for the join. geopandas.sjoin backed by the built-in R-tree assigns each polygon to its tile in near-linear time, avoiding an O(n·m) bounding-box scan.
• Stream results to a columnar store. Append each tile’s output to a Parquet or GeoPackage rather than concatenating GeoDataFrames in memory, so peak memory stays bounded by the largest tile.

These outputs rarely terminate here: extracted canopy metrics commonly combine with multispectral indices to validate forest health, so a tuned overlay is the on-ramp to Vegetation Index Calculation in Python, and the polygon layer itself often originates from a defensible Spatial Plot Sampling Design.

Frequently Asked Questions

Related

• Coordinate Reference Systems for Forestry — reconcile datums and projections before any overlay.
• Merging Shapefiles and Rasters in Python — grid-topology matching and rasterisation for the same integration problem.
• Vegetation Index Calculation in Python — turn extracted metrics into NDVI/EVI forest-health indicators.
• Spatial Plot Sampling Design — produce the polygon layer the overlay summarises.

Up: Ecological GIS Data Foundations in Python