Python LISA 局部空间自相关实战：找出谁跟谁扎堆

from esda.moran import Moran_Local

# 假设 w 是上篇构建的 KNN 权重矩阵，gdf 是 GeoDataFrame
lisa = Moran_Local(gdf['volatility'].values, w)

# 提取结果
gdf['lisa_I'] = lisa.Is          # 每个城市的局部 Moran's I
gdf['lisa_p'] = lisa.p_sim       # p 值（排列检验）
gdf['lisa_q'] = lisa.q           # 象限（1=HH, 2=LH, 3=LL, 4=HL）

# 分类
gdf['cluster'] = '不显著'
gdf.loc[(gdf['lisa_p'] < 0.05) & (gdf['lisa_q'] == 1), 'cluster'] = 'HH'
gdf.loc[(gdf['lisa_p'] < 0.05) & (gdf['lisa_q'] == 2), 'cluster'] = 'LH'
gdf.loc[(gdf['lisa_p'] < 0.05) & (gdf['lisa_q'] == 3), 'cluster'] = 'LL'
gdf.loc[(gdf['lisa_p'] < 0.05) & (gdf['lisa_q'] == 4), 'cluster'] = 'HL'

# 统计
print(gdf['cluster'].value_counts())

不显著    62
HH        4
HL        2
LH        1
LL        1

sig = gdf[gdf['lisa_p'] < 0.05][['CITY', 'volatility', 'cum_return', 'cluster', 'lisa_p']]
print(sig.sort_values('cluster').to_string(index=False))

 CITY  volatility  cum_return cluster  lisa_p
  秦皇岛    0.70%      51.4%      LL   0.027
  南宁      0.79%      61.2%      LH   0.020
  丹东      0.87%      43.7%      HL   0.045
  重庆      0.82%     104.9%      HL   0.045
  三亚      2.04%     186.8%      HH   0.002
  北海      1.30%      84.8%      HH   0.004
  海口      1.82%     184.6%      HH   0.002
  湛江      0.87%      69.1%      HH   0.002

import geopandas as gpd
import matplotlib.pyplot as plt

# 加载行政区划底图（需要一个中国地级市 shapefile）
admin_gdf = gpd.read_file('2023年地级.shp')

fig, ax = plt.subplots(1, 1, figsize=(16, 12))

# 画底图
admin_gdf.plot(ax=ax, color='#fafafa', edgecolor='#dddddd', linewidth=0.3)

# 颜色方案
LISA_COLORS = {'不显著': '#cccccc', 'HH': '#e15759', 'LH': '#76b7b2', 'LL': '#4e79a7', 'HL': '#f28e2b'}

for cluster, color in LISA_COLORS.items():
    mask = gdf['cluster'] == cluster
    if mask.any():
        subset = gdf[mask]
        subset.plot(ax=ax, color=color, markersize=100, alpha=0.85,
                    label=cluster, edgecolor='black' if cluster != '不显著' else None,
                    linewidth=0.5 if cluster != '不显著' else 0)

# 标注显著城市
for _, row in gdf[gdf['lisa_p'] < 0.05].iterrows():
    ax.annotate(row['CITY'], (row.geometry.x, row.geometry.y),
                fontsize=10, fontweight='bold', xytext=(10, 10),
                textcoords='offset points',
                bbox=dict(boxstyle='round,pad=0.2', facecolor='white', alpha=0.8, edgecolor='gray'))

ax.set_xlim(73, 136)
ax.set_ylim(17, 54)
ax.set_title('波动率 LISA 聚类图', fontsize=16, fontweight='bold')
ax.legend(loc='lower left', fontsize=10, framealpha=0.9)
ax.grid(True, alpha=0.2)
plt.savefig('lisa_volatility.png', dpi=150, bbox_inches='tight')

fig, ax = plt.subplots(1, 1, figsize=(16, 12))
admin_gdf.plot(ax=ax, color='#fafafa', edgecolor='#dddddd', linewidth=0.3)

# 气泡大小 = 波动率，颜色 = 波动率
gdf.plot(ax=ax, column='volatility', cmap='YlOrRd',
         markersize=gdf['volatility'] * 40, alpha=0.8,
         legend=True, edgecolor='black', linewidth=0.3,
         legend_kwds={'label': '波动率', 'orientation': 'horizontal', 'shrink': 0.6})

# 标注 Top 5
for _, row in gdf.nlargest(5, 'volatility').iterrows():
    ax.annotate(f\"{row['CITY']}\\n{row['volatility']*100:.0f}%\",
                (row.geometry.x, row.geometry.y), fontsize=9, fontweight='bold',
                xytext=(10, 10), textcoords='offset points',
                bbox=dict(boxstyle='round,pad=0.2', facecolor='white', alpha=0.8))

ax.set_xlim(73, 136)
ax.set_ylim(17, 54)
ax.set_title('70城房价波动率空间分布', fontsize=16, fontweight='bold')
plt.savefig('volatility_map.png', dpi=150, bbox_inches='tight')

for k in [3, 5, 7]:
    w_k = lps.weights.KNN.from_dataframe(gdf, k=k)
    w_k.transform = 'r'
    lisa_k = Moran_Local(gdf['volatility'].values, w_k)
    sig_mask = lisa_k.p_sim < 0.05
    hh_mask = sig_mask & (lisa_k.q == 1)
    hh_cities = gdf[hh_mask]['CITY'].tolist()
    print(f'k={k}: {sig_mask.sum()} 个显著, HH={hh_cities}')

k=3: 6 个显著, HH=['三亚', '海口', '湛江']
k=5: 8 个显著, HH=['三亚', '北海', '海口', '湛江']
k=7: 11 个显著, HH=['三亚', '北海', '广州', '海口', '深圳', '湛江']