锐多宝
锐多宝的朋友圈
锐多宝

网格矢量如何计算莫兰指数

引言

遇到一个问题，计算矢量网格的莫兰指数。

概念解释

莫兰指数

莫兰指数（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()
  1. 使用 ogr 库打开矢量数据文件（ESRI Shapefile），以读写模式打开。
  2. 获取矢量数据的图层。

提取属性值和坐标

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)
  1. 遍历图层中的每个要素（feature）。
  2. 获取要素的几何体（geometry），并计算其质心坐标。
  3. 将质心坐标添加到 coords 列表中。
  4. 将指定字段（'singlearea'）的属性值添加到 values 列表中。
  5. 将属性值和坐标转换为 NumPy 数组。

创建权重矩阵

knn = pysal.lib.weights.KNN(coords, k=4)
knn.transform = 'r'
  1. 使用 pysal 库的 KNN 函数创建 k 最近邻权重矩阵，设置 k=4。
  2. 对权重矩阵进行行标准化。

计算局部莫兰指数

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)
  1. 使用 pysal 库的 Moran_Local 函数计算每个网格的局部莫兰指数。
  2. 打印计算得到的局部莫兰指数。

将局部莫兰指数添加到矢量数据属性表

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)</code></pre><ol><li>创建一个新的字段（'LISA_I'）来存储局部莫兰指数。</li><li>重新打开矢量数据集并获取图层。</li><li>遍历图层中的每个要素。</li><li>使用 <code>layer.GetFeature(i)</code> 获取要素，并将对应的局部莫兰指数赋值给新字段。</li><li>更新要素的属性表。</li></ol><h3>关闭数据集并销毁数据源</h3><pre><code>dataset.Destroy()
dataset = None
print("局部莫兰指数已成功添加到矢量数据属性表中。")
  1. 关闭矢量数据集。
  2. 销毁数据源以释放资源。
  3. 打印提示信息，表示局部莫兰指数已成功添加到矢量数据的属性表中。

完整代码

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完成莫兰指数的计算，然后把计算结果写入到属性表里，

2024年04月09日
锐多宝 .