Forestry GIS Python · Forestry · Ecology

Merging Shapefiles and Rasters in Python with rasterio.features.rasterize and sample_gen

Operationalizing spatial data integration requires a clear conceptual shift: in a forestry workflow the term “merge” rarely denotes direct concatenation. Instead, it describes a controlled overlay where polygon boundaries are rasterized to a common grid, raster values are sampled into vector attributes, or gridded surfaces are clipped to a vector extent. Vectors and rasters inhabit fundamentally different mathematical domains — continuous Euclidean geometry versus a discrete matrix indexed by an affine transform. This page covers the narrow, high-frequency sub-task of binding a stand-boundary shapefile to a LiDAR-derived raster (canopy height, a vegetation index, or a terrain layer) without corrupting topology or measurement units. It sits within the Raster-Vector Overlay Techniques methodology, which in turn is part of the broader Ecological GIS Data Foundations in Python framework for reproducible conservation analysis.

When to Use Each Overlay Strategy

There is no single “merge” operation — the correct tool depends on whether your output should be a raster, an enriched vector, or a clipped subset. Choose deliberately before writing code, because the wrong strategy silently produces statistically invalid results rather than an error.

Goal Strategy Tool Output
Burn stand IDs / class labels onto the raster grid Rasterize vector → raster rasterio.features.rasterize() Integer/float raster aligned to the source grid
Attach a single pixel value to each plot Point sampling at centroids rasterio.sample.sample_gen Vector with a new attribute column
Aggregate pixels within each polygon (mean, percentile) Zonal statistics rasterstats.zonal_stats() Vector with per-polygon summary columns
Restrict a raster to a polygon extent Clip / mask rasterio.mask.mask() Cropped raster + updated transform

If you only need to clip a raster to a boundary, do not rasterize — use mask(). If you need a per-tree or per-plot value, sampling at centroids is cheap; if you need a distribution across each stand, you must aggregate with zonal statistics. The CRS alignment step below is a non-negotiable prerequisite for all four strategies, and depends on the same projection discipline covered in how to fix CRS mismatches in geopandas.

Choosing an overlay strategy by desired output A stand-boundary shapefile and an aligned raster feed a decision: if the output must be a raster, rasterize the vector; for one value per plot, sample at centroids; for a per-polygon distribution, run zonal statistics; to restrict the raster to a boundary, clip with mask. Each branch names its rasterio or rasterstats function. Stand shapefile (vector boundaries) Aligned raster (same target CRS) What output do you need? A raster of stand IDs features.rasterize() One value per plot sample.sample_gen Per-polygon distribution rasterstats.zonal_stats() Raster clipped to extent mask.mask() Pick the strategy before writing code — the wrong one fails silently, not loudly.

Coordinate Reference System Alignment

The most frequent failure point in overlay workflows is unaligned coordinate systems. When a forest inventory plot shapefile uses EPSG:32615 (WGS 84 / UTM zone 15N) and a LiDAR-derived canopy height raster uses EPSG:4269 (NAD83 geographic), direct overlay either silently produces NaN-filled output or raises rasterio.errors.CRSError. Before any spatial intersection, transform both datasets to a single projected system whose linear units (metres) preserve area for ecological metrics.

Reproject the vector with to_crs() and warp the raster with an explicit transformation pipeline:

import geopandas as gpd
import rasterio
from rasterio.warp import calculate_default_transform, reproject, Resampling

target_crs = "EPSG:32615"  # UTM Zone 15N for area-preserving metrics

# Vector CRS alignment
plots_gdf = gpd.read_file("inventory_plots.shp").to_crs(target_crs)

# Raster CRS alignment
with rasterio.open("canopy_height.tif") as src:
    transform, width, height = calculate_default_transform(
        src.crs, target_crs, src.width, src.height, *src.bounds
    )
    kwargs = src.meta.copy()
    kwargs.update({
        "crs": target_crs,
        "transform": transform,
        "width": width,
        "height": height,
    })
    with rasterio.open("canopy_height_aligned.tif", "w", **kwargs) as dst:
        for i in range(1, src.count + 1):
            reproject(
                source=rasterio.band(src, i),
                destination=rasterio.band(dst, i),
                src_transform=src.transform,
                src_crs=src.crs,
                dst_transform=transform,
                dst_crs=target_crs,
                resampling=Resampling.bilinear,
            )

Use Resampling.bilinear for continuous surfaces (height, NDVI) and Resampling.nearest for categorical rasters (land-cover classes) so resampling never invents fractional class codes. Official reprojection guidance is maintained in the GeoPandas CRS documentation.

Minimal Reproducible Example

The complete, copy-pasteable workflow below rasterizes stand polygons onto the aligned raster grid, then samples the raster back into each plot’s attributes. Both halves assume plots_gdf and the raster already share target_crs from the step above.

import numpy as np
import rasterio
import geopandas as gpd
from rasterio.features import rasterize
from rasterio.sample import sample_gen

ALIGNED = "canopy_height_aligned.tif"

