您好，登录后才能下订单哦！
密码登录
登录注册
点击 登录注册 即表示同意《亿速云用户服务条款》
# 如何利用Python计算全局莫兰指数
## 一、空间自相关与莫兰指数简介
### 1.1 空间自相关的概念
空间自相关（Spatial Autocorrelation）是地理信息系统和空间统计学中的核心概念，用于描述地理空间上邻近位置的观测值之间的相关性。当相邻区域具有相似属性值时，表现为正空间自相关；当相邻区域属性差异显著时，表现为负空间自相关；随机分布则表现为无空间自相关。
### 1.2 全局莫兰指数的定义
全局莫兰指数（Global Moran's I）由Patrick Alfred Pierce Moran于1950年提出，是最常用的空间自相关度量指标之一。其计算公式为：
$$
I = \frac{n}{\sum_{i=1}^{n}\sum_{j=1}^{n}w_{ij}} \cdot \frac{\sum_{i=1}^{n}\sum_{j=1}^{n}w_{ij}(x_i-\bar{x})(x_j-\bar{x})}{\sum_{i=1}^{n}(x_i-\bar{x})^2}
$$
其中：
- $n$：空间单元数量
- $x_i, x_j$：位置$i$和$j$的观测值
- $\bar{x}$：所有观测值的均值
- $w_{ij}$：空间权重矩阵元素
### 1.3 指数取值范围与解释
| 莫兰指数值 | 空间自相关类型 | 地理分布特征 |
|------------|----------------|------------------------|
| I ≈ 1 | 强正相关 | 相似值聚集（高-高/低-低）|
| I ≈ 0 | 无相关性 | 随机分布 |
| I ≈ -1 | 强负相关 | 相异值聚集（高-低相邻） |
## 二、Python计算环境准备
### 2.1 必需库及其功能
```python
# 核心计算库
import numpy as np
import pandas as pd
import libpysal as lps # 空间权重矩阵构建
from esda.moran import Moran # 莫兰指数计算
import geopandas as gpd # 空间数据处理
# 可视化库
import matplotlib.pyplot as plt
import seaborn as sns
from splot.esda import plot_moran # 莫兰指数可视化
# 使用conda安装（推荐）
conda install -c conda-forge libpysal esda geopandas
# 使用pip安装
pip install libpysal esda geopandas matplotlib seaborn
# 加载示例数据集（以美国各州人口数据为例）
us_states = gpd.read_file(lps.examples.get_path('us48.shp'))
us_income = pd.read_csv(lps.examples.get_path('usjoin.csv'))
# 数据合并
merged_data = us_states.merge(us_income, left_on='STATE_NAME', right_on='Name')
| 权重类型 | 适用场景 | 构建方法 | 特点 |
|---|---|---|---|
| 邻接权重 | 行政区划数据 | Queen/Rook相邻规则 | 简单直观，边界定义明确 |
| 距离权重 | 点数据或连续空间 | 固定距离/反距离 | 反映实际空间相互作用 |
| K最近邻权重 | 点密度不均匀分布 | 每个点的K个最近邻居 | 保证每个单元有相同邻居数 |
| 核密度权重 | 复杂空间交互模式 | 高斯核/二次核函数 | 平滑过渡，参数敏感 |
# 1. Queen邻接权重
w_queen = lps.weights.Queen.from_dataframe(merged_data)
# 2. 反距离权重（阈值100km）
points = merged_data.geometry.centroid
coords = [(x,y) for x,y in zip(points.x, points.y)]
w_dist = lps.weights.DistanceBand.from_array(coords, threshold=100, binary=False)
# 3. K最近邻权重（K=4）
w_knn = lps.weights.KNN.from_array(coords, k=4)
# 权重矩阵标准化（行标准化）
w_queen.transform = 'R'
# 绘制权重连接图
fig, ax = plt.subplots(figsize=(12,10))
merged_data.plot(ax=ax, edgecolor='grey', facecolor='white')
w_queen.plot(merged_data, ax=ax,
edge_kws=dict(linewidth=1.5, color='red', alpha=0.6),
node_kws=dict(marker=''))
plt.title("Queen Contiguity Spatial Weights")
plt.show()
# 选择分析变量（这里使用州收入中位数）
y = merged_data['median_income']
# 计算莫兰指数
moran = Moran(y, w_queen)
# 输出结果
print(f"Moran's I 值: {moran.I:.4f}")
print(f"标准化统计量: {moran.z_rand:.4f}")
print(f"P值: {moran.p_rand:.4f}")
# 结果解读
if moran.p_rand < 0.05:
if moran.I > 0:
print("在0.05水平下显著正相关")
else:
print("在0.05水平下显著负相关")
else:
print("未表现出显著空间自相关")
# 莫兰散点图
from splot.esda import moran_scatterplot
fig, ax = moran_scatterplot(moran, aspect_equal=True)
plt.title(f"Moran Scatterplot (I = {moran.I:.3f})")
plt.show()
# 空间自相关分布图
from splot.esda import plot_local_autocorrelation
plot_local_autocorrelation(moran, merged_data, 'median_income')
plt.show()
# 蒙特卡洛模拟（999次 permutation）
from esda.moran import Moran
moran_mc = Moran(y, w_queen, permutations=999)
# 绘制模拟分布
plt.figure(figsize=(10,6))
plt.hist(moran_mc.sim, bins=50, density=True)
plt.axvline(x=moran_mc.I, color='red', linestyle='--')
plt.title(f'Moran I Distribution (p={moran_mc.p_sim:.4f})')
plt.xlabel("Moran's I Value")
plt.ylabel("Density")
plt.show()
# 加载中国城市数据
china_cities = gpd.read_file('china_cities.shp')
gdp_data = pd.read_excel('city_gdp_2020.xlsx')
# 数据合并与清洗
china_cities['city'] = china_cities['NAME'].str.replace('市','')
merged = china_cities.merge(gdp_data, on='city')
merged = merged.dropna(subset=['gdp'])
# 坐标参考系统转换
merged = merged.to_crs('EPSG:4326')
# 使用KNN+反距离组合权重
coords = [(x,y) for x,y in zip(merged.geometry.centroid.x,
merged.geometry.centroid.y)]
# 先构建KNN基础连接
w_knn = lps.weights.KNN.from_array(coords, k=8)
# 转换为反距离权重
distances = lps.weights.DistanceBand.from_array(coords, threshold=500)
w_comb = lps.weights.adaptive_kernelW(w_knn, bandwidth=5)
moran_gdp = Moran(merged['gdp'], w_comb, permutations=9999)
print(f"""中国城市GDP空间自相关分析结果：
Moran's I: {moran_gdp.I:.3f}
Z-score: {moran_gdp.z_rand:.3f}
P-value: {moran_gdp.p_sim:.5f}
""")
# 绘制热力图
fig, ax = plt.subplots(1, 2, figsize=(16,6))
merged.plot(column='gdp', cmap='OrRd', scheme='quantiles',
legend=True, ax=ax[0])
ax[0].set_title('中国城市GDP分布')
moran_scatterplot(moran_gdp, ax=ax[1])
ax[1].set_title('Moran Scatterplot')
plt.tight_layout()
NaN值错误：
# 检查并处理缺失值
print(y.isnull().sum())
y = y.fillna(y.mean()) # 或使用插值法
权重矩阵不匹配：
# 确保数据顺序一致
w = lps.weights.util.reorder(w, merged.index.tolist())
投影系统警告：
# 统一坐标参考系统
data = data.to_crs('EPSG:3395') # 世界墨卡托投影
from scipy.sparse import csr_matrix
w_sparse = csr_matrix(w.full()[0])
敏感性分析：
# 测试不同权重矩阵的影响
weights_list = [w_queen, w_knn, w_dist]
for w in weights_list:
moran = Moran(y, w)
print(f"{w.__class__.__name__}: I={moran.I:.3f}")
空间滞后可视化：
from splot.esda import lisa_cluster
lisa_cluster(moran, merged)
plt.show()
from esda.moran import Moran_Local
# 计算局部莫兰指数
lisa = Moran_Local(y, w_queen, permutations=9999)
# 识别显著热点
merged['hotspot'] = lisa.q # 聚类类型
merged['p_value'] = lisa.p_sim
# 绘制LISA聚类图
fig, ax = plt.subplots(figsize=(12,10))
merged.plot(column='hotspot', categorical=True,
legend=True, ax=ax)
plt.title('LISA Clusters Map')
# 使用splm库进行面板数据分析
from splm import spreg
# 构建时空权重
w_time = lps.weights.block_weights(time_periods)
w_spacetime = lps.weights.w_spatial_temporal(w_queen, w_time)
# 时空回归模型
model = spreg.GM_Error_Het(panel_data, w=w_spacetime)
# 使用numba加速计算
from numba import jit
@jit(nopython=True)
def moran_matrix_calc(x, w, n):
# 高性能矩阵运算实现
pass
# 并行计算多个变量
from joblib import Parallel, delayed
results = Parallel(n_jobs=4)(
delayed(Moran)(merged[col], w_queen)
for col in ['gdp', 'pop', 'income'])
本文系统介绍了利用Python计算全局莫兰指数的完整流程，从理论基础到实战应用，涵盖了数据准备、权重矩阵构建、指数计算、结果可视化等关键环节。掌握空间自相关分析方法，能够帮助研究者发现地理数据中隐藏的空间模式，为城市规划、区域经济、环境科学等领域研究提供有力工具。
随着空间大数据时代的到来，莫兰指数等空间统计方法将与机器学习、时空建模等技术深度融合。建议读者后续可深入学习： 1. 空间计量经济学模型 2. 点模式分析方法 3. 空间交互网络建模 4. 时空数据挖掘技术
附录：推荐学习资源 - 《Spatial Statistical Data Analysis》- Luc Anselin - PySAL官方文档（https://pysal.org） - GeoDa空间分析教程 - ESRI空间统计白皮书 “`
免责声明：本站发布的内容（图片、视频和文字）以原创、转载和分享为主，文章观点不代表本网站立场，如果涉及侵权请联系站长邮箱：is@yisu.com进行举报，并提供相关证据，一经查实，将立刻删除涉嫌侵权内容。