How to Batch Clip Many Layers to a Single Boundary in Python

You have a study area — a catchment, a city limit, a national park — and a folder full of roads, rivers, parcels, and land-use layers that all sprawl far beyond it. You need every layer trimmed to that one boundary so downstream analysis only ever touches the area you care about. Clicking through a desktop tool once per layer is slow and inconsistent, and the layers rarely share the boundary's CRS. This guide shows you how to load the boundary once, march through the folder, reproject each layer to match, clip it with gpd.clip(), skip anything that falls outside, drop empty results, and write clean output — without one bad file killing the run.

Problem statement

You have one boundary polygon and a directory of vector layers, and you are hitting some mix of these symptoms:

The same clip repeated by hand — you open each layer, load the mask, run clip, export, and repeat, which is tedious and easy to do inconsistently across a dozen files.
CRS mismatch stops the clip — gpd.clip() needs the layer and the mask in the same CRS, and your folder is a jumble of EPSG:4326, a national grid, and Web Mercator.
Layers that never touch the boundary — some datasets sit entirely outside the study area, so clipping them wastes time and returns an empty frame you then try to write.
Empty or invalid outputs — a clip that returns zero features still gets written as a useless file, and one invalid geometry raises an error that aborts the whole batch.
No consistent output — clipped files end up scattered with inconsistent names and formats, so you cannot tell what was trimmed from what.

The goal: one repeatable script that reads the boundary a single time, clips every layer to it in the correct CRS, quietly skips the layers that miss, drops the empty results, and writes a tidy set of clipped outputs with a summary of what happened.

Quick answer

Read the boundary once and note its CRS. Loop over the folder with pathlib, read each layer, reproject it to the boundary CRS with to_crs(), then clip with gpd.clip(layer, boundary). Wrap each file in its own try/except so one bad layer does not abort the batch, skip any layer whose bounds do not overlap the boundary, drop clips that come back empty, and write the rest to GeoPackage.

Several vector layers each clipped to one shared study-area boundary, producing trimmed output layers. — Load the boundary once, then reproject and clip every layer in the folder to that single mask.

import geopandas as gpd
from pathlib import Path

src_dir = Path("data/layers")
out_dir = Path("data/clipped")
out_dir.mkdir(parents=True, exist_ok=True)

# Load the study-area boundary ONCE and dissolve it to a single mask.
boundary = gpd.read_file("data/boundary.gpkg")
boundary = boundary.dissolve()          # one polygon covering the study area
mask_crs = boundary.crs

PATTERNS = ("*.shp", "*.geojson", "*.gpkg")

summary = []
for path in sorted(p for pat in PATTERNS for p in src_dir.glob(pat)):
    try:
        layer = gpd.read_file(path)
        if layer.crs != mask_crs:
            layer = layer.to_crs(mask_crs)          # clip needs matching CRS

        clipped = gpd.clip(layer, boundary)
        if clipped.empty:
            summary.append((path.name, 0, "skipped: no overlap"))
            continue

        out_path = out_dir / f"{path.stem}_clip.gpkg"
        clipped.to_file(out_path, driver="GPKG")
        summary.append((path.name, len(clipped), "clipped"))
    except Exception as exc:
        summary.append((path.name, 0, f"error: {exc}"))

for name, n, status in summary:
    print(f"{name:28} {n:6} features  | {status}")

That is the whole job. The rest of this article explains each decision — dissolving the mask, matching CRS, skipping misses, dropping empties — so you can adapt it safely instead of running it blind. For the mechanics of a single clip, see how to clip spatial data in Python with GeoPandas; for the general folder-loop pattern, see the batch-processing cornerstone guide.

Step-by-step solution

Load the boundary once

The boundary is fixed for the whole run, so read it a single time before the loop — never re-read it per file. If the boundary file has several rows (multiple polygons, or a country split into regions), collapse it to one geometry with dissolve() so gpd.clip() treats the whole study area as a single mask.

import geopandas as gpd

boundary = gpd.read_file("data/boundary.gpkg")
boundary = boundary.dissolve()      # merge all rows into one mask polygon

