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?