Counting points in overlapping polygons using GeoPandas

Question

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 pandas import read_csv
from geopandas import GeoDataFrame, read_file, points_from_xy

#import buffer polygon layer
gBuffer = read_file('London Buffer.zip')
df1 = gBuffer.head()

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

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

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

#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))
    
#Add the points
gBuffer['number of Crime points'] = gpd.GeoSeries(pts_in_polys) 

How can I solve this problem?

Answer 1

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 = read_file('London Buffer.zip')

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.