print(boundary.crs)                 # remember this CRS
print(len(boundary), "mask feature")

Reading the mask once keeps the loop fast and guarantees every layer is clipped against exactly the same shape.

Match the CRS before every clip

gpd.clip() compares coordinates directly, so the layer and the boundary must share a CRS. In GeoPandas 1.x a mismatch raises an error rather than guessing. Reproject each layer to the boundary's CRS with to_crs() inside the loop, and guard against layers that have no CRS at all.

import geopandas as gpd

boundary = gpd.read_file("data/boundary.gpkg").dissolve()
layer = gpd.read_file("data/layers/roads.shp")

if layer.crs is None:
    raise ValueError("Layer has no CRS; set it before clipping")

if layer.crs != boundary.crs:
    layer = layer.to_crs(boundary.crs)   # move layer onto the mask's CRS

It does not matter whether you reproject the layer to the mask or the mask to the layer, but reprojecting each layer to one fixed mask CRS keeps every output in the same frame. See how to reproject spatial data in GeoPandas if the CRS handling is new to you.

Clip the layer to the boundary

With CRS matched, clipping is one call. gpd.clip(layer, boundary) keeps only the parts of layer that fall inside the mask and cuts geometries that cross the edge — a road that leaves the study area is trimmed at the boundary line, not dropped whole.

import geopandas as gpd

clipped = gpd.clip(layer, boundary)

print(len(clipped), "features after clip")
print(clipped.total_bounds)          # now within the boundary extent

gpd.clip() preserves the input layer's attributes and CRS, so clipped is ready to write straight to disk.

Skip layers that never intersect

A layer sitting entirely outside the study area clips to nothing. You can let gpd.clip() run and return an empty frame, but it is faster to check the bounding boxes first and skip the layer before doing any geometric work.

import geopandas as gpd

# Quick reject: do the extents even overlap?
lx0, ly0, lx1, ly1 = layer.total_bounds
bx0, by0, bx1, by1 = boundary.total_bounds

overlaps = not (lx1 < bx0 or lx0 > bx1 or ly1 < by0 or ly0 > by1)
if not overlaps:
    print("Layer is outside the study area — skipping")

Bounding-box overlap is a cheap necessary condition: if the extents miss, the geometries cannot intersect. If they do overlap, the real clip decides what actually falls inside.

Drop empty clips and write output

Even when extents overlap, a clip can still return zero features (the layer skirts the boundary without entering it). Never write an empty result. Check clipped.empty, skip it if true, and write the rest to GeoPackage with a consistent _clip suffix.

import geopandas as gpd
from pathlib import Path

out_dir = Path("data/clipped")
out_dir.mkdir(parents=True, exist_ok=True)

clipped = gpd.clip(layer, boundary)
if clipped.empty:
    print("Nothing inside the boundary — not writing")
else:
    clipped.to_file(out_dir / "roads_clip.gpkg", driver="GPKG")

GeoPackage is the better output format here: one file per layer, the CRS stored internally, no field-name truncation, and no sidecar files to lose.

Handle errors per file

Wrap the read-reproject-clip-write cycle for each file in try/except so a corrupt geometry, a missing CRS, or an unreadable file produces a logged error instead of a stack trace that stops the whole batch. Collect one summary row per file as you go.

import geopandas as gpd
from pathlib import Path

src_dir = Path("data/layers")
out_dir = Path("data/clipped")
out_dir.mkdir(parents=True, exist_ok=True)
boundary = gpd.read_file("data/boundary.gpkg").dissolve()

for path in sorted(src_dir.glob("*.shp")):
    try:
        layer = gpd.read_file(path).to_crs(boundary.crs)
        clipped = gpd.clip(layer, boundary)
        if clipped.empty:
            print(f"SKIP {path.name}: no overlap")
            continue
        clipped.to_file(out_dir / f"{path.stem}_clip.gpkg", driver="GPKG")
        print(f"OK   {path.name}: {len(clipped)} features")
    except Exception as exc:
        print(f"FAIL {path.name}: {exc}")

Code examples

Example 1: Clip one layer as the building block

Get a single clip working before you wrap it in a loop. Read the layer, read and dissolve the boundary, match the CRS, then clip.

