How to Clip Spatial Data in Python with GeoPandas

On this page

Problem statement

Clipping cuts one layer to the outline of another: roads to a city boundary, land use to a project area. geopandas.clip() does it in one call, and most of what goes wrong with it is silent. Layers in different coordinate systems clip to nothing with only a warning; lines cut by the boundary change geometry type; and the usual advice for speeding it up — prefilter, dissolve the mask first — repeats work clip already does.

The examples clip OpenStreetMap's roads and land use for Bristol, in longitude and latitude, to ONS boundaries in British National Grid, and measure each of those.

Quick answer

Clipping trims an input layer to a boundary, cutting features at the boundary edge and dropping everything outside.

Put both layers in one CRS, then clip by the polygon layer:

import geopandas as gpd

target = gpd.read_file("data/roads.shp")
boundary = gpd.read_file("data/city_boundary.shp")

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

clipped = gpd.clip(target, boundary)
clipped.to_file("output/roads_clipped.shp")
print(f"{len(target):,} roads in the extract -> {len(clipped):,} inside the city boundary")

9,698 roads in the extract -> 9,351 inside the city boundary

Features wholly outside go, features crossing the boundary are cut, and the attributes of the clipped layer come through unchanged.

Load, align CRSs, clip to a polygon, save

The first and third catch the empty-result case that puzzles everyone once.

Load the input spatial layers

First read the dataset you want to clip and the polygon layer that will act as the clip boundary.

import geopandas as gpd

roads = gpd.read_file("data/roads.shp")
city_boundary = gpd.read_file("data/city_boundary.shp")

Check that both layers use the same CRS

Before clipping, both layers should be in the same coordinate reference system.

print(roads.crs)
print(city_boundary.crs)

EPSG:4326
EPSG:27700

Here they differ — OpenStreetMap arrives in longitude and latitude, the ONS boundary in British National Grid — and clipping anyway does not fail:

import warnings

with warnings.catch_warnings(record=True) as caught:
    warnings.simplefilter("always")
    mismatched = gpd.clip(roads, city_boundary)
print(len(mismatched), "roads;", [str(w.message)[:60] for w in caught])

0 roads; ['CRS mismatch between the CRS of left geometries and the CRS ']

An empty result and a warning that is easy to lose in a log. The degrees and the metres do not overlap, so GeoPandas checks the two bounding boxes, finds no intersection and returns nothing. If they differ, reproject one layer to match the other:

if roads.crs != city_boundary.crs:
    city_boundary = city_boundary.to_crs(roads.crs)

Confirm the clip layer is polygon-based

The clip boundary should usually be a polygon or multipolygon layer.

You can inspect geometry types like this:

print(city_boundary.geom_type.unique().tolist())

['MultiPolygon']

Typical valid clip geometry types are:

Polygon
MultiPolygon

If your boundary layer is made of points or lines, it is not the right input for a standard vector clip workflow.

Run the clip operation with GeoPandas

Use gpd.clip() to clip the input data by the boundary layer.

roads_clipped = gpd.clip(roads, city_boundary)
print(len(roads), "->", len(roads_clipped), "|", roads_clipped.geom_type.value_counts().to_dict())

9698 -> 9351 | {'LineString': 9326, 'MultiLineString': 25}

This returns:

features fully inside the boundary unchanged
line and polygon features crossing the boundary cut to the boundary edge
point features outside the boundary removed
features fully outside removed

Save the clipped result to a new file

Write the output to a new dataset so the original files stay unchanged.

Save as a shapefile:

roads_clipped.to_file("output/roads_clipped.shp")

Or save as GeoJSON:

roads_clipped.to_file("output/roads_clipped.geojson", driver="GeoJSON")
print(sorted(p.name for p in Path("output").glob("roads_clipped.*")))

['roads_clipped.cpg', 'roads_clipped.dbf', 'roads_clipped.geojson', 'roads_clipped.prj', 'roads_clipped.shp', 'roads_clipped.shx']

Bristol roads clipped by boundaries

The data used on this page

Bristol. The layers to clip come from Geofabrik's OpenStreetMap extract for the city (ODbL), read from the .osm.pbf and written in the formats the examples use: main and residential roads as a shapefile, land-use polygons as both a shapefile and GeoJSON. The boundaries come from the ONS Open Geography Portal (Open Government Licence) in British National Grid: the city's full-resolution boundary, and its 34 wards with the ward name in a NAME column. One ward, Hotwells and Harbourside, stands in for a project area.

from pathlib import Path

import geopandas as gpd
import requests

Path("data").mkdir(exist_ok=True)
Path("output").mkdir(exist_ok=True)
PBF = Path("bristol.osm.pbf")
if not PBF.exists():
    PBF.write_bytes(requests.get("https://download.geofabrik.de/europe/united-kingdom/england/"
                                 "bristol-260101.osm.pbf", timeout=300).content)
ONS = "https://services1.arcgis.com/ESMARspQHYMw9BZ9/arcgis/rest/services/"

