How to Select Features by Location in GeoPandas

On this page

Problem statement

Selecting by location keeps the features of one layer that stand in some relation to another — inside, touching, crossing, within a distance. In practice:

selecting points inside administrative boundaries
selecting parcels within a study area
selecting roads that cross water
filtering buildings within a buffer distance of a river

In GeoPandas, this means comparing geometries in one GeoDataFrame to geometries in another and keeping only the matching features. The key requirements are:

both layers must use the same CRS
you must choose the correct spatial predicate, such as within, intersects, or contains

Every example below runs on real layers for the county of Rutland: its 15 wards from the ONS, 21,253 Land Registry parcels, and pubs, school grounds, roads, water and rivers from OpenStreetMap.

Quick answer

The usual GeoPandas workflow for selecting features by location is:

read both layers
make sure both use the same CRS
choose the correct spatial predicate
use geopandas.sjoin() when matching one full layer to another

Example:

pubs = gpd.read_file(f"{OSM}gis_osm_pois_free_1.shp", where="fclass = 'pub'")
wards = gpd.read_file("data/wards.geojson")

pubs = pubs.to_crs(wards.crs)

selected = gpd.sjoin(pubs, wards, predicate="within", how="inner")
print(len(pubs), "pubs in the extract,", len(selected), "inside a Rutland ward")
print(selected["WD24NM"].value_counts().head(3).to_dict())

33 pubs in the extract, 32 inside a Rutland ward
{'Braunston & Martinsthorpe': 5, 'Uppingham': 4, 'Ryhall & Casterton': 4}

This keeps only points that fall inside the polygons — and shows that Geofabrik's county extract reaches a little beyond the county.

Selecting Rutland features by location

Real data sits on boundaries constantly — pick the predicate deliberately.

Load the layers you want to compare

Read source files and check geometry types

Start by loading the datasets, and confirm what kinds of features you are working with:

schools = gpd.read_file(f"{OSM}gis_osm_pois_a_free_1.shp", where="fclass = 'school'")
parcels = gpd.read_file(PARCELS, columns=["INSPIREID"], DOWNLOAD_SCHEMA="NO")
for name, layer in (("pubs", pubs), ("schools", schools), ("wards", wards), ("parcels", parcels)):
    print(f"{name:8} {len(layer):>6} {layer.geom_type.value_counts().to_dict()}")

pubs         33 {'Point': 33}
schools      26 {'Polygon': 25, 'MultiPolygon': 1}
wards        15 {'Polygon': 15}
parcels   21253 {'Polygon': 21253}

Points against polygons, polygons against polygons — the predicates behave differently for each, as the next sections show.

Check for invalid geometries

Invalid geometries can cause spatial operations to fail or return unexpected results:

for name, layer in (("schools", schools), ("wards", wards), ("parcels", parcels)):
    print(f"{name:8} invalid: {int((~layer.is_valid).sum())}")

schools  invalid: 0
wards    invalid: 0
parcels  invalid: 0

None here. If there were, repair with make_valid() and check the result's geometry types; buffer(0) also repairs some polygons but can silently drop parts of a self-intersecting one.

Make sure both layers use the same CRS

Check the CRS of each layer

print("schools:", schools.crs.to_epsg(), "| wards:", wards.crs.to_epsg(), "| parcels:", parcels.crs.to_epsg())

schools: 4326 | wards: 27700 | parcels: 27700

Why a CRS mismatch causes wrong results

OpenStreetMap layers arrive in longitude and latitude; the ONS and Land Registry layers are in British National Grid metres. Compare them unconverted and nothing matches:

with warnings.catch_warnings(record=True) as caught:
    warnings.simplefilter("always")
    unmatched = gpd.sjoin(schools, wards, predicate="intersects")
print(len(unmatched), "matches |", str(caught[0].message).split(".")[0] if caught else "no warning")

schools = schools.to_crs(wards.crs)
print(len(gpd.sjoin(schools, wards, predicate="intersects")), "matches after schools.to_crs(wards.crs)")

