Geopandas创建带唯一值统计的点密度栅格 求类似ESRI PointDensity_sa函数

一、彻底优化原循环逻辑 #

你的双重循环逐个像元查询的方式效率极低，核心问题是每次查询都要遍历整个数据集。可以通过空间索引+批量处理大幅提速：

1. 预生成所有像元中心的点集 #

先把所有像元的中心坐标转换成GeoDataFrame，这样可以利用Geopandas的空间索引批量查询：

import geopandas as gpd
import numpy as np
from shapely.geometry import Point

# 生成像元中心坐标
x_coords = np.arange(bbox[0], bbox[2], cellsize) + cellsize/2
y_coords = np.arange(bbox[1], bbox[3], cellsize) + cellsize/2
xx, yy = np.meshgrid(x_coords, y_coords)
centers = [Point(x, y) for x, y in zip(xx.flatten(), yy.flatten())]
grid_centers = gpd.GeoDataFrame({'geometry': centers}, crs=bb.crs)

2. 用空间索引快速查询邻域 #

利用GeoDataFrame的sindex（R-tree空间索引），批量查询每个像元中心缓冲区范围内的点，再统计唯一flkey数量：

# 构建原始数据的空间索引
sindex = bb.sindex

# 初始化结果数组
mappa = np.zeros((len(y_coords), len(x_coords)))

for idx, center in enumerate(grid_centers.geometry):
    # 创建查询缓冲区
    buffer = center.buffer(radius)
    # 用空间索引获取候选点（大幅缩小查询范围）
    possible_matches_idx = list(sindex.intersection(buffer.bounds))
    possible_matches = bb.iloc[possible_matches_idx]
    # 精确筛选在缓冲区内的点
    actual_matches = possible_matches[possible_matches.geometry.within(buffer)]
    # 统计唯一flkey数量
    unique_count = actual_matches['flkey'].nunique()
    # 映射回栅格数组
    row = idx // len(x_coords)
    col = idx % len(x_coords)
    mappa[row, col] = unique_count / period

这种方式通过空间索引先过滤掉绝大多数无关点，再做精确的空间判断，比原循环快几十到上百倍。

二、类似ESRI PointDensity_sa的替代工具 #

Python生态里有几个可以实现类似点密度分析的工具：

• PySAL库：pysal.esda.KernelDensity可以生成核密度栅格，支持自定义带宽（对应你的radius），可基于它的逻辑修改为统计唯一flkey数量。
• Scipy KDTree：如果是小范围平面坐标（经纬度需转投影），可以用scipy.spatial.KDTree快速查找每个像元中心半径内的点，再统计唯一值：

from scipy.spatial import KDTree

# 转换为投影坐标（经纬度转平面，避免球面距离误差）
bb_proj = bb.to_crs(epsg=3857)
grid_centers_proj = grid_centers.to_crs(epsg=3857)

# 构建KDTree
tree = KDTree(np.array([bb_proj.geometry.x, bb_proj.geometry.y]).T)
# 批量查询每个中心半径内的点索引
indices = tree.query_ball_point(np.array([grid_centers_proj.geometry.x, grid_centers_proj.geometry.y]).T, r=radius)

# 统计唯一flkey
for idx, idx_list in enumerate(indices):
    if idx_list:
        unique_count = bb_proj.iloc[idx_list]['flkey'].nunique()
    else:
        unique_count = 0
    row = idx // len(x_coords)
    col = idx % len(x_coords)
    mappa[row, col] = unique_count / period

三、只处理高数据密度区域的思路 #

1. 先做低分辨率密度概查：用10倍于目标cellsize的尺寸生成粗略密度栅格，筛选出密度大于阈值的区域。
2. 裁剪目标栅格范围：只对高密区域对应的目标分辨率像元进行详细计算，跳过空白或低密度区域。
3. 代码示例：

# 第一步：生成低分辨率密度栅格
low_cellsize = cellsize * 10
low_x = np.arange(bbox[0], bbox[2], low_cellsize) + low_cellsize/2
low_y = np.arange(bbox[1], bbox[3], low_cellsize) + low_cellsize/2
low_xx, low_yy = np.meshgrid(low_x, low_y)
low_centers = [Point(x,y) for x,y in zip(low_xx.flatten(), low_yy.flatten())]
low_grid = gpd.GeoDataFrame({'geometry': low_centers}, crs=bb.crs)

# 计算低分辨率密度
low_density = np.zeros((len(low_y), len(low_x)))
sindex = bb.sindex
for idx, center in enumerate(low_grid.geometry):
    buffer = center.buffer(radius)
    possible_matches = bb.iloc[list(sindex.intersection(buffer.bounds))]
    actual_matches = possible_matches[possible_matches.geometry.within(buffer)]
    low_density[idx//len(low_x), idx%len(low_x)] = len(actual_matches)

# 第二步：筛选高密区域（阈值自定义，比如大于5个点）
high_density_mask = low_density > 5

# 第三步：只处理高密区域对应的目标像元
for low_row in range(len(low_y)):
    for low_col in range(len(low_x)):
        if high_density_mask[low_row, low_col]:
            # 计算当前低分辨率像元对应的目标像元范围
            start_x = bbox[0] + low_col * low_cellsize
            end_x = start_x + low_cellsize
            start_y = bbox[1] + low_row * low_cellsize
            end_y = start_y + low_cellsize
            # 遍历该范围内的目标像元
            target_x_coords = np.arange(start_x, end_x, cellsize) + cellsize/2
            target_y_coords = np.arange(start_y, end_y, cellsize) + cellsize/2
            for tx in target_x_coords:
                for ty in target_y_coords:
                    col = int((tx - bbox[0])/cellsize)
                    row = int((ty - bbox[1])/cellsize)
                    if 0<=row<len(y_coords) and 0<=col<len(x_coords):
                        center = Point(tx, ty)
                        buffer = center.buffer(radius)
                        possible_matches = bb.iloc[list(sindex.intersection(buffer.bounds))]
                        actual_matches = possible_matches[possible_matches.geometry.within(buffer)]
                        mappa[row, col] = actual_matches['flkey'].nunique() / period

这样可以跳过大部分无数据或低密度区域，进一步减少计算量。

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