gpd.read_file(PBF, layer="lines", columns=["osm_id", "name", "highway"],
              where="highway IN ('motorway','trunk','primary','secondary','tertiary','residential')"
              ).to_file("data/roads.shp")
landuse = gpd.read_file(PBF, layer="multipolygons", where="landuse IS NOT NULL", columns=["osm_id", "landuse"])
landuse.to_file("data/land_use.shp")
landuse.to_file("data/parcels.geojson")
gpd.read_file(ONS + "Local_Authority_Districts_May_2024_Boundaries_UK_BFC/FeatureServer/0/query"
              "?where=LAD24CD%3D%27E06000023%27&outFields=LAD24NM&outSR=27700&f=geojson"
              ).to_file("data/city_boundary.shp")
wards = gpd.read_file(ONS + "Wards_December_2024_Boundaries_UK_BGC/FeatureServer/0/query"
                      "?where=LAD24CD%3D%27E06000023%27&outFields=WD24CD,WD24NM&outSR=27700&f=geojson"
                      ).rename(columns={"WD24NM": "NAME"})
wards.to_file("data/districts.shp")
wards[wards["NAME"] == "Hotwells and Harbourside"].to_file("data/project_area.geojson")
for f in ("roads.shp", "land_use.shp", "city_boundary.shp", "districts.shp", "project_area.geojson"):
    g = gpd.read_file(f"data/{f}")
    print(f"{f:22} {len(g):6,} features, {g.crs.to_epsg()}")

roads.shp               9,698 features, 4326
land_use.shp            5,289 features, 4326
city_boundary.shp           1 features, 27700
districts.shp              34 features, 27700
project_area.geojson        1 features, 27700

Example 1: Clip a shapefile by an administrative boundary

This example clips a roads shapefile to a city boundary shapefile.

import geopandas as gpd

roads = gpd.read_file("data/roads.shp")
city_boundary = gpd.read_file("data/city_boundary.shp")

if roads.crs != city_boundary.crs:
    city_boundary = city_boundary.to_crs(roads.crs)

roads_clipped = gpd.clip(roads, city_boundary)

roads_clipped.to_file("output/roads_within_city.shp")
print(f"{roads_clipped.to_crs(27700).length.sum() / 1000:,.0f} km of road inside Bristol")

1,091 km of road inside Bristol

This is a typical workflow when you need road features only for one municipality or study area.

Example 2: Clip GeoJSON data with a polygon layer

The same process works for GeoJSON.

import geopandas as gpd

parcels = gpd.read_file("data/parcels.geojson")
project_area = gpd.read_file("data/project_area.geojson")

if parcels.crs != project_area.crs:
    project_area = project_area.to_crs(parcels.crs)

parcels_clipped = gpd.clip(parcels, project_area)

parcels_clipped.to_file("output/parcels_project_area.geojson", driver="GeoJSON")
print(len(parcels_clipped), "land-use polygons in the project area:",
      parcels_clipped["landuse"].value_counts().head(4).to_dict())

77 land-use polygons in the project area: {'grass': 45, 'residential': 13, 'flowerbed': 5, 'retail': 4}

Example 3: Clip using a single polygon from a larger boundary dataset

Sometimes your boundary file contains many areas and you only want one of them.

import geopandas as gpd

land_use = gpd.read_file("data/land_use.shp")
districts = gpd.read_file("data/districts.shp")

selected_district = districts[districts["NAME"] == "Central"]

if land_use.crs != selected_district.crs:
    selected_district = selected_district.to_crs(land_use.crs)

land_use_clipped = gpd.clip(land_use, selected_district)

land_use_clipped.to_file("output/land_use_central_district.shp")
print(len(selected_district), "district;", len(land_use_clipped), "land-use polygons clipped to it")

1 district; 187 land-use polygons clipped to it

This is useful when clipping by one administrative unit, project area, or management zone from a larger boundary dataset.

Clip against bounding boxes, joins and overlays

The first cause accounts for most of them, and it never raises an exception.

clip is a cut, not a selection: it keeps the part of each feature that lies inside the mask, so a road crossing the city boundary loses the length outside it, and one that crosses twice comes back as a MultiLineString — 25 of Bristol's did. That is what separates it from the operations it is often confused with.

Clip vs bounding box filter

A bounding box filter selects whole features that touch a rectangle. It does not cut features to the true polygon boundary. Three ways to ask for "roads in a 300 m square of Bristol", in British National Grid:

from shapely.geometry import box

r = roads.to_crs(27700)
square = box(358_500, 174_000, 358_800, 174_300)
by_bbox = r.iloc[r.sindex.query(square)]                 # bounding boxes overlap
by_cx = r.cx[358_500:358_800, 174_000:174_300]           # geometry intersects
by_clip = gpd.clip(r, square)                            # geometry cut to the square
print(f"bbox {len(by_bbox)}, .cx {len(by_cx)}, clip {len(by_clip)}")
print(f"length: .cx {by_cx.length.sum() / 1000:.1f} km, clip {by_clip.length.sum() / 1000:.1f} km")