import geopandas as gpd

layer = gpd.read_file("data/layers/roads.shp")
boundary = gpd.read_file("data/boundary.gpkg").dissolve()

if layer.crs != boundary.crs:
    layer = layer.to_crs(boundary.crs)

clipped = gpd.clip(layer, boundary)
clipped.to_file("data/clipped/roads_clip.gpkg", driver="GPKG")

print(f"{len(layer)} -> {len(clipped)} features after clipping")

Once this runs cleanly on one file, the folder loop is just this block repeated with error handling around it.

Example 2: The folder loop — reproject then clip

This is the core batch pattern: load the mask once, then reproject and clip every layer in the folder to it.

from pathlib import Path
import geopandas as gpd

src_dir = Path("data/layers")
out_dir = Path("data/clipped")
out_dir.mkdir(parents=True, exist_ok=True)

boundary = gpd.read_file("data/boundary.gpkg").dissolve()
PATTERNS = ("*.shp", "*.geojson", "*.gpkg")

for path in sorted(p for pat in PATTERNS for p in src_dir.glob(pat)):
    layer = gpd.read_file(path)

    if layer.crs is None:
        print(f"Skipping {path.name}: no CRS")
        continue
    if layer.crs != boundary.crs:
        layer = layer.to_crs(boundary.crs)

    clipped = gpd.clip(layer, boundary)
    clipped.to_file(out_dir / f"{path.stem}_clip.gpkg", driver="GPKG")
    print(f"{path.name}: {len(clipped)} features")

Example 3: Skip non-intersecting layers via bounds and `intersects`

For big folders, reject layers that cannot possibly overlap before spending time on the clip. A cheap bounding-box test filters most misses; intersects against the mask geometry is the exact check.

from pathlib import Path
import geopandas as gpd

boundary = gpd.read_file("data/boundary.gpkg").dissolve()
mask_geom = boundary.geometry.iloc[0]
bx0, by0, bx1, by1 = boundary.total_bounds

for path in sorted(Path("data/layers").glob("*.shp")):
    layer = gpd.read_file(path).to_crs(boundary.crs)
    lx0, ly0, lx1, ly1 = layer.total_bounds

    # Cheap extent test first, exact geometry test second.
    if lx1 < bx0 or lx0 > bx1 or ly1 < by0 or ly0 > by1:
        print(f"{path.name}: extent outside study area — skip")
        continue
    if not layer.intersects(mask_geom).any():
        print(f"{path.name}: no feature touches the boundary — skip")
        continue

    clipped = gpd.clip(layer, boundary)
    print(f"{path.name}: {len(clipped)} features clipped")

Example 4: Clip to a bounding box instead of a polygon mask

Sometimes you only want a rectangular window, not the exact boundary shape. gpd.clip() accepts a (xmin, ymin, xmax, ymax) tuple directly, and shapely.box() turns those same numbers into a polygon if you prefer an explicit geometry.

import geopandas as gpd
from shapely.geometry import box

layer = gpd.read_file("data/layers/roads.shp")

# Option A: pass the bounds tuple straight to clip (must be in the layer CRS).
bbox = (-2.60, 51.44, -2.55, 51.48)
clipped_a = gpd.clip(layer, bbox)

# Option B: build an explicit box polygon in a GeoDataFrame.
window = gpd.GeoDataFrame(geometry=[box(*bbox)], crs=layer.crs)
clipped_b = gpd.clip(layer, window)

print(len(clipped_a), len(clipped_b))   # same result, two styles

A bounding box is faster and fine for a quick rectangular crop; use the real polygon mask when the study area is an irregular shape and you must not keep anything outside it.

Example 5: Robust loop that drops empty clips and logs

The production version: skip missing CRS, skip empty clips, catch per-file errors, and print a summary you can audit.

from pathlib import Path
import geopandas as gpd

src_dir = Path("data/layers")
out_dir = Path("data/clipped")
out_dir.mkdir(parents=True, exist_ok=True)

boundary = gpd.read_file("data/boundary.gpkg").dissolve()
PATTERNS = ("*.shp", "*.geojson", "*.gpkg")

