matlab求莫兰指数程序,python计算莫兰指数(Moran's I)并绘制地区热力图——以中国各省pm2.5为例...
[TOC]
程序简介
计算省的pm2.5平均值作为观测矩阵,省会的距离的倒数作为空间权重矩阵,计算全局莫兰指数为0.49,显著性检验p值为3.75>1.96,得出中国地区的pm2.5存在空间正相关
输入:中国各地区pm2.5值 输出:莫兰散点图、莫兰指数、莫兰检验数、地区热力图
莫兰指数(Moran’s I)是空间自相关系数的一种,其值分布在[-1,1],用于判别空间是否存在自相关。比如:一个小区内有钱人会集聚住在豪宅别墅区,平民们(比如我- -)则会集聚住在普通的楼房区,这里的个人财富就是一种观测值,呈现高高聚集,低低聚集的现象。
程序/数据集下载
点击进入下载地址
代码分析
导入模块,路径
# -*- coding: utf-8 -*-
from Module.CalPosition import calDistance
from Module.MoranI import moranI
import pandas as pd
import numpy as np
import json
import os
from pyecharts import Map
import re
#路径目录
baseDir = ''#当前目录
staticDir = os.path.join(baseDir,'Static')#静态文件目录
resultDir = os.path.join(baseDir,'Result')#结果文件目录
整理省-pm2.5表格作为观测值矩阵,将省内所有市的pm2.5的算术平均值作为该省的pm2.5,查看前5行
#读取省市字典
with open(staticDir+'/中国省市字典.json','r') as f:
china = json.loads(f.read())
#省-pm2.5字典
provincePm = {'province':[],'pm2.5':[]}
#读取各市pm2.5,求平均值作为整个省的pm2.5值
for province in china:
pm = []#该省所有市的pm2.5列表
for city in china[province]:
#文件名没有市这个字
city = city.replace('市','')
path = staticDir+'/中国各市2个月的pm2.5/%s.csv'%city
if not os.path.exists(path):
continue
cityPm = pd.read_csv(path,engine="python")['pm25']
pm.extend(cityPm.values.tolist())
#列表不为空才添加
if pm:
provincePm['province'].append(province)
provincePm['pm2.5'].append(np.mean(pm))
provincePm = pd.DataFrame(provincePm).set_index(['province'])
#观测值矩阵
provincePm = provincePm.drop(['西藏藏族自治区','海南省'],axis=0)
provincePm.head()
pm2.5 | |
---|---|
province | |
北京市 | 37.807018 |
天津市 | 45.245614 |
上海市 | 24.228070 |
重庆市 | 26.368421 |
河南省 | 37.522478 |
构建空间权重矩阵,该矩阵反应了地区间空间联系,一般距离越近,权重越大,但是对角线为0 这里要调用calDistance函数计算两地距离,该函数的详情可以查看这篇博客
程序在得到距离后,取其负3次幂作为两地的权重,查看矩阵前5行5列
#构建空间权重矩阵
if not os.path.exists(resultDir+'/距离矩阵.xlsx'):
spaceMatrix = pd.DataFrame({},index=provincePm.index,columns=provincePm.index)
for province1 in spaceMatrix.index:
for province2 in spaceMatrix.columns:
if not np.isnan(spaceMatrix.loc[province1,province2]):
continue
#两地距离
try:
distance = calDistance(province1,province2)
spaceMatrix.loc[province1,province2] = distance
spaceMatrix.loc[province2,province1] = distance
except:
spaceMatrix.loc[province1,province2] = np.nan
spaceMatrix.loc[province2,province1] = np.nan
continue
#地点相同,距离取无穷大,不然后面的倒数会报错
if not distance:
spaceMatrix.loc[province1,province2] = 1e+10
#删除空行和列
spaceMatrix = spaceMatrix.drop(['西藏藏族自治区','海南省'],axis=0)
spaceMatrix = spaceMatrix.drop(['西藏藏族自治区','海南省'],axis=1)
spaceMatrix.to_excel(resultDir+'/距离矩阵.xlsx')
else:
spaceMatrix = pd.read_excel(resultDir+'/距离矩阵.xlsx',index_col=0)
#取反距离,次幂可去1~3
spaceMatrix = spaceMatrix**-3
spaceMatrix.iloc[0:5,0:5]
北京市 | 天津市 | 上海市 | 重庆市 | 河南省 | |
---|---|---|---|---|---|
province | |||||
北京市 | 1.000000e-30 | 6.807960e-16 | 8.255561e-19 | 3.228904e-19 | 4.260133e-18 |
天津市 | 6.807960e-16 | 1.000000e-30 | 1.150491e-18 | 3.355669e-19 | 5.412329e-18 |
上海市 | 8.255561e-19 | 1.150491e-18 | 1.000000e-30 | 3.316709e-19 | 1.808123e-18 |
重庆市 | 3.228904e-19 | 3.355669e-19 | 3.316709e-19 | 1.000000e-30 | 1.415375e-18 |
河南省 | 4.260133e-18 | 5.412329e-18 | 1.808123e-18 | 1.415375e-18 | 1.000000e-30 |
现在有了省-PM2.5表格作为观测值矩阵,空间权重矩阵,输入到moranI函数中可以得到输出结论并自动绘制莫兰散点图,其中moranI函数位置和代码在下文给出。
莫兰散点图,横纵坐标一般是进行分析的数据的离差值,莫兰散点图经常用来研究局部空间不稳定性,其四个象限分别对应于区域单元与其邻居之间四种类型的局部空间联系形式。
#计算莫兰指数
result = moranI(spaceMatrix,provincePm)
print('全局莫兰指数为',np.round(result['I']['value'],2),'说明地区之间的pm2.5正相关')
print('全局莫兰检验数为',np.round(result['ZI_N']['value'],2),'>1.96,在0.05的水平上显著,认为pm2.5存在空间相关性')
全局莫兰指数为 0.49 说明地区之间的pm2.5正相关
全局莫兰检验数为 3.75 >1.96,在0.05的水平上显著,认为pm2.5存在空间相关性
上文中用到的moranI函数位置:Module/MoranI.py,该代码可以直接运行进行测试
# -*- coding: utf-8 -*-
import matplotlib.pyplot as plt
import numpy as np
import os
def moranI(W,X):
'''
W:空间权重矩阵
X:观测值矩阵
归一化空间权重矩阵后进行moran检验,实例https://bbs.pinggu.org/thread-3568074-1-1.html
'''
W = np.array(W)
X = np.array(X)
X = X.reshape(1,-1)
W = W/W.sum(axis=1)#归一化
n = W.shape[0]#空间单元数
Z = X - X.mean()#离差阵
S0 = W.sum()
S1 = 0
for i in range(n):
for j in range(n):
S1 += 0.5*(W[i,j]+W[j,i])**2
S2 = 0
for i in range(n):
S2 += (W[i,:].sum()+W[:,i].sum())**2
#全局moran指数
I = np.dot(Z,W)
I = np.dot(I,Z.T)
I = n/S0*I/np.dot(Z,Z.T)
#在正太分布假设下的检验数
EI_N = -1/(n-1)
VARI_N = (n**2*S1-n*S2+3*S0**2)/(S0**2*(n**2-1))-EI_N**2
ZI_N = (I-EI_N)/(VARI_N**0.5)
#在随机分布假设下检验数
EI_R = -1/(n-1)
b2 = 0
for i in range(n):
b2 += n*Z[0,i]**4
b2 = b2/((Z*Z).sum()**2)
VARI_R = n*((n**2-3*n+3)*S1-n*S2+3*S0**2)-b2*((n**2-n)*S1-2*n*S2+6*S0**2)
VARI_R = VARI_R/(S0**2*(n-1)*(n-2)*(n-3))-EI_R**2
ZI_R = (I-EI_R)/(VARI_R**0.5)
#计算局部moran指数
Ii = list()
for i in range(n):
Ii_ = n*Z[0,i]
Ii__ = 0
for j in range(n):
Ii__ += W[i,j]*Z[0,j]
Ii_ = Ii_*Ii__/((Z*Z).sum())
Ii.append(Ii_)
Ii = np.array(Ii)
#局部检验数
ZIi = list()
EIi = Ii.mean()
VARIi = Ii.var()
for i in range(n):
ZIi_ = (Ii[i]-EIi)/(VARIi**0.5)
ZIi.append(ZIi_)
ZIi = np.array(ZIi)
#moran散点图
#用来正常显示中文标签
plt.rcParams['font.sans-serif']=['SimHei']
#用来正常显示负号
plt.rcParams['axes.unicode_minus']=False
fig = plt.figure()
ax = fig.add_subplot(1,1,1)
ax.spines['top'].set_color('none')
ax.spines['right'].set_color('none')
ax.xaxis.set_ticks_position('bottom')
ax.spines['bottom'].set_position(('data', 0))
ax.yaxis.set_ticks_position('left')
ax.spines['left'].set_position(('data', 0))
WZ = np.dot(Z,W)
ax.scatter(Z,WZ,c='k')
x1 = range(int(Z.min()),int(Z.max()+1))
y1 = range(int(Z.min()),int(Z.max()+1))
ax.plot(x1,y1,'k--',label='x=y')
x2 = list(range(int(Z.min()),int(Z.max()+1)))
y2 = np.array(x2)*I[0][0]
ax.plot(x2,y2,'k-',label='I*x=y')
ax.legend(loc='upper right')
imgPath = os.path.join(os.getcwd(),'莫兰散点图.png')
ax.set_title('莫兰散点图')
fig.savefig(imgPath)
return {
'I':{'value':I[0,0],'desc':'全局moran指数'},
'ZI_N':{'value':ZI_N[0,0],'desc':'正太分布假设下检验数'},
'ZI_R':{'value':ZI_R[0,0],'desc':'随机分布假设下检验数'},
'Ii':{'value':Ii,'desc':'局部moran指数'},
'ZIi':{'value':ZIi,'desc':'局部检验数'},
'img':{'path':imgPath,'desc':'莫兰散点图路径'}
}
if __name__ == "__main__":
w = [
[0,1,1,0,0],
[1,0,1,1,0],
[1,1,0,1,0],
[0,1,1,0,1],
[0,0,0,1,0]
]
w = np.array(w)
x = [
[8,6,6,3,2]
]
x = np.array(x)
print(moranI(w,x))
最后根据上文得到的局部莫兰指数,绘制中国各省的热力图
#热力图
provice = list(re.sub(r'市|省|.族自治区|维吾尔自治区|自治区','',index) for index in spaceMatrix.index)
values = list(result['Ii']['value'])
map = Map("pm2.5莫兰指数热力图", '', width=1200, height=600)
map.add("", provice, values, visual_range=[-1, 1], maptype='china', is_visualmap=True,
visual_text_color='#000')
map.render(path=resultDir+"/pm2.5莫兰指数热力图.html")
matlab求莫兰指数程序,python计算莫兰指数(Moran's I)并绘制地区热力图——以中国各省pm2.5为例...相关推荐
- 全局莫兰指数_R:计算莫兰系数 R code for calculating Moran#x27;s I
封图来自网络,侵删 莫兰指数(Moran's I)可能是测量空间上集聚的一个重要的指数,当方差归一化之后,它的值会落到[-1, 1]的区间上. 当 的时候,我们可以认为某个指标在空间上 ,也就是说空间 ...
- matlab求点介数程序,matlab_bgl 一个很有用的计算网络中每个节点介数的程序,对 分析 Cloud Computing 云 266万源代码下载- www.pudn.com...
文件名称: matlab_bgl下载 收藏√ [ 5 4 3 2 1 ] 开发工具: Others 文件大小: 2098 KB 上传时间: 2016-10-26 下载次数: 0 提 供 者 ...
- python利用以下公式求π的值_使用Python计算 π 值
π是一个无数人追随的真正的神奇数字.我不是很清楚一个永远重复的无理数的迷人之处.在我看来,我乐于计算π,也就是计算π的值.因为π是一个无理数,它 是无限的.这就意味着任何对π的计算都仅仅是个近似值.如 ...
- matlab求根的原程序,MATLAB求根程序求帮助
我有一个函数g=m*(m1/m)*(m2/m)*(r^2)*rm*Wc1*cos(2*u),我需要一个程序求这个函数在在(-2*pi,2*pi)上所有的根. 其中u是自变量,其他所有参数都是已知.已知 ...
- 求瑞年的java程序,java 计算瑞年的方法
任何语言都有可能计算某一年是否为瑞年的方法,也就是说一年有 366 天,每隔4 年就出现一次.最基本的算法如下:if year is divisible by 400 then is_leap_yea ...
- c语言求平均数double,编写程序以计算浮点值的平均值
定义一个函数,用于计算任意数量的浮点值的平均值.double类型值的数组在数组参数中传递给函数.读取从键盘输入的任意数量的值并输出平均值. 实现代码 #define __STDC_WANT_LIB_E ...
- matlab求矩阵行最简形,计算矩阵行最简行的命令
初等行变换之互换两行 Public Sub Matrix_Specify_Tow_Row_Exchange(Row_A_Index As Integer, Row_B_Index As Integer ...
- matlab求cos角,科学网—MATLAB求太阳高度角的小程序 - 张乐乐的博文
参考链接:http://bbs.06climate.com/forum.php?mod=viewthread&tid=36366 代码部分: function HSI=calHSI(year, ...
- 用python计算身体BMI指数
tall=float(input("请输入身高/m:")) kg=float(input("请输入体重/kg:")) BMI=kg/tall/tall prin ...
最新文章
- Java入门教程五(数字和日期处理)
- php中插入表格 标签,PHP_HTML中的表格元素,一,table标签。tablegt - phpStudy
- Activity的四种加载模式(转载)
- UVA 818 Cutting Chains 切断圆环链 (暴力dfs)
- Google的特殊功能
- 一个牛人给的java九点建议
- AlphaBlend失败,错误码87
- 王之泰201771010131《面向对象程序设计(java)》第十三周学习总结
- 红蓝对抗——蓝军(CheckList)总结
- [技术分享 – FCS 篇] 驭龙五式3之飞龙在天:安装 FCS 服务器
- 求最大公约数,最小公倍数
- 用Python画樱花树的代码
- Akita与脉冲云的关系
- 自定义Firefox、IE收藏夹
- 自然语言处理——词性标注、词干提取、词形还原
- 计算身体质量指数BMI
- 生物信息学软件_生物信息学视频教程大赛
- fiddler模拟进行接口测试
- JAVA毕业设计计算机实验中心网站计算机源码+lw文档+系统+调试部署+数据库
- ie退出全屏快捷键_视频剪辑快捷键大全
热门文章
- [应用相关]NUCLEO-H7A3ZI-Q板子学习
- HTML5-iframe-frameset
- ProjectDay04
- tankbot 机器人_优必选首款履带式Jimu机器人 TankBot 登陆Apple Store零售店
- int3断点指令的原理和示例
- 腾讯云:聚焦“双十一”背后 不容忽视的电商风控与安全
- 利用端口,进程,文件,服务和日志信息来排查系统安全
- Redis实战:如何构建类微博的亿级社交平台
- python 爬虫小试牛刀(request,BeautifulSoup库的实战)
- 解决Linux服务器中TCP的FIN_WAIT2,CLOSE_WAIT状态连接过多的问题