# --- 1. Rasterize stand boundaries onto the source grid -------------------
with rasterio.open(ALIGNED) as src:
    out_shape = src.shape           # (rows, cols)
    out_transform = src.transform   # affine mapping pixel <-> CRS coords
    nodata = src.nodata

    shapes = ((geom, val) for geom, val in
              zip(plots_gdf.geometry, plots_gdf["stand_id"]))

    stand_grid = rasterize(
        shapes=shapes,
        out_shape=out_shape,
        transform=out_transform,
        fill=0,
        all_touched=True,         # keep thin riparian / boundary pixels
        dtype="int32",
    )

# --- 2. Sample the raster value at each plot centroid ---------------------
def sample_raster_at_centroids(raster_path, vector_gdf, column="sampled_value"):
    gdf = vector_gdf.copy()
    coords = [(pt.x, pt.y) for pt in gdf.geometry.centroid]
    with rasterio.open(raster_path) as src:
        if src.crs.to_epsg() != gdf.crs.to_epsg():
            raise ValueError("Vector and raster CRS differ — reproject first.")
        sampled = list(sample_gen(src, coords))
        nodata = src.nodata
    gdf[column] = [
        float(v[0]) if (v[0] != nodata and np.isfinite(v[0])) else np.nan
        for v in sampled
    ]
    return gdf

plots_enriched = sample_raster_at_centroids(ALIGNED, plots_gdf, "canopy_ht_m")

sample_gen streams (x, y) coordinate tuples and reads only the pixels it needs, so it never loads the full raster into RAM — the correct choice for per-plot extraction across a regional inventory. For polygon-level distributions (mean, 95th-percentile height) rather than a single centroid value, swap in rasterstats.zonal_stats(), which performs windowed I/O internally.

Parameter Reference

The defaults of rasterize() and sample_gen are tuned for general use, not for narrow ecological boundaries. The values below are the ones that matter when stand edges, riparian buffers, or nodata gaps must be preserved exactly.

Parameter Type Default Recommended Ecological rationale
all_touched bool False True for thin features False only burns pixels whose centre falls inside a polygon, dropping narrow riparian corridors and stand slivers entirely
fill scalar 0 0 or src.nodata Background value for unburned pixels; reuse the raster’s nodata so empty area is masked, not treated as stand 0
dtype str float64 match source (int32 for IDs, float32 for indices) A float32 index silently truncates if the output array is int32; an ID field overflows int8
out_shape / transform tuple / Affine — always pass src.shape + src.transform Guarantees identical grid topology; mismatched values raise shapes not aligned
resampling (warp) enum nearest bilinear for continuous, nearest for categorical Bilinear smooths height/NDVI; nearest preserves discrete class codes
merge_alg enum MergeAlg.replace MergeAlg.add for overlap counts replace lets later polygons overwrite earlier ones; add accumulates where stands overlap

Expected Output and Verification

A correct rasterization yields an array of out_shape whose unique values are exactly the stand_id set plus the fill value, with non-background pixels covering the union of the polygons. A correct sampling pass yields one finite value per plot, with NaN only where a centroid falls on nodata. Assert these invariants before trusting downstream metrics:

# Grid topology matches the source raster exactly
assert stand_grid.shape == out_shape
assert set(np.unique(stand_grid)) <= set(plots_gdf["stand_id"]).union({0})

# Burned area is plausible (non-empty, not the whole scene)
burned = np.count_nonzero(stand_grid != 0)
assert 0 < burned < stand_grid.size

# Every plot received a value; flag suspicious nodata hits
n_nan = plots_enriched["canopy_ht_m"].isna().sum()
print(f"{n_nan} of {len(plots_enriched)} plots sampled nodata")
assert plots_enriched["canopy_ht_m"].dropna().between(0, 80).all()  # m

For a visual check, write stand_grid to a GeoTIFF with the source profile and overlay it on the original raster in QGIS — the burned polygons must register pixel-for-pixel with the underlying surface. The 0–80 m bound on canopy height is a sanity gate; tighten it to your forest type to catch unit errors (feet vs. metres) early.

How all_touched changes which pixels a polygon burns Two identical raster grids with the same stand polygon and a thin riparian corridor overlaid. With all_touched=False only pixels whose centre falls inside the polygon are burned, and the narrow corridor drops out entirely. With all_touched=True every pixel the polygon intersects is burned, so the thin corridor survives. all_touched=False burns pixels whose centre is inside all_touched=True burns every intersected pixel Thin riparian strip lost Thin riparian strip preserved Shaded cells are burned; the inner line traces the narrow corridor only one pixel wide.

Common Pitfalls

  • CRS mismatch read as empty output. When the vector and raster disagree on projection, sampling returns nodata for every centroid and rasterization burns nothing — verify gdf.total_bounds overlaps src.bounds before blaming the data.
  • ValueError: shapes not aligned. The out_shape you passed does not match the affine transform. Always read both src.shape and src.transform from the same opened raster and pass them together.
  • Silent dtype truncation. Rasterizing a float32 vegetation index into the default int32 floor every value to an integer with no warning. Set dtype explicitly to match the source band.
  • Nodata counted as a real value. sample_gen returns the raster’s raw nodata sentinel (often -9999 or 0); compare each sample against src.nodata and convert to np.nan so means and percentiles are not poisoned.

When integrating these outputs into conservation policy mapping, confirm the merged dataset preserves original measurement units and topological integrity, so downstream work — such as the canopy and biomass metrics built on Raster-Vector Overlay Techniques — stays audit-ready.