下面代码为第二问所做的工作，请基于问题二的工作，“考虑第2问所反映的从“暴雨”到“灾害”中上述三方面因素的角色及其交互作用，确定暴雨成灾的临界条件；”
下面是原问题参考：
2.	近年来，以暴雨为代表的极端天气事件对人类的生产生活造成了越来越难以忽视的影响。请结合附件中所给的数据，建立数学模型，说明地形-气候相互作用在极端天气形成过程中的作用。
3.	降雨、地形和土地利用对于暴雨等极端天气灾害的形成都具有不可忽视的影响。这其中，降雨的时空变异性和不可控性都最强；土地利用作为自然条件和人类活动的综合结果，虽然也随时空演化，但具有一定可控性；地形是最为稳定、不易改变的因素。请考虑第2问所反映的从“暴雨”到“灾害”中上述三方面因素的角色及其交互作用，确定暴雨成灾的临界条件；并结合第1问中降雨量和土地利用/土地覆被变化的历史时空演化特征，对2025~2035年间中国境内应对暴雨灾害能力最为脆弱的地区进行预测。请以地图的形式呈现你们的预测结果。
# %% [markdown]
# ## 数据准备 相关性

# %%
import pickle
from mgwr.gwr import MGWR
import statsmodels.api as sm
from tqdm import tqdm
import geopandas
import numpy.ma as ma
from matplotlib.colors import Normalize
import cartopy.feature as cfeature
import cartopy.crs as ccrs
from scipy.ndimage import zoom
from splot.esda import plot_moran, moran_scatterplot, lisa_cluster
from mgwr.sel_bw import Sel_BW
from mgwr.gwr import GWR
from shapely.geometry import Point
import geopandas as gpd
from pysal.explore import esda
from pysal.lib import weights
import pysal
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

# 设置中文字体支持
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)

# 2. 读取气温数据并计算平均气温
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))])

sum_data = None
count_data = 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 sum_data is None:
                sum_data = np.zeros(
                    (temp_height, temp_width), dtype=np.float32)
                count_data = np.zeros(
                    (temp_height, temp_width), dtype=np.int32)

            valid_mask = ~np.isnan(temp_data)
            sum_data[valid_mask] += temp_data[valid_mask]
            count_data[valid_mask] += 1

mean_temp_data = np.full_like(sum_data, np.nan, dtype=np.float32)
valid_mask = count_data > 0
mean_temp_data[valid_mask] = sum_data[valid_mask] / count_data[valid_mask]

# 3. 将高程数据重采样到气温数据的网格
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
)

# 4. 读取并处理降水数据
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) | (precip == 0), np.nan, precip)

# 计算平均降水量
mean_precip = np.nanmean(precip, axis=0)  # (lat, lon)

# 关闭数据集
dataset.close()

# 检查平均降水量是否有有效数据
if np.isnan(mean_precip).all():
    print("平均降水量数据全部为空，请检查缺失值处理。")
else:
    print(f"降水平均值最大值: {np.nanmax(mean_precip)}, 最小值: {np.nanmin(mean_precip)}")

# 5. 将降水数据重采样到气温数据的网格
# 创建目标网格的经纬度坐标
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 = mean_precip.flatten()
mask = ~np.isnan(values)
points = points[mask]
values = values[mask]
target_points = np.column_stack((lon_grid.flatten(), lat_grid.flatten()))
mean_precip_resampled = griddata(
    points, values, target_points, method='linear')
mean_precip_resampled = mean_precip_resampled.reshape(lat_grid.shape)

# 6. 进行Spearman相关性分析
# 将数据展开并对齐
elev_flat = dem_resampled.flatten()
temp_flat = mean_temp_data.flatten()
precip_flat = mean_precip_resampled.flatten()

data_stack = np.column_stack((elev_flat, temp_flat, precip_flat))
mask = ~np.isnan(data_stack).any(axis=1)
data_stack = data_stack[mask]

elev_data = data_stack[:, 0]
temp_data = data_stack[:, 1]
precip_data = data_stack[:, 2]

