统计多面shapefile内不同点shapefile点数的高效实现方法

你原有的实现方案每次处理一类点都需要单独做空间连接、分组统计、全表合并，重复操作多、冗余开销大，处理数十个大型点文件时效率会非常低。以下是针对性优化后的实现方案：

核心优化思路 #

• 仅对多边形做一次预处理，提前构建空间索引、设置唯一主键，避免重复计算
• 批量循环处理所有点文件，每次统计完成后直接追加列到多边形表，避免多次全表合并
• 读取点文件时直接过滤范围外的点，减少无效计算

完整实现代码 #

import pandas as pd
import numpy as np
import geopandas as gpd
from shapely.geometry import Polygon
import os

# ---------------------- 1. 预处理多边形数据（仅执行1次） ----------------------
# 读取你的多面shapefile，替换为实际文件路径
# polygons = gpd.read_file(r"你的多边形文件路径.shp")

# 示例测试数据，实际使用时替换为上面的读取代码
p1 = Polygon( [(0, 0), (0, 1), (0.5, 0.5)])
p2 = Polygon([(.75, 0), (1, .5), (.75, 1), (1.75, 0)])
p3 = Polygon([(2, 0), (2, 1), (3, 0)])
polygons = gpd.GeoDataFrame([p1, p2, p3], columns={'geometry'})
polygons['name'] = ["a", "b", "c"]

# 步骤1：设置多边形唯一标识列为索引（如果有FID/OBJECTID也可以用，确保每个多边形唯一即可）
polygons = polygons.set_index('name', drop=False)
# 步骤2：提前构建空间索引，大幅提升后续空间连接速度
polygons.sindex

# ---------------------- 2. 配置待处理的点文件列表 ----------------------
# 每个元素格式为 (点文件路径, 统计结果列名)，直接把50个点文件依次添加即可
point_config = [
    # (r"E:\gis_data\red_points.shp", "red点数量"),
    # (r"E:\gis_data\green_points.shp", "green点数量"),
    # 剩余点文件按上面格式补充
]

# 示例用内存点数据测试，实际使用时删除下面这行，用上面的文件路径配置
redPoints =  gpd.GeoDataFrame(gpd.points_from_xy( np.random.uniform(0, 3, 50), np.random.uniform(0, 1, 50)), columns = ['geometry'])
greenPoints =  gpd.GeoDataFrame(gpd.points_from_xy( np.random.uniform(1, 3, 50), np.random.uniform(0, 1, 50)), columns = ['geometry'])
point_config = [(redPoints, "red点数量"), (greenPoints, "green点数量")]

# ---------------------- 3. 批量统计点数量 ----------------------
for point_source, count_col in point_config:
    # 读取点文件，mask参数直接过滤掉不在多边形范围内的点，大幅减少计算量
    # 实际使用时取消下面这行注释，替换内存读取逻辑
    # point_gdf = gpd.read_file(point_source, mask=polygons)
    point_gdf = point_source
    
    # 空间连接，用left join保留所有多边形
    sjoin_result = gpd.sjoin(polygons, point_gdf, how='left', predicate='intersects')
    # 按多边形索引统计点数量，生成统计Series
    count_series = sjoin_result.groupby(level=0).size().rename(count_col)
    # 直接追加统计列到多边形表，无需全表merge
    polygons = polygons.join(count_series)
    # 空值（没有点的多边形）填充为0
    polygons[count_col] = polygons[count_col].fillna(0).astype(int)

# 最终结果导出为shapefile
# polygons.to_file(r"E:\gis_data\统计结果.shp", encoding='utf-8')

额外提速建议 #

• 所有数据统一坐标系：处理前确保点和多边形的坐标系一致，不一致先调用.to_crs()转换，避免计算错误
• 优化存储格式：如果点文件过大，可先将shp转为geoparquet格式存储，读取速度比shp快10倍以上，大幅降低IO开销
• 并行处理：因为不同点类型的统计完全独立，可以用concurrent.futures多进程并行处理多个点文件，适合几十上百个点文件的场景
• 复杂多面拆分：如果存在面数极多的复杂Multipolygon，可先调用polygons = polygons.explode(index_parts=True)拆为单面对象再统计，计算效率更高

内容的提问来源于stack exchange，提问作者Gerard