Python与R空间权重矩阵对比：3种邻接规则与2种距离矩阵的代码实现

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)

关键差异对比：

提示：在实际项目中，建议使用相同的shapefile作为两种语言的输入，确保比较的公平性。地理坐标系应统一为WGS84（EPSG: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(ax=ax, edgecolor='gray') 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], color='red', linewidth=0.5) plt.title('Queen邻接关系可视化') plt.show()

R实现：

# 创建Queen邻接 queen_nb <- poly2nb(china_sf, queen=TRUE) print(paste("广东省邻居数:", length(queen_nb[[which(china_sf$NAME=="广东省")]]))) # 可视化 plot(st_geometry(china_sf), col='lightblue') plot(queen_nb, coords=st_centroid(st_geometry(china_sf)), add=TRUE, col='red', pch=19, cex=0.6) title(main="Queen邻接关系")

2.2 Rook邻接（仅共边）

Python实现：

rook_w = lps.weights.Rook.from_dataframe(gdf) print("Rook模式下广东省邻居数:", len(rook_w.neighbors['广东省']))

R实现：

rook_nb <- poly2nb(china_sf, queen=FALSE) 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)))

性能对比测试（省级数据，单位：秒）：

3. 距离矩阵构建实战

距离矩阵考虑了空间单元之间的连续关系，比邻接矩阵更精细。

3.1 地理距离矩阵

Python实现：

# 计算质心间地理距离 centroids = gdf.geometry.centroid coords = list(zip(centroids.x, centroids.y)) dist_w = lps.weights.DistanceBand.from_array(coords, threshold=500000) # 500km阈值 # 自定义距离衰减函数 def inverse_distance(dist, power=1): return 1 / (dist ** power) dist_w.transform = inverse_distance

R实现：

# 计算地理距离 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, glist=dist_weights)

3.2 经济距离矩阵

经济距离矩阵结合了地理距离和经济指标差异：

Python实现：

# 假设有GDP数据 gdf['gdp'] = np.random.uniform(100, 1000, len(gdf)) # 模拟GDP数据 def economic_weight(gdf, alpha=0.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, FUN=function(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, r=500000)

R解决方案：

# 使用spdep的dnearneigh优化 optimized_nb <- dnearneigh(coords, 0, 500000, longlat=TRUE)

4.2 矩阵标准化处理

空间权重矩阵常需要行标准化：

Python实现：

queen_w.transform = 'R' # 行标准化

R实现：

queen_listw <- nb2listw(queen_nb, style="W") # 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(f'queen/{k}', data=np.array(v))

R的RDS存储：

saveRDS(queen_listw, "queen_weights.rds")

5. 技术选型建议

根据实际项目需求，可以参考以下决策矩阵：

我在实际城市分析项目中发现，对于千万级空间单元，Python的geopandas+libpysal组合配合Dask并行化，相比R有显著性能优势。但在空间计量经济学模型拟合时，R的spatialreg包提供了更完整的模型诊断工具。