# 创建数据框
df = pd.DataFrame({
    '海拔': elev_data,
    '气温': temp_data,
    '降水': precip_data
})

# 计算Spearman相关矩阵
corr = df.corr(method='spearman')

# 7. 绘制相关矩阵热力图
plt.figure(figsize=(6, 5))
sns.heatmap(corr, annot=True, cmap='coolwarm', square=True)
plt.title('海拔、气温和降水的Spearman相关性矩阵', fontsize=16)
plt.show()


# %% [markdown]
# ## 莫兰 lisa

# %%
# 导入必要的库

# 新增的库
# 空间分析库

# 导入用于绘图的函数

# 设置中文字体支持
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)
w.transform = 'r'  # 权重归一化

# 对高温日频率进行全局莫兰指数分析
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}")

# 绘制莫兰散点图
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()

# **计算高温日频率的局部莫兰指数（LISA）**
high_temp_lisa = esda.Moran_Local(gdf['high_temp_freq'], w)

# **绘制高温日频率的LISA聚类图**
fig, ax = plt.subplots(1, figsize=(10, 6))
lisa_cluster(high_temp_lisa, gdf, p=0.05, ax=ax)
ax.set_title('高温日频率的LISA聚类图')
plt.show()

# **计算暴雨日频率的局部莫兰指数（LISA）**
heavy_rain_lisa = esda.Moran_Local(gdf['heavy_rain_freq'], w)

# **绘制暴雨日频率的LISA聚类图**
fig, ax = plt.subplots(1, figsize=(10, 6))
lisa_cluster(heavy_rain_lisa, gdf, p=0.05, ax=ax)
ax.set_title('暴雨日频率的LISA聚类图')
plt.show()

# %% [markdown]
# ## 数据准备

# %%

# 假设 dem_data 是一个大尺寸的数组
# 计算缩放因子
zoom_factor = (lat_grid.shape[0] / dem_data.shape[0],
               lat_grid.shape[1] / dem_data.shape[1])
dem_data_resampled = zoom(dem_data, zoom_factor)

# 打印每个数组的长度
print(f'lon length: {len(lon_grid.flatten())}')
print(f'lat length: {len(lat_grid.flatten())}')
print(f'high_count length: {len(high_temp_count.flatten())}')
print(f'heavy_count length: {len(heavy_rain_count_resampled.flatten())}')
print(f'elev length: {len(dem_data_resampled.flatten())}')
print(f'temp_mean length: {len(avg_temp_data.flatten())}')
print(f'temp_max length: {len(max_temp_data.flatten())}')
print(f'temp_min length: {len(min_temp_data.flatten())}')
print(f'precip length: {len(heavy_rain_count_resampled.flatten())}')

# 创建 DataFrame，加入气温和降水数据
data = pd.DataFrame({
    'lon': lon_grid.flatten(),
    'lat': lat_grid.flatten(),
    'high_count': high_temp_count.flatten(),
    'heavy_count': heavy_rain_count_resampled.flatten(),
    'elev': dem_data_resampled.flatten(),
    'temp_mean': avg_temp_data.flatten(),
    'temp_max': max_temp_data.flatten(),
    'temp_min': min_temp_data.flatten(),
    'precip': heavy_rain_count_resampled.flatten()  # 降水数据
})

# 去除缺失值
data = data.dropna(subset=['high_count', 'heavy_count', 'elev', 'temp_mean',
                   'temp_max', 'temp_min', 'precip']).reset_index(drop=True)


# %%

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

# --- 读取并处理高程数据 ---
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
total_temp_data = None
min_temp_data = None
max_temp_data = None
temp_transform = None
temp_crs = None
temp_width = None
temp_height = None
temp_bounds = None
temp_res = None
day_count = 0

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:
        day_count += 1
        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)
                total_temp_data = np.zeros(
                    (temp_height, temp_width), dtype=np.float32)
                min_temp_data = temp_data.copy()
                max_temp_data = temp_data.copy()
            else:
                # 更新最小值和最大值
                min_temp_data = np.nanmin(
                    np.array([min_temp_data, temp_data]), axis=0)
                max_temp_data = np.nanmax(
                    np.array([max_temp_data, temp_data]), axis=0)

            # 标记高温日（假设日平均气温超过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_temp_data = np.nansum([total_temp_data, temp_data], axis=0)

