Spatial Joins, Dissolve, and Overlay in GeoPandas

The previous lesson made reprojecting a habit: check the CRS of every layer before doing anything with it, and reproject whichever one does not match. That habit pays off directly here, because a spatial join is really just a question about distance and containment, and both of those go wrong silently on data in the wrong CRS. A join between bus stops in EPSG:4326 and census tracts in EPSG:3857 will not raise an error. It will just return nonsense, because "within" means something different when one layer measures in degrees and the other in meters.

With both layers on the same CRS, the question a spatial join answers is one that comes up constantly: given a table of points and a table of polygons, which polygon does each point sit inside? A city has thousands of addresses and a few hundred zoning districts. A delivery service has thousands of drop-off pins and a few dozen service zones. An analyst working by hand in QGIS would run Vector Analysis, Join Attributes by Location, and click through a dialog. The same operation in GeoPandas is one function call, and because it is code, it can run again tomorrow on a fresh export without anyone touching a mouse.

Tagging points with the polygon they fall in

GeoPandas calls this operation sjoin, short for spatial join, and it lives in the geopandas module itself rather than on a single GeoDataFrame, because it needs two of them.

import geopandas as gpd

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

# both layers must already share a CRS, or the join is meaningless
stops = stops.to_crs(tracts.crs)

joined = gpd.sjoin(stops, tracts, predicate="within")
print(joined[["stop_name", "GEOID"]].head(3))

predicate is the argument that names the spatial relationship being tested. "within" asks, for every row in the left GeoDataFrame (stops), whether its geometry sits entirely inside a geometry in the right GeoDataFrame (tracts). Other predicates exist, "intersects" and "contains" among them, but "within" is the natural choice for points against polygons, since a point cannot partly overlap anything. It is either inside a tract or it is not.

The result, joined, is a new GeoDataFrame with every column from stops plus every column from tracts for the matching row, including an index_right column that records which row of tracts supplied the match. A stop that falls inside tract 001200 now carries that tract's GEOID (the Census Bureau's geographic identifier) sitting right next to its own name and coordinates, without the tract's polygon geometry itself, which GeoPandas drops from the right side to avoid carrying two geometry columns at once.

Keeping every row or only the matches

sjoin takes a how argument that decides what happens to a row on the left side that finds no match at all, a bus stop sitting just outside every tract boundary, say, because of a sliver of bad digitizing or a stop technically in the water.

