How to Count Points in Polygons with GeoPandas

On this page

Problem statement

Counting points per polygon is a spatial join and a groupby — and on Bristol's street-level crimes, the join matches nothing until the two layers share a CRS, and most wards come back missing rather than zero. The same pattern counts:

counting addresses inside parcels
counting stores inside districts
counting incident locations inside census tracts
counting GPS observations inside management zones

The expected output is usually a polygon layer with a new count field, such as point_count, where each polygon stores the number of matching points.

In GeoPandas, the standard workflow is:

load the point and polygon layers
make sure both layers use the same CRS
spatially match points to polygons
aggregate the matches by polygon
join the counts back to the polygon layer

Quick answer

A spatial join counts the points inside each polygon — polygons with none keep a count of zero.

To count points in polygons with GeoPandas:

use geopandas.sjoin() to match point features to polygons
group the joined records by a polygon ID
count matched points per polygon
merge the counts back into the polygon GeoDataFrame

Both layers must use the same CRS before the spatial join, and the polygon layer should have a stable unique ID such as polygon_id.

import geopandas as gpd

points = gpd.read_file("data/crimes.gpkg")
polygons = gpd.read_file("data/districts.gpkg")

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

polygons = polygons.reset_index(drop=True)
polygons["polygon_id"] = polygons.index

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

counts = (
    joined.dropna(subset=["polygon_id"])
    .groupby("polygon_id")
    .size()
    .reset_index(name="point_count")
)

result = polygons.merge(counts, on="polygon_id", how="left")
result["point_count"] = result["point_count"].fillna(0).astype(int)
print(result.nlargest(3, "point_count")[["WD24NM", "point_count"]].to_string(index=False))
print((result["point_count"] == 0).sum(), "of", len(result), "wards with none")

       WD24NM  point_count
      Central          894
Lawrence Hill          259
       Ashley          213
25 of 34 wards with none

Align, join, count and join the counts back

Without the final merge and fillna, empty polygons silently disappear.

Load the point and polygon layers

Read both datasets with gpd.read_file(). GeoPandas supports shapefiles, GeoJSON, and GeoPackage files.

import geopandas as gpd

points = gpd.read_file("data/crimes.gpkg")
polygons = gpd.read_file("data/districts.gpkg")

print(points[["category", "month", "location.street.name"]].head(3))
print(polygons[["WD24CD", "WD24NM"]].head(3))
print(points.geom_type.unique().tolist(), polygons.geom_type.unique().tolist())

                category    month             location.street.name
0  anti-social-behaviour  2026-06             On or near Nightclub
1  anti-social-behaviour  2026-06          On or near Parking Area
2  anti-social-behaviour  2026-06  On or near Theatre/concert Hall
      WD24CD                         WD24NM
0  E05010885                         Ashley
1  E05010886  Avonmouth and Lawrence Weston
2  E05010887                     Bedminster
['Point'] ['MultiPolygon']

Check that:

the point layer contains point geometries
the polygon layer contains polygon geometries
the polygon layer has a usable identifier field, or you can create one

Check and align coordinate reference systems

Point-in-polygon operations require both layers to use the same CRS.

print("Points CRS:", points.crs)
print("Polygons CRS:", polygons.crs)

print("before reprojecting:", len(gpd.sjoin(points, polygons, predicate="within")), "matches")

if points.crs != polygons.crs:
    points = points.to_crs(polygons.crs)
print("after:", len(gpd.sjoin(points, polygons, predicate="within")), "matches")

Points CRS: EPSG:4326
Polygons CRS: EPSG:27700
before reprojecting: 0 matches
after: 1691 matches

The crimes come in longitude and latitude and the wards in British National Grid: joined as they are, nothing matches, and the only sign is a UserWarning about the CRS. The join does not fail.

Make sure polygons have a unique ID

You should group counts by a stable unique field, not by a name field that may repeat.

If the polygon layer does not already have a unique ID, create one:

if "polygon_id" not in polygons.columns:
    polygons = polygons.reset_index(drop=True)
    polygons["polygon_id"] = polygons.index + 1