# 计算平均气温
avg_temp_data = total_temp_data / day_count

# --- 读取并处理降水数据，识别暴雨日，计算最小值，最大值，平均值 ---
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', -99.9)  # 设定默认缺失值为-99.9

# 处理缺失值
precip = np.where((precip == missing_value) | (precip == 0), 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)

# 计算平均降水、最小值、最大值
avg_precip_data = np.nanmean(precip, axis=0)
min_precip_data = np.nanmin(precip, axis=0)
max_precip_data = np.nanmax(precip, axis=0)

# 关闭数据集
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()))
target_points = np.column_stack((lon_grid.flatten(), lat_grid.flatten()))

# 要重采样的变量
variables_to_resample = {
    'heavy_rain_count': heavy_rain_count,
    'avg_precip_data': avg_precip_data,
    'min_precip_data': min_precip_data,
    'max_precip_data': max_precip_data
}

resampled_variables = {}

for var_name, var_data in variables_to_resample.items():
    values = var_data.flatten()
    mask = ~np.isnan(values)
    points_var = points[mask]
    values_var = values[mask]
    resampled = griddata(points_var, values_var,
                         target_points, method='linear')
    resampled = resampled.reshape(lat_grid.shape)
    resampled_variables[var_name] = resampled

# --- 将高程数据重采样到气温数据的网格 ---
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
)

# 3. 创建 DataFrame，加入气温和降水数据
data = pd.DataFrame({
    'longitude': lon_grid.flatten(),
    'latitude': lat_grid.flatten(),
    'high_temp_count': high_temp_count.flatten(),
    'heavy_rain_count': resampled_variables['heavy_rain_count'].flatten(),
    'elevation': dem_resampled.flatten(),
    'avg_temp': avg_temp_data.flatten(),
    'min_temp': min_temp_data.flatten(),
    'max_temp': max_temp_data.flatten(),
    'avg_precip': resampled_variables['avg_precip_data'].flatten(),
    'min_precip': resampled_variables['min_precip_data'].flatten(),
    'max_precip': resampled_variables['max_precip_data'].flatten()
})

# 去除高温日和暴雨日均为0的格点
data = data[~((data['high_temp_count'] == 0) &
              (data['heavy_rain_count'] == 0))]

# 去除缺失值
data = data.dropna(subset=['high_temp_count', 'heavy_rain_count', 'elevation', 'avg_temp',
                   'min_temp', 'max_temp', 'avg_precip', 'min_precip', 'max_precip']).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')

# 输出 gdf 的样本
print(gdf.head())

# 保存 gdf 以便后续使用
gdf.to_file("gdf_with_temp_precip.shp")

# 4. 绘制高温日与暴雨日的频数地图
# 绘制高温日频数地图
plt.figure(figsize=(10, 6))
plt.pcolormesh(lon_grid, lat_grid, high_temp_count, cmap='hot', shading='auto')
plt.colorbar(label='高温日天数')
plt.title('高温日频数分布图')
plt.xlabel('经度')
plt.ylabel('纬度')
plt.show()

# 绘制暴雨日频数地图
plt.figure(figsize=(10, 6))
plt.pcolormesh(lon_grid, lat_grid,
               resampled_variables['heavy_rain_count'], cmap='Blues', shading='auto')
plt.colorbar(label='暴雨日天数')
plt.title('暴雨日频数分布图')
plt.xlabel('经度')
plt.ylabel('纬度')
plt.show()

# 5. 导出频数数据到 CSV 文件
# 导出高温日频数数据
high_temp_df = data[['longitude', 'latitude', 'high_temp_count']]
high_temp_df.to_csv('high_temp_frequency.csv', index=False)

# 导出暴雨日频数数据
heavy_rain_df = data[['longitude', 'latitude', 'heavy_rain_count']]
heavy_rain_df.to_csv('heavy_rain_frequency.csv', index=False)

