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