import esda
import pandas as pd
import geopandas as gpd
from geopandas import GeoDataFrame
import libpysal as lps
import numpy as np
import matplotlib.pyplot as plt
from shapely.geometry import Point
import contextily as ctx
from pylab import figure, scatter, show
import warnings
warnings.filterwarnings('ignore')
%matplotlib inline
关于geopandas和contextily库的安装，参考我的另一篇博客https://www.jianshu.com/p/a1496e316bb4
读取数据
一个是地图信息，一个是需要研究的指标。
# 除中华人民共和国.shp外其余三个文件必须也在，只有一个则读取出错
# 中华人民共和国.dbf
# 中华人民共和国.prj
# 中华人民共和国.shx
gdf_data = gpd.read_file('data/中华人民共和国.shp',encoding='utf_8_sig')
feature_df = pd.read_excel('data/清洁低碳.xlsx')
构建邻接矩阵
wq = lps.weights.Queen.from_dataframe(gdf_data) #使用Quuen式邻接矩阵
wq.transform = 'r' #标准化矩阵
# 两个省没有直接相邻的省（34行数据不清楚是什么，下载地图shp文件中自带的）
print(gdf_data.loc[[20,31],'name'])
'''
('WARNING: ', 20, ' is an island (no neighbors)')
('WARNING: ', 31, ' is an island (no neighbors)')
('WARNING: ', 34, ' is an island (no neighbors)')
20 海南省
31 台湾省
Name: name, dtype: object
'''
centroids = gdf_data.geometry.centroid # 计算多边形几何中心
# 本图将邻接矩阵覆盖在整个地图上，如果只需要邻接矩阵，删除ax = gdf_data.plot(figsize=(10, 10),cmap="Blues")
ax = gdf_data.plot(figsize=(10, 10),cmap="Blues")
plt.plot(centroids.x, centroids.y, '.')
for k, neighs in wq.neighbors.items():
origin = centroids[k]
for neigh in neighs:
segment = centroids[[k,neigh]]
plt.plot(segment.x, segment.y, '-')
plt.title('Queen Neighbor Graph')
plt.axis('off')
plt.show()
plt.savefig("save/{}.jpg".format('wq'),dpi=300) # 图片模糊问题，需要指定图片的dpi大小
局部中国地图.jpg
批量计算2010-2020年的全局莫兰指数
合并数据，自动去掉不研究的数据（不带指标的省份）
# 实际数据只研究30个省份，不包含香港澳门，台湾和西藏
gdf_drop = gdf_data.merge(feature_df,left_on='name',right_on='地区')
wq = lps.weights.Queen.from_dataframe(gdf_drop)
wq.transform = 'r' # 标准化矩阵
centroids = gdf_drop.geometry.centroid # 计算多边形几何中心
wr = lps.weights.Rook.from_dataframe(gdf_drop)# 使用Rook式邻接矩阵
Moran_I_list = []
random_Z_list = []
random_Zp_list = []
norm_Z_list = []
norm_Zp_list = []
for i in range(12, 23): # 按照列名位置
y_name = gdf_drop.columns[i]
y = gdf_drop[y_name]
mi = esda.moran.Moran(y, wq)
Moran_I_list.append(mi.I)
random_Z_list.append(mi.z_rand)
random_Zp_list.append(mi.p_rand)
norm_Z_list.append(mi.z_norm)
norm_Zp_list.append(mi.p_norm)
Moran_I_df = pd.DataFrame()
Moran_I_df['year'], Moran_I_df['Moran_I'], Moran_I_df['random_Z'], Moran_I_df['random_p'], Moran_I_df['norm_Z'], Moran_I_df['norm_p'] = \
(range(2010, 2021), Moran_I_list, random_Z_list, random_Zp_list, norm_Z_list, norm_Zp_list)
print(Moran_I_df)
'''
year Moran_I random_Z random_p norm_Z norm_p
0 2010 0.4xxxxx 4.xxxxx 0.0xxxxx 4.1xxxxx 0.0xxxxx
1 2011 0.4xxxxx 3.xxxxx 0.0xxxxx 3.7xxxxx 0.0xxxxx
2 2012 0.4xxxxx 3.xxxxx 0.0xxxxx 3.8xxxxx 0.0xxxxx
3 2013 0.4xxxxx 3.xxxxx 0.0xxxxx 3.9xxxxx 0.0xxxxx
4 2014 0.4xxxxx 4.xxxxx 0.0xxxxx 4.0xxxxx 0.0xxxxx
（省略）
'''
以2020年为例展示计算的莫兰指数和检验结果
y = gdf_drop["2020y"]
mi = esda.moran.Moran(y, wq)
print("Moran's I 值为：",mi.I)
print("随机分布假设下Z检验值为：",mi.z_rand)
print("随机分布假设下Z检验的P值为：",mi.p_rand)
print("正态分布假设下Z检验值为：",mi.z_norm)
print("正态分布假设下Z检验的P值为：",mi.p_norm)
'''
Moran's I 值为： 0.4542425813757728
随机分布假设下Z检验值为： 4.xxxxxx
随机分布假设下Z检验的P值为： 6.xxxxxxe-05
正态分布假设下Z检验值为： 4.xxxxxx
正态分布假设下Z检验的P值为： 5.xxxxxx-05
'''
绘制2020年的莫兰散点图
from splot.esda import plot_moran
fig, ax = plot_moran(mi, zstandard=True, figsize=(11,5),aspect_equal=False)
fig.show()
fig.savefig("save/{}.jpg".format('清洁低碳全局莫兰散点图'),dpi=300) # 图片模糊问题，需要指定图片的dpi大小
清洁低碳全局莫兰散点图.jpg
2020年全局莫兰指数检验通过，继续做局部莫兰指数(全局不通过，不用做局部)
from splot.esda import moran_scatterplot
from esda.moran import Moran_Local
from splot.esda import lisa_cluster
from splot.esda import plot_local_autocorrelation
y = gdf_drop['2020y'].values
w = lps.weights.distance.Kernel.from_dataframe(gdf_drop, fixed=False, k=15)
w.transform = 'r'
moran_loc = Moran_Local(y, w)
# 局部莫兰指数
loc_moran_df = pd.DataFrame({'name':gdf_drop.name, 'loc_moran':moran_loc.Is}) # moran_loc.Is是各省的局部莫兰指数
print(loc_moran_df)
'''
name loc_moran
0 北京市 -0.033960
1 天津市 0.319125
2 河北省 -0.086364
3 山西省 0.000686
4 内蒙古自治区 0.068562
5 辽宁省 -0.063257
6 吉林省 0.002306
7 黑龙江省 -0.012012
8 上海市 0.001538
9 江苏省 0.025877
10 浙江省 0.307162
(省略)
'''
局部莫兰散点图
fig, ax = moran_scatterplot(moran_loc, p=0.1,aspect_equal=True)
#ax.set_xlabel('FSTGDPRATE')
#ax.set_ylabel('Spatial Lag of FSTGDPRATE')
#plt.show()
fig.savefig("save/{}.jpg".format('清洁低碳局部莫兰散点图'),dpi=300) # 图片模糊问题，需要指定图片的dpi大小
清洁低碳局部莫兰散点图.jpg
空间分布图(不展示组合图)
强调：没有计算西藏地区，所以中国大陆地图中缺失西藏地区
fig, ax = lisa_cluster(moran_loc, gdf_drop, p=0.1, figsize = (9,9))
fig.savefig("save/{}.jpg".format('清洁低碳聚集区的空间分布图'),dpi=300) # 图片模糊问题，需要指定图片的dpi大小
#绘制组合图(不展示组合图)
#fig, ax = plot_local_autocorrelation(moran_loc, gdf_drop, '2020y',p=0.1,figsize = (27,6))
#fig.savefig("save/{}.jpg".format('清洁低碳组合图'),dpi=300) # 图片模糊问题，需要指定图片的dpi大小
能源高质量发展聚集区的空间分布图.jpg
参考：
https://www.yisu.com/zixun/503529.html
https://www.yisu.com/zixun/504927.html
实操问题
做完全局莫兰指数后，计算局部莫兰指数的函数 Moran_Local 一直出错，其中提到 numba 库的问题。我查询报错，发现有个提升 numba 版本的说法，然后我升级到最新版本。随后直接导致 esda 库和 libpysal 库导入失败，导致全局莫兰指数部分无法运行。随后我卸载了numba，esda 和 libpysal，再次安装 numba，esda 和 libpysal，依然无法导入 esda 和 libpysal。尝试多种方法后，我再一次卸载了 numba，esda 和 libpysal，然后只安装 esda 库，神奇的是我可以导入esda，甚至可以直接导入 libpysal，连之前出错的 Moran_Local 函数也可以运行了。我推测是Moran_Local 依赖的 numba 库版本和我的 numba 库版本不一致，在我将 numba 库版本升级成最新版时，又与 esda 和 libpysal 不匹配。最后我删除 numba 后，只安装最新的 esda，会自己下载合适的 numba，然后就可以运行 Moran_Local 函数（esda.moran.Moran_Local）。
此处展示我的环境：
- Python 3.6.5
- numba 0.38.0
- esda 2.4.1
- libpysal 4.5.0
