1、是不是有什么限制了纬度？为什么红点到了大约北纬33就全没了？切割的还特别整齐
2、最后的绘图，红点使用25-35年的预测的“暴雨成灾”频数来作深浅变化。
3、问题改完后输出全部代码
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import os
from sklearn.model_selection import train_test_split, GridSearchCV, StratifiedKFold
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import classification_report, confusion_matrix, accuracy_score
from sklearn.preprocessing import StandardScaler, MinMaxScaler
import warnings

warnings.filterwarnings("ignore")

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

# 导入 LSTM 所需的库
from tensorflow.keras.models import Sequential
from tensorflow.keras.layers import LSTM, Dense, Dropout

# 需要安装 geopandas 和 contextily 库，如果未安装，请执行以下命令：
# pip install geopandas contextily

import geopandas as gpd
import contextily as ctx

# 1. 数据读取和预处理
# -----------------------

# 定义数据文件路径
data_file = r'F:\HuaweiCup\Q3\disaster_index.csv'

# 读取数据
data = pd.read_csv(data_file)

# 显示数据的基本信息
print("数据集基本信息：")
print(data.info())
print("\n数据集的前5行：")
print(data.head())

# 将 'disaster_index' 转化为二分类标签
# 设置阈值，例如，根据数据的分布设定阈值
threshold = data['disaster_index'].quantile(0.75)  # 使用第75百分位数作为阈值
data['disaster'] = data['disaster_index'].apply(lambda x: 1 if x > threshold else 0)

# 删除不需要的列，例如 'disaster_index'
data.drop(['disaster_index'], axis=1, inplace=True)

# 将 'land_use' 转换为分类变量（如果不是数值型）
if data['land_use'].dtype == 'object' or str(data['land_use'].dtype) == 'category':
    data['land_use'] = data['land_use'].astype('category').cat.codes

# 检查是否存在缺失值
print("\n缺失值统计：")
print(data.isnull().sum())

# 如果有缺失值，可以选择填充，这里我们使用均值填充
data.fillna(data.mean(), inplace=True)

# 处理类别不平衡问题（可选）
# 查看类别分布
print("\n类别分布：")
print(data['disaster'].value_counts())

# 如果类别不平衡严重，可以考虑使用下采样或过采样方法平衡数据
from imblearn.over_sampling import SMOTE

X = data.drop('disaster', axis=1)
y = data['disaster']

# 使用 SMOTE 进行过采样
smote = SMOTE(random_state=42)
X_res, y_res = smote.fit_resample(X, y)

# 2. 特征选择和数据集划分
# ----------------------------

# 特征列和目标列
features = ['longitude', 'latitude', 'elevation', 'slope', 'aspect',
            'max_temperature', 'max_precipitation', 'land_use',
            'population', 'GDP']

target = 'disaster'

# 提取特征和目标变量
X = X_res[features]
y = y_res

# 特征标准化
scaler = StandardScaler()
X[features] = scaler.fit_transform(X[features])

# 将数据集划分为训练集和测试集，比例为 80% : 20%
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, stratify=y, random_state=42)

# 3. 模型训练
# ----------------------------

# 定义随机森林模型，参数可以根据需要调整
rf_model = RandomForestClassifier(random_state=42)

# 设定参数网格，用于网格搜索
param_grid = {
    'n_estimators': [50, 100, 150],            # 决策树的个数
    'max_features': ['sqrt', 'log2'],          # 特征选择的最大数量
    'max_depth': [5, 10, 15],                  # 决策树的最大深度
    'min_samples_split': [2, 5, 10],           # 内部节点再划分所需最小样本数
    'min_samples_leaf': [1, 2, 4],             # 叶子节点最小样本数
    'criterion': ['gini', 'entropy']           # 划分标准
}

# 使用交叉验证寻找最佳参数
cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)
grid_search = GridSearchCV(estimator=rf_model,
                           param_grid=param_grid,
                           cv=cv,
                           n_jobs=-1,
                           scoring='accuracy')
# 训练模型
grid_search.fit(X_train, y_train)

# 输出最佳参数
print("\n最佳参数：")
print(grid_search.best_params_)

# 使用最佳参数建立模型
best_rf_model = grid_search.best_estimator_

# 输出特征重要性
importances = best_rf_model.feature_importances_
feature_importance = pd.Series(
    importances, index=features).sort_values(ascending=False)
print("\n特征重要性：")
print(feature_importance)

# 可视化特征重要性
plt.figure(figsize=(10, 6))
feature_importance.plot(kind='bar')
plt.title('特征重要性')
plt.ylabel('重要性')
plt.tight_layout()
plt.show()

# 4. 模型评估
# ----------------------------

# 在测试集上进行预测
y_pred = best_rf_model.predict(X_test)

# 输出分类报告
print("\n模型评估报告：")
print(classification_report(y_test, y_pred))