how="inner", the default, drops any stop that matches no tract. The result only contains rows that found a home. how="left" keeps every stop regardless, filling GEOID and the other tract columns with NaN (pandas's marker for a missing value) wherever no match turned up.

inner_joined = gpd.sjoin(stops, tracts, how="inner", predicate="within")
left_joined = gpd.sjoin(stops, tracts, how="left", predicate="within")

print(len(stops), len(inner_joined), len(left_joined))

Three fewer rows in the inner join than in the original 412 stops means three stops matched no tract at all. The left join still reports all 412, because it keeps them with blank tract columns instead of dropping them. Which one to reach for depends on the question being asked. Counting stops per tract calls for inner, since a stop with no tract contributes nothing to any tract's count anyway. Auditing data quality, finding which stops need a second look, calls for left, because the whole point is to see the unmatched rows, not lose them.

Turning a join into a count or an average

A join by itself answers "which polygon does this point belong to." The more common question is one step further: how many points landed in each polygon, or what is the average of some point attribute per polygon. That second step is an ordinary pandas groupby, run on the result of the join exactly as it would run on any other table, because once the join has attached a GEOID to every stop, the geometry has already done its job.

counts = joined.groupby("GEOID").size().reset_index(name="stop_count")
print(counts.sort_values("stop_count", ascending=False).head(4))

groupby("GEOID") gathers every row that shares the same tract identifier into one bucket. .size() counts the rows in each bucket, which here means the number of bus stops in that tract, and .reset_index(name="stop_count") turns the result back into an ordinary two-column table instead of pandas's grouped-index format. The same pattern with .mean() in place of .size() would average a numeric column instead, average ridership per tract, say, if the stops carried a ridership figure. Nothing about this step is spatial anymore. The spatial work already happened in the join; the aggregation is plain tabular arithmetic riding on top of it.

Merging features that share an attribute

A spatial join attaches a category to individual features without changing their shapes. dissolve does the opposite kind of merge: it takes many polygons that share a value in some column and fuses them into one polygon per distinct value, unioning their geometries together.

tracts = gpd.read_file("census_tracts.gpkg")
regions = tracts.dissolve(by="region", aggfunc={"population": "sum"})
print(regions[["population"]])

Passing by="region" tells dissolve which column names the groups. Every tract sharing the same region value gets its geometry unioned into a single multipart or single-part polygon, and the returned GeoDataFrame has one row per distinct region instead of one row per tract. The aggfunc argument decides what happens to the other columns during that collapse, since population cannot simply be picked from one of the merged tracts, it has to be added up across all of them. "sum" does that addition; other columns could use "mean", "max", or a custom function depending on what the column actually measures. This is the scripted version of drawing a county map by grouping and dissolving a parcel layer by hand in QGIS, except it runs on a schedule and never drifts from a documented rule about how the aggregation works.

Cutting one polygon layer by another

overlay solves a different problem again: combining two polygon layers by their geometry, not by a shared attribute column. Flood zones and parcels rarely share any identifier at all, but an analyst still needs to know exactly which piece of each parcel sits inside the flood zone, because a parcel that is only ten percent inside the zone is a very different risk than one that is fully inside it.

flood_zones = gpd.read_file("flood_zones.gpkg")
parcels = gpd.read_file("parcels.gpkg").to_crs(flood_zones.crs)

flooded_parts = gpd.overlay(parcels, flood_zones, how="intersection")
flooded_parts["flood_area_m2"] = flooded_parts.geometry.area
print(flooded_parts[["parcel_id", "flood_area_m2"]].head(2))

how="intersection" keeps only the geometry where the two input layers actually overlap, and it clips each parcel down to exactly that overlapping shape, which is why flooded_parts["parcel_id"] can now carry a brand new, smaller area than the original parcel had. Other overlay modes exist. "union" keeps every piece from both layers, overlapping and not. "difference" keeps only the part of the first layer that does not overlap the second, the parts of each parcel safely outside the flood zone, in this example. overlay is the scripted twin of running Vector Overlay, Intersection in a desktop GIS, with the advantage that the exact CRS, predicate, and column handling are written down in a file instead of remembered by whoever ran the tool last.

A worked example: bus stops per census tract

Putting the pieces from this lesson together produces a small, repeatable report: how many bus stops fall in each census tract, sorted from the busiest tract down.

import geopandas as gpd

stops = gpd.read_file("bus_stops.gpkg")
tracts = gpd.read_file("census_tracts.gpkg")
stops = stops.to_crs(tracts.crs)

joined = gpd.sjoin(stops, tracts, how="inner", predicate="within")
report = (
    joined.groupby("GEOID")
    .size()
    .reset_index(name="stop_count")
    .merge(tracts[["GEOID", "NAME"]], on="GEOID")
    .sort_values("stop_count", ascending=False)
)
print(report[["NAME", "stop_count"]].head(5))

Every line in that script maps to something covered above. Reprojecting brings both layers onto one CRS. sjoin with how="inner" and predicate="within" attaches each stop's tract to it and drops any stop that landed outside every tract. groupby and size turn that attachment into a count per tract. The final merge is an ordinary pandas join, not a spatial one, that pulls in the tract's human-readable name so the report reads as "Downtown Core" instead of a fourteen-digit GEOID. A city planner asking "which neighborhoods are best served by transit" gets an answer built entirely from open data and four lines of real work, and the exact same script runs again next year on a refreshed stops file without anyone re-clicking a single dialog box.

The next lesson turns to rasters instead of vectors, reading a satellite image band with rasterio and pulling the pixel value under a specific point. The join, dissolve, and overlay covered here only ever combine vector features with other vector features. Rasters need their own tools, because a grid of cells has no attribute table and no polygon boundaries to test containment against, just a value sitting at every location on a regular lattice.