Python | 绘制黑底的水平空间分布图

news2025/4/13 13:11:37

写在前面

记录一下之前为了做PPT汇报画的一张图，虽然最后也没怎么用上。为了方面以后再需要，这里把代码和数据整理放到GitHub上。有兴趣的也可以玩玩

需要的数据

风场数据可以从ERA5的官网下载

https://cds.climate.copernicus.eu/datasets/reanalysis-era5-single-levels?tab=overview

GPCP降水为日平均数据，下载地址为：

https://www.ncei.noaa.gov/data/global-precipitation-climatology-project-gpcp-daily/access/

ETOPO1地形数据的下载地址为：

https://www.ngdc.noaa.gov/mgg/global/relief/ETOPO1/data/bedrock/grid_registered/netcdf/

地形数据的介绍为：

http://page.oceanmax.cn/Download_data_ETOP.aspx

以上风场变量选择为850hpa，时间为2004年06月22日，降水时间相同

特别说明

此次绘图的脚本(.ipynb)，包括之前的绘图的脚本，都是测试跑通过的，都会上传到GitHub上。由于此次使用的数据过大，远远超过25Mb（GitHub上单次上传最大文件容量）。所以测试的数据就没有放上去了，但是从提高的数据下载链接，是可以得到所有使用到的数据。

本质上，代码本身比数据更重要一点，因为我使用的数据基本上都是netcdf格式，你看看我写的逻辑基本上换成自己的数据都可以同理实现。

绘制内容

850hpa水平风场
降水
地形

总体上逻辑很简单，就是三个图的叠加，套娃一样。

再此基础上，把背景版设置为黑色，凸显出地形和降水，这里需要注意的一个是，如果想让降水在地形上显示，必须手动设置降水的范围大于10，否则就得设置地形的透明度。

这里绘制的顺序是先绘制地形，后绘制降水，最后是风场和海岸线

而且，很奇怪的点是：单独画风场叠加降水，以及地形图都很快；但是一起画，速度就贼慢，大约需要90s。这也导致我当时后面再去微调的心思都没有了。

以下最终呈现的结果，其实最初的想法是上面是原始场，下面是滤波场，显示滤波后的扰动在原始场中造成的降水的异常。体现高频波动对于降水的调制作用，但是由于调一次图就得等2min，给我等烦了，也就只有下面的1.0版本了，将就着看吧。

python代码

# ======================== 配置文件参数 ========================
DATA_PATH = 'O:/JAS/H-drive-16350906/ERA5_200406/'
TERRAIN_PATH = 'M:/old_desktop/test_pycode/ETOPO1_Bed_g_gmt4.grd'
SAVE_PATH = './mean.png'

# ======================== 数据预处理函数 ========================
def adjust_longitude(ds):
    
    lon_name = 'longitude'  # whatever name is in the data
    ds['longitude_adjusted'] = xr.where(ds[lon_name] < 0, ds[lon_name]%360,\
                          ds[lon_name])
    ds = (
       ds
       .swap_dims({lon_name: 'longitude_adjusted'})
       .sel(**{'longitude_adjusted': sorted(ds.longitude_adjusted)})
       .drop_vars(lon_name))
    ds = ds.rename({'longitude_adjusted': lon_name})
    return ds

def load_data(file, adjust_lon=False):
    """加载并处理数据"""
    ds = xr.open_dataset(DATA_PATH + file).sortby('latitude')
    return adjust_longitude(ds) if adjust_lon else ds

# ======================== 主数据处理 ========================
# 加载数据集
pre_ds = load_data('GPCP-daily_2004.nc')
u_ds = load_data('uwnd.200406.nc', adjust_lon=True)
v_ds = load_data('vwnd.200406.nc', adjust_lon=True)

# 公共选择参数
common_sel = {
    'latitude': slice(-20, 26),
    'longitude': slice(100, 200),
    'time': slice('2004-06-22', '2004-06-28')
}

# 提取数据
pre = pre_ds.precip.sel(**common_sel)[::2]
lon = pre.longitude.data
lat = pre.latitude.data
x, y = np.meshgrid(lon, lat)
u = u_ds.u.sel(level=850, **common_sel)[::8]
v = v_ds.v.sel(level=850, **common_sel)[::8]
ulon=u.longitude.data
ulat=u.latitude.data
x_u,y_u=np.meshgrid(ulon,ulat)
# ======================== 地形数据处理 ========================
terrain = xr.open_dataset(TERRAIN_PATH).sel(
    x=slice(110, 180), y=slice(-10, 25)).z
tlon = terrain.x.data
tlat = terrain.y.data
# ======================== 绘图配置 ========================
# 颜色映射配置
colors = np.vstack([
    plt.cm.terrain(np.linspace(0, 0.15, 250)),
    plt.cm.terrain(np.linspace(0.23, 0.85, 250))
])
terrain_cmap = LinearSegmentedColormap.from_list('hls', colors)

norm = cm.colors.Normalize(vmin=-8000, vmax=8000)
# 绘图参数
PLOT_PARAMS = {
    'pre_clevs': np.linspace(0, 200, 21),
    'terrain_levels': np.arange(-8000, 8000, 500),
    'quiver_step': 2,
    'title_labels': ['2004.06.24', '2004.06.26']
}

# ======================== 绘图函数 ========================

def create_map(ax, u, v, pre,terrain,title):
    extent = [110, 180,-10, 25]
    ax.set_extent(extent, crs=ccrs.PlateCarree())
    cf = ax.contourf(terrain.x.data,terrain.y.data,terrain,levels=np.arange(-8000,8000,500),
                     transform=ccrs.PlateCarree(),cmap=terrain_cmap,norm=norm,
                      extend='both')
    ax.contourf(pre,  [15,20,30,40,50,60,70,80,90,100],  cmap=cmaps.cmocean_ice_r,
                extent=[110, 200, -20, 26], 
                        transform=ccrs.PlateCarree(),)
    ax.set_title(title,color='white',loc='right',fontsize=15)
    ax.set_title('Precipitation',color='white',loc='left',fontsize=15)
    gl = ax.gridlines(draw_labels=True, linewidth=0.5, color='gray', alpha=0.5, linestyle='--')
    gl.xlabels_top = False
    gl.ylabels_right = False
    gl.xlocator = mticker.FixedLocator([105,120, 135, 150, 165, 180, 185])  # 调整经度刻度位置
    gl.xlabel_style = {'color': 'white','size': 15}  # 设置经度标签颜色为白色
    gl.ylabel_style = {'color': 'white','size': 15}  # 设置纬度标签颜色为白色
    ax.coastlines('10m')
    step = 2
    Q = ax.quiver(x_u[::step, ::step], y_u[::step, ::step], u[::step, ::step], v[::step, ::step],
                  transform=ccrs.PlateCarree())

# ======================== 主绘图程序 ========================
fig, axs = plt.subplots(2, 1, figsize=(10, 12), dpi=300,
                       subplot_kw={'projection': ccrs.PlateCarree(180)})
fig.patch.set_facecolor('black')

for i, ax in enumerate(axs):
    print(i)
    create_map(ax, u.data[i+1], v.data[i+1], pre[i+1], terrain,PLOT_PARAMS['title_labels'][i])

plt.tight_layout()
# fig.savefig(SAVE_PATH, dpi=300, bbox_inches='tight')
plt.show()