HF.086 I Python组合图绘制

GEE遥感训练营 · 公众号 · · 2024-08-05 12:30

正文

请到「今天看啥」查看全文

在科研领域，数据可视化是传达复杂信息的关键手段之一。组合图作为一种高效的数据可视化方法，通过将多个子图和主图结合在一张图表中，能够展现多维度的信息，从而帮助研究人员更直观地理解数据趋势和关系。

在 HF.084 中，我们已经重点讲述了组合图的意义及分类，本文以全球水文站点日尺度径流数据（GRDC）为例，筛选出时间序列大于60年的站点，对每个站点的数据进行Mann-Kendall、Levene及KPSS检验，分析全球各水文站点年最大径流的一致性，结果通过组合图进行展示。

绘图逻辑

首先，使用 热图(heatmap) 作为二维图例，以Levene及KPSS检验得到的p值作为两个维度的值，对筛选后的水文站点进行一致性分析。

然后，通过绘制全球地图以及组合 放大图 ，显示局部细节区域，包括美国、欧洲、南美洲、非洲南部和澳大利亚。

接着，根据Mann-Kendall检验的统计值，对非一致性的水文站点进行趋势判断，并替换放大子图中的显示值。

最后，使用饼图统计三种检验中各类水文站点的占比，简要展示各个检验的结果。

实现过程

首先，加载并验证数据，确保数据包含必要的列，如经纬度、Levene Test p-value、KPSS Test p-value、Mann-Kendall Test Statistic等。同时，处理数据中的极端值，以避免极端值对可视化效果的影响。

import pandas as pdimport geopandas as gpdimport matplotlib.pyplot as pltimport numpy as npimport matplotlib.colors as mcolorsfrom sklearn.preprocessing import StandardScaler

def load_station_data(file_path):    return pd.read_excel(file_path)

def validate_columns(df, required_columns):    missing_columns = [col for col in required_columns if col not in df.columns]    if missing_columns:        raise ValueError(f"The data file must contain the following columns: {', '.join(missing_columns)}")

def filter_extreme_values(df, col, lower_quantile=0.05, upper_quantile=0.95):    lower_bound = df[col].quantile(lower_quantile)    upper_bound = df[col].quantile(upper_quantile)    return df[(df[col] >= lower_bound) & (df[col] <= upper_bound)]

将经纬度数据转换为地理数据框，以便进行空间绘图。

def create_geodataframe(df, lon_col='Longitude', lat_col='Latitude'):    gdf = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df[lon_col], df[lat_col]))    gdf.set_crs(epsg=4326, inplace=True)    return gdf

创建二维颜色映射，根据Levene Test p-value和KPSS Test p-value的值创建二维颜色映射，以便在主图中展示站点的稳定性分析。

def create_2d_color_map(df, stat_col, pval_col):    stat_min, stat_max = df[stat_col].min(), df[stat_col].max()    pval_min, pval_max = df[pval_col].min(), df[pval_col].max()

    stat_diff_max = abs(stat_max - 0.05)    stat_diff_min = abs(stat_min - 0.05)    pval_diff_max = abs(pval_max - 0.05)    pval_diff_min = abs(pval_min - 0.05)

    stat_extend = max(stat_diff_max, stat_diff_min)    pval_extend = max(pval_diff_max, pval_diff_min)

    stat_bins = np.array([0.05 - stat_extend, 0.05, 0.05 + stat_extend])    pval_bins = np.array([0.05 - pval_extend, 0.05, 0.05 + pval_extend])

    color_matrix = np.zeros((3, 3, 4))    cmap = plt.get_cmap('viridis')    for i in range(3):        for j in range(3):            color_matrix[j, i] = cmap((i * 3 + j) / 9)

    return stat_bins, pval_bins, color_matrix

绘制主图和子图，主图展示全球水文站点的稳定性分析，子图分别放大展示美国、欧洲、南美洲、非洲南部和澳大利亚的局部细节，并用颜色表示Mann-Kendall Test Statistic的值。

