Spatial Joins: Combining Data by Location

The previous lesson, overlay operations, took two layers and cut their actual shapes into a new geometry. An intersection kept only the ground two layers share. A union merged every shape from both. A clip trimmed one layer down to the outline of another. In every case the output was new geometry, built from the input shapes themselves. A spatial join does something different. It never moves a single coordinate. Every geometry stays exactly as drawn, and instead the other layer's attributes get copied over, matched to whichever features actually sit near, inside, or against them.

The key question versus the location question

The foundations module covered the attribute join: attaching one table to another by a column both tables share, a FIPS code, a customer number, a parcel ID. Two rows join because the values in that column are identical strings. Location never enters the decision. A county's population figure lands on the right polygon because both tables happen to carry the same five digit code, not because anyone checked where the polygon sits on the Earth.

A spatial join answers the same basic question, whose information belongs with which feature, but drops the shared column entirely. It tests the geometries themselves. Instead of asking whether this row's code equals that row's code, it asks whether this shape actually touches, contains, or sits near that shape. No key is required, and no shared column is required either. A city's crime reports and a state's neighborhood boundaries can come from two agencies that have never coordinated an ID scheme in their lives, and a spatial join still lines them up correctly, because it only cares where each shape is.

That is the whole reason this operation gets its own name instead of just being another kind of merge. A join by key can only ever be as good as the ID system behind it. A join by location needs no ID system at all.

Point in polygon: tagging a report with its neighborhood

The clearest case for a spatial join is a point that needs the properties of whatever polygon it falls inside. A city's crime reports arrive as points: a latitude, a longitude, an incident type, nothing about which neighborhood the crime happened in. A layer of neighborhood polygons exists separately, each one carrying a name and a handful of demographic fields. Nothing links the two tables by key, because the reports were never tagged with a neighborhood ID in the first place. The only fact available is where each point sits, and that is exactly enough.

import geopandas as gpd

crimes = gpd.read_file("crime_reports.gpkg")
neighborhoods = gpd.read_file("neighborhoods.gpkg")

tagged = gpd.sjoin(crimes, neighborhoods, predicate="within", how="left")
print(tagged[["incident_id", "neighborhood_name"]].head(3))

gpd.sjoin is GeoPandas' spatial join function, and it takes the same left and right shape as the attribute join's .merge(), a left table and a right table, but instead of an on= column it takes a predicate= argument. predicate="within" tests each crime point against every neighborhood polygon and keeps the match where the point sits inside that polygon's interior. how="left" plays the same role it played in the attribute join lesson: every crime report survives the join even if a stray point falls outside all mapped neighborhoods, with its neighborhood fields left null instead of the whole row disappearing.

Choosing the predicate that fits the question

The predicate argument is doing the real matching work, the same way intersects, within, and contains did back in the lesson on spatial relationships. within tests whether each point sits inside a polygon's interior, so a report that lands exactly on a shared boundary line can slip past it. intersects is looser, catching any point that touches a polygon anywhere, including its edge, which is the safer default when a report should never be silently dropped over a boundary technicality. contains reads the identical relationship from the other direction and is the natural choice when the polygon layer sits on the left side of the join rather than the point layer. All three of these still require actual overlap between the two shapes.

Sometimes overlap is not the question being asked at all. A delivery company wants to know which warehouse is closest to an address, even though the address point will never sit inside a warehouse's small footprint. That calls for a different kind of match entirely, one based on distance rather than containment, and GeoPandas gives it a separate function:

import geopandas as gpd

addresses = gpd.read_file("addresses.gpkg").to_crs("EPSG:32618")
warehouses = gpd.read_file("warehouses.gpkg").to_crs("EPSG:32618")

nearest = gpd.sjoin_nearest(addresses, warehouses, distance_col="dist_m")
print(nearest[["address_id", "warehouse_name", "dist_m"]].head(3))

sjoin_nearest exists as its own function rather than another predicate string, because nearest is not a yes-or-no overlap test at all. It measures the distance from each address to every candidate warehouse and keeps whichever one comes back smallest, overlap or no overlap, and distance_col names the new column that carries that measured distance along with the match.

One match becomes many, and a groupby follows

Point-in-polygon joins rarely stop at one match per polygon. A single neighborhood polygon usually contains dozens or hundreds of crime points, which means after the earlier join, that neighborhood's name appears dozens or hundreds of times over in the result, once per point tagged to it. This is the natural shape of the operation, not a mistake that needs fixing. It commits the analyst to a second step almost every time: collapsing many rows per polygon back down into one summary row per polygon before anything gets mapped, charted, or handed to someone else.

