程序简介
计算省的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")