0 matches | CRS mismatch between the CRS of left geometries and the CRS of right geometries
28 matches after schools.to_crs(wards.crs)

GeoPandas warned about the mismatch, and still returned an empty result rather than an error. Reproject one layer to the other's CRS — to_crs(other.crs) — before every spatial comparison.

Select features by spatial relationship

There are two main patterns:

compare a layer to one merged geometry
compare every feature in one layer to every feature in another with sjoin()

Use the first pattern for simple filtering against a single area. Use sjoin() for feature-to-feature matching between full layers.

`within`, `intersects` and `touches` on the same boundary

The predicate decides what happens at the edge, and real data sits on edges constantly. Here are Rutland's parcels against the Oakham North East ward:

ward = wards.loc[wards["WD24NM"] == "Oakham North East", "geometry"].iloc[0]
counts = {predicate: int(getattr(parcels, predicate)(ward).sum())
          for predicate in ("within", "intersects", "touches", "overlaps")}
print(counts)
print(f"crossing the ward boundary: {counts['intersects'] - counts['within'] - counts['touches']}")

{'within': 1570, 'intersects': 1609, 'touches': 0, 'overlaps': 39}
crossing the ward boundary: 39

within keeps only parcels entirely inside the ward; intersects adds the 39 that cross its boundary. touches found none at all — not because no parcel runs along the ward's edge, but because the ONS and the Land Registry draw their boundaries independently, so a parcel on the edge overlaps the ward by a sliver instead of sharing the line exactly. On real data, "touches" is almost always "overlaps a little". Pick the predicate deliberately: "parcels in the ward" means within for a count of whole parcels, and intersects for every parcel that the ward affects.

The predicate comparison above, counted.

Use `intersects` to select overlapping or touching features

Roads that cross open water, for example:

roads = gpd.read_file(f"{OSM}gis_osm_roads_free_1.shp").to_crs(wards.crs)
water = gpd.read_file(f"{OSM}gis_osm_water_a_free_1.shp").to_crs(wards.crs)

water_geom = water.union_all()
roads_on_water = roads[roads.intersects(water_geom)]
print(len(roads_on_water), "road segments touch or cross mapped water:", roads_on_water["fclass"].value_counts().to_dict())
crossing = roads_on_water[roads_on_water.crosses(water_geom)]
print(len(crossing), "of them cross it rather than only touching")

24 road segments touch or cross mapped water: {'service': 13, 'footway': 8, 'track': 2, 'track_grade1': 1}
11 of them cross it rather than only touching

None of the 24 is a public road — they are service roads, footways and tracks at the edges of Rutland Water and smaller ponds — and only 11 cross water rather than touching it. intersects counts a line that only touches a water polygon's edge as well as one that runs across it; crosses keeps only the second kind. Which one you want is a question about the data, not about GeoPandas.

Use `contains` when the reference geometry must contain another

contains is the reverse of within. Which ward contains the site of Uppingham School?

school = schools.loc[schools["name"] == "Uppingham School", "geometry"]
print("wards containing the whole site:", wards.loc[wards.contains(school.iloc[0]), "WD24NM"].tolist())
print("wards containing a point inside it:",
      wards.loc[wards.contains(school.representative_point().iloc[0]), "WD24NM"].tolist())

wards containing the whole site: ['Uppingham']
wards containing a point inside it: ['Uppingham']

Choose the right predicate

Use:

within for features completely inside another geometry
contains for geometries that completely contain another
intersects for any overlap or touch
touches when only shared boundaries matter
crosses for line-crossing cases

Use spatial join for layer-to-layer selection

When `geopandas.sjoin()` is the simplest option

If you want to compare every feature in one layer against every feature in another and keep matching rows, sjoin() is usually the cleanest approach — and it brings the other layer's attributes along:

joined = gpd.sjoin(schools, wards, how="inner", predicate="intersects")
print(len(schools), "school sites ->", len(joined), "rows")
repeated = joined[joined.index.duplicated(keep=False)]
print(repeated.groupby("name")["WD24NM"].apply(list).to_dict())

