把下面的莫兰指数代码加一个LISA聚类图（输出完整代码）
# 导入必要的库
import os
import numpy as np
import rasterio
from rasterio.warp import reproject, Resampling
import matplotlib.pyplot as plt
import netCDF4 as nc
from scipy.interpolate import griddata
import pandas as pd
import seaborn as sns

# 新增的库
# 空间分析库
import pysal
from pysal.lib import weights
from pysal.explore import esda
import geopandas as gpd
from shapely.geometry import Point
from mgwr.gwr import GWR
from mgwr.sel_bw import Sel_BW

# **导入 plot_moran**
from splot.esda import plot_moran

# 设置中文字体支持
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

# 1. 定义高温日和暴雨事件
# 我们将基于之前的代码，修改和扩展以实现新的任务

# --- 读取并处理高程数据（与之前相同） ---
dem_file_path = r"F:\HuaweiCup\DATA\data1\Geo\TIFF\chinadem_geo.tif"

with rasterio.open(dem_file_path) as dem_src:
    dem_data = dem_src.read(1).astype(np.float32)
    dem_transform = dem_src.transform
    dem_crs = dem_src.crs
    dem_nodata = dem_src.nodata
    dem_bounds = dem_src.bounds
    dem_res = dem_src.res

    # 处理缺失值
    dem_data = np.where((dem_data == dem_nodata) | (dem_data == -99.9) | (dem_data == 0), np.nan, dem_data)

# --- 读取气温数据，识别高温日 ---
data_root_dir = r"F:\HuaweiCup\DATA\data2\day_avg"
year_dirs = sorted([os.path.join(data_root_dir, d) for d in os.listdir(data_root_dir) if os.path.isdir(os.path.join(data_root_dir, d))])

# 初始化高温日计数
high_temp_count = None
temp_transform = None  # 初始化
temp_crs = None
temp_width = None
temp_height = None
temp_bounds = None
temp_res = None

for year_dir in year_dirs:
    tif_files = sorted([os.path.join(year_dir, f) for f in os.listdir(year_dir) if f.endswith('_avg.tif')])

    for tif_file in tif_files:
        with rasterio.open(tif_file) as temp_src:
            temp_data = temp_src.read(1).astype(np.float32)
            temp_nodata = temp_src.nodata

            # 在第一次迭代时，获取仿射变换和其他属性
            if temp_transform is None:
                temp_transform = temp_src.transform
                temp_crs = temp_src.crs
                temp_width = temp_src.width
                temp_height = temp_src.height
                temp_bounds = temp_src.bounds
                temp_res = temp_src.res

            # 处理缺失值
            temp_data = np.where((temp_data == temp_nodata) | (temp_data == -99.9) | (temp_data == 0), np.nan, temp_data)

            # 初始化高温日计数矩阵
            if high_temp_count is None:
                high_temp_count = np.zeros((temp_height, temp_width), dtype=np.int32)

            # 标记高温日（假设日平均气温超过35℃）
            high_temp_mask = temp_data > 35.0
            high_temp_mask = np.where(np.isnan(temp_data), False, high_temp_mask)

            # 更新高温日计数
            high_temp_count[high_temp_mask] += 1

# 计算高温日的总天数
total_days = len(year_dirs) * 365  # 简单假设每年365天
# 或者实际计算天数
# total_days = high_temp_count.max()

# 计算高温日频率
high_temp_frequency = high_temp_count / total_days

# --- 读取并处理降水数据，识别暴雨日 ---
precip_file_path = r"F:\HuaweiCup\DATA\data3\CHM_PRE_0.25dg_19612022.nc\CHM_PRE_0.25dg_19612022.nc"

# 打开 NetCDF 文件
dataset = nc.Dataset(precip_file_path, mode='r')

# 获取数据
precip_var = dataset.variables['pre']
precip = precip_var[:]  # (time, lat, lon)
latitude = dataset.variables['latitude'][:]
longitude = dataset.variables['longitude'][:]

# 获取缺失值标识符
missing_value = getattr(precip_var, 'missing_value', None)
if missing_value is None:
    missing_value = getattr(precip_var, '_FillValue', None)
if missing_value is None:
    missing_value = -99.9  # 默认缺失值
print(f"降水数据使用的缺失值标识符: {missing_value}")

# 确保降水数据为浮点型
precip = precip.astype(np.float32)

