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?