Python 空间解释性机器学习实战:基于 RF 与 GeoShapley 的全流程代码解析

# =========================================================================================# ====================================== 1. 环境设置 =======================================# =========================================================================================import osimport numpy as npimport pandas as pdimport matplotlibimport matplotlib.pyplot as pltimport matplotlib.colors as mcolors

# =========================================================================================# ======================================2.颜色库=======================================# =========================================================================================COLOR_SCHEMES = {    1: ['#086c76', '#3dbbc7', '#ffd5cf', '#f19c83', '#cc2a2e'],}scheme_id = 5#设置要使用的配色方案colors = COLOR_SCHEMES[scheme_id] #提取配色方案

注意：所有内容均为我的个人理解，可能存在错误或不足，使用时还要阅读原文

# =========================================================================================# ====================================== 3.GeoShapley值空间绘图函数===========================# =========================================================================================def plot_spatial_shap_china(df_plot, plot_values, china_boundary, feature_name, save_dir, plot_type):    #创建画布    fig, ax = plt.subplots(figsize=(12, 8), dpi=120)    # 将shp绘制在图上    china_boundary.plot(ax=ax,  #坐标轴                        facecolor='none',  #填充色                        edgecolor='black',  #边界线颜色                        linewidth=1.0,  #边界线宽度                        zorder=1)  #层    idx_max = np.argmax(plot_values) #数值最大位置    idx_min = np.argmin(plot_values) #数值最小位置    minx, miny, maxx, maxy = china_boundary.total_bounds #四至范围，最小经/纬度，最大经/纬度    buffer = 2.0 #边缘的缓冲距离    #标题    ax.set_title(f'{feature_name} ({plot_type})', #文本                 fontsize=22, #大小                 pad=15, #标题与图的间距                 fontweight='bold') #加粗    ax.grid(True, linestyle='-', color='#e0e0e0', alpha=0.7, zorder=0) #背景网格线    #设置边框线    for spine in ax.spines.values():        spine.set_edgecolor('black')        spine.set_linewidth(2.0)    ax.tick_params(axis='both', which='major', length=8, width=2.0, labelsize=14) #主刻度样式    color_labels = [f"{bounds[i]:.4f} ~ {bounds[i + 1]:.4f}" for i in range(5)]    #颜色图例矩形    color_legend_elements = [mpatches.Patch(color=colors[i], label=color_labels[i]) for i in range(5)]    png_path = os.path.join(save_dir, f'spatial_{plot_type}_{feature_name}_{scheme_id}.png')    pdf_path = os.path.join(save_dir, f'spatial_{plot_type}_{feature_name}_{scheme_id}.pdf')    plt.savefig(png_path, dpi=300, bbox_inches='tight')    plt.savefig(pdf_path, bbox_inches='tight')    plt.close() #关闭

# =========================================================================================# ====================================== 4.SVC值空间绘图函数===========================# =========================================================================================def plot_spatial_svc_china(df_plot, plot_values, china_boundary, feature_name, save_dir, plot_type):    #创建画布    fig, ax = plt.subplots(figsize=(12, 8), dpi=120)    # 将shp绘制在图上    china_boundary.plot(ax=ax,  #坐标轴                        facecolor='none',  #填充色                        edgecolor='black',  #边界线颜色                        linewidth=1.0,  #边界线宽度                        zorder=1)  #层    minx, miny, maxx, maxy = china_boundary.total_bounds #四至    ax.grid(True,linestyle='-',color='#e0e0e0',alpha=0.7,zorder=0)    #边框线    for spine in ax.spines.values():        spine.set_edgecolor('black')        spine.set_linewidth(2.0)    #设置刻度样式    ax.tick_params(axis='both', #双轴                   which='major', #主刻度                   length=8, #刻度线长度                   width=2.0, #刻度线粗细                   labelsize=14) #标签字号    png_path = os.path.join(save_dir, f'spatial_{plot_type}_{feature_name}_{scheme_id}.png')    pdf_path = os.path.join(save_dir, f'spatial_{plot_type}_{feature_name}_{scheme_id}.pdf')    plt.savefig(png_path, dpi=300, bbox_inches='tight')    plt.savefig(pdf_path, bbox_inches='tight')    plt.close() #关闭

# =========================================================================================# ====================================== 5.GeoShapley蜂巢图绘图函数===========================# =========================================================================================def plot_custom_summary(geo_results, save_dir):    ax.tick_params(axis='both',which='major',labelsize=14,width=1.8,length=6)    #设置轴刻度标注大小    for label in (ax.get_xticklabels() + ax.get_yticklabels()):        label.set_fontweight('bold')    ax.xaxis.label.set_size(16) #X轴标题大小    ax.xaxis.label.set_weight('bold') #X轴标题加粗    ax.yaxis.label.set_size(16) #Y轴标题大小    ax.yaxis.label.set_weight('bold') #Y轴标题加粗    #保存    plt.savefig(os.path.join(save_dir, f"summary_plot_{scheme_id}.png"), dpi=300, bbox_inches='tight')    plt.savefig(os.path.join(save_dir, f"summary_plot_{scheme_id}.pdf"),bbox_inches='tight')    plt.close() #关闭

