ArcGIS: calculate the number of points in a polygon unit with Python GeoPandas batch processing

created at 07-01-2021 views: 10

In ArcMap, the operation is very simple. Use the Spatial Join tool directly. Pay attention to two points: In the output options, control the attribute fields to be included in the output (all by default, including First, Sum, Mode, Count, etc.) ), in the matching conditions, you can choose options such as intersection, location, inclusion, boundary contact, full inclusion, center within the element range, etc. For the selection of these two important parameters, please refer to ArcGIS Spatial Join.

In Python, Geopandas is selected for processing. The difference is that there is no default summary function in ArcMap for the generated DataFrame, but the summary function of Pandas can be used. The code is as follows:

import geopandas as gpd

input_path = "E:\\Data\\shp\\"
output_path = "E:\\Data\\join\\"
# Read the surface data, adjust the name of the latitude and longitude, and set the coordinate system
gdf_polys = gpd.read_file(input_path+"poly.shp")
gdf_polys = gdf_polys.rename({'lng': 'lng_c', 'lat': 'lat_c'})
gdf_polys.to_crs(crs="EPSG:4326", inplace=True)
# Read point data, if necessary, set the coordinate system
gdf_points = gpd.read_file(input_path+"point.shp")

# Spatial join
gdfsjoin = gpd.sjoin(gdf_points, gdf_polys, how='left', op='within', lsuffix='_pt', rsuffix='_py')
# aggregate
df = gdfsjoin.drop('geometry', axis=1).groupby(by=['poly_index']).aggregate({'point_id': 'count'})
# output
df.to_csv(output_path+'join.csv', index=False, header=True, encoding='gbk')
created at:07-01-2021
edited at: 07-02-2021: