空间分布模式分析——核密度估计

[ZhangAilan]

核密度估计

1）读取地震shp数据，将坐标系从（WGS84，经纬度）转为Web Mercator（EPSG：3857，米），并提取每个点的平面坐标xy
2）根据地震点的范围创建分析网格用于KDE插值，增加50km缓冲区避免边界裁切
3）核密度估计：from sklearn.neighbors import KernelDensity，使用高斯核密度模型对地震点分布进行估计，bandwidth（控制平滑程度，默认选500km），对网格点坐标进行估值输出每个位置的密度值

def create_grid(coords, buffer_km=50, grid_size=100):
    """
    创建规则网格
    
    Args:
        coords: 点坐标数组
        buffer_km: 缓冲区大小（公里）
        grid_size: 网格大小
        
    Returns:
        X_grid, Y_grid: 网格坐标
        grid_points: 网格点坐标数组
        x_min, x_max, y_min, y_max: 边界坐标
    """
    
    # 获取数据边界
    x_min, y_min = coords.min(axis=0)
    x_max, y_max = coords.max(axis=0)
    
    # 扩展边界，增加缓冲区
    buffer = buffer_km * 1000  # 转换为米
    x_min -= buffer
    x_max += buffer
    y_min -= buffer
    y_max += buffer
    
    # 构造栅格网格
    x_grid = np.linspace(x_min, x_max, grid_size)
    y_grid = np.linspace(y_min, y_max, grid_size)
    X_grid, Y_grid = np.meshgrid(x_grid, y_grid)
    
    # 将网格点转换为二维数组
    grid_points = np.vstack([X_grid.ravel(), Y_grid.ravel()]).T
    
    print(f"缓冲区大小: {buffer_km} km")
    print(f"网格分辨率: {grid_size}x{grid_size}")
    print(f"网格点数量: {grid_points.shape[0]}")
    print(f"分析区域X范围: {x_min:.2f} - {x_max:.2f} 米")
    print(f"分析区域Y范围: {y_min:.2f} - {y_max:.2f} 米")
    
    return X_grid, Y_grid, grid_points, x_min, x_max, y_min, y_max


def perform_kde_analysis(coords, grid_points, X_grid, bandwidth_km=500):
    """
    执行核密度估计分析
    
    Args:
        coords: 点坐标数组
        grid_points: 网格点坐标数组
        X_grid: X方向网格
        bandwidth_km: 带宽（公里）
        
    Returns:
        density_grid: 密度网格
        density: 密度值数组
        kde: KDE模型
    """
    
    # 设置带宽（转换为米）
    bandwidth = bandwidth_km * 1000
    print(f"使用的带宽: {bandwidth_km} km ({bandwidth:.0f} 米)")
    
    # 创建KDE模型
    kde = KernelDensity(kernel='gaussian', bandwidth=bandwidth)
    kde.fit(coords)
    
    # 计算网格点的密度值
    log_density = kde.score_samples(grid_points)
    density = np.exp(log_density)
    
    # 将密度值重新整形为网格形状
    density_grid = density.reshape(X_grid.shape)
    
    print(f"密度值范围: {density.min():.2e} - {density.max():.2e}")
    print("核密度估计计算完成")
    
    return density_grid, density, kde