Counting points in overlapping polygons using GeoPandas

I have converted a 200m x 200m point grid of Greater London into a multypolygon 500m radius buffer layer for each point in the grid. What this means is that I have over 100,000 overlapping polygons.

I also have a years worth of crime data as a point layer with lat longs (over 1.1million crimes x 12 columns of data)

I am trying to find the most efficient way to count the number of crime points in each polygon buffer. As the polygon buffers are overlapping the crime points will overlap too for all of the buffers.

The spatial join in GeoPandas doesn't seem to work, maybe because the polygons are overlapping?

If I use "inner" join I just get a blank dataframe back. If I use "left" join then I just get all the crime rows (1.1million) with the buffer polygon columns to the right all as "nan". And vice versa if I use "right" join - just the buffer rows (100,000) with crime columns as nan. See the code below:

``````import pandas as pd
import geopandas as gpd
from geopandas import GeoDataFrame, read_file, points_from_xy

#import buffer polygon layer

#import crime csv

#drop nan rows from coords
crime2 = crime[crime['Longitude'].notna()]

#geocode crime points
gCrime = GeoDataFrame(crime2, geometry=points_from_xy(crime2['Longitude'], crime2['Latitude']))

#set equal crs
gCrime.crs = gBuffer.crs

#spatial join data
BufferCrime = gpd.sjoin(gCrime, gBuffer, how="inner")
``````

The other solution is to iterate over each polygon and count the number of points but this will take forever given that it has to do 100,000 x 1,100,000 iterations.

``````# Loop over polygons with index i.
for i, poly in gBuffer.iterrows():

#list of points in this poly
pts_in_this_poly = []

#loop over all points
for j, pt in gCrime.iterrows():
if poly.geometry.contains(pt.geometry):
# Add it to the list
pts_in_this_poly.append(pt.geometry)

pts_in_polys.append(len(pts_in_this_poly))

gBuffer['number of Crime points'] = gpd.GeoSeries(pts_in_polys)
``````

How can I solve this problem?

You can play it a bit smart and use some numerical calculations to make your code more efficient.

It is better to iterate over polygons. Although it is 100,000, you only need to do it once so your PC can handle it.

lets import the data

``````crime = read_csv('2020-2021 London Crime.csv')

#drop nan rows from coords
crime2 = crime[crime['Longitude'].notna()]

#geocode crime points
gCrime = GeoDataFrame(crime2, geometry=points_from_xy(crime2['Longitude'], crime2['Latitude']))
``````

In order to simplify the `sjoin` operation, we need the bounds of each feature of your buffers:

``````#import buffer polygon layer

gBuffer [['minx','miny','maxx','maxy']] = gBuffer.geometry.bounds
``````

Now we are going to iterate over buffers. We can use the buffer bounds to filter the points and then use the spatial join. This will reduce the cross-match between all of the points and the polygons.

``````joined_result = gpd.GeoDataFrame()
for i, poly in gBuffer.iterrows():
selected_points = gCrime[
(gCrime.geometry.x < poly["maxx"]) &
(gCrime.geometry.x > poly["minx"]) &
(gCrime.geometry.y < poly["maxy"]) &
(gCrime.geometry.y > poly["miny"])
]

row = pd.DataFrame(poly).T
row = gpd.GeoDataFrame(row, geometry = 'geometry') # converting row to GDF for sjoin

temp_result = gpd.jsoin(selected_points, row, op='intersects', how='inner')

joined_result = joined_result.append(temp_result, ignore_index=True)
``````

The `joined_result` is your final result. Since you are using a grid, you can use the `selected_points` GDF as the join result if you don't mind the buffer with straight corners:

``````joined_result = gpd.GeoDataFrame()
for i, poly in gBuffer.iterrows():
selected_points = gCrime[
(gCrime.geometry.x < poly["maxx"]) &
(gCrime.geometry.x > poly["minx"]) &
(gCrime.geometry.y < poly["maxy"]) &
(gCrime.geometry.y > poly["miny"])
]

selected_points['buffer_ID'] = poly[<ID VALUE KEY>]
joined_result = joined_result.append(selected_points)
``````

These solutions might take a bit of time, but they will be much less power- and time-consuming Since we are using numerical operations to optimize the spatial join.