def




    
 plot_stations_on_world_map(gdf, shapefile_path, stat_bins, pval_bins, color_matrix, stat_col, pval_col):    world = gpd.read_file(shapefile_path)

    fig, ax_main = plt.subplots(figsize=(36, 28))    world.plot(ax=ax_main, color='lightgray', edgecolor='black')

    for idx, row in gdf.iterrows():        stat_bin = np.digitize([row[stat_col]], stat_bins)[0] - 1        pval_bin = np.digitize([row[pval_col]], pval_bins)[0] - 1        color = color_matrix[pval_bin, stat_bin]        ax_main.scatter(row['geometry'].x, row['geometry'].y, color=color, s=50)

    ax_main.grid(True, which='both', linestyle='--', linewidth=1)    ax_main.set_xticks(np.arange(-180, 181, 30))    ax_main.set_yticks(np.arange(-90, 91, 30))    ax_main.set_xticklabels([f'{x}°' for x in np.arange(-180, 181, 30)])    ax_main.set_yticklabels([f'{y}°' for y in np.arange(-90, 91, 30)])    ax_main.set_xlabel('Longitude')    ax_main.set_ylabel('Latitude')

    insets = {        'USA': [-130, -60, 20, 55, [0.15, 0.71, 0.30, 0.30]],        'Europe': [-10, 30, 38, 65, [0.58, 0.74, 0.26, 0.26]],        'South America': [-80, -30, -55, -5, [0.1, 0.02, 0.28, 0.21]],        'Southern Africa': [20, 33, -33, -23, [0.40, 0.02, 0.22, 0.22]],        'Australia': [130, 155, -45, -20, [0.65, 0.02, 0.20, 0.24]]    }

    filtered_gdf = filter_extreme_values(gdf, 'Mann-Kendall Test Statistic')    scaler = StandardScaler()    filtered_gdf['Mann-Kendall Test Statistic'] = scaler.fit_transform(filtered_gdf[['Mann-Kendall Test Statistic']])    max_abs_value = max(abs(filtered_gdf['Mann-Kendall Test Statistic'].min()), abs(filtered_gdf['Mann-Kendall Test Statistic'].max()))    norm = plt.Normalize(vmin=-max_abs_value, vmax=max_abs_value)    cmap = plt.get_cmap('coolwarm')

    for region, (xmin, xmax, ymin, ymax, bbox) in insets.items():        ax_inset = fig.add_axes(bbox)        world.plot(ax=ax_inset, color='lightgray', edgecolor='black')        for idx, row in filtered_gdf.iterrows():            if xmin <= row['geometry'].x <= xmax and ymin <= row['geometry'].y <= ymax:                color = cmap(norm(row['Mann-Kendall Test Statistic']))                ax_inset.scatter(row['geometry'].x, row['geometry'].y, color=color, s=50)

        ax_inset.set_xlim(xmin, xmax)        ax_inset.set_ylim(ymin, ymax)        ax_inset.grid(True, which='both', linestyle='--', linewidth=0.5)        ax_inset.set_xticks([])        ax_inset.set_yticks([])

    cax_main = fig.add_axes([0.18, 0.40, 0.06, 0.08])    cax_main.imshow(color_matrix, origin='lower', aspect='auto', extent=[stat_bins.min(), stat_bins.max(), pval_bins.min(), pval_bins.max()])    cax_main.set_xticks([stat_bins[0], stat_bins[1], stat_bins[2]])    cax_main.set_yticks([pval_bins[0], pval_bins[1], pval_bins[2]])    cax_main.set_xticklabels([f'{stat_bins[0]:.2f}', '0.05', f'{stat_bins[2]:.2f}'])    cax_main.set_yticklabels([f'{pval_bins[0]:.2f}', '0.05', f'{pval_bins[2]:.2f}'])    cax_main.set_xlabel('Levene Test')    cax_main.set_ylabel('KPSS Test')    cax_main.set_title('Stationarity Analysis')
    cax_inset = fig.add_axes([0.18, 0.35, 0.06, 0.02])    sm_inset = plt.cm.ScalarMappable(cmap=cmap, norm=norm)    sm_inset.set_array([])    fig.colorbar(sm_inset, cax=cax_inset, orientation='horizontal', label='Mann-Kendall Test Statistic')
    # 画饼图    ax_pie1 = fig.add_axes([0.12, 0.25, 0.12, 0.12])    ax_pie2 = fig.add_axes([0.35, 0.25, 0.12, 0.12])    ax_pie3 = fig.add_axes([0.58, 0.25, 0.12, 0.12])
    def plot_pie_chart(ax, sizes, labels, title, colors):        ax.pie(sizes, labels=labels, colors=colors, autopct='%1.1f%%', startangle=90)        ax.axis('equal')        ax.set_title(title)
    mk_stat_counts = [sum(gdf['Mann-Kendall Test Statistic'] > 0), sum(gdf['Mann-Kendall Test Statistic'] <= 0)]    plot_pie_chart(ax_pie1, mk_stat_counts, ['Increase', 'Decrease'], 'Mann-Kendall Test Statistic', ['#66c2a5', '#fc8d62'])
    levene_pval_counts = [sum(gdf['Levene Test p-value'] > 0.05), sum(gdf['Levene Test p-value'] <= 0.05)]    plot_pie_chart(ax_pie2, levene_pval_counts, ['Not Significant', 'Significant'], 'Levene Test p-value', ['#8da0cb', '#e78ac3'])
    kpss_pval_counts = [sum(gdf['KPSS Test p-value'] > 0.05), sum(gdf['KPSS Test p-value'] <= 0.05)]    plot_pie_chart(ax_pie3, kpss_pval_counts, ['Not Significant', 'Significant'], 'KPSS Test p-value', ['#a6d854', '#ffd92f'])
    plt.show()

结果呈现

全球水文站点一致性分析图（ 二维热图图例 ）

全球水文站点一致性分析图（ 局部放大图 ）

全球水文站点一致性分析图（ 双图例图 ）

全球水文站点一致性分析图（饼图）

绘图总结

组合图是科研数据可视化中非常重要的一种工具，通过多个子图和主图的结合，能够在一张图表中展示多维度的信息。本文通过实际案例详细介绍了组合图的绘制过程，包括数据加载与验证、数据预处理、创建GeoDataFrame、生成颜色映射以及绘制地图等步骤。最终生成的组合图不仅能够展示全球水文站点的分布，还能直观地反映各站点的变化趋势以及整体统计，为科研人员提供了丰富的信息支持。

一图胜千言！水文图绘改版后致力于分享水文相关的精美图表，为读者提供作图思路和经验，帮助大家制作更漂亮丰富的图表。同时欢迎留言咨询绘图难点，我们会针对性地分享相关绘制经验。另外也期待读者踊跃来稿，分享更好的构图思维和技巧，稿件可发送至邮箱[email protected], 或者联系微信17339888901投稿。

编辑: 徐淑高马孟良｜校稿：Hydro90编委团