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

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

引言

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

概念解释

莫兰指数

莫兰指数(Moran’s Index)是一种空间自相关指标,用于衡量空间数据的相似性和聚集程度。它可以用来描述一个区域与其邻近区域之间的属性值的相关性。莫兰指数的取值范围通常在-1到1之间。

  • 当莫兰指数接近1时,表示空间数据呈现出正相关,即相似的值倾向于聚集在一起。
  • 当莫兰指数接近-1时,表示空间数据呈现出负相关,即不同的值倾向于聚集在一起。
  • 当莫兰指数接近0时,表示空间数据呈现出随机分布,没有明显的空间自相关性。

knearst=4?

knearst=4矩阵是一种空间权重矩阵,用于定义空间数据中每个观测点的邻域。在这种矩阵中,每个观测点的邻域由其最近的4个点组成。

示意图,这个用距离小时

解决思路

计算矢量数据中每个要素(网格)的局部莫兰指数,并将计算结果添加到矢量数据的属性表中。我做了一个示意矢量,如图所示:

因为需要涉及到矢量数据的操作,这里我们使用gdal

还涉及到莫兰指数,我们使用pysal,这个包用于空间权重矩阵的构建、空间自相关指标的计算、空间回归模型的估计等。

初始化和读取矢量数据

import numpy as np
import pysal
from osgeo import ogrdriver = 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)
  1. 创建一个新的字段(‘LISA_I’)来存储局部莫兰指数。
  2. 重新打开矢量数据集并获取图层。
  3. 遍历图层中的每个要素。
  4. 使用 layer.GetFeature(i) 获取要素,并将对应的局部莫兰指数赋值给新字段。
  5. 更新要素的属性表。

关闭数据集并销毁数据源

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 = Noneprint("标准化后的局部莫兰指数已成功添加到矢量数据属性表中。")

效果展示

运行完代码,效果为:

总结

使用gdal负责空间数据处理,使用pysal完成莫兰指数的计算,然后把计算结果写入到属性表里,

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.mzph.cn/news/803418.shtml

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

MySQL之函数:字符串函数、日期函数、数值函数、流程函数