# 导出包含所有信息的 CSV 文件
data.to_csv('temperature_precipitation_data.csv', index=False)

# %%

# 假设 high_temp_count 和 heavy_rain_count 有 NaN 值
# 将 NaN 值屏蔽
heavy_rain_count_masked = ma.masked_invalid(
    resampled_variables['heavy_rain_count'])
high_temp_count_masked = ma.masked_invalid(high_temp_count)

# 定义经纬度范围，适当缩小范围
lon_min, lon_max = 77, 136  # 设定经度范围
lat_min, lat_max = 18, 53   # 设定纬度范围

# 创建图形和坐标轴
fig, ax = plt.subplots(figsize=(10, 6), subplot_kw={
                       'projection': ccrs.PlateCarree()})

# 设置地图显示的经纬度范围
ax.set_extent([lon_min, lon_max, lat_min, lat_max], crs=ccrs.PlateCarree())

# 添加地理特征
ax.add_feature(cfeature.LAND.with_scale('50m'), facecolor='lightgray')
ax.add_feature(cfeature.OCEAN.with_scale('50m'), facecolor='lightblue')
ax.add_feature(cfeature.COASTLINE.with_scale('50m'), linewidth=0.5)
ax.add_feature(cfeature.BORDERS.with_scale(
    '50m'), linestyle=':', linewidth=0.5)
ax.add_feature(cfeature.LAKES.with_scale('50m'),
               facecolor='lightblue', alpha=0.5)
ax.add_feature(cfeature.RIVERS.with_scale('50m'), linewidth=0.5)

# 设置最大值截断，增强频数图的对比度
vmax_value_rain = 1000  # 根据实际数据设置一个合理的最大值
norm_rain = Normalize(vmin=0, vmax=vmax_value_rain)

# 绘制暴雨日频数地图，解决奇怪的框框问题
c1 = ax.pcolormesh(lon_grid, lat_grid, heavy_rain_count_masked, cmap='Blues',
                   shading='auto', norm=norm_rain, transform=ccrs.PlateCarree())
plt.colorbar(c1, ax=ax, label='暴雨日天数')

# 设置暴雨日图的标题
ax.set_title('暴雨日频数分布图')
ax.set_xlabel('经度')
ax.set_ylabel('纬度')

plt.show()

# 创建新的图形和坐标轴用于高温日频数地图
fig, ax = plt.subplots(figsize=(10, 6), subplot_kw={
                       'projection': ccrs.PlateCarree()})

# 设置地图显示的经纬度范围
ax.set_extent([lon_min, lon_max, lat_min, lat_max], crs=ccrs.PlateCarree())

# 添加地理特征
ax.add_feature(cfeature.LAND.with_scale('50m'), facecolor='lightgray')
ax.add_feature(cfeature.OCEAN.with_scale('50m'), facecolor='lightblue')
ax.add_feature(cfeature.COASTLINE.with_scale('50m'), linewidth=0.5)
ax.add_feature(cfeature.BORDERS.with_scale(
    '50m'), linestyle=':', linewidth=0.5)
ax.add_feature(cfeature.LAKES.with_scale('50m'),
               facecolor='lightblue', alpha=0.5)
ax.add_feature(cfeature.RIVERS.with_scale('50m'), linewidth=0.5)

# 设置最大值截断，增强高温频数图的对比度
vmax_value_temp = 5  # 根据实际数据设置一个合理的最大值
norm_temp = Normalize(vmin=0, vmax=vmax_value_temp)

# 绘制高温日频数地图
c2 = ax.pcolormesh(lon_grid, lat_grid, high_temp_count_masked, cmap='hot_r',
                   shading='auto', norm=norm_temp, transform=ccrs.PlateCarree())
plt.colorbar(c2, ax=ax, label='高温日天数')

# 设置高温日图的标题
ax.set_title('高温日频数分布图')
ax.set_xlabel('经度')
ax.set_ylabel('纬度')

plt.show()


# %%