# 处理缺失值
precip = np.where((precip == missing_value), np.nan, precip)

# 初始化暴雨日计数矩阵
precip_height, precip_width = len(latitude), len(longitude)
heavy_rain_count = np.zeros((precip_height, precip_width), dtype=np.int32)

# 标记暴雨日（降水量超过50mm/day）
heavy_rain_mask = precip > 50.0
heavy_rain_mask = np.where(np.isnan(precip), False, heavy_rain_mask)

# 按时间维度累加暴雨日计数
heavy_rain_count = np.nansum(heavy_rain_mask, axis=0)

# 计算暴雨日的总天数
total_precip_days = precip.shape[0]  # 时间维度的长度

# 计算暴雨日频率
heavy_rain_frequency = heavy_rain_count / total_precip_days

# 关闭数据集
dataset.close()

# --- 将降水数据重采样到气温数据的网格 ---
# 创建目标网格的经纬度坐标
lon_min, lat_min, lon_max, lat_max = temp_bounds.left, temp_bounds.bottom, temp_bounds.right, temp_bounds.top
lon_array = np.linspace(lon_min + temp_res[0]/2, lon_max - temp_res[0]/2, temp_width)
lat_array = np.linspace(lat_max - temp_res[1]/2, lat_min + temp_res[1]/2, temp_height)
lon_grid, lat_grid = np.meshgrid(lon_array, lat_array)

# 降水数据的经纬度坐标
precip_lon_grid, precip_lat_grid = np.meshgrid(longitude, latitude)

# 插值暴雨日频率
points = np.column_stack((precip_lon_grid.flatten(), precip_lat_grid.flatten()))
values = heavy_rain_frequency.flatten()
mask = ~np.isnan(values)
points = points[mask]
values = values[mask]
target_points = np.column_stack((lon_grid.flatten(), lat_grid.flatten()))
heavy_rain_freq_resampled = griddata(points, values, target_points, method='linear')
heavy_rain_freq_resampled = heavy_rain_freq_resampled.reshape(lat_grid.shape)

# --- 将高程数据重采样到气温数据的网格（与之前相同） ---
dem_resampled = np.zeros((temp_height, temp_width), dtype=np.float32)

reproject(
    source=dem_data,
    destination=dem_resampled,
    src_transform=dem_transform,
    src_crs=dem_crs,
    dst_transform=temp_transform,
    dst_crs=temp_crs,
    resampling=Resampling.bilinear
)

# 2. 对高温、暴雨天气进行莫兰指数分析

# 创建 DataFrame
data = pd.DataFrame({
    'longitude': lon_grid.flatten(),
    'latitude': lat_grid.flatten(),
    'high_temp_freq': high_temp_frequency.flatten(),
    'heavy_rain_freq': heavy_rain_freq_resampled.flatten(),
    'elevation': dem_resampled.flatten()
})

# 去除缺失值
data = data.dropna(subset=['high_temp_freq', 'heavy_rain_freq', 'elevation']).reset_index(drop=True)

# 在删除缺失值后，创建 geometry 列
geometry = [Point(xy) for xy in zip(data['longitude'], data['latitude'])]

# 转换为 GeoDataFrame
gdf = gpd.GeoDataFrame(data, geometry=geometry, crs='EPSG:4326')

# 生成空间权重矩阵（例如，使用 Queen 邻接方式）
w = weights.contiguity.Queen.from_dataframe(gdf)

# 对高温日频率进行全局莫兰指数分析
high_temp_moran = esda.Moran(gdf['high_temp_freq'], w)
print(f"高温日频率的莫兰指数：{high_temp_moran.I}, p-value: {high_temp_moran.p_sim}")

# 对暴雨日频率进行全局莫兰指数分析
heavy_rain_moran = esda.Moran(gdf['heavy_rain_freq'], w)
print(f"暴雨日频率的莫兰指数：{heavy_rain_moran.I}, p-value: {heavy_rain_moran.p_sim}")

# **绘制莫兰散点图**
from splot.esda import moran_scatterplot

fig, axes = plt.subplots(1, 2, figsize=(12, 5))
moran_scatterplot(high_temp_moran, ax=axes[0], scatter_kwds={'s':5})
axes[0].set_title('高温日频率的莫兰散点图')
moran_scatterplot(heavy_rain_moran, ax=axes[1], scatter_kwds={'s':5})
axes[1].set_title('暴雨日频率的莫兰散点图')
plt.show()