This is the same one-to-many situation the foundations module covered for attribute joins, where many business license rows shared one FIPS code and needed counting before they could join cleanly onto a single polygon per county. The fix here follows the identical shape, run with a groupby after the spatial match instead of before an attribute match:

tagged = gpd.sjoin(crimes, neighborhoods, predicate="within")
counts = tagged.groupby("neighborhood_name").size()
print(counts.head(3))

tagged.groupby("neighborhood_name").size() gathers every row that shares the same neighborhood name into one group and counts how many rows landed in it, turning a table with one row per crime back into a table with one row per neighborhood. Counting is not the only aggregation a groupby can run here. Swapping .size() for .groupby("neighborhood_name")["assessed_value"].mean() would compute the average assessed value of every property point that fell inside each neighborhood instead of a plain count, same shape, different question about the points that landed there.

Worked example: bus stops per 1000 residents

Putting the whole pipeline together end to end. The data is a GeoPackage of bus stop points and a GeoPackage of census tract polygons, the tracts already carrying a population field from an earlier attribute join. The task: for every tract, count how many bus stops fall inside it, then express that count as a rate per 1000 residents so a tract of 40,000 people is not unfairly compared against one of 4,000 on raw stop count alone.

import geopandas as gpd

stops = gpd.read_file("bus_stops.gpkg")
tracts = gpd.read_file("tracts.gpkg")

print(stops.crs == tracts.crs)
stops = stops.to_crs(tracts.crs)

joined = gpd.sjoin(stops, tracts, predicate="within")
stop_counts = joined.groupby("tract_id").size().reset_index(name="stop_count")

tracts = tracts.merge(stop_counts, on="tract_id", how="left")
tracts["stop_count"] = tracts["stop_count"].fillna(0)
tracts["stops_per_1000"] = tracts["stop_count"] / tracts["population"] * 1000

print(tracts[["tract_id", "population", "stop_count", "stops_per_1000"]].head(3))

stops.crs == tracts.crs prints False here, which is worth catching before the join rather than after. A spatial join compares raw coordinates between the two layers, and if one layer sits in a different coordinate reference system than the other, GeoPandas will warn about the mismatch but still run the comparison, silently lining up numbers that do not describe the same points in space. stops.to_crs(tracts.crs) fixes that by reprojecting the stops to match the tracts exactly, so the predicate that follows is comparing geometries that actually agree on what a coordinate means.

gpd.sjoin(stops, tracts, predicate="within") matches each stop to the tract polygon that contains it, same operation as the crime reports example. joined.groupby("tract_id").size() collapses that stop-level table into one row per tract, and .reset_index(name="stop_count") turns that grouped result back into an ordinary table with tract_id and stop_count as real columns instead of an index, which the merge on the next line needs in order to find tract_id at all.

The merge that follows is worth pausing on, because it is not a spatial join at all. tracts.merge(stop_counts, on="tract_id", how="left") matches the summarized counts back onto the tract polygons using tract_id, a shared key both tables now carry. The spatial join's job already finished the moment stops were matched to tracts. Reattaching the summarized result onto the tract geometry is the plain attribute join from the foundations module, doing exactly what it was built for. A real workflow like this one almost always chains both kinds of join: one to match by location, one to reattach a summary by key.

tracts["stop_count"].fillna(0) covers a gap that a .groupby() cannot fill on its own. A tract with zero bus stops never produced a matching row in joined at all, so it never shows up in stop_counts, and the merge leaves it null rather than zero, the same absence-versus-zero trap a raw groupby summary always carries into a merge. Filling that null with 0 before dividing turns a missing value into the correct number for a tract genuinely served by no stops. The last line, stop_count / population * 1000, is the rate itself: stops per resident, scaled up to a per-1000 figure that reads naturally and lets a tract of 40,000 people compare fairly against one of 4,000.

Every join in this lesson followed the same shape underneath: pick a predicate, or a distance for nearest, run it across every candidate pair of geometries, and expect a one-to-many result that needs a groupby before it goes anywhere near a map or a report. That shape is why a spatial join belongs in the same family as the attribute join rather than being some unrelated skill. Both attach outside information onto a layer's rows, one by key and one by location, and a workflow like the bus stop rate above usually needs both in sequence. The next lesson in this module leaves discrete points and polygons behind and turns to data stored as a grid of cells, where distance, overlap, and a join mean something a little different once nothing has a clean edge to test against.