# 读取保存的 GeoDataFrame
gdf = geopandas.read_file("gdf_with_temp_precip.shp")

# 查看数据
print(gdf.columns)
print(gdf.head())


# %%
# 保存为 CSV
gdf.to_csv("gdf_with_temp_precip.csv", index=False)

# 查看前几行数据
print(gdf.head(100))


# %%

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

# 假设 data 已经存在，包含您的数据
# 请确保在运行此代码之前定义了 data DataFrame
# data = pd.DataFrame({...})

# 将经度和纬度转换为几何对象
data['geometry'] = [Point(xy)
                    for xy in zip(data['longitude'], data['latitude'])]
gdf = gpd.GeoDataFrame(data, geometry='geometry')

# 1. 优化数据类型
gdf['elevation'] = gdf['elevation'].astype('float32')
gdf['avg_temp'] = gdf['avg_temp'].astype('float32')
gdf['min_temp'] = gdf['min_temp'].astype('float32')
gdf['max_temp'] = gdf['max_temp'].astype('float32')
gdf['avg_precip'] = gdf['avg_precip'].astype('float32')
gdf['min_precip'] = gdf['min_precip'].astype('float32')
gdf['max_precip'] = gdf['max_precip'].astype('float32')
gdf['high_temp_count'] = gdf['high_temp_count'].astype('float32')
gdf['heavy_rain_count'] = gdf['heavy_rain_count'].astype('float32')

# 2. 减少数据精度
gdf['longitude'] = gdf['longitude'].round(4)
gdf['latitude'] = gdf['latitude'].round(4)
gdf['elevation'] = gdf['elevation'].round(2)
gdf['avg_temp'] = gdf['avg_temp'].round(2)
gdf['min_temp'] = gdf['min_temp'].round(2)
gdf['max_temp'] = gdf['max_temp'].round(2)
gdf['avg_precip'] = gdf['avg_precip'].round(2)
gdf['min_precip'] = gdf['min_precip'].round(2)
gdf['max_precip'] = gdf['max_precip'].round(2)
gdf['high_temp_count'] = gdf['high_temp_count'].round(2)
gdf['heavy_rain_count'] = gdf['heavy_rain_count'].round(2)

# 打印去除缺失值前的数据行数
print(f"去除缺失值前数据行数: {gdf.shape[0]}")

# 去除缺失值
gdf = gdf.dropna(subset=['high_temp_count', 'heavy_rain_count', 'elevation', 'avg_temp',
                 'min_temp', 'max_temp', 'avg_precip', 'min_precip', 'max_precip']).reset_index(drop=True)

# 去除缺失值后的数据行数
print(f"去除缺失值后数据行数: {gdf.shape[0]}")

# 检查去除缺失值后数据是否为空
if gdf.shape[0] == 0:
    raise ValueError("去除缺失值后数据为空，请检查数据。")

# 3. 使用部分数据进行采样
gdf_sampled = gdf.sample(frac=0.00228, random_state=42).reset_index(drop=True)

# 提取坐标
coords = list(zip(gdf_sampled.geometry.x, gdf_sampled.geometry.y))

# 提取自变量和因变量（高温模型）
# 包含“降水最大值”“降水最小值”
X_high_temp = gdf_sampled[['elevation',
                           'avg_precip', 'min_precip', 'max_precip']]
X_high_temp = sm.add_constant(X_high_temp)  # 添加常数项
Y_high_temp = gdf_sampled['high_temp_count'].values.reshape(-1, 1)  # 调整维度

# 提取自变量和因变量（降水模型）
# 包含“温度最大值”“温度最小值”
X_heavy_rain = gdf_sampled[['elevation', 'avg_temp', 'min_temp', 'max_temp']]
X_heavy_rain = sm.add_constant(X_heavy_rain)  # 添加常数项
Y_heavy_rain = gdf_sampled['heavy_rain_count'].values.reshape(-1, 1)  # 调整维度

