Performing Left Joins with GeoPandas sjoin

A left spatial join attaches attributes from a second layer to a primary layer while guaranteeing that every primary geometry survives — unmatched rows keep their geometry and receive NaN for the appended columns. This guide is for anyone enriching a point or polygon layer in GeoPandas who cannot afford to silently drop records. It sits under Spatial Joins & Merging in Geospatial Data Ingestion & Processing Workflows.

Why This Approach / What Goes Wrong

The default gpd.sjoin is an inner join: any left geometry that matches no right feature is dropped from the output. For enrichment that is almost never what you want — you asked "which district contains each sensor?", not "delete the sensors outside every district." Passing how="left" keeps the full left layer and fills the right-side columns with NaN where no spatial match exists, so the count of your records is preserved and the misses are visible rather than lost.

Two things silently defeat the join. First, the coordinate reference system. sjoin evaluates its predicate on raw coordinates, so if the two layers disagree on their coordinate reference system every predicate is false and you get an all-NaN result — a mismatch you fix by aligning projections, ideally validated against the definitions in Coordinate Systems with PyProj. Second, the predicate itself: intersects (the default) treats a point lying exactly on a polygon boundary as a match to both neighbouring polygons, duplicating the left row. The predicate is a Shapely geometry operation evaluated pairwise against candidates from an R-tree index, so choosing within or contains when you mean strict containment removes that boundary ambiguity. Getting the join type right is a separate decision from getting the mechanism right — see Spatial Join vs Attribute Join in GeoPandas if you are unsure whether a spatial join is even the correct tool.

Left spatial join data flow with GeoPandas sjoin A sensors point layer and a districts polygon layer, both in EPSG:4326, feed gpd.sjoin with how="left" and predicate="within". The output table keeps all three sensor rows: sensor 1 maps to Riverside, sensor 2 to Hillcrest, and sensor 3, whose point falls in no polygon, is retained with NaN in the district column. how="left" keeps every left row — misses become NaN sensors (points) 123 EPSG:4326 districts (polygons) 101 102 EPSG:4326 gpd.sjoin( … ) how="left" predicate="within" sensor_id district_name 1Riverside 2Hillcrest 3 NaN point in no district 3 rows in → 3 rows out · unmatched columns = NaN
Left spatial join data flow with GeoPandas sjoin

Prerequisites

conda install -c conda-forge "geopandas=0.14.*" "shapely=2.0.*" "pandas=2.0.*"

Three version differences change the answer on this page. GeoPandas 0.13 dropped the optional pygeos/rtree index backends, so from that release on the index is always Shapely's STRtree and there is no backend flag left to tune. GeoPandas 1.0 renamed the match column: it is index_right in 0.x, but when the right frame carries a named index that name is used instead, so hard-coding the string "index_right" in downstream code breaks on data you did not control. And the dwithin predicate — "within a fixed distance of", without materializing a buffer — only exists when Shapely 2.0 is linked against GEOS 3.10 or newer; check shapely.geos_version before writing it into a pipeline that has to run on someone else's machine.

Step-by-Step Implementation

1. Build (or load) the two layers. Here a point layer of air-quality sensors will inherit the administrative districts polygon each one falls inside. Real pipelines would gpd.read_file(...); synthetic geometries keep the block runnable in isolation.

import geopandas as gpd
from shapely.geometry import Point, Polygon

sensors = gpd.GeoDataFrame(
    {"sensor_id": [1, 2, 3], "reading": [18.4, 22.1, 9.7]},
    geometry=[Point(0, 0), Point(5, 5), Point(10, 10)],
    crs="EPSG:4326",
)

districts = gpd.GeoDataFrame(
    {"district_id": [101, 102], "district_name": ["Riverside", "Hillcrest"]},
    geometry=[
        Polygon([(-1, -1), (1, -1), (1, 1), (-1, 1)]),
        Polygon([(4, 4), (6, 4), (6, 6), (4, 6)]),
    ],
    crs="EPSG:4326",
)

2. Align the CRS before the join — never rely on the join to warn you. Recent GeoPandas raises on a mismatch, but calling .to_crs() explicitly is the safe habit.

if sensors.crs != districts.crs:
    districts = districts.to_crs(sensors.crs)

3. Perform the left join and specify the predicate explicitly. how="left" keeps every sensor; naming predicate="within" makes point-in-polygon containment unambiguous rather than trusting the intersects default.

