I'm working towards plotting a large GIS dataset of which I've shown a sample above of about 1/6 of the data. I'm happy with how quickly the data loads in, and bokeh
renders the html nearly instantaneously. However, I've encountered a pretty hot loop in my code that is not scaling well as I increase the 1) number of rows and 2) resolution of the polygons. I'm just getting killed in the #count points
loop and am wondering if there isn't a better way of doing this?
I found the suggestion for a loop off a GIS readthedoc.io and was happy with its performance for a few thousand points a couple months ago. But now the project needs to process a GeoDataFrame
with >730000 rows. Is there a better method I'm suppose to be using to count the number of points in each polygon? I'm on a modern desktop to do the computation but the project has access to Azure resources so maybe that's most people professionally do this sort of computation? I'd prefer to do the computation locally but it means my desktop might have to sit at max cpu cycles overnight or longer which isn't a thrilling prospect. I'm using Python 3.8.2 & Conda 4.3.2.
from shapely.geometry import Polygon
import pysal.viz.mapclassify as mc
import geopandas as gpd
def count_points(main_df, geo_grid, levels=5):
"""
outputs a gdf of polygons with a columns of classifiers to be used for color mapping
"""
pts = gpd.GeoDataFrame(main_df["geometry"]).copy()
#counts points
pts_in_polys = []
for i, poly in geo_grid.iterrows():
pts_in_this_poly = []
for j, pt in pts.iterrows():
if poly.geometry.contains(pt.geometry):
pts_in_this_poly.append(pt.geometry)
pts = pts.drop([j])
nums = len(pts_in_this_poly)
pts_in_polys.append(nums)
geo_grid['number of points'] = gpd.GeoSeries(pts_in_polys) #Adds number of points in each polygon
# Adds Quantiles column
classifier = mc.Quantiles.make(k=levels)
geo_grid["class"] = geo_grid[["number of points"]].apply(classifier)
# Adds Polygon grid points to new geodataframe
geo_grid["x"] = geo_grid.apply(getPolyCoords, geom="geometry", coord_type="x", axis=1)
geo_grid["y"] = geo_grid.apply(getPolyCoords, geom="geometry", coord_type="y", axis=1)
polygons = geo_grid.drop("geometry", axis=1).copy()
return polygons