网格矢量如何计算莫兰指数
引言
遇到一个问题，计算矢量网格的莫兰指数。
概念解释
莫兰指数
莫兰指数（Moran’s Index）是一种空间自相关指标，用于衡量空间数据的相似性和聚集程度。它可以用来描述一个区域与其邻近区域之间的属性值的相关性。莫兰指数的取值范围通常在-1到1之间。
- 当莫兰指数接近1时，表示空间数据呈现出正相关，即相似的值倾向于聚集在一起。
- 当莫兰指数接近-1时，表示空间数据呈现出负相关，即不同的值倾向于聚集在一起。
- 当莫兰指数接近0时，表示空间数据呈现出随机分布，没有明显的空间自相关性。
knearst=4?
knearst=4矩阵是一种空间权重矩阵，用于定义空间数据中每个观测点的邻域。在这种矩阵中，每个观测点的邻域由其最近的4个点组成。
解决思路
计算矢量数据中每个要素（网格）的局部莫兰指数，并将计算结果添加到矢量数据的属性表中。我做了一个示意矢量，如图所示：
因为需要涉及到矢量数据的操作，这里我们使用gdal。
还涉及到莫兰指数，我们使用pysal，这个包用于空间权重矩阵的构建、空间自相关指标的计算、空间回归模型的估计等。
初始化和读取矢量数据
import numpy as np
import pysal
from osgeo import ogr
driver = ogr.GetDriverByName('ESRI Shapefile')
SHP_PATH = r"矢量数据.shp"
dataset = driver.Open(SHP_PATH, 1)
layer = dataset.GetLayer()
- 使用
ogr库打开矢量数据文件（ESRI Shapefile），以读写模式打开。 - 获取矢量数据的图层。
提取属性值和坐标
values = []
coords = []
for feature in layer:
geom = feature.GetGeometryRef()
centroid = geom.Centroid()
coords.append([centroid.GetX(), centroid.GetY()])
values.append(feature.GetField('singlearea'))
values = np.array(values)
coords = np.array(coords)
- 遍历图层中的每个要素（feature）。
- 获取要素的几何体（geometry），并计算其质心坐标。
- 将质心坐标添加到
coords列表中。 - 将指定字段（‘singlearea’）的属性值添加到
values列表中。 - 将属性值和坐标转换为 NumPy 数组。
创建权重矩阵
knn = pysal.lib.weights.KNN(coords, k=4)
knn.transform = 'r'
- 使用
pysal库的KNN函数创建 k 最近邻权重矩阵，设置k=4。 - 对权重矩阵进行行标准化。
计算局部莫兰指数
local_moran = pysal.explore.esda.Moran_Local(values, knn)
print("局部莫兰指数：", local_moran.Is)
# 标准化局部莫兰指数
min_value = np.min(local_moran.Is)
max_value = np.max(local_moran.Is)
normalized_local_moran = (local_moran.Is - min_value) / (max_value - min_value) * 2 - 1
print("标准化后的局部莫兰指数：", normalized_local_moran)
- 使用
pysal库的Moran_Local函数计算每个网格的局部莫兰指数。 - 打印计算得到的局部莫兰指数。
将局部莫兰指数添加到矢量数据属性表
lisa_field = ogr.FieldDefn('LISA_I', ogr.OFTReal)
layer.CreateField(lisa_field)
dataset = None
dataset = driver.Open(SHP_PATH, 1)
layer = dataset.GetLayer()
for i in range(layer.GetFeatureCount()):
feature = layer.GetFeature(i)
feature.SetField('LISA_I', float(local_moran.Is[i]))
layer.SetFeature(feature)
- 创建一个新的字段（‘LISA_I’）来存储局部莫兰指数。
- 重新打开矢量数据集并获取图层。
- 遍历图层中的每个要素。
- 使用
layer.GetFeature(i)获取要素，并将对应的局部莫兰指数赋值给新字段。 - 更新要素的属性表。
关闭数据集并销毁数据源
dataset.Destroy()
dataset = None
print("局部莫兰指数已成功添加到矢量数据属性表中。")
- 关闭矢量数据集。
- 销毁数据源以释放资源。
- 打印提示信息，表示局部莫兰指数已成功添加到矢量数据的属性表中。
完整代码
import numpy as np
import pysal
from osgeo import ogr
# 打开矢量数据文件（以读写模式打开）
driver = ogr.GetDriverByName('ESRI Shapefile')
SHP_PATH = r"矢量数据 - 副本.shp"
dataset = driver.Open(SHP_PATH, 1)
layer = dataset.GetLayer()
# 提取属性值和坐标
values = []
coords = []
for feature in layer:
geom = feature.GetGeometryRef()
centroid = geom.Centroid()
coords.append([centroid.GetX(), centroid.GetY()])
values.append(feature.GetField('cenlan'))
# 将属性值和坐标转换为NumPy数组
values = np.array(values)
coords = np.array(coords)
# 创建k最近邻权重矩阵（knearst=4）
knn = pysal.lib.weights.KNN(coords, k=4)
# 行标准化权重矩阵
knn.transform = 'r'
# 计算每个网格的局部莫兰指数
local_moran = pysal.explore.esda.Moran_Local(values, knn)
print("局部莫兰指数：", local_moran.Is)
# 标准化局部莫兰指数
min_value = np.min(local_moran.Is)
max_value = np.max(local_moran.Is)
normalized_local_moran = (local_moran.Is - min_value) / (max_value - min_value) * 2 - 1
print("标准化后的局部莫兰指数：", normalized_local_moran)
# 将标准化后的局部莫兰指数添加到矢量数据属性表，使用有效的字段名称
lisa_field = ogr.FieldDefn('LISA_I', ogr.OFTReal)
layer.CreateField(lisa_field)
# 重新打开数据集并获取图层
dataset = None
dataset = driver.Open(SHP_PATH, 1)
layer = dataset.GetLayer()
# 使用 layer.GetFeature(i) 获取要素并更新，使用更新后的字段名称
for i in range(layer.GetFeatureCount()):
feature = layer.GetFeature(i)
feature.SetField('LISA_I', float(normalized_local_moran[i]))
layer.SetFeature(feature)
# 关闭数据集并销毁数据源
dataset.Destroy()
dataset = None
print("标准化后的局部莫兰指数已成功添加到矢量数据属性表中。")
效果展示
运行完代码，效果为：
总结
使用gdal负责空间数据处理，使用pysal完成莫兰指数的计算，然后把计算结果写入到属性表里，
