python计算局部空间自相关——问题记录

 问题描述

        在用esda包计算空间自相关时，按照官方示例代码：

    gdf = gpd.read_file(shp)
    # 使用Quuen式邻接矩阵
    w = lps.weights.Queen.from_dataframe(gdf)  
    # 标准化矩阵
    w.transform = 'r'  
    y = gdf['col']
    # 全局空间自相关
    m = esda.moran.Moran(y, w)
    # 局部空间自相关
    ml = esda.moran.Moran_Local(y,w)

        全局空间自相关能获得结果，但是计算局部空间自相关时会报如下错误：

ValueError: cannot assign slice from input of different size

        大概知道造成问题的直接原因应该是大小不匹配。但是，始终不知道是哪两个大小不匹配，以及怎么样让它们匹配。

---

解决思路

        后续在用GeoDa软件计算空间自相关时候，注意到了结果上面一行字：

        于是就想到，是不是把gdf里面没有邻居的对象删了就没问题了。

        然后，对代码做了如下更改：

    gdf = gpd.read_file(shp)
    # 使用Quuen式邻接矩阵
    w = lps.weights.Queen.from_dataframe(gdf)  
    # 标准化矩阵
    w.transform = 'r'

    # 在计算完邻接矩阵后，按照无邻接对象的唯一标识删除gdf里的对象
    # 注意！！！shp文件里面需要有一列记录FID（这里是FID_1），根据这个删除
    gdf = gdf[~gdf['FID_1'].isin(w.islands)]
    # 再重新计算邻接矩阵，这样就没有无邻接对象了
    w = lps.weights.Queen.from_dataframe(gdf)
    w.transform = 'r'

    y = gdf['col']
    # 全局空间自相关
    m = esda.moran.Moran(y, w)
    # 局部空间自相关
    ml = esda.moran.Moran_Local(y,w)

        注意！！！shp文件里面需要添加一列记录FID（这里是FID_1），根据这个删除！！！

        结果证明，删去之后就能正常运行了。

---

其他——导出局部空间自相关的聚类

        在解决完不能正常运行的问题后，我又想把局部空间自相关的聚类结果导出来，方便自己制图。于是，研究起了esda库里的代码。esda库能够绘制聚类地图，我就从绘制地图的函数（lisa_cluster）开始溯源，看看能不能找到储存聚类结果的变量。找到的最关键代码如下：

gdf.assign(cl=labels).plot(column='cl', categorical=True,
                           k=2, cmap=hmap, linewidth=0.1, ax=ax,
                           edgecolor='white', legend=legend,
                           legend_kwds=legend_kwds, **kwargs)

        原函数是把'labels'添加给了gdf，然后基于'labels'进行制图。可以推得，'labels'储存的就是聚类结果，且大小与gdf相同（事实也确实如此，详情请看_viz_utils.py文件中的moran_hot_cold_spots函数）。所以，只要获得'labels'变量就能得到聚类结果。而'labels'变量是从mask_local_auto函数得到的，但是函数返回的'labels'是文本列表（'ns','HH','HL','LL','LH'），我后续想要转栅格，所以就需要再做一个反变换，把文本转为数值。综上，可以整理出如何获得聚类结果的方法：

from splot._viz_utils import mask_local_auto

csv_path = 'your_csv_path'
_,_,_,labels = mask_local_auto(ml,p=0.05)
cluster_labels_dict = {'ns':0, 'HH':1, 'LH':2, 'LL':3, 'HL':4}
label_cl = [cluster_labels_dict[cl] for cl in labels]
# 这里按需求修改，我想根据Id连接表，就保留了'Join_Id'和'cl'列
gdf_with_cl = gdf.assign(cl = label_cl)[['Join_Id','cl']]
gdf_with_cl.to_csv(csv_path)

注意： 这里的'Join_Id'需要与上文的'FID_1'进行区分，两者不是一个。也就是说原来的shp文件里至少有两个字段：

FID_1：与shp文件的FID相同，起到删除对象索引的作用

Join_Id：这个是用来连接表的索引。我的文件两个ID并不匹配，如果匹配可以用一个，我是因为有特殊要求。

---

结果展示：

GeoDa聚类地图：

导出csv，然后连接表，利用Arcgis pro 制图：