# =========================================================================================# ====================================== 6.部分依赖图(PDP)绘图函数===========================# =========================================================================================def plot_custom_pdp(geo_results, save_dir):    #调用内置函数绘图    geo_results.partial_dependence_plots(gam_curve=True, #平滑GAM曲线拟合                                         edgecolors='white', #边颜色                                         lw=0.8) #边线宽    #画布    fig = plt.gcf()    fig.set_size_inches(15, 10) #调整子图尺寸    # 手动控制所有子图离边框左边缘的留白比例    fig.subplots_adjust(left=0.08, #左边缘                        right=0.88, #右边缘                        top=0.92, #上边缘                        bottom=0.08, #下边缘                        wspace=0.35, #水平间隙                        hspace=0.45) #垂直间隙        ax.yaxis.label.set_weight('bold') #Y标题加粗        path_col = ax.collections[0] #取出第一组点的渲染层        offsets = path_col.get_offsets() #点的XY坐标        target_vals = offsets[:, 1] #提取每个点的Y轴坐标数值        norm = mcolors.Normalize(vmin=target_vals.min(), vmax=target_vals.max()) #映射[min, max]至[0,1]        face_colors = cmap(norm(target_vals)) #使用规范化后的数值从之前定义的色带中提取出对应的一组颜色        path_col.set_facecolors(face_colors) #覆盖默认颜色    for label in cbar.ax.get_yticklabels():        label.set_fontweight('bold')    cbar.outline.set_linewidth(1.5)    #保存    plt.savefig(os.path.join(save_dir, f"partial_dependence_plots_{scheme_id}.png"), dpi=300, bbox_inches='tight')    plt.savefig(os.path.join(save_dir, f"partial_dependence_plots_{scheme_id}.pdf"),bbox_inches='tight')    plt.close() #关闭

# =========================================================================================# ====================================== 7.绘制特征全局贡献柱状图绘图函数===========================# =========================================================================================def plot_custom_contribution(geo_results, save_dir):    for i, patch in enumerate(ax.patches):        if i < half: #如果当前条柱属于前半部分，代表非空间解释部分            patch.set_facecolor(colors[1]) #颜色        else: # 如果属于后半部分，空间解释相关            patch.set_facecolor(colors[-2])    #创建图例    non_geo_patch = mpatches.Patch(color=colors[1], label='Non-Geo')    geo_patch = mpatches.Patch(color=colors[-2], label='Geo')    #添加图例    ax.legend(handles=[non_geo_patch, geo_patch], loc='lower right', fontsize=14, frameon=True)    ax.xaxis.label.set_size(16) #X轴标题大小    ax.xaxis.label.set_weight('bold') #x轴标题加粗    ax.yaxis.label.set_size(16) #标题字体大小    ax.yaxis.label.set_weight('bold') #标题加粗    plt.tight_layout() #自动调整    #保存    plt.savefig(os.path.join(save_dir, f"contribution_bar_plot_{scheme_id}.png"), dpi=300, bbox_inches='tight')    plt.savefig(os.path.join(save_dir, f"contribution_bar_plot_{scheme_id}.pdf"),bbox_inches='tight')    plt.close() #关闭

# =========================================================================================# ====================================== 8.GeoShapley值空间绘图函数===========================# =========================================================================================if __name__ == "__main__":    BASE_OUTPUT_DIR = r"haley的地理空间可解释性分析与可视化绘图" #保存路径    GEOSHAPLEY_DIR = os.path.join(BASE_OUTPUT_DIR, "GeoShapley_Plots") #GeoShapley图路径    SVC_DIR = os.path.join(BASE_OUTPUT_DIR, "SVC_Plots") #SVC图保存路径    OTHER_DIR = os.path.join(BASE_OUTPUT_DIR, "Other_Analysis") #非空见图保存路径    #不存在就创建    os.makedirs(GEOSHAPLEY_DIR, exist_ok=True)    os.makedirs(SVC_DIR, exist_ok=True)    os.makedirs(OTHER_DIR, exist_ok=True)    LOCAL_SHP_PATH = r"2023年省级.shp" #shp文件路径    LOCAL_DATA_PATH = os.path.join(BASE_OUTPUT_DIR, "synthetic_china_data.xlsx") #模型数据路径    china_boundary = gpd.read_file(LOCAL_SHP_PATH) # 读取shp数据    df = pd.read_excel(LOCAL_DATA_PATH) # 读取模型数据    features = ['LUE', 'GDPE', 'POPE', 'COHESION', 'FRAC_CV', 'IJI', 'ROAD', 'SLOPE', 'CR'] #特征    geo_features = features + ['LON', 'LAT'] #补充经纬度    X = df[geo_features] #提取特征数据    y = df['Target_ER'] #提取目标数据    #划分数据集    X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)    #设置网格参数    param_grid = {        'n_estimators': [50, 100,200,500],        'max_depth': [6, 7,8,9],        'min_samples_split': [2, 5]    }    #实例化网格搜索优化器对象    grid_search = GridSearchCV(estimator=RandomForestRegressor(random_state=42),param_grid=param_grid,cv=3,scoring='r2',n_jobs=-1,verbose=2)

scheme_id = 1#设置要使用的配色方案

2.设置保存的路径，执行部分：

BASE_OUTPUT_DIR = r"基于RF与Geoshaley的地理空间可解释性分析与可视化绘图" #保存路径

3.设置shp文件的保存地址，执行部分：

LOCAL_SHP_PATH = r"2023年省级.shp" #shp文件路径

LOCAL_DATA_PATH = os.path.join(BASE_OUTPUT_DIR, "synthetic_china_data.xlsx") #模型数据路径

5.设置常规特征，执行部分：

features = ['LUE', 'GDPE', 'POPE', 'COHESION', 'FRAC_CV', 'IJI', 'ROAD', 'SLOPE', 'CR'] #特征

geo_features = features + ['LON', 'LAT'] #补充经纬度

7.设置目标变量，执行部分：

y = df['Target_ER'] #提取目标数据

8.设置超参数，执行部分：

param_grid = {    'n_estimators': [50, 100,200,500],    'max_depth': [6, 7,8,9],    'min_samples_split': [2, 5]}

9.设置是否进行批量绘图，执行部分：

plot_all = False