How to Aggregate Spatial Data by Region in GeoPandas

On this page

Problem statement

Summarising points by region — collisions per district, injuries per ward, docks per neighbourhood — is a spatial join followed by an ordinary groupby. The join is the easy part. The hard part is the rows it does not match, because they leave the groupby without a trace: a total that is a few per cent short looks exactly like a correct one.

This page aggregates a month of New York City's reported motor-vehicle collisions by community district, and accounts for every collision that did not land in one. One in 23 did not, and none of the reasons was geometry.

Quick answer

Join the points to the regions, count what did not match, then group by the region key:

import geopandas as gpd

points = gpd.read_file("data/incidents.gpkg")
regions = gpd.read_file("data/districts.gpkg")

if points.crs != regions.crs:
    points = points.to_crs(regions.crs)

joined = gpd.sjoin(
    points,
    regions[["district_id", "geometry"]],
    how="left",
    predicate="within"
)

summary = (
    joined.groupby("district_id")
    .agg(
        point_count=("incident_id", "count"),
        total_injured=("injured", "sum"),
        avg_injured=("injured", "mean")
    )
    .reset_index()
)

result = regions.merge(summary, on="district_id", how="left")
result["point_count"] = result["point_count"].fillna(0).astype(int)

print(f"{len(points):,} collisions, {joined['district_id'].isna().sum():,} matched no district; "
      f"{len(regions)} districts, {(result.point_count == 0).sum()} with none")
print(result.nlargest(3, "point_count")[["district_id", "point_count", "total_injured"]].to_string(index=False))

7,483 collisions, 325 matched no district; 71 districts, 1 with none
 district_id  point_count  total_injured
         412          249          173.0
         413          244          140.0
         305          228          165.0

Aggregation pipeline

The aggregate-by-region pipeline: assign, then summarize.

Join, account for the unmatched, aggregate, merge back

A spatial join labels each point; groupby does the arithmetic.

1. Load the layers and join them

import geopandas as gpd

points = gpd.read_file("data/incidents.gpkg")
districts = gpd.read_file("data/districts.gpkg")
points = points.to_crs(districts.crs)                 # a no-op here; never assume it is

joined = gpd.sjoin(points, districts[["district_id", "geometry"]], how="left", predicate="within")
print(f"{len(points):,} collisions, {len(districts)} district polygons, {len(joined):,} joined rows")

7,483 collisions, 71 district polygons, 7,483 joined rows

how="left" keeps every collision, matched or not, which is what makes the next step possible. Group by the district code rather than a name, which can be spelled two ways.

2. Account for every collision the join did not place

unmatched = joined["district_id"].isna()
empty = points.geometry.is_empty
zero = ~empty & (points.geometry.x == 0) & (points.geometry.y == 0)
print(f"{unmatched.sum()} unmatched: {empty.sum()} with no location, {zero.sum()} at 0°, 0°, "
      f"{(unmatched & ~empty & ~zero).sum()} elsewhere")

both = gpd.sjoin(points[~empty], districts[["district_id", "geometry"]], predicate="intersects")
print(f"intersects: {both.index.duplicated().sum()} collisions matched two districts")

jia = joined["district_id"] % 100 > 18          # 164, 226, 480…: parks and airports, not districts
print(f"{jia.sum()} collisions in {joined.loc[jia, 'district_id'].nunique()} of the "
      f"{(districts.district_id % 100 > 18).sum()} joint interest areas among {len(districts)} polygons")

325 unmatched: 90 with no location, 181 at 0°, 0°, 54 elsewhere
intersects: 0 collisions matched two districts
140 collisions in 11 of the 12 joint interest areas among 71 polygons

One collision in 23 matches no district, and only 54 of the 325 are about geometry — points just outside every district polygon. The other 271 are data: 90 published with no location, and 181 at latitude 0, longitude 0, a placeholder rather than a place. With how="left" all of them stay in joined with an empty district_id, and step 3's groupby drops them without a word. On these coordinates intersects placed no collision twice. The other catch is in the region layer: 12 of its 71 polygons are joint interest areas — large parks and airports, with codes such as 164 and 480 — and 140 collisions landed in 11 of them. Whether they belong in the total is a decision about the question, not about the join.

3. Count, sum and average per district

summary = (joined.groupby("district_id")
           .agg(collisions=("incident_id", "size"),
                injured=("injured", "sum"),
                injured_per_collision=("injured", "mean"))
           .reset_index()
           .astype({"district_id": int}))