results = []
for path in sorted(p for pat in PATTERNS for p in src_dir.glob(pat)):
    try:
        layer = gpd.read_file(path)
        if layer.crs is None:
            results.append((path.name, 0, "skip: no CRS"))
            continue
        if layer.crs != boundary.crs:
            layer = layer.to_crs(boundary.crs)

        clipped = gpd.clip(layer, boundary)
        if clipped.empty:
            results.append((path.name, 0, "skip: no overlap"))
            continue

        clipped.to_file(out_dir / f"{path.stem}_clip.gpkg", driver="GPKG")
        results.append((path.name, len(clipped), "ok"))
    except Exception as exc:
        results.append((path.name, 0, f"error: {exc}"))

print(f"\nProcessed {len(results)} files")
for name, n, status in results:
    print(f"  {name:28} {n:6}  {status}")

Explanation

The whole workflow rests on three facts about gpd.clip().

First, clip is an intersection that cuts geometries at the mask edge. Unlike a spatial filter that keeps or drops whole features, gpd.clip(layer, mask) returns the geometric intersection: a polygon straddling the boundary comes back trimmed to the part inside, a line is cut where it crosses. Attributes and CRS carry through unchanged, so the output is a normal GeoDataFrame ready to write.

Second, clip requires a shared CRS. Coordinates only mean anything relative to their CRS. [400000, 100000] in a national grid and [-2.5, 51.4] in WGS84 cannot be compared until they are in the same frame, so clipping across mismatched CRSs is meaningless — GeoPandas 1.x raises rather than silently returning garbage. Reprojecting each layer to the fixed mask CRS with to_crs() is what makes the comparison valid, and it keeps every output in one consistent frame.

Third, a clip can legitimately return nothing. A layer entirely outside the study area, or one that only skirts its edge, produces an empty GeoDataFrame. That is a valid answer, not an error — which is why the loop checks clipped.empty and skips instead of writing a zero-feature file. Testing bounding-box overlap first is an optimisation on top of this: extents that miss guarantee an empty clip, so you can skip the expensive geometry work entirely.

Put together, the pattern is: read the mask once, and for each file reproject to the mask CRS, optionally reject on non-overlapping extents, clip, drop empties, and write — all inside a try/except so no single file can stop the batch.

Edge cases or notes

Invalid geometries break the clip

gpd.clip() runs a geometric overlay, and overlays fail on invalid geometry — a self-intersecting polygon or a bowtie ring will raise a GEOSException mid-clip. If a layer or the boundary itself throws during clipping, repair the geometry first with make_valid() (or buffer(0) for simple cases) before retrying. See how to fix invalid geometries in GeoPandas for the full repair recipe.

A multi-part boundary needs dissolving

If your boundary file has several rows — separate polygons, or one row per administrative unit — pass it to gpd.clip() as-is only if you want the union of all of them. Calling dissolve() first merges every row into a single mask so the clip treats the whole study area as one shape and does not accidentally clip against just the first feature.

Bounding box versus polygon mask

Clipping to a bounding box keeps everything inside a rectangle; clipping to a polygon keeps only what falls inside the actual shape. A bbox is faster and good enough for a rough crop or a preview, but for an irregular study area it will retain features in the corners of the rectangle that lie outside the real boundary. Use the polygon mask whenever "outside the area" must genuinely mean outside.

Empty results are valid, not failures

Do not treat a zero-feature clip as an error. A layer can be perfectly good and simply have nothing inside your study area. Log it as skipped and move on; writing an empty file just litters the output folder with layers that break later steps expecting geometry.

The mask CRS drives everything downstream

Whatever CRS the boundary is in becomes the CRS of every clipped output, because you reproject each layer to match it. If you need outputs in a specific CRS — metres for area calculations, or EPSG:4326 for sharing — set the boundary to that CRS once at the top with to_crs() and let the whole batch inherit it.

Very large layers and memory

gpd.read_file() loads an entire layer into memory before clipping. For a folder of large files, process one at a time inside the loop (as shown) so only one layer is resident at once. If a single layer is too big to fit, read just the boundary's extent with the bbox= argument of gpd.read_file() to pre-filter features at read time, then clip the smaller result.

