Python/Geopandas实现空间密度（空间平滑）方法咨询

嘿，这个需求我太熟悉了——把离散的1/0存在/缺失点转换成指定半径内的密度分布，本质就是**核密度估计(Kernel Density Estimation, KDE)**或者基于固定半径的点密度计算，用Geopandas配合几个空间工具就能高效搞定，给你整理几个实用方案：

方案一：Geopandas + Scipy 核密度估计（灵活易上手） #

这个方案适合中小规模数据集，步骤清晰，可视化方便。核心思路是先把地理坐标转成投影坐标系（保证距离计算准确），再用Scipy的gaussian_kde计算密度，最后转回地理坐标可视化。

代码示例： #

import geopandas as gpd
import numpy as np
from scipy.stats import gaussian_kde
import matplotlib.pyplot as plt

# 1. 加载你的点数据（假设是包含geometry列的GeoDataFrame，值为1/0）
gdf = gpd.read_file("your_points.shp")  # 或者从CSV加载后转换为GeoDataFrame

# 2. 关键：转换为投影坐标系（比如UTM，根据数据所在区域选对应EPSG）
# 示例：假设数据在WGS84（EPSG:4326），转成UTM Zone 10N（EPSG:32610）
gdf_proj = gdf.to_crs(epsg=32610)

# 3. 提取点的坐标数组，只保留存在的点（值为1的）
points = np.array([gdf_proj.geometry.x, gdf_proj.geometry.y]).T
existing_points = points[gdf_proj['your_value_column'] == 1]

# 4. 设置核密度的带宽（对应你要的指定半径，单位是投影坐标系的单位，比如米）
bandwidth = 1000  # 比如1000米

# 5. 计算核密度
kde = gaussian_kde(existing_points.T, bw_method=bandwidth / np.std(existing_points))

# 6. 创建网格用于绘制密度图（可以用原始数据的范围生成）
xmin, ymin, xmax, ymax = gdf_proj.total_bounds
x_grid = np.linspace(xmin, xmax, 200)
y_grid = np.linspace(ymin, ymax, 200)
X, Y = np.meshgrid(x_grid, y_grid)
positions = np.vstack([X.ravel(), Y.ravel()])
density = kde(positions).reshape(X.shape)

# 7. 转回地理坐标系并可视化
# 把网格转成GeoDataFrame
grid_gdf = gpd.GeoDataFrame(
    geometry=gpd.points_from_xy(X.ravel(), Y.ravel()),
    crs=gdf_proj.crs
).to_crs(gdf.crs)
grid_gdf['density'] = density.ravel()

# 绘制密度图
fig, ax = plt.subplots(figsize=(10, 8))
gdf.plot(ax=ax, color='gray', markersize=5, alpha=0.5)
grid_gdf.plot(ax=ax, column='density', cmap='viridis', markersize=10, alpha=0.8, legend=True)
plt.title("存在点的核密度分布（半径1000米）")
plt.show()

方案二：Geopandas + Rasterio 栅格化密度（适合大数据） #

如果你的数据集很大（比如几十万上百万个点），栅格化的方式会更高效，直接把每个栅格单元内的存在点数除以单元面积（或者按指定半径统计邻域点数）。

代码示例（固定半径邻域密度）： #

import geopandas as gpd
import rasterio
from rasterio.features import rasterize
from scipy.ndimage import uniform_filter

# 1. 加载数据并转投影
gdf = gpd.read_file("your_points.shp")
gdf_proj = gdf.to_crs(epsg=32610)
existing_points = gdf_proj[gdf_proj['your_value_column'] == 1]

# 2. 设置栅格参数（分辨率比如100米，对应你需要的精度）
resolution = 100
xmin, ymin, xmax, ymax = gdf_proj.total_bounds
width = int((xmax - xmin) / resolution)
height = int((ymax - ymin) / resolution)

# 3. 把存在点栅格化（每个点所在栅格设为1，其他0）
transform = rasterio.transform.from_origin(xmin, ymax, resolution, resolution)
raster = rasterize(
    [(geom, 1) for geom in existing_points.geometry],
    out_shape=(height, width),
    transform=transform,
    fill=0,
    dtype=np.float32
)

# 4. 用均匀滤波计算指定半径内的密度（窗口大小=半径/分辨率）
radius = 1000
window_size = int(radius / resolution)
density_raster = uniform_filter(raster, size=window_size) / (np.pi * radius**2)  # 单位：点/平方米

# 5. 保存或可视化栅格
with rasterio.open(
    'density_raster.tif',
    'w',
    driver='GTiff',
    height=height,
    width=width,
    count=1,
    dtype=density_raster.dtype,
    crs=gdf_proj.crs,
    transform=transform,
) as dst:
    dst.write(density_raster, 1)

# 可视化的话可以用rasterio.plot.show()

方案三：PySal + Geopandas 专业空间统计（高效且功能丰富） #

PySal是专门做空间数据分析的库，它的KernelDensity类可以直接处理地理数据，还支持多种核函数和带宽选择方法，适合需要更严谨空间统计的场景。

代码示例： #

import geopandas as gpd
from pysal.explore.esda import KernelDensity

# 1. 加载数据（PySal支持直接用GeoDataFrame）
gdf = gpd.read_file("your_points.shp")
existing_points = gdf[gdf['your_value_column'] == 1]

# 2. 设置参数：带宽（半径，单位米，注意这里需要数据是投影坐标系）
gdf_proj = existing_points.to_crs(epsg=32610)
coords = list(zip(gdf_proj.geometry.x, gdf_proj.geometry.y))
kd = KernelDensity(coords, bandwidth=1000, fixed=True)

# 3. 计算每个点的密度值（或者生成网格密度）
# 如果要生成网格密度，可以用PySal的grid模块
from pysal.lib import grid
bbox = gdf_proj.total_bounds
grid_df = grid.Grid(bbox, 100)  # 100米分辨率的网格
grid_coords = list(zip(grid_df.centroid.x, grid_df.centroid.y))
grid_df['density'] = kd.pdf(grid_coords)

# 4. 转回地理坐标系可视化
grid_df = grid_df.to_crs(gdf.crs)
fig, ax = plt.subplots(figsize=(10,8))
gdf.plot(ax=ax, color='gray', markersize=5)
grid_df.plot(ax=ax, column='density', cmap='viridis', alpha=0.7, legend=True)
plt.title("PySal核密度分布")
plt.show()

关键注意事项 #

• 投影转换：一定要把经纬度（EPSG:4326）转换成投影坐标系（比如UTM），因为经纬度的单位是度，计算距离会严重失真，所有距离相关的参数（半径、带宽）都要对应投影坐标系的单位（通常是米）。
• 带宽/半径选择：如果不确定合适的半径，可以尝试用数据的平均最近邻距离作为参考，或者用交叉验证来优化带宽。
• 性能优化：大数据集优先用方案二（栅格化）或方案三（PySal，底层用C加速），方案一的Scipy KDE在数据量超过10万时会变慢。

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