bbox 26, .cx 20, clip 20
length: .cx 3.3 km, clip 1.7 km

.cx is stricter than a bounding-box test — six roads whose bounding boxes overlap the square never enter it, and .cx leaves them out — but it keeps the whole of each road, so its total length includes everything outside. clip has the same features and only the parts inside.

Clip vs spatial join

A spatial join adds attributes based on spatial relationships. It does not trim geometry.

Clip vs overlay intersection

overlay(..., how="intersection") is more general and combines attributes from both layers. clip() is simpler when you only want to trim one layer by a boundary.

Invalid masks, multipart masks and speed

Invalid geometries can break the clip process

Self-intersections or other invalid geometry problems can cause clipping errors.

Check validity:

print(boundary.is_valid.value_counts().to_dict(), "|", districts.is_valid.value_counts().to_dict(),
      "|", land_use.is_valid.value_counts().to_dict())

{True: 1} | {True: 34} | {True: 5287, False: 2}

The ONS boundaries are valid; the raw OSM land use is not entirely. Repair with make_valid, preferably in structure mode, rather than the old buffer(0) trick, which quietly drops parts of some shapes:

land_use["geometry"] = land_use.geometry.make_valid(method="structure", keep_collapsed=False)
print(land_use.is_valid.all())

True

Multi-part boundaries may affect output size

A clip layer made of many polygons does not cut features once per polygon. gpd.clip unions the mask before it cuts, so clipping by Bristol's 34 wards and by the wards dissolved into one shape gives the same answer:

import time

t0 = time.perf_counter(); by_wards = gpd.clip(target.to_crs(27700), districts); t1 = time.perf_counter()
by_dissolved = gpd.clip(target.to_crs(27700), districts.dissolve()); t2 = time.perf_counter()
print(len(by_wards), len(by_dissolved), "features; same geometry:",
      bool(by_wards.sort_index().geom_equals(by_dissolved.sort_index()).all()),
      f"| {t1 - t0:.2f} s vs {t2 - t1:.2f} s")

9348 9348 features; same geometry: True | 0.18 s vs 0.18 s

Fragmented output comes from the boundary's shape — a coastline, an island, a hole — not from how many rows it has. To cut features at every internal boundary and keep the pieces separate, that is overlay(..., how="intersection"), not clip.

Why two familiar speed-ups change nothing, and why .cx is not a clip.

Attribute data is preserved from the input layer

The clipped output keeps attributes from the layer being clipped. It does not automatically add fields from the boundary layer.

If you need attributes from both layers, use an overlay or spatial join instead.

Performance considerations for large datasets

Clipping very large layers can be slow, especially with complex polygons.

Practical ways to improve performance:

keep only necessary columns before clipping
simplify a very detailed boundary if the result allows it — the cost is in the cutting
use a projected CRS for regional workflows

Two common tips do nothing. clip already queries the input's spatial index for candidates before cutting, so prefiltering by extent repeats work it does itself, and it already unions a many-part boundary, so dissolving first saves nothing.

Internal links

FAQ

How do I clip a shapefile in Python with GeoPandas?

Read both layers with gpd.read_file(), make sure they share the same CRS, then run gpd.clip().

clipped = gpd.clip(input_layer, boundary_layer)

Save the result with to_file().

Why does `geopandas.clip()` return an empty result?

The most common reasons are:

the layers use different CRS values
the datasets do not actually overlap
the clip layer is not polygon-based
one of the layers has invalid geometry

Check CRS, geometry type, and validity first.

Do both layers need the same CRS before clipping?

Yes. In practice, they should be in the same CRS before running gpd.clip(). If not, results may be empty or spatially wrong.

What is the difference between clip and overlay in GeoPandas?

clip() trims one layer by a polygon boundary and keeps attributes from the input layer. overlay() performs more general spatial operations like intersection and can combine attributes from both layers.

Can I clip a GeoDataFrame using a bounding box instead of a polygon?

Yes. Pass a tuple of (minx, miny, maxx, maxy) in the layer's CRS to gpd.clip(), for example gpd.clip(roads, (minx, miny, maxx, maxy)). This trims to a rectangle rather than to a true polygon edge.

How is `gpd.clip()` different from filtering with `.cx`?

The .cx indexer selects whole features whose geometry intersects a coordinate rectangle — the geometry itself, not its bounding box — so it never cuts anything. gpd.clip() trims crossing features at the boundary edge, which is the right choice when you need geometry that stops at the boundary.

Does `gpd.clip()` keep the same index as the input layer?

Yes, clipped features retain the original index of the input layer, which lets you join attributes back afterward. Pass keep_geom_type=True if you want to drop geometry types that differ from the input — the lines and points produced where a polygon only touches the boundary.

Why are some of my polygons sliced into many small pieces after clipping?

Because the boundary's shape cuts them: an inlet, an island or a hole in the mask splits any feature that crosses it into parts. It is not the number of boundary rows — clip unions a many-part mask before cutting, so dissolving it first changes nothing.