print(polygons["polygon_id"].is_unique, polygons["WD24NM"].is_unique)

True True

Using district_name or tract_name alone can fail if multiple polygons share the same name.

Run a spatial join to match points to polygons

For most point counting workflows, it is easiest to join points to polygons. Each matched point gets the polygon ID it falls inside.

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

print(joined[["category", "polygon_id"]].head(3))

                category  polygon_id
0  anti-social-behaviour           8
1  anti-social-behaviour           8
2  anti-social-behaviour           8

This keeps all point records and adds polygon attributes where a match exists.

within is a good default when your rule is to count points strictly inside polygons. If points on boundaries should also count, use a different predicate such as intersects.

Count points per polygon

Now group the joined point records by polygon ID and count the number of matched rows.

counts = (
    joined.dropna(subset=["polygon_id"])
    .groupby("polygon_id")
    .size()
    .reset_index(name="point_count")
)

print(counts.head(), f"\n{len(counts)} of {len(polygons)} polygons appear")

   polygon_id  point_count
0           1          213
1           8          894
2           9            5
3          10            5
4          11           47 
9 of 34 polygons appear

Polygons with no matching points do not appear in this table yet — here most of the city's wards, because the crimes were fetched for one mile around the centre.

Join the counts back to the polygons

Merge the aggregated count table back into the polygon GeoDataFrame.

result = polygons.merge(counts, on="polygon_id", how="left")
result["point_count"] = result["point_count"].fillna(0).astype(int)

print(result[["polygon_id", "WD24NM", "point_count"]].head())

   polygon_id                         WD24NM  point_count
0           1                         Ashley          213
1           2  Avonmouth and Lawrence Weston            0
2           3                     Bedminster            0
3           4     Bishopston and Ashley Down            0
4           5                   Bishopsworth            0

Using a left merge from the polygon layer ensures that polygons with no points are kept in the final output.

Counted by crime ID and by distinct location, from Example 3.

Save the result

Write the final polygons with counts to a new file.

result.to_file("output/census_tracts_point_counts.gpkg", driver="GPKG")

You can also export to shapefile or GeoJSON:

result.to_file("output/census_tracts_point_counts.shp")
result.to_file("output/census_tracts_point_counts.geojson", driver="GeoJSON")
print(list(gpd.read_file("output/census_tracts_point_counts.shp").columns))

['WD24CD', 'WD24NM', 'polygon_id', 'point_coun', 'geometry']

GeoPackage is often the safer choice because it handles field names and data types better than shapefiles — the shapefile above came back with point_count cut to point_coun, the format's ten-character limit.

Bristol crimes counted by area

The data used on this page

Street-level crime from data.police.uk (Open Government Licence) for June 2026: every crime recorded within one mile of the centre of Bristol. The locations are anonymised — each crime is moved to the nearest of a fixed set of map points on streets and public places — so many crimes share a point. The polygons are Bristol's 34 wards from the ONS (full resolution, British National Grid).

from pathlib import Path

import geopandas as gpd
import pandas as pd
import requests

Path("data").mkdir(exist_ok=True)
Path("output").mkdir(exist_ok=True)
if not Path("data/crimes.gpkg").exists():
    rows = requests.get("https://data.police.uk/api/crimes-street/all-crime",
                        params={"lat": 51.4545, "lng": -2.5879, "date": "2026-06"}, timeout=120).json()
    df = pd.json_normalize(rows)
    gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df["location.longitude"].astype(float),
                                                     df["location.latitude"].astype(float)),
                     crs="EPSG:4326").to_file("data/crimes.gpkg")
    gpd.read_file("https://services1.arcgis.com/ESMARspQHYMw9BZ9/arcgis/rest/services/"
                  "Wards_December_2024_Boundaries_UK_BFC/FeatureServer/0/query"
                  "?where=LAD24CD%3D%27E06000023%27&outFields=WD24CD,WD24NM&outSR=27700&f=geojson"
                  ).to_file("data/districts.gpkg")