# 输出混淆矩阵
confusion = confusion_matrix(y_test, y_pred)
print("\n混淆矩阵：")
print(confusion)

# 计算准确率
accuracy = accuracy_score(y_test, y_pred)
print("\n模型准确率： {:.2f}%".format(accuracy * 100))

# 绘制学习曲线
from sklearn.model_selection import learning_curve

train_sizes, train_scores, test_scores = learning_curve(
    estimator=best_rf_model,
    X=X_train,
    y=y_train,
    cv=cv,
    scoring='accuracy',
    n_jobs=-1,
    train_sizes=np.linspace(0.1, 1.0, 10),
    shuffle=True,
    random_state=42
)

train_scores_mean = np.mean(train_scores, axis=1)
test_scores_mean = np.mean(test_scores, axis=1)

plt.figure(figsize=(10, 6))
plt.plot(train_sizes, train_scores_mean, 'o-', color='r', label='训练集分数')
plt.plot(train_sizes, test_scores_mean, 'o-', color='g', label='交叉验证集分数')
plt.xlabel('训练样本数')
plt.ylabel('得分')
plt.title('学习曲线')
plt.legend(loc='best')
plt.show()

# 5. 未来数据预测
# ----------------------------

# 使用 LSTM 预测未来特征值

# 定义需要预测的特征列表（时间序列特征）
time_series_features = ['max_temperature', 'max_precipitation',
                        'population', 'GDP']

# 将年份转换为日期格式，方便 LSTM 使用
data['date'] = pd.to_datetime(data['year'], format='%Y')

# 创建一个字典，存储每个特征的预测结果
future_feature_predictions = {}

# 定义 LSTM 模型
def create_lstm_model(input_shape):
    model = Sequential()
    model.add(LSTM(50, return_sequences=True, input_shape=input_shape))
    model.add(Dropout(0.2))
    model.add(LSTM(50))
    model.add(Dropout(0.2))
    model.add(Dense(1))
    model.compile(optimizer='adam', loss='mean_squared_error')
    return model

for feature in time_series_features:
    print(f"\n正在预测特征：{feature}")

    # 准备数据
    df_feature = data.groupby('date')[feature].mean().reset_index()
    df_feature.set_index('date', inplace=True)

    # 数据归一化
    lstm_scaler = MinMaxScaler(feature_range=(0, 1))
    scaled_data = lstm_scaler.fit_transform(df_feature)

    # 定义时间步长
    time_step = 10  # 可以调整时间步长

    # 创建输入输出数据集
    X = []
    y = []
    for i in range(len(scaled_data) - time_step):
        X.append(scaled_data[i:i + time_step, 0])
        y.append(scaled_data[i + time_step, 0])

    X = np.array(X)
    y = np.array(y)

    # 调整输入形状
    X = X.reshape(X.shape[0], X.shape[1], 1)

    # 拆分训练集
    train_size = int(len(X) * 0.8)
    X_train, y_train = X[:train_size], y[:train_size]

    # 创建和训练LSTM模型
    model = create_lstm_model((X_train.shape[1], 1))
    model.fit(X_train, y_train, epochs=1250, batch_size=32, verbose=0)

    # 预测未来数据
    n_future = 2035 - data['year'].max()
    future_predictions = []

    # 获取最后的 time_step 输入
    last_sequence = scaled_data[-time_step:].reshape(1, time_step, 1)
    # 预测未来年份
    for _ in range(n_future):
        next_pred = model.predict(last_sequence)[0, 0]
        future_predictions.append(next_pred)
        # 更新序列
        last_sequence = np.append(last_sequence[:,1:,:], [[[next_pred]]], axis=1)

    # 反归一化
    future_predictions = lstm_scaler.inverse_transform(np.array(future_predictions).reshape(-1, 1))

    # 生成未来年份
    last_year = data['year'].max()
    future_years = np.arange(last_year + 1, 2036)

    # 创建预测结果 DataFrame
    forecast_feature = pd.DataFrame({'year': future_years, 'yhat': future_predictions.flatten()})

    # 将预测结果存入字典
    future_feature_predictions[feature] = forecast_feature

    # 可视化预测结果（可选）
    # 绘制历史数据和未来预测数据
    df_feature_reset = df_feature.reset_index()
    plt.figure(figsize=(12, 6))
    plt.plot(df_feature_reset['date'], df_feature_reset[feature], label='历史数据')
    plt.plot(pd.to_datetime(forecast_feature['year'], format='%Y'), forecast_feature['yhat'], label='预测数据', linestyle='--')
    plt.title(f"{feature} 的预测结果")
    plt.xlabel('年份')
    plt.ylabel(feature)
    plt.legend()
    plt.show()

print("\n预测完成。")

# 构建未来数据集
# 首先，获取数据集中的经纬度范围，生成未来年份的经纬度点
min_lon, max_lon = data['longitude'].min(), data['longitude'].max()
min_lat, max_lat = data['latitude'].min(), data['latitude'].max()