# 检查数据形状，确保没有维度问题
print(f"X_high_temp shape: {X_high_temp.shape}")  # 自变量矩阵
print(f"Y_high_temp shape: {Y_high_temp.shape}")  # 因变量
print(f"X_heavy_rain shape: {X_heavy_rain.shape}")  # 自变量矩阵
print(f"Y_heavy_rain shape: {Y_heavy_rain.shape}")  # 因变量

# 检查是否有缺失值
if np.any(np.isnan(X_high_temp)) or np.any(np.isnan(Y_high_temp)):
    raise ValueError("高温模型输入数据中存在缺失值，请处理缺失值。")

if np.any(np.isnan(X_heavy_rain)) or np.any(np.isnan(Y_heavy_rain)):
    raise ValueError("降水模型输入数据中存在缺失值，请处理缺失值。")

# 确保行数大于列数
if X_high_temp.shape[0] <= X_high_temp.shape[1]:
    raise ValueError("高温模型自变量矩阵的行数应大于列数，请检查输入数据。")

if X_heavy_rain.shape[0] <= X_heavy_rain.shape[1]:
    raise ValueError("降水模型自变量矩阵的行数应大于列数，请检查输入数据。")

# 获取样本数量
n_samples = gdf_sampled.shape[0]

# 设置带宽的最小值和最大值
bw_min_ht = int(n_samples * 0.05)  # 最小带宽为样本数的 5%
bw_max_ht = int(n_samples * 0.3)   # 最大带宽为样本数的 30%

bw_min_hr = int(n_samples * 0.05)  # 最小带宽为样本数的 5%
bw_max_hr = int(n_samples * 0.3)   # 最大带宽为样本数的 30%

# 4. 选择带宽（高温模型）
print("选择高温日频率带宽中...")
selector_ht = Sel_BW(coords, Y_high_temp, X_high_temp, multi=True)
bws_ht = selector_ht.search(bw_min=bw_min_ht, bw_max=bw_max_ht)
print(f"高温日频率的 MGWR 模型最佳带宽：{bws_ht}")

# 5. 构建 MGWR 模型（高温模型）
print("构建高温日频率 MGWR 模型中...")
mgwr_model_ht = MGWR(coords, Y_high_temp, X_high_temp,
                     selector_ht, n_jobs=-1)  # 启用并行计算
mgwr_results_ht = mgwr_model_ht.fit()

print(mgwr_results_ht.summary())

# 6. 提取 MGWR 结果（高温模型）
gdf_sampled['mgwr_const_coef_ht'] = mgwr_results_ht.params[:, 0]
gdf_sampled['mgwr_elevation_coef_ht'] = mgwr_results_ht.params[:, 1]
gdf_sampled['mgwr_avg_precip_coef_ht'] = mgwr_results_ht.params[:, 2]
gdf_sampled['mgwr_min_precip_coef_ht'] = mgwr_results_ht.params[:, 3]
gdf_sampled['mgwr_max_precip_coef_ht'] = mgwr_results_ht.params[:, 4]

# 7. 绘制高温日频率 MGWR 模型的系数分布
# 高程系数
fig, ax = plt.subplots(1, 1, figsize=(8, 6))
gdf_sampled.plot(column='mgwr_elevation_coef_ht',
                 ax=ax, cmap='coolwarm', legend=True)
ax.set_title('高温日频率 MGWR 模型的高程系数分布')
plt.xlabel('经度')
plt.ylabel('纬度')
plt.show()

# 平均降水系数
fig, ax = plt.subplots(1, 1, figsize=(8, 6))
gdf_sampled.plot(column='mgwr_avg_precip_coef_ht',
                 ax=ax, cmap='coolwarm', legend=True)
ax.set_title('高温日频率 MGWR 模型的平均降水系数分布')
plt.xlabel('经度')
plt.ylabel('纬度')
plt.show()

# 最小降水系数
fig, ax = plt.subplots(1, 1, figsize=(8, 6))
gdf_sampled.plot(column='mgwr_min_precip_coef_ht',
                 ax=ax, cmap='coolwarm', legend=True)
ax.set_title('高温日频率 MGWR 模型的最小降水系数分布')
plt.xlabel('经度')
plt.ylabel('纬度')
plt.show()