26 school sites -> 28 rows
{'Catmose College': ['Oakham North West', 'Barleythorpe'], 'Oakham School': ['Oakham South', 'Oakham North East']}

Two school sites — Catmose College and Oakham School — produced two rows each, because each intersects two wards. An inner join keeps a row per match, so counting rows is not counting schools. Deduplicate on the original index — or assign each site to one ward by its representative point.

Select features that are NOT inside a boundary

Run a left join, then keep the rows that found no match:

joined = gpd.sjoin(pubs, wards, how="left", predicate="within")
outside = joined[joined["index_right"].isna()]
print(len(outside), "pubs outside every Rutland ward, e.g.", outside["name"].dropna().head(3).tolist())

1 pubs outside every Rutland ward, e.g. ['The Castle Inn']

Wards, parcels and distance-based selections

The data used on this page

Rutland's wards from the ONS December 2024 boundaries (OGL); HM Land Registry's INSPIRE parcels for the county (© Crown copyright and database rights; OGL); and pubs, school grounds, roads, water and waterways from Geofabrik's OpenStreetMap extract of 1 January 2026 (© OpenStreetMap contributors, ODbL).

import warnings
from pathlib import Path

import geopandas as gpd
import requests

Path("data").mkdir(exist_ok=True)
Path("output").mkdir(exist_ok=True)
downloads = {
    "data/rutland.shp.zip": "https://download.geofabrik.de/europe/united-kingdom/england/rutland-260101-free.shp.zip",
    "data/Rutland.zip": "https://use-land-property-data.service.gov.uk/datasets/inspire/download/Rutland.zip",
}
for target, url in downloads.items():
    if not Path(target).exists():
        Path(target).write_bytes(requests.get(url, timeout=300).content)
if not Path("data/wards.geojson").exists():
    Path("data/wards.geojson").write_bytes(requests.get(
        "https://services1.arcgis.com/ESMARspQHYMw9BZ9/arcgis/rest/services/"
        "Wards_December_2024_Boundaries_UK_BFC/FeatureServer/0/query",
        params={"where": "LAD24CD='E06000017'", "outFields": "WD24NM", "outSR": 27700, "f": "geojson"},
        timeout=120).content)
OSM = "/vsizip/data/rutland.shp.zip/"
PARCELS = "/vsizip/data/Rutland.zip/Land_Registry_Cadastral_Parcels.gml"
print(sorted(p.name for p in Path("data").iterdir()))

['Rutland.zip', 'rutland.shp.zip', 'wards.geojson']

Example 1: Select points within administrative boundaries, one merged area

county = wards.union_all()
inside = pubs[pubs.within(county)]
print(f"pubs inside Rutland: {len(inside)} of {len(pubs)}")

pubs inside Rutland: 32 of 33

Merging the wards into one geometry answers "inside the county?" without producing a row per ward.

Example 2: Select parcels inside a study area and save them

selected_parcels = parcels[parcels.within(ward)]
selected_parcels.to_file("output/oakham_ne_parcels.geojson", driver="GeoJSON")
selected_parcels.to_file("output/oakham_ne_parcels.shp")
back = gpd.read_file("output/oakham_ne_parcels.geojson")
print(len(selected_parcels), "parcels saved;", "GeoJSON read back in EPSG:", back.crs.to_epsg())

1570 parcels saved; GeoJSON read back in EPSG: 27700

A GeoJSON written by GDAL from British National Grid keeps its CRS for GDAL, but the GeoJSON standard expects longitude and latitude; convert with to_crs(4326) before handing it to a web map.

Example 3: Distance-based selection

For distance-based selection, work in a projected CRS so the distance is in metres rather than degrees:

rivers = gpd.read_file(f"{OSM}gis_osm_waterways_free_1.shp", where="fclass = 'river'").to_crs(wards.crs)
buildings = gpd.read_file(f"{OSM}gis_osm_buildings_a_free_1.shp").to_crs(wards.crs)