Choosing the how and predicate arguments for a spatial join A three-question decision staircase. First question: must every left row survive the join? Answering no leaves the default how equals inner, which drops unmatched left features silently. Answering yes moves to how equals left and the second question: do the geometries overlap at all? If they do not, use sjoin_nearest with a max_distance on a metric CRS. If they do, the third question asks whether a point on a shared border should match both polygons. Yes gives predicate equals intersects and duplicate left rows to collapse; no gives predicate equals within, strict containment with one row per left feature. Choosing how= and predicate= before you run the join Must every left row survive the join? Do the two geometries overlap at all? Should a point on a shared border match both polygons? how="inner" (the default) left features with no match are dropped silently sjoin_nearest(..., max_distance=500) on a metric CRS, so 500 means 500 metres, not degrees predicate="intersects" border points match both — expect duplicate left rows predicate="within" strict containment · one row per left feature no no yes yes · how="left" yes no Both defaults — inner and intersects — are the wrong end of every branch when the goal is enrichment.
Two independent decisions live in one call: how controls which left rows survive, predicate controls what counts as a match.
sensors_in_district = gpd.sjoin(
    sensors,
    districts[["district_id", "district_name", "geometry"]],
    how="left",
    predicate="within",
)

print(sensors_in_district[["sensor_id", "reading", "district_id", "district_name"]])
# Expected output:
#    sensor_id  reading  district_id district_name
# 0          1     18.4        101.0     Riverside
# 1          2     22.1        102.0     Hillcrest
# 2          3      9.7          NaN           NaN   <- Point(10, 10) is in no district

The third sensor matched no polygon: its row is retained, and the right-side columns are NaN. GeoPandas also adds an index_right column carrying the matched right-layer index, which doubles as a null-check for unmatched rows.

4. Repair the dtypes the NaN fill silently changed. Look closely at the printed frame: district_id came back as 101.0, not 101. NumPy's int64 has no null value, so the moment a left join introduces one missing row pandas promotes the whole column to float64. That promotion is invisible until an id round-trips through str() and arrives at a database as "101.0", or until two ids beyond 2⁵³ collide because float precision ran out. Cast the key columns to a nullable integer dtype immediately after the join, while you still know which ones were integers.

# NaN forces int64 -> float64. Int64 (capital I) is pandas' nullable integer.
sensors_in_district["district_id"] = sensors_in_district["district_id"].astype("Int64")

print(sensors_in_district["district_id"].tolist())
# [101, 102, <NA>]      not [101.0, 102.0, nan]

5. Collapse one-to-many matches deterministically. A left join is only row-preserving when each left geometry matches at most one right feature. Overlapping zoning polygons, a point on a shared border under intersects, or a polygon-to-polygon join all break that assumption and repeat the left row once per match. drop_duplicates fixes the row count but picks an arbitrary winner; when the layers are polygons the defensible rule is largest shared area, which needs the actual intersection rather than the join output.

import geopandas as gpd

# parcels and flood_zones are both polygons in EPSG:25832 (metres)
overlaps = gpd.overlay(parcels, flood_zones, how="intersection", keep_geom_type=True)
overlaps["shared_m2"] = overlaps.area

# One winner per parcel: the zone it shares the most ground with
winners = (
    overlaps.sort_values("shared_m2")
    .groupby("parcel_id", as_index=False)
    .tail(1)[["parcel_id", "zone_code", "shared_m2"]]
)

# Merge the winner back onto the FULL parcel layer to keep the left cardinality
parcels_zoned = parcels.merge(winners, on="parcel_id", how="left")
assert len(parcels_zoned) == len(parcels)

overlay drops parcels that touch no zone, which is exactly why the winners are merged back onto parcels rather than used directly — the merge restores the unmatched rows the left join promised to keep. Use keep_geom_type=True so a parcel that merely shares an edge with a zone contributes a zero-area line to the intersection instead of a spurious sliver polygon.

6. Recover the matched right geometry when you need it. sjoin returns the left geometry column and discards the right one; the joined frame knows which district a sensor fell in but not where that district is. index_right is the handle for getting it back — it holds the right frame's index label, so a plain positional lookup reattaches the polygon under a second, explicitly named geometry column.

matched = sensors_in_district["index_right"].dropna().astype("Int64")

sensors_in_district.loc[matched.index, "district_geom"] = (
    districts.geometry.loc[matched].to_numpy()
)

# Distance from each sensor to its own district boundary, in CRS units
edge_gap = gpd.GeoSeries(sensors_in_district["district_geom"], crs=districts.crs)
sensors_in_district["edge_m"] = sensors_in_district.geometry.distance(edge_gap.boundary)

Only one geometry column can be active at a time. sensors_in_district.geometry still refers to the sensor points, and set_geometry("district_geom") switches the frame over if you would rather export the polygons — a distinction that matters the moment you write the result to a file, because most drivers persist a single geometry column.

Verification

A left join must never change the number of left rows, and the misses should be exactly the geometries outside every polygon.

Row counts for the same three sensors under three join settings Three cards compare the output of the same three-sensor join. With how equals inner and predicate within, only two rows come back and sensor 3 is dropped without a message. With how equals left and predicate within, all three rows survive and sensor 3 carries NaN. With how equals left and predicate intersects, sensor 2 sits on a shared border, matches two polygons and is emitted twice, so four rows come back from three inputs. Only an explicit assertion on the row count catches both the silent drop and the silent duplicate. Three sensors in — two, three, or four rows out how="inner" predicate="within" how="left" predicate="within" how="left" predicate="intersects" sensor 1 · Riverside sensor 2 · Hillcrest sensor 3 · dropped, no warning 2 rows out sensor 1 · Riverside sensor 2 · Hillcrest sensor 3 · NaN, kept and visible 3 rows out · count preserved sensor 1 · Riverside sensor 2 · Hillcrest sensor 2 · second match on border sensor 3 · NaN 4 rows out assert len(joined) == len(sensors) catches both the silent drop and the silent duplicate.
The row count is the fastest signal a left join gives you: below the input count means dropped features, above it means a one-to-many expansion.
# Row count is preserved — the defining property of how="left"
assert len(sensors_in_district) == len(sensors), (
    "Row count changed — a point matched multiple polygons (border/overlap) "
    "or how != 'left'"
)

# Unmatched rows are visible, not dropped
unmatched = sensors_in_district["index_right"].isna().sum()
print(f"Sensors outside every district: {unmatched}")   # Sensors outside every district: 1

# A wholesale all-NaN result is the classic CRS-mismatch signature
matched_share = sensors_in_district["district_id"].notna().mean()
assert matched_share > 0, "Zero matches — inputs are almost certainly in different CRSs"

Row count alone is a weak test, because a dropped row and an extra row cancel out: two sensors silently duplicated on a border while two others fell outside every district still totals three. Assert on the index instead. A left join preserves the left index verbatim, so comparing the index sets catches both directions at once, and comparing the index for uniqueness catches the duplication independently of the count.

import pandas as pd

# Identity, not just cardinality: the same labels come back, in the same order
pd.testing.assert_index_equal(sensors_in_district.index, sensors.index)

# Uniqueness is the real one-to-one guarantee
assert sensors_in_district.index.is_unique, "a left feature matched several polygons"

# The match rate is the number worth logging, not asserting — it drifts with data
rate = sensors_in_district["district_id"].notna().mean()
print(f"matched {rate:.1%} of sensors")   # matched 66.7%

A useful habit in a scheduled pipeline is to assert a floor on the match rate rather than a fixed value — assert rate > 0.95 on a layer that historically matches 99% turns a bad upstream delivery into a failed run instead of a quietly emptier map. Pin the floor to what the data has actually done, and treat a sudden drop as a CRS or extent problem until proven otherwise.

Edge Cases & Debugging

Frequently Asked Questions

Should I put the bigger layer on the left or the right? Correctness decides this, not performance: how="left" protects whichever layer is on the left, so the layer whose rows must all survive goes there. The cost asymmetry is real but small — GeoPandas builds the R-tree over the right layer and queries it once per left geometry, so a tall left layer against a short right layer is the cheaper arrangement, and it is usually also the semantically correct one (many sensors, few districts). If you genuinely need every right row preserved instead, use how="right" rather than swapping the frames, because swapping also inverts which geometry column the result carries.

Why does how="left" return fewer rows than the left layer on some GeoPandas versions? It does not, and if it appears to, the row loss happened before the join. The usual culprits are a dropna() earlier in the chain, an explode() that changed the index, or null geometries in the left layer — a None geometry participates in no predicate and, depending on version, is either kept with NaN matches or raises. Assert sensors.geometry.notna().all() upstream so the question never arises.

Is sjoin(how="left") the same as buffering and intersecting? No, and the difference is both semantic and expensive. A buffer creates new geometry — an approximated circle with a fixed vertex count — and then joins against that approximation, so the result depends on the buffer's resolution and inflates memory by storing a polygon per left feature. When the intent is "within 500 m of", predicate="dwithin" with a distance argument evaluates the true distance with no intermediate geometry, and sjoin_nearest(..., max_distance=500) answers the related "closest one within 500 m" question. Reach for a buffer only when you actually want the buffered polygon as an output.

When does a left join stop scaling, and what replaces it? The output, not the input, sets the ceiling: sjoin materializes one row per matching pair, so a left layer of five million points against densely overlapping polygons can produce far more rows than either input and exhaust memory during concatenation rather than during the index query. Watch the match multiplicity on a sample first (joined.index.value_counts().max()). Below a few million left rows a single-process join is fine; past that, chunk the left layer and join each chunk against the same right layer, or move the whole operation to Parallel Spatial Joins with Dask-GeoPandas, which partitions both sides spatially so each worker compares only nearby features.

How do I keep the unmatched rows and still learn why they were unmatched? Separate "outside every polygon" from "matched nothing because the data is broken". After the join, take the unmatched subset and test it against the right layer's total footprint: misses.within(districts.union_all()) returns False for genuinely exterior features and True for ones that sit inside the coverage but fell through a gap between polygons — a topology defect to fix upstream via Topology Validation & Repair rather than a join setting to tweak.

Related