Clipping a Raster with a GeoPackage Layer

Name the layer, filter at read time, reproject, then clip:

import fiona, geopandas as gpd, rasterio
from rasterio.mask import mask

print(fiona.listlayers("aoi.gpkg"))                          # never assume the first layer
aoi = gpd.read_file("aoi.gpkg", layer="districts", where="region = 'Western'")
with rasterio.open("scene.tif") as src:
    out, transform = mask(src, aoi.to_crs(src.crs).geometry, crop=True, filled=False)

GeoPackage is the usual way AOIs arrive now, and its multi-layer nature is the first trap. This page belongs to automated image clipping and cropping in Satellite Processing Workflows & Index Pipelines.

One File, Several Layers

Each layer also carries its own CRS, so two layers from the same file may need different transformations. Reprojecting to the raster CRS after reading, every time, removes the question.

Environment & Setup

pip install "geopandas>=0.14" "pyogrio>=0.7" "rasterio>=1.3.0" "shapely>=2.0"

Complete Working Example

from pathlib import Path

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


def read_aoi(gpkg: str, layer: str, *, where: str | None = None,
             bbox: tuple | None = None) -> gpd.GeoDataFrame:
    layers = [name for name, _ in pyogrio.list_layers(gpkg)]
    if layer not in layers:
        raise KeyError(f"layer {layer!r} not in {layers}")
    gdf = gpd.read_file(gpkg, layer=layer, where=where, bbox=bbox, engine="pyogrio")
    gdf["geometry"] = gdf.geometry.make_valid()
    return gdf[~gdf.geometry.is_empty]


def clip_union(raster: str, aoi: gpd.GeoDataFrame, out_path: str) -> None:
    with rasterio.open(raster) as src:
        geom = aoi.to_crs(src.crs).union_all()               # one dissolved AOI
        arr, transform = mask(src, [geom], crop=True, filled=True, nodata=src.nodata or 0)
        profile = src.profile | {"height": arr.shape[1], "width": arr.shape[2],
                                 "transform": transform, "nodata": src.nodata or 0}
    with rasterio.open(out_path, "w", **profile) as dst:
        dst.write(arr)


def clip_per_feature(raster: str, aoi: gpd.GeoDataFrame, id_field: str, out_dir: str) -> int:
    Path(out_dir).mkdir(parents=True, exist_ok=True)
    written = 0
    with rasterio.open(raster) as src:
        feats = aoi.to_crs(src.crs)
        for fid, geom in zip(feats[id_field], feats.geometry):
            if not geom.intersects(box(*src.bounds)):
                continue                                  # feature outside this scene
            arr, transform = mask(src, [geom], crop=True, filled=True, nodata=src.nodata or 0)
            profile = src.profile | {"height": arr.shape[1], "width": arr.shape[2],
                                     "transform": transform, "nodata": src.nodata or 0}
            with rasterio.open(Path(out_dir) / f"{fid}.tif", "w", **profile) as dst:
                dst.write(arr)
            written += 1
    return written

Passing where and bbox to the reader pushes the filtering into SQLite, so a national layer of a million features is never fully loaded when only one region is needed. make_valid before clipping prevents the topology errors hand-digitised boundaries so often contain; the clipping mechanics for irregular shapes are in how to clip rasters to irregular polygon boundaries.

Union or Per Feature

The choice follows the unit of analysis. A study area delivered as one product wants the union. Per-district statistics, per-field chips or per-client deliveries want one output per feature, named by a stable identifier from the layer rather than by row order. For very many small features, writing thousands of tiny GeoTIFFs is slow; computing statistics directly, as in computing zonal statistics with rasterstats, is often what was really wanted.

Keeping AOIs Consistent Across Runs

GeoPackages are edited, and an AOI that silently changes between runs makes results incomparable. Recording the layer name, the filter, the feature count and a hash of the dissolved geometry in each output’s tags turns that into something checkable: if a later run’s hash differs, the AOI changed, and the difference in results has a known cause rather than a mysterious one. The tagging mechanics are in reading and writing GDAL tags and band descriptions.

Clipping Many Scenes to One Layer

A common production pattern is the reverse of the examples above: one AOI layer, many scenes. Reading and reprojecting the layer once per scene is wasteful when every scene shares a CRS, and wrong when they do not, because each scene may sit in a different UTM zone.

The efficient arrangement is to read the layer once, group scenes by CRS, reproject the layer once per distinct CRS, and cache the result. Pre-filtering the features to each scene’s footprint with a spatial index before masking avoids calling mask on features that cannot intersect, which for a large parcel layer is most of them. Both steps are small, and together they turn a job dominated by repeated vector work into one dominated by raster reads, which is where the time should go. The catalog-driven variant of this loop is covered in cropping a STAC item to an AOI without full download.

Verification

import numpy as np
import rasterio

with rasterio.open("clip.tif") as out:
    arr = out.read(1, masked=True)
    px = abs(out.transform.a * out.transform.e)
valid_area = arr.count() * px
aoi_area = aoi.to_crs(out.crs).union_all().area
print(f"valid {valid_area/1e6:.2f} km², AOI {aoi_area/1e6:.2f} km²")
assert valid_area <= aoi_area * 1.05, "clip larger than the AOI"

Common Errors

The clip used the wrong polygons

No layer was named, so the first layer was read. Always pass layer=.

ValueError: Input shapes do not overlap raster

The AOI was not reprojected, or the filter selected features outside the scene. Reproject and check intersection first.

Per-feature outputs overwrite each other

Files were named by a non-unique attribute. Use a stable unique identifier.

Reading the GeoPackage is slow

The whole layer was loaded. Filter with where or bbox at read time.

Frequently Asked Questions

Q: Why does reading my GeoPackage return the wrong features? A GeoPackage can hold many layers, and reading without naming one returns the first. List the layers and pass the layer name explicitly every time.

Q: Should I clip to each feature or to their union? To the union when the features together form one study area, and per feature when each is a separate unit of analysis such as a field or a district. The choice decides whether you get one output or many.

Q: Can I read only the features I need? Yes. GeoPackage is SQLite underneath, so a where clause or a bounding-box filter at read time avoids loading the whole layer, which matters for national parcel datasets.

Q: What about shapefiles and GeoJSON? The same code reads them — they simply have one layer each. GeoPackage is preferable for new work because it holds several layers, keeps full attribute names and stores the CRS unambiguously.

Q: Should the clipped output be cropped or kept at full scene extent? Cropped, almost always: crop=True shrinks the file to the AOI’s bounding box and saves storage and every later read. Keep the full extent only when outputs must stack exactly with other full-scene products.

Related

• Automated Image Clipping and Cropping — the parent topic.
• Handling Partial Overlap between AOI and Scene — when the AOI extends beyond the image.
• Cropping a STAC Item to an AOI without Full Download — the remote-data variant.
• Rasterizing Vector Polygons onto a Raster Grid — the burn rules underneath masking.