字符串函数:用于对文本数据进行操作和处理 CONCAT:用于将多个字符串连接成一个字符串。 SELECT CONCAT(Hello, , World); -- 输出: Hello WorldSUBSTRING:用于截取字符串的子串。包前不包后,从1开始 SELECT SUBSTRING(MySQL, 1, …

这4大流氓软件,以后看见记得绕着走。

第一个,2345浏览器。时不时弹出广告,还会自动修改你的默认主页,并且很难修改回去。 第二个,搜狗输入法。别以为它打字很方便,实际上弹广告也很方便。 第三个,快压。解压不仅给你弹广告,还能让…

TCP-IP详解卷一:协议——阅读总结

该内容适合程序员查看 第1章 概述 1.1 引言 WAN全称是 Wide Area Network,中文名为广域网。 LAN全称是 Local Area Network,中文名为局域网。 1.2分层 ICP/IP协议族通常被认为是一个四层协议系统 分层协议应用层Telnet、FTP和e-mail运输层TCP和UDP网…

在Terraform中定义模块依赖关系

在Terraform中定义模块依赖关系 是指在Terraform配置文件中明确声明模块之间的依赖关系,以确保正确的部署和管理基础设施。 Terraform是一个开源的基础设施即代码工具,它允许开发人员使用简单的声明性语言来定义和管理基础设施资源。模块是Terraform中…

EVAL-21489-EZLITE原理图下载方法

1,进入官网,搜索“21489”: 2,下拉找到EVAL-21489-EZLITE,点击“文件”: 3,点击“电路板设计数据库”的“查看全部”: 4,点击下载即可: 5,下载完成…

系统设计之订单系统中如何防止商品超卖

一、策略 在订单系统中,防止商品超卖是一个至关重要的问题,涉及到多个方面的策略和技术。以下是一些防止商品超卖的主要措施: 库存实时同步与检查: 确保订单系统和库存系统之间的数据实时同步,每次订单生成或取消时&…

Python 全栈体系【四阶】(二十五)

第五章 深度学习 二、计算机视觉基本理论 11. 图像梯度处理 11.1 什么是图像梯度 图像梯度计算的是图像变化的速度。对于图像的边缘部分,其灰度值变化较大,梯度值也较大;相反,对于图像中比较平滑的部分,其灰度值变化…

【鸿蒙开发】ArkTS和组件

1. 初识ArkTS语言 ArkTS是HarmonyOS优选的主力应用开发语言。ArkTS围绕应用开发在TypeScript生态基础上做了进一步扩展,继承了TS的所有特性。 当前,ArkTS在TS的基础上主要扩展了如下能力: 基本语法:ArkTS定义了声明式UI描述、自…

一招搞定vcruntime140_1.dll无法继续执行代码的解决方法

在我们日常频繁地与计算机互动、依赖其高效处理各类任务的过程中,偶尔会遭遇一些突发的技术问题,导致原本顺畅的操作流程被迫中断。其中一种常见的困扰便是系统弹出一则明确且令人颇感困惑的错误提示:“由于找不到vcruntime140_1.dll文件&…

C语言的显式类型转换和隐式类型转换详细讲解

目录 一、类型转换 1、显式类型转换 2、隐式类型转换 二、算术转换 三、总结 每个编译器都会对表达式做两件事情,一是判断表达式中操作符的优先级和结合性,二是判断表达式中的操作数类型是否一致,如果不一致则需要进行类型转换。第一点在…

机器学习(五) -- 监督学习(2) -- k近邻

系列文章目录及链接 目录 前言 一、K近邻通俗理解及定义 二、原理理解及公式 1、距离度量 四、接口实现 1、鸢尾花数据集介绍 2、API 3、流程 3.1、获取数据 3.2、数据预处理 3.3、特征工程 3.4、knn模型训练 3.5、模型评估 3.6、结果预测 4、超参数搜索-网格搜…

从零开始学RSA:已知n,e,d求p,q和私钥文件修复

(8)已知n,e,d求p,q 一看这个标题你就应该有个觉悟,n一定无法直接分解得到p和q。 题目: 10-存货5 题目给出了两个文件,一个是加密脚本chall.py,一个是加密后输出的内容output.txt。 分析一下加密脚本: from gmpy2 import invertf…

相机模型浅析

相机模型 文章目录 相机模型四个坐标系针孔相机模型世界坐标系到相机坐标系相机坐标系到图像坐标系图像坐标到像素坐标 四个坐标系 ①世界坐标系:是客观三维世界的绝对坐标系,也称客观坐标系。因为数码相机安放在三维空间中,我们需要世界坐标…

Python3 replace()函数使用详解:字符串的艺术转换

博主猫头虎的技术世界 🌟 欢迎来到猫头虎的博客 — 探索技术的无限可能! 专栏链接: 🔗 精选专栏: 《面试题大全》 — 面试准备的宝典!《IDEA开发秘籍》 — 提升你的IDEA技能!《100天精通鸿蒙》 …

JavaScript(1)神秘的编程技巧

大家都感兴趣的箭头函数 箭头函数在许多场景中都可以发挥作用,尤其适用于简化函数声明和提高代码的可读性。以下是箭头函数可以使用的一些常见方面: (1)回调函数: 箭头函数特别适合作为回调函数,例如在事…

RuntimeError: Library cublas64_12.dll is not found or cannot be loaded

运行guillaumekln/faster-whisper-large-v2模型进行语音识别的时候报错了 RuntimeError: Library cublas64_12.dll is not found or cannot be loaded 代码: from faster_whisper import WhisperModelmodel WhisperModel("H:\\model\\guillaumekln\\faster…

Linux系统安装内网穿透实现固定公网地址访问本地MinIO服务

文章目录 前言1. 创建Buckets和Access Keys2. Linux 安装Cpolar3. 创建连接MinIO服务公网地址4. 远程调用MinIO服务小结5. 固定连接TCP公网地址6. 固定地址连接测试 正文开始前给大家推荐个网站,前些天发现了一个巨牛的 人工智能学习网站, 通俗易懂&am…

你不知道的JavaScript---深入理解 JavaScript 作用域

你好,我是小白Coding日志,一个热爱技术的程序员。在这里,我分享自己在编程和技术世界中的学习心得和体会。希望我的文章能够给你带来一些灵感和帮助。欢迎来到我的博客,一起在技术的世界里探索前行吧! 1. 什么是作用域…

MY-Java高级面试题

1. jdk1.7 到 jdk1.8 Map 发生了什么变化 ( 底层 )? 1.8 之后 hashMap 的数据结构发生了变化,从之前的单纯的数组 链表结构变成数组 链 表 红黑树。也就是说在 JVM 存储 hashMap 的 K-V 时仅仅通过 key 来决定每一个 entry 的存 储槽位&…

网络安全:重要性与应对措施

1. 网络安全的重要性 随着互联网的普及和信息技术的快速发展,网络安全问题已经变得日益突出。网络攻击者可以通过各种手段窃取个人信息、破坏系统、传播病毒等,给个人和社会带来巨大的损失。因此,网络安全已经成为信息化时代的重要问题之一。…