# 最大降水系数
fig, ax = plt.subplots(1, 1, figsize=(8, 6))
gdf_sampled.plot(column='mgwr_max_precip_coef_ht',
                 ax=ax, cmap='coolwarm', legend=True)
ax.set_title('高温日频率 MGWR 模型的最大降水系数分布')
plt.xlabel('经度')
plt.ylabel('纬度')
plt.show()

# 8. 选择带宽（降水模型）
print("选择暴雨日频率带宽中...")
selector_hr = Sel_BW(coords, Y_heavy_rain, X_heavy_rain, multi=True)
bws_hr = selector_hr.search(bw_min=bw_min_hr, bw_max=bw_max_hr)
print(f"暴雨日频率的 MGWR 模型最佳带宽：{bws_hr}")

# 9. 构建 MGWR 模型（降水模型）
print("构建暴雨日频率 MGWR 模型中...")
mgwr_model_hr = MGWR(coords, Y_heavy_rain, X_heavy_rain,
                     selector_hr, n_jobs=-1)  # 启用并行计算
mgwr_results_hr = mgwr_model_hr.fit()

print(mgwr_results_hr.summary())

# 10. 提取 MGWR 结果（降水模型）
gdf_sampled['mgwr_const_coef_hr'] = mgwr_results_hr.params[:, 0]
gdf_sampled['mgwr_elevation_coef_hr'] = mgwr_results_hr.params[:, 1]
gdf_sampled['mgwr_avg_temp_coef_hr'] = mgwr_results_hr.params[:, 2]
gdf_sampled['mgwr_min_temp_coef_hr'] = mgwr_results_hr.params[:, 3]
gdf_sampled['mgwr_max_temp_coef_hr'] = mgwr_results_hr.params[:, 4]

# 11. 绘制暴雨日频率 MGWR 模型的系数分布
# 高程系数
fig, ax = plt.subplots(1, 1, figsize=(8, 6))
gdf_sampled.plot(column='mgwr_elevation_coef_hr',
                 ax=ax, cmap='coolwarm', legend=True)
ax.set_title('暴雨日频率 MGWR 模型的高程系数分布')
plt.xlabel('经度')
plt.ylabel('纬度')
plt.show()

# 平均温度系数
fig, ax = plt.subplots(1, 1, figsize=(8, 6))
gdf_sampled.plot(column='mgwr_avg_temp_coef_hr',
                 ax=ax, cmap='coolwarm', legend=True)
ax.set_title('暴雨日频率 MGWR 模型的平均温度系数分布')
plt.xlabel('经度')
plt.ylabel('纬度')
plt.show()

# 最小温度系数
fig, ax = plt.subplots(1, 1, figsize=(8, 6))
gdf_sampled.plot(column='mgwr_min_temp_coef_hr',
                 ax=ax, cmap='coolwarm', legend=True)
ax.set_title('暴雨日频率 MGWR 模型的最小温度系数分布')
plt.xlabel('经度')
plt.ylabel('纬度')
plt.show()

# 最大温度系数
fig, ax = plt.subplots(1, 1, figsize=(8, 6))
gdf_sampled.plot(column='mgwr_max_temp_coef_hr',
                 ax=ax, cmap='coolwarm', legend=True)
ax.set_title('暴雨日频率 MGWR 模型的最大温度系数分布')
plt.xlabel('经度')
plt.ylabel('纬度')
plt.show()

# %% [markdown]
# ## 保存参数

# %%
# # 保存模型
# with open('mgwr_model.pkl', 'wb') as f:
#     pickle.dump(mgwr_results, f)

# 如果希望同时保存自变量、因变量以及带宽信息，可以一起保存
with open('mgwr_temp_model_and_data.pkl', 'wb') as f:
    pickle.dump((mgwr_model_ht, mgwr_results_ht,
                X, Y_high_temp, selector_ht), f)

with open('mgwr_rain_model_and_data.pkl', 'wb') as f:
    pickle.dump((mgwr_model_hr, mgwr_results_hr,
                X, Y_heavy_rain, selector_hr), f)