如何在GeoPandas中按指定规模聚合正方形网格面要素
2026-7-3
高效滑动聚合规则网格的GeoPandas方案
问题描述
我用GeoPandas处理空间数据，现有一个大行数的GeoDataFrame，每行是EPSG:4326投影下的等尺寸地面正方形，构成特定区域的规则网格。需要根据输入值将网格按滑动窗口式的2×2（4个一组）或3×3（9个一组）聚合：
- 输入值为2时，以每个正方形为窗口左下角，合并右侧、上方及右上角的正方形，仅保留窗口内4个正方形都存在的聚合结果，最终得到7个新正方形
- 输入值为3时，合并3×3范围内的9个正方形，仅保留所有9个正方形都存在的窗口，最终得到2个新正方形（边缘无法凑组的窗口舍弃）
示例网格代码：
import shapely import geopandas as gpd a0 = shapely.geometry.Point([0, 0]).buffer(0.5, cap_style = 'square') a1 = shapely.affinity.translate(a0, xoff=1) a2 = shapely.affinity.translate(a0, xoff=2) a3 = shapely.affinity.translate(a0, xoff=3) a4 = shapely.affinity.translate(a0, xoff=4) b0 = shapely.affinity.translate(a0, yoff=1) b1 = shapely.affinity.translate(b0, xoff=1) b2 = shapely.affinity.translate(b0, xoff=2) b3 = shapely.affinity.translate(b0, xoff=3) b4 = shapely.affinity.translate(b0, xoff=4) c0 = shapely.affinity.translate(b0, yoff=1) c1 = shapely.affinity.translate(c0, xoff=1) c2 = shapely.affinity.translate(c0, xoff=2) c3 = shapely.affinity.translate(c0, xoff=3) liste_geo = [a0, a1, a2, a3, a4, b0, b1, b2, b3, b4, c0, c1, c2, c3] column = [1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4] row = [0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 2, 2, 2, 2] gdf = gpd.GeoDataFrame({'column' : column, 'row' : row, 'geometry' : liste_geo}, crs = 'epsg:4326', index = range(len(liste_geo))) gdf.plot(facecolor = 'none')
我自己实现的函数在处理10万行左右的大型GeoDataFrame时速度极慢，现有代码：
def define_spot_grid(gdf_base_grid, size_base, size_swath): """ gdf_base_grid : our initial grid size_base : size of the initial square size_swath : size of our bigger square """ rapport = size_swath/size_base unique_column = gdf_base_grid['column'].unique() unique_row = gdf_base_grid['row'].unique() liste_bigger_spot = [] for k, v in gdf_base_grid.sort_values(by = ['row', 'column']).iterrows(): # actually for every square i have a row and column information current_row, current_column = v['row'], v['column'] bound_row, bound_column = v['row'] + rapport - 1, v['column'] + rapport - 1 spot_bl = gdf_base_grid[(gdf_base_grid['row'] == current_row) & (gdf_base_grid['column'] == current_column)] spot_tl = gdf_base_grid[(gdf_base_grid['row'] == bound_row) & (gdf_base_grid['column'] == current_column)] spot_tr = gdf_base_grid[(gdf_base_grid['row'] == bound_row) & (gdf_base_grid['column'] == bound_column)] spot_br = gdf_base_grid[(gdf_base_grid['row'] == current_row) & (gdf_base_grid['column'] == bound_column)] # checking if all these squares exist if (not spot_bl.empty) and (not spot_tl.empty) and (not spot_tr.empty) and (not spot_br.empty): point_bl = spot_bl.unary_union.centroid point_tl = spot_tl.unary_union.centroid point_tr = spot_tr.unary_union.centroid point_br = spot_br.unary_union.centroid # we define a spot with 4 points bigger_spot = shapely.geometry.Polygon([point_bl, point_tl, point_tr, point_br]) gdf_bigger_spot = gpd.GeoDataFrame({'geometry' : bigger_spot}, crs = 'epsg:4326', index = [0]) # we check how many spot our bigger spot intersects using a spatial join sjoin = gpd.sjoin(gdf_bigger_spot, gdf_base_grid, how = 'right', predicate = 'intersects').dropna() # example : if rapport = 3 we need to get 9 squares, if it is less we do not keep if len(sjoin) == rapport**2: liste_bigger_spot.append(sjoin.unary_union) gdf_bigger = gpd.GeoDataFrame({'geometry' : liste_bigger_spot}, crs = 'epsg:4326', index = range(len(liste_bigger_spot))) return gdf_bigger
求更高效的解决方案。
高效解决方案
原方法性能瓶颈在于：逐行循环所有网格单元、多次重复的空间查询、不必要的空间连接操作。利用网格的规则行列属性，可以完全规避这些耗时操作，将时间复杂度从O(N)（N为单元数）降低到O(W)（W为有效窗口数，远小于N）。
核心思路
- 预存网格单元索引：将所有
(row, column)对存入集合，实现O(1)时间复杂度的存在性检查 - 生成有效窗口范围：仅遍历可能形成完整聚合窗口的起始行列，避免无效循环
- 批量筛选合并几何：利用布尔索引快速定位窗口内的所有单元，合并几何形状
实现代码
import geopandas as gpd import shapely def define_spot_grid_fast(gdf_base_grid, size_base, size_swath): """ 高效滑动聚合规则网格的函数 gdf_base_grid : 初始网格GeoDataFrame，需包含row、column列 size_base : 基础正方形的尺寸 size_swath : 聚合后大正方形的尺寸 """ rapport = int(size_swath / size_base) if rapport not in (2, 3): raise ValueError("仅支持2×2或3×3的滑动聚合") # 复制原数据避免修改输入 gdf = gdf_base_grid.copy() # 创建(row, column)的集合，用于快速检查单元是否存在 cell_coords = set(gdf[['row', 'column']].itertuples(index=False, name=None)) # 获取网格的最大行列值，确定窗口的有效起始范围 max_row = gdf['row'].max() max_col = gdf['column'].max() valid_start_rows = [r for r in gdf['row'].unique() if r + rapport - 1 <= max_row] valid_start_cols = [c for c in gdf['column'].unique() if c + rapport - 1 <= max_col] aggregated_geoms = [] # 遍历所有可能的有效窗口起始点 for start_r in valid_start_rows: end_r = start_r + rapport - 1 window_rows = range(start_r, end_r + 1) for start_c in valid_start_cols: end_c = start_c + rapport - 1 window_cols = range(start_c, end_c + 1) # 检查窗口内所有单元是否都存在 all_cells_exist = all((r, c) in cell_coords for r in window_rows for c in window_cols) if all_cells_exist: # 筛选窗口内的所有几何并合并 mask = (gdf['row'].between(start_r, end_r)) & (gdf['column'].between(start_c, end_c)) merged_geometry = shapely.ops.unary_union(gdf.loc[mask, 'geometry']) aggregated_geoms.append(merged_geometry) # 生成结果GeoDataFrame并设置CRS result_gdf = gpd.GeoDataFrame({'geometry': aggregated_geoms}, crs=gdf_base_grid.crs) return result_gdf
性能优化点
- 减少循环次数：从遍历所有单元改为遍历有效窗口起始点，循环次数大幅减少
- 快速存在性检查：集合查询替代多次布尔索引，速度提升数十倍
- 避免空间操作：用行列属性筛选替代空间连接，完全规避耗时的空间计算
验证示例
用提供的示例数据测试：
# 测试2×2滑动聚合 gdf_2x2 = define_spot_grid_fast(gdf, size_base=1, size_swath=2) print(f"2×2聚合结果行数：{len(gdf_2x2)}") # 输出7，符合预期 gdf_2x2.plot(facecolor='none', edgecolor='red') # 测试3×3滑动聚合 gdf_3x3 = define_spot_grid_fast(gdf, size_base=1, size_swath=3) print(f"3×3聚合结果行数：{len(gdf_3x3)}") # 输出2，符合预期 gdf_3x3.plot(facecolor='none', edgecolor='blue')
内容的提问来源于stack exchange，提问作者Magic Mushroom
相关产品推荐
相关产品推荐