Internal links

FAQ

How do I clip many layers to one boundary in Python?

Read the boundary once and dissolve it to a single mask, then loop over the folder with pathlib. For each layer, reproject it to the boundary's CRS with to_crs(), call gpd.clip(layer, boundary), skip the result if it is empty, and write the rest to GeoPackage. Wrap each file in try/except so one failure does not stop the batch.

Why does `gpd.clip()` need the same CRS for the layer and the mask?

Clipping compares coordinates directly to find the intersection, and coordinates only have meaning within their CRS. If the layer and mask are in different CRSs the comparison is meaningless, so GeoPandas 1.x raises an error instead of guessing. Reproject the layer to the boundary CRS with to_crs() before every clip.

How do I skip layers that fall outside the boundary?

Compare bounding boxes first: if the layer's total_bounds does not overlap the boundary's total_bounds, the geometries cannot intersect, so continue past it. For an exact test, use layer.intersects(mask_geom).any(). Either way, also check clipped.empty after the clip, because a layer can overlap the extent yet still leave nothing inside the real polygon.

What is the difference between clipping to a bounding box and a polygon mask?

A bounding box keeps everything inside a rectangle defined by (xmin, ymin, xmax, ymax); a polygon mask keeps only what falls inside the actual boundary shape. gpd.clip() accepts a bounds tuple or a shapely.box() for the rectangle, and a GeoDataFrame for the polygon. Use the polygon when the study area is irregular and you must exclude everything genuinely outside it.

Why does my clip return an empty GeoDataFrame?

The layer has nothing inside the boundary — either it sits entirely outside the study area, or it only touches the edge without crossing in. That is a valid result, not a bug. Check clipped.empty and skip writing the file. If you expected features, confirm the layer and mask are in the same CRS; a mismatch left uncorrected can place the layer far from the boundary.

My clip fails with a geometry error — what do I do?

gpd.clip() runs a geometric overlay that fails on invalid geometry, such as self-intersecting polygons. Repair the offending layer (or the boundary) with make_valid() before clipping, and put the clip inside a try/except so a single bad file is logged and skipped rather than aborting the run. See the guide on fixing invalid geometries for the details.

Should I clip in the boundary's CRS or the layer's CRS?

Clip in the boundary's CRS and reproject each layer to match it. This keeps every clipped output in one consistent frame, which matters if you later merge or compare them. If you need a specific output CRS, reproject the boundary to it once at the top and the whole batch inherits it.

How do I combine all the clipped layers into one file afterwards?

Clipping and merging are separate steps. Once every layer is clipped to the boundary, concatenate the results into a single layer or write them as layers of one GeoPackage. See how to merge many shapefiles into one file for the merge patterns; keep every clipped output in the same CRS first so the coordinates line up.

How to Batch Clip Many Layers to a Single Boundary in Python #

Problem statement #

Quick answer #

Step-by-step solution #

Load the boundary once #

Match the CRS before every clip #

Clip the layer to the boundary #

Skip layers that never intersect #

Drop empty clips and write output #

Handle errors per file #

Code examples #

Example 1: Clip one layer as the building block #

Example 2: The folder loop — reproject then clip #

Example 3: Skip non-intersecting layers via bounds and intersects #

Example 4: Clip to a bounding box instead of a polygon mask #

Example 5: Robust loop that drops empty clips and logs #

Explanation #

Edge cases or notes #

Invalid geometries break the clip #

A multi-part boundary needs dissolving #

Bounding box versus polygon mask #

Empty results are valid, not failures #

The mask CRS drives everything downstream #

Very large layers and memory #

Internal links #

FAQ #

How do I clip many layers to one boundary in Python? #

Why does gpd.clip() need the same CRS for the layer and the mask? #

How do I skip layers that fall outside the boundary? #

What is the difference between clipping to a bounding box and a polygon mask? #

Why does my clip return an empty GeoDataFrame? #

My clip fails with a geometry error — what do I do? #

Should I clip in the boundary's CRS or the layer's CRS? #

How do I combine all the clipped layers into one file afterwards? #