Python | 基于栅格数据绘制空间分布、频率分布和剖面图

#示例数据 https://pan.baidu.com/s/1kCV24MVFyXtKfnVtds_iKQ?pwd=dec3 提取码: dec3import matplotlibmatplotlib.use('TkAgg')import numpy as npimport rasterioimport matplotlib.pyplot as pltimport cartopy.crs as ccrsimport cartopy.feature as cfeaturefrom cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
# 配置字体plt.rcParams['font.sans-serif'] = ['Times New Roman']plt.rcParams['axes.unicode_minus'] = False
# 栅格数据路径tif_path = r"D:\PML_WUE等\出图\空间分布\T_ET_MEAN.tif"
with rasterio.open(tif_path) as src:    data = src.read(1)  # 读取第一个波段数据    lon_min, lat_min, lon_max, lat_max = src.bounds  # 获取地理范围    transform = src.transform  # 获取仿射变换参数
data = np.flipud(data)
# 去除异常值data[data > 1] = np.nandata[data <= 0] = np.nan
# 检查数据范围print(f"数据范围：最小值={np.nanmin(data)}, 最大值={np.nanmax(data)}")
# 生成经纬度网格lon = np.linspace(lon_min, lon_max, data.shape[1])lat = np.linspace(lat_min, lat_max, data.shape[0])
# 裁剪纬度范围为 -60 到 90lat_mask = (lat >= -60) & (lat <= 90)lat = lat[lat_mask]data = data[lat_mask, :]
lon, lat = np.meshgrid(lon, lat)
# 计算随纬度变化的平均值lat_mean = np.nanmean(data, axis=1)
# 绘图函数# 绘图函数def plot_combined():    fig = plt.figure(figsize=(12, 6))
    # 空间分布图    ax_map = fig.add_axes([0.05, 0.1, 0.6, 0.6], projection=ccrs.PlateCarree())    ax_map.set_title('Annual average of T/ET', loc='left', fontsize=16)
    # 设置显示范围（-60到90纬度）    ax_map.set_extent([lon_min, lon_max, -60, 90], crs=ccrs.PlateCarree())
    # 添加海岸线并加粗边框    ax_map.add_feature(cfeature.COASTLINE.with_scale('110m'), linewidth=0.9)
    # 设置经纬度刻度    ax_map.set_xticks(np.arange(-180, 181, 60), crs=ccrs.PlateCarree())    ax_map.set_yticks(np.arange(-60, 91, 30), crs=ccrs.PlateCarree())    ax_map.xaxis.set_major_formatter(LongitudeFormatter())    ax_map.yaxis.set_major_formatter(LatitudeFormatter())    ax_map.tick_params(labelsize=10)
    # 绘制等值线图    levels = [0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1]    contour = ax_map.contourf(lon, lat, data, levels=levels, cmap='RdYlGn', extend='both', transform=ccrs.PlateCarree())
    # 设置图例    cbar = plt.colorbar(contour, ax=ax_map, orientation='horizontal', pad=0.08, aspect=50, format='%.2f')    cbar.set_label('T/ET', fontsize=18)    cbar.ax.tick_params(labelsize=14)    cbar.set_ticks([0,0.2,0.4,0.6,0.8,1])
    # 剖面图    ax_profile = fig.add_axes([0.65, 0.237, 0.08, 0.463])  # 确保与空间分布图对齐    ax_profile.plot(lat_mean, lat, color='black', linewidth=0.8)  # 颜色改为黑色
    # 设置剖面图样式    ax_profile.set_xlim(0, 1)    ax_profile.set_ylim(-60, 90)    ax_profile.grid(True, linestyle='--', alpha=0.5)
    # 去除纵坐标刻度和注释    ax_profile.yaxis.set_visible(False)
    # 去除标题    ax_profile.set_title('')
    # 横坐标保留刻度但设置字体大小为10号    ax_profile.tick_params(axis='x', labelsize=10)
    # 频数分布直方图    ax_hist = fig.add_axes([0.08, 0.25, 0.1, 0.25])  # 左侧直方图位置，高度为 1/4，宽度为 1/3    flattened_data = data.flatten()    flattened_data = flattened_data[~np.isnan(flattened_data)]  # 去除 NaN
    # 计算频数和相对百分比    counts, bin_edges = np.histogram(flattened_data, bins=levels)    percentages = counts / counts.sum() * 100  # 百分比
    # 计算固定宽度    bar_width = (levels[1] - levels[0])  # 每个柱子的宽度相等
    # 绘制柱子并上色    for i in range(len(levels) - 1):        ax_hist.bar(            levels[i] + bar_width / 2, percentages[i], width=bar_width,            color=plt.cm.RdYlGn((levels[i] + levels[i + 1]) / 2 / max(levels)),            edgecolor='black',            align='center'  # 确保柱子居中对齐        )        # 在柱子顶部垂直标注百分比        ax_hist.text(            levels[i] + bar_width / 2, percentages[i] + 0.5,  # 偏移量调整位置            f"{percentages[i]:.1f}%", ha='center', va='bottom', fontsize=8,            rotation=90  # 文字垂直显示        )
    # 设置直方图样式    ax_hist.set_xlim(levels[0], levels[-1])    ax_hist.axis('off')  # 去掉坐标轴、刻度和边框
    # 保存图片    # 保存图片    plt.savefig(        r'D:\PML_WUE等\出图\空间分布\gzh_t_et.jpg',        dpi=300, bbox_inches='tight', pad_inches=0.01    )    plt.show()
plot_combined()