# 定义未来年份
future_years = np.arange(data['year'].max() + 1, 2036)

# 生成经纬度网格点（根据实际需求调整分辨率）
grid_res = 1.0  # 1度分辨率
lons = np.arange(min_lon, max_lon + grid_res, grid_res)
lats = np.arange(min_lat, max_lat + grid_res, grid_res)
lon_grid, lat_grid = np.meshgrid(lons, lats)
grid_points = np.vstack([lon_grid.ravel(), lat_grid.ravel()]).T

# 构建未来数据集
future_data_list = []
for year in future_years:
    for lon, lat in grid_points:
        future_data_list.append([lon, lat, year])

future_df = pd.DataFrame(future_data_list, columns=['longitude', 'latitude', 'year'])

# 为未来数据添加其他特征
# 对于非时间序列特征，仍然采用平均值或最近年份的数据
static_features = ['elevation', 'slope', 'aspect', 'land_use']

# 获取静态特征的平均值
average_static_features = data.groupby(['longitude', 'latitude'])[static_features].mean().reset_index()

# 将静态特征合并到未来数据集中
future_df = pd.merge(future_df, average_static_features, on=['longitude', 'latitude'], how='left')

# 将时间序列预测的特征添加到未来数据集中
for feature in time_series_features:
    # 获取该特征的预测结果
    forecast_feature = future_feature_predictions[feature]
    # 将预测结果映射到相应的年份
    forecast_feature['year'] = forecast_feature['year'].astype(int)
    # 将预测特征与未来数据集合并
    future_df = pd.merge(future_df, forecast_feature[['year', 'yhat']], on='year', how='left')
    # 重命名列
    future_df.rename(columns={'yhat': feature}, inplace=True)

# 剔除缺失值
future_df.dropna(inplace=True)

# 提取特征
X_future = future_df[features]

# 对特征进行标准化
X_future[features] = scaler.transform(X_future[features])

# 进行预测
future_predictions = best_rf_model.predict(X_future)

# 将预测结果添加到数据集中
future_df['disaster_prediction'] = future_predictions

# 筛选预测为高风险（暴雨灾害）的地区
vulnerable_areas = future_df[future_df['disaster_prediction'] == 1]

# 输出结果
print("\n2025~2035 年预测的暴雨灾害高风险地区：")
print(vulnerable_areas[['longitude', 'latitude', 'year']])

# 保存结果到 CSV 文件
output_file = r'F:\HuaweiCup\Q3\predicted_vulnerable_areas_2025_2035.csv'
vulnerable_areas.to_csv(output_file, index=False)
print(f"\n预测结果已保存到 {output_file}")

# 6. 可视化高风险地区，加入中国地图
# -----------------------------

import cartopy.feature as cfeature  
import cartopy.crs as ccrs  
import geopandas as gpd  
import matplotlib.pyplot as plt  

# 将未来预测的高风险地区转换为 GeoDataFrame  
geometry = gpd.points_from_xy(vulnerable_areas['longitude'], vulnerable_areas['latitude'])  
geo_vulnerable_areas = gpd.GeoDataFrame(vulnerable_areas, geometry=geometry)  

# 设置坐标系为 WGS84  
geo_vulnerable_areas.set_crs(epsg=4326, inplace=True)  

# 绘制地图  
fig, ax = plt.subplots(figsize=(10, 8), subplot_kw={'projection': ccrs.PlateCarree()})  

# 添加地理特征  
ax.add_feature(cfeature.COASTLINE, linewidth=0.5)  
ax.add_feature(cfeature.LAND, edgecolor='black', facecolor='lightgray')  
ax.add_feature(cfeature.OCEAN, facecolor='lightblue')  
ax.add_feature(cfeature.LAKES, alpha=0.5, facecolor='lightblue')  
ax.add_feature(cfeature.RIVERS, linewidth=0.5)  

# 将地图的范围设置为中国区域  
ax.set_extent([73, 135, 18, 54], crs=ccrs.PlateCarree())  

# 绘制高风险地区  
geo_vulnerable_areas.plot(ax=ax, color='red', markersize=5, label='高风险地区', alpha=0.5, transform=ccrs.PlateCarree())  

# 添加标题和标签  
plt.title('2025~2035 年暴雨灾害高风险地区预测', fontsize=20)  
plt.xlabel('经度', fontsize=14)  
plt.ylabel('纬度', fontsize=14)  
plt.legend()  

# 显示经纬度网格线和标签  
gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,  
                  linewidth=0.2, color='gray', alpha=0.5)  
gl.top_labels = False  
gl.right_labels = False  
gl.xlabel_style = {'size': 10, 'color': 'gray'}  
gl.ylabel_style = {'size': 10, 'color': 'gray'}  

# 调整布局  
plt.tight_layout()  

# 显示地图  
plt.show()