crimes = gpd.read_file("data/crimes.gpkg")
print(len(crimes), "crimes at", crimes.geometry.to_wkt().nunique(), "distinct points;",
      len(gpd.read_file("data/districts.gpkg")), "wards")

1691 crimes at 425 distinct points; 34 wards

Example 1: Basic point count in polygons

import geopandas as gpd

points = gpd.read_file("data/crimes.gpkg")
polygons = gpd.read_file("data/districts.gpkg")

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

if "polygon_id" not in polygons.columns:
    polygons = polygons.reset_index(drop=True)
    polygons["polygon_id"] = polygons.index + 1

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

counts = (
    joined.dropna(subset=["polygon_id"])
    .groupby("polygon_id")
    .size()
    .reset_index(name="point_count")
)

result = polygons.merge(counts, on="polygon_id", how="left")
result["point_count"] = result["point_count"].fillna(0).astype(int)
print(result["point_count"].sum(), "crimes counted in", (result["point_count"] > 0).sum(), "wards")

1691 crimes counted in 9 wards

Example 2: Count points using `intersects` instead of `within`

joined_i = gpd.sjoin(
    points,
    polygons[["polygon_id", "geometry"]],
    how="left",
    predicate="intersects"
)
print(len(joined), "rows with within,", len(joined_i), "with intersects")

1691 rows with within, 1691 with intersects

Use this when you want to count points inside polygons as well as points that fall exactly on polygon boundaries. On the crime points it changes nothing: floating-point locations almost never land exactly on a boundary line. It matters when points are derived from the boundaries — survey markers, vertices, points snapped to the same linework:

from shapely.geometry import MultiPoint

shared = polygons.geometry.iloc[0].boundary.intersection(polygons.geometry.iloc[1:].union_all().boundary)
markers = gpd.GeoDataFrame(geometry=[shared.representative_point()], crs=polygons.crs)
print("within:    ", len(gpd.sjoin(markers, polygons, predicate="within")), "matches")
print("intersects:", len(gpd.sjoin(markers, polygons, predicate="intersects")), "matches")

within:     0 matches
intersects: 2 matches

A point on the line two wards share is in neither ward by within and in both by intersects — counted twice. Neither is right by default; the rule has to be chosen, and deduplicated if it double-counts.

Example 3: Count unique point IDs instead of all joined rows

joined["location_key"] = joined.geometry.to_wkt()
per_ward = joined.dropna(subset=["polygon_id"]).groupby("polygon_id").agg(
    crimes=("id", "nunique"), locations=("location_key", "nunique"))
print(per_ward.sum().to_dict())
print(per_ward.sort_values("crimes", ascending=False).head(3))

{'crimes': 1691, 'locations': 425}
            crimes  locations
polygon_id                   
8              894        159
23             259         98
1              213         71

Use this if your point layer may contain duplicate records and you only want to count unique features — and decide what "unique" means first. Counted by crime ID, every one of the crimes counts; counted by distinct location, a few hundred points do, because police.uk snaps many crimes to the same anonymised point. Neither is wrong; they answer different questions.

Example 4: Export counted polygons to file

result.to_file("output/district_store_counts.gpkg", driver="GPKG")
print(gpd.list_layers("output/district_store_counts.gpkg").to_dict("records"))

[{'name': 'district_store_counts', 'geometry_type': 'MultiPolygon'}]

This output can be opened directly in QGIS for styling, labeling, or further analysis.

Counting starts with a spatial join

gpd.sjoin() is used because counting points in polygons starts with a spatial relationship: which points fall inside which polygons. The spatial join attaches polygon attributes to matching point features.

The first two cover nearly every case, and neither raises an error.

The count is calculated after the join because the join creates the row-level matches first. Once each point has a polygon ID, standard pandas grouping can count how many points belong to each polygon.

That is why the workflow is:

spatial join first
aggregation second

There is also a difference between:

counting joined rows with .size()
counting unique point IDs with .nunique()

If every point feature should count once, .size() is usually correct. If your point data may contain duplicates and you want unique features only, count a point ID field instead.

Predicate choice matters:

use within when boundary points should not count
use intersects when boundary points should count too

For large datasets, performance depends on:

number of points and polygons
geometry complexity
spatial index availability

GeoPandas uses a spatial index when available, which is much faster than checking every point against every polygon manually.

Boundary points, empty polygons and multipart areas

Points on polygon boundaries

Boundary cases depend on the predicate:

within usually excludes points exactly on the boundary
intersects includes points inside polygons and points touching polygon boundaries

Choose the predicate that matches your counting rule.

Polygons with no points

If you only look at the grouped count table, polygons with no matches are missing. To keep them, always merge the counts back to the full polygon layer with how="left" and fill null values with zero.

Duplicate or multipart polygons

If multiple polygons share the same polygon_id, their counts will be combined. Make sure the ID is truly unique.

Multipart polygons are supported. A point inside any part of the multipart feature will count toward that polygon.

Invalid geometries

Broken polygon geometries can cause unexpected join results. If results look wrong, check validity:

invalid = polygons[~polygons.is_valid]
print(len(invalid), "invalid polygons")

0 invalid polygons

If invalid features exist, repair them before running the join.

CRS issues

A mismatched CRS is one of the most common reasons for empty results. Always compare points.crs and polygons.crs before running sjoin().

Internal links

If you want the broader pattern behind this workflow, see spatial joins in GeoPandas.

For related tasks, see:

If your join returns no matches, check why GeoPandas spatial join is returning empty results.

FAQ

How do I count points in polygons with GeoPandas?

Use geopandas.sjoin() to match points to polygons, then group by a polygon ID and count the matched records. Merge the counts back into the polygon layer.

What is the difference between `within` and `intersects` when counting points?

within counts points strictly inside polygons. intersects counts points inside polygons and points that touch polygon boundaries. Boundary points can change your totals.

Why are some polygons missing from my output after counting points?

They are usually missing because only matched polygons appear in the grouped count table. Merge the counts back to the full polygon layer with a left join and fill missing values with zero.

How do I keep polygons that contain no points?

Use a left merge from the polygon GeoDataFrame to the count table:

result = polygons.merge(counts, on="polygon_id", how="left")
result["point_count"] = result["point_count"].fillna(0).astype(int)

This keeps all polygons and assigns 0 where no points were found.

Can I count points in polygons without creating a polygon_id column?

Yes. Reset the polygon index and group on it directly: after joined = gpd.sjoin(points, polygons, predicate="within"), the polygon index appears as index_right, so joined.groupby("index_right").size() gives counts you can join back by index. A stable polygon_id is still clearer and safer when polygon names repeat.

How do I sum a point attribute per polygon instead of just counting?

After the spatial join, group by the polygon ID and aggregate the value column instead of using .size(). For example joined.groupby("polygon_id")["sales"].sum().reset_index(name="total_sales"), then merge back with how="left" and fill missing polygons with zero.

Is sjoin or a contains-based loop faster for counting points in polygons?

gpd.sjoin() is faster because it uses a spatial index, so it avoids testing every point against every polygon:

import time

t = time.perf_counter()
loop = {w.polygon_id: int(sum(w.geometry.contains(p) for p in points.geometry)) for w in polygons.itertuples()}
t_loop = time.perf_counter() - t
t = time.perf_counter()
fast = gpd.sjoin(points, polygons, predicate="within").groupby("polygon_id").size()
t_join = time.perf_counter() - t
print(f"loop {t_loop:.2f} s, sjoin {t_join:.3f} s, same counts: {fast.to_dict() == {k: v for k, v in loop.items() if v}}")

loop 0.15 s, sjoin 0.004 s, same counts: True

Why do some points get counted in two polygons at once?

This happens when polygons overlap or share boundaries and you used predicate="intersects", so a boundary point matches more than one polygon. Use predicate="within" to count strictly interior points, or deduplicate the join on the point ID before counting.

How can I speed up counting when I have millions of points?

Make sure both layers share a CRS and let sjoin use the spatial index it builds automatically. pygeos is no longer needed — its vectorised engine became Shapely 2, which current GeoPandas requires. For very large jobs, join in chunks or use dask-geopandas to manage memory.