print(summary.sort_values("injured", ascending=False).head(5).round(2).to_string(index=False))
print(f"placed in districts: {summary.collisions.sum():,} of {len(points):,} collisions")

 district_id  collisions  injured  injured_per_collision
         412         249      173                   0.69
         305         228      165                   0.72
         413         244      140                   0.57
         318         177      137                   0.77
         407         210      117                   0.56
placed in districts: 7,158 of 7,483 collisions

The district codes are a borough digit and a district number, so the most injuries were in Queens district 12, 173 from 249 collisions, and Brooklyn district 5, 165 from 228. The last line is the check that matters: 7,158 collisions placed out of 7,483, and the 325 missing are exactly the ones step 2 accounted for.

4. Put the totals back on the map, with zeros kept apart from gaps

out = districts.merge(summary, on="district_id", how="left")
print("districts with no collisions:", int(out["collisions"].isna().sum()))
out["collisions"] = out["collisions"].fillna(0).astype(int)    # none happened: a real zero
out["injured"] = out["injured"].fillna(0)
# injured_per_collision stays empty where there were no collisions: there is no average of nothing
out.to_file("output/district_collisions.gpkg")

districts with no collisions: 1

One of the 71 polygons had no collision at all in June 2025. Its count and injury total are real zeros and are filled as such; its injuries per collision stays empty, because there is no average of nothing, and a zero there would drag down any mean taken across districts.

Bike-share stations summarised by district

The data used on this page

Real layers for New York City. The incidents are every motor-vehicle collision reported to the NYPD in June 2025, from NYC Open Data's Motor Vehicle Collisions – Crashes table, with the number of people injured and killed. The districts are the Department of City Planning's community districts, whose BoroCD code becomes district_id. The stations are Citi Bike's live GBFS station list with each station's dock capacity. The collisions are written exactly as published — including the rows with no location.

import json
import os

import geopandas as gpd
import pandas as pd
import requests

os.makedirs("data", exist_ok=True)
os.makedirs("output", exist_ok=True)
if not os.path.exists("data/incidents.gpkg"):
    reply = requests.get("https://data.cityofnewyork.us/resource/h9gi-nx95.csv", timeout=300, params={
        "$limit": 50000, "$where": "crash_date between '2025-06-01T00:00:00' and '2025-06-30T23:59:59'",
        "$select": "collision_id,latitude,longitude,number_of_persons_injured,number_of_persons_killed"})
    reply.raise_for_status()
    open("data/collisions.csv", "wb").write(reply.content)
    c = pd.read_csv("data/collisions.csv")
    gpd.GeoDataFrame({"incident_id": c.collision_id, "injured": c.number_of_persons_injured,
                      "killed": c.number_of_persons_killed},
                     geometry=gpd.points_from_xy(c.longitude, c.latitude), crs=4326
                     ).to_file("data/incidents.gpkg")

    cd = gpd.read_file("https://services5.arcgis.com/GfwWNkhOj9bNBqoJ/arcgis/rest/services/"
                       "NYC_Community_Districts/FeatureServer/0/query?where=1%3D1&outFields=BoroCD&outSR=4326&f=geojson")
    cd.rename(columns={"BoroCD": "district_id"})[["district_id", "geometry"]].to_file("data/districts.gpkg")

    st = requests.get("https://gbfs.citibikenyc.com/gbfs/en/station_information.json", timeout=60).json()
    st = pd.DataFrame(st["data"]["stations"])
    gpd.GeoDataFrame(st[["station_id", "name", "capacity"]], geometry=gpd.points_from_xy(st.lon, st.lat),
                     crs=4326).to_file("data/stations.geojson", driver="GeoJSON")

for f in ("incidents.gpkg", "districts.gpkg", "stations.geojson"):
    print(f, len(gpd.read_file(f"data/{f}")))

incidents.gpkg 7483
districts.gpkg 71
stations.geojson 2520

This example counts Citi Bike stations by district and summarizes their capacity — the number of docks.

import geopandas as gpd

stations = gpd.read_file("data/stations.geojson")
districts = gpd.read_file("data/districts.gpkg")

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

joined = gpd.sjoin(
    stations,
    districts[["district_id", "geometry"]],
    how="left",
    predicate="within"
)

district_docks = (
    joined.groupby("district_id")
    .agg(
        station_count=("station_id", "count"),
        docks_sum=("capacity", "sum"),
        docks_mean=("capacity", "mean")
    )
    .reset_index()
    .astype({"district_id": int})
)

districts_out = districts.merge(district_docks, on="district_id", how="left")
districts_out["station_count"] = districts_out["station_count"].fillna(0).astype(int)
districts_out["docks_sum"] = districts_out["docks_sum"].fillna(0)
print(district_docks.sort_values("docks_sum", ascending=False).head(3).round(1).to_string(index=False))
print(f"{(districts_out.station_count == 0).sum()} of {len(districts_out)} polygons have no station; "
      f"{joined.district_id.isna().sum()} stations fall in none")

 district_id  station_count  docks_sum  docks_mean
         105             74       4992        67.5
         103             79       4118        52.1
         301            119       3769        31.7
27 of 71 polygons have no station; 109 stations fall in none

Manhattan's districts 5 and 3 hold the most docks, 4,992 and 4,118, and 27 of the 71 polygons have no station at all. The other count is the reason to look: 109 stations fall in no district. They are Citi Bike's stations in Jersey City and Hoboken — every one lies west of the Hudson, outside New York's boroughs — so the system's total and the city's total are two different numbers.

Export the summarized polygons

Save the result for use in QGIS, a web map, or another script.

districts_out.to_file("output/district_docks.gpkg", driver="GPKG")

Why the missing rows leave no trace

All three "aggregate by region"; only one gives you polygons ready to map.

An aggregation by region is two operations: a spatial join that writes a district code onto each point, and a groupby on that code. Each is correct on its own, and the gap between them is where totals go short. A collision the join could not place gets an empty code, and groupby drops empty keys by default — so it is not counted as "unknown", it is not counted at all. Nothing warns, and the sum of the district totals is simply smaller than the number of collisions.

The reasons a point goes unplaced are mostly not spatial. Here they were a missing location, a 0°, 0° placeholder, and a handful of points off the edge of the city's polygons. Each calls for a different decision, which is why the count has to come before the groupby, broken down by reason. The same goes for the region layer: New York's community-district file includes parks and airports as polygons of their own, and whether a collision on an airport road belongs in the total is a question about the analysis, not about the join.

Overlapping regions, lines and polygons, and averages

Overlapping or touching regions can place a point twice. On New York's districts intersects duplicated nothing; on a layer with overlaps, check joined.index.duplicated().sum() before summing.
Lines and polygons cross boundaries. A road or parcel spanning two districts joins to both; for totals shared in proportion to length or area, overlay first and aggregate the pieces.
size counts rows, count counts non-null values in one column. After a left join they differ for exactly the unmatched rows.
A mean has no zero. Fill missing counts and sums with 0 after merging back; leave averages empty.

Internal links

For the broader concept behind this workflow, see Spatial Joins in GeoPandas: How They Work and When to Use Them.

If your main task is just point counting, see How to Count Points Within Polygons in GeoPandas.

If you need to merge polygons and summarize their attributes, see How to Dissolve Polygons and Aggregate Attributes in GeoPandas.

If your join returns no matches or obviously wrong results, check How to Fix CRS Mismatch and Empty Spatial Join Results in GeoPandas.

How to aggregate events to route segments — the same question along a line

FAQ

How do I find the most common category per region instead of a sum or mean?

Use a named aggregation with a lambda that returns the mode, for example top_class=("class", lambda s: s.mode().iloc[0]) inside groupby().agg(). This gives the most frequent value in each region, with ties resolved by the first sorted value.

How can I aggregate by two columns, such as region and year, at the same time?

Pass a list to groupby, like joined.groupby(["district_id", "year"]).agg(...). The result has one row per region-year combination, which you can later pivot or merge back to the region polygons as needed.

How do I compute area-weighted statistics when features span multiple regions?

A plain sjoin() cannot split a feature across regions, so use gpd.overlay(features, regions, how="intersection") to cut each feature at region boundaries first. Recompute each piece's area in a projected CRS, then weight the values by that area before summing per region.

Can I sum the polygon area inside each region instead of summarizing an attribute?

Yes. After assigning features to regions, compute area on a projected layer with joined["piece_area"] = joined.geometry.area, then groupby("district_id")["piece_area"].sum(). Use overlay() rather than sjoin() if features cross region boundaries, so only the part inside each region is measured.