river_buffer = rivers.buffer(100).union_all()
nearby_buildings = buildings[buildings.intersects(river_buffer)]
print(f"buildings within 100 m of a river: {len(nearby_buildings)} of {len(buildings)}")

near = gpd.sjoin(buildings, rivers[["name", "geometry"]], predicate="dwithin", distance=100)
print("same answer with predicate='dwithin':", near.index.nunique())

buildings within 100 m of a river: 189 of 5513
same answer with predicate='dwithin': 189

predicate="dwithin" with a distance does the same selection without building a buffer — and without the buffer's tiny polygon approximation of a circle.

Boolean predicates against a spatial join

There are two practical ways to select features by location in GeoPandas.

The first two account for nearly every surprising selection result.

The first is to use boolean geometry methods such as within() and intersects() directly on a GeoSeries. This is a good fit when you want to compare one layer against a single geometry, such as a study area created by merging polygons with union_all().

The second is to use gpd.sjoin() when you want feature-to-feature comparison between two full layers. This is usually the better option when each feature in one dataset may match one or more features in another — and it is where duplicated rows come from.

The main decision is:

use geometry methods for one geometry or one merged area
use sjoin() for full layer-to-layer matching

Predicate choice also matters:

within excludes anything that reaches the boundary from outside
intersects includes touching and overlapping cases
contains depends on geometry direction

This pattern is useful in automation because it can be reused in scripts: load layers, align CRS, run the spatial filter, and export the result.

Empty results from CRS mismatches, and duplicates

CRS mismatch is the most common reason for empty results: GeoPandas warns and returns nothing. Always compare gdf1.crs and gdf2.crs.
Invalid geometries can break spatial matching or give unreliable output. Check is_valid before joining or filtering.
Boundary-touching features may not be returned by within. Use intersects if touching the boundary should count.
Duplicate matches occur with sjoin() when one feature matches several polygons, as a school site spanning two wards did above.
Large datasets are usually faster with sjoin() than manual loops, but performance still depends on geometry complexity and data size.
Null geometries should be removed or handled before spatial operations.

Internal links

For the broader workflow, see Spatial joins in GeoPandas.

Related tasks:

For exporting filtered results, see How to Export GeoJSON in Python with GeoPandas.

If your query returns nothing, see Why GeoPandas Spatial Join Returns No Results.

FAQ

How do I select points inside polygons in GeoPandas?

Use within() against a polygon geometry or use gpd.sjoin() with predicate="within".

selected = gpd.sjoin(points, polygons, how="inner", predicate="within")

Why does my spatial selection return no rows?

Common causes are:

CRS mismatch
invalid geometries
wrong predicate choice
null or empty geometries

Can I export selected features to a new file?

Yes. After filtering, save the result with to_file().

selected.to_file("output/selected_features.geojson", driver="GeoJSON")

Is `sjoin()` better than geometry methods?

Use geometry methods for simple comparisons against one geometry. Use sjoin() when comparing two layers feature by feature.

How do I select features within a distance of another layer?

Buffer the reference layer in a projected CRS, then test intersection. For example, reproject to metres, build target.buffer(500).union_all(), and keep features where gdf.intersects(buffered) is true. The buffer distance is in the CRS units, so use a metric CRS such as EPSG:32633.

How do I keep the polygon attributes when selecting points with `sjoin()`?

gpd.sjoin() already carries the right-layer attributes into the result, so a within join of points to polygons returns each point row with its containing polygon's columns attached. Drop the leftover index_right column afterward if you do not need it.

How do I select features that are NOT inside a boundary?

Run a left join with predicate="within", then keep the unmatched rows where the joined index is null:

joined = gpd.sjoin(points, area, how="left", predicate="within")
outside = joined[joined["index_right"].isna()]

Does the order of layers matter in `gpd.sjoin()`?

Yes. The first argument is the left layer whose rows and geometry are preserved, and the predicate is read left-to-right, so sjoin(points, polygons, predicate="within") tests whether each point is within a polygon. Swapping the layers changes both the kept geometry and the meaning of the predicate.