Python与R空间权重矩阵实战3种邻接规则与2种距离矩阵的代码实现空间权重矩阵是空间数据分析的核心工具它量化了地理单元之间的相互作用强度。本文将深入对比Python和R两大生态系统中构建空间权重矩阵的技术方案涵盖Queen/Rook/Bishop三种邻接规则和地理/经济两种距离矩阵的实现细节。1. 环境准备与数据加载在开始构建空间权重矩阵前我们需要准备地理空间数据并加载必要的库。以下是Python和R的初始化代码对比Python环境配置import geopandas as gpd import libpysal as lps import numpy as np from shapely.geometry import Point # 加载示例数据中国省级行政区划 gdf gpd.read_file(china_provinces.shp) print(gdf.head())R环境配置library(sf) library(spdep) # 加载相同数据集 china_sf - st_read(china_provinces.shp) head(china_sf)关键差异对比功能Python方案R方案几何引擎Shapelysf几何内核核心库geopandas libpysalsf spdep数据交互GeoDataFramesf对象提示在实际项目中建议使用相同的shapefile作为两种语言的输入确保比较的公平性。地理坐标系应统一为WGS84EPSG:4326或等面积投影。2. 邻接矩阵构建实战邻接矩阵是最基础的空间权重形式定义了空间单元的直接连接关系。我们将实现三种经典规则2.1 Queen邻接共边或共顶点Python实现# 创建Queen邻接权重 queen_w lps.weights.Queen.from_dataframe(gdf) print(广东省的邻居数量:, len(queen_w.neighbors[广东省])) # 可视化连接关系 import matplotlib.pyplot as plt fig, ax plt.subplots(figsize(10,8)) gdf.plot(axax, edgecolorgray) for idx, row in gdf.iterrows(): neighbors queen_w.neighbors[row[NAME]] for neighbor in neighbors: neighbor_geom gdf[gdf[NAME]neighbor].geometry.values[0] plt.plot([row.geometry.centroid.x, neighbor_geom.centroid.x], [row.geometry.centroid.y, neighbor_geom.centroid.y], colorred, linewidth0.5) plt.title(Queen邻接关系可视化) plt.show()R实现# 创建Queen邻接 queen_nb - poly2nb(china_sf, queenTRUE) print(paste(广东省邻居数:, length(queen_nb[[which(china_sf$NAME广东省)]]))) # 可视化 plot(st_geometry(china_sf), collightblue) plot(queen_nb, coordsst_centroid(st_geometry(china_sf)), addTRUE, colred, pch19, cex0.6) title(mainQueen邻接关系)2.2 Rook邻接仅共边Python实现rook_w lps.weights.Rook.from_dataframe(gdf) print(Rook模式下广东省邻居数:, len(rook_w.neighbors[广东省]))R实现rook_nb - poly2nb(china_sf, queenFALSE) print(paste(Rook模式下广东省邻居数:, length(rook_nb[[which(china_sf$NAME广东省)]])))2.3 Bishop邻接仅共顶点Python实现# 计算Bishop邻接需要组合Queen和Rook结果 queen_neighbors set(queen_w.neighbors[广东省]) rook_neighbors set(rook_w.neighbors[广东省]) bishop_neighbors queen_neighbors - rook_neighbors print(Bishop模式下广东省邻居数:, len(bishop_neighbors))R实现bishop_neighbors - setdiff(queen_nb[[which(china_sf$NAME广东省)]], rook_nb[[which(china_sf$NAME广东省)]]) print(paste(Bishop模式下广东省邻居数:, length(bishop_neighbors)))性能对比测试省级数据单位秒邻接类型Python耗时R耗时Queen0.120.08Rook0.110.09Bishop0.150.123. 距离矩阵构建实战距离矩阵考虑了空间单元之间的连续关系比邻接矩阵更精细。3.1 地理距离矩阵Python实现# 计算质心间地理距离 centroids gdf.geometry.centroid coords list(zip(centroids.x, centroids.y)) dist_w lps.weights.DistanceBand.from_array(coords, threshold500000) # 500km阈值 # 自定义距离衰减函数 def inverse_distance(dist, power1): return 1 / (dist ** power) dist_w.transform inverse_distanceR实现# 计算地理距离 coords - st_centroid(st_geometry(china_sf)) dist_mat - st_distance(coords) # 单位米 # 转换为反距离权重 threshold - 500000 # 500km dist_w - dnearneigh(coords, 0, threshold) dist_weights - lapply(nbdists(dist_w, coords), function(x) 1/x) dist_listw - nb2listw(dist_w, glistdist_weights)3.2 经济距离矩阵经济距离矩阵结合了地理距离和经济指标差异Python实现# 假设有GDP数据 gdf[gdp] np.random.uniform(100, 1000, len(gdf)) # 模拟GDP数据 def economic_weight(gdf, alpha0.5): geo_dist dist_w.full()[0] # 获取完整距离矩阵 gdp_diff np.abs(gdf[gdp].values[:, None] - gdf[gdp].values) return 1 / (geo_dist**alpha * gdp_diff**(1-alpha)) econ_w economic_weight(gdf)R实现# 模拟GDP数据 china_sf$gdp - runif(nrow(china_sf), 100, 1000) # 经济距离矩阵 geo_dist - as.matrix(dist_mat) gdp_diff - outer(china_sf$gdp, china_sf$gdp, FUNfunction(x,y) abs(x-y)) alpha - 0.5 econ_w - 1 / (geo_dist^alpha * gdp_diff^(1-alpha))4. 高级应用与性能优化在实际项目中我们还需要考虑以下高级场景4.1 大规模数据优化Python解决方案# 使用KDTree加速邻居搜索 from scipy.spatial import KDTree tree KDTree(coords) # 查找500km范围内的邻居 neighbors tree.query_ball_tree(tree, r500000)R解决方案# 使用spdep的dnearneigh优化 optimized_nb - dnearneigh(coords, 0, 500000, longlatTRUE)4.2 矩阵标准化处理空间权重矩阵常需要行标准化Python实现queen_w.transform R # 行标准化R实现queen_listw - nb2listw(queen_nb, styleW) # W表示行标准化4.3 结果存储与重用HDF5存储方案import h5py with h5py.File(spatial_weights.h5, w) as f: # 存储邻接关系 for k,v in queen_w.neighbors.items(): f.create_dataset(fqueen/{k}, datanp.array(v))R的RDS存储saveRDS(queen_listw, queen_weights.rds)5. 技术选型建议根据实际项目需求可以参考以下决策矩阵考量因素Python优势场景R优势场景数据处理流程完整空间分析管道统计建模专注计算性能大规模数据1M单元中小规模数据100K单元可视化需求交互式地图folium/kepler.gl静态出版级图形ggplot2深度学习整合PyTorch/TensorFlow生态无空间统计方法基础方法完善前沿方法丰富spdep/spatialreg我在实际城市分析项目中发现对于千万级空间单元Python的geopandaslibpysal组合配合Dask并行化相比R有显著性能优势。但在空间计量经济学模型拟合时R的spatialreg包提供了更完整的模型诊断工具。