Python绘制任意两点剖面图

 作者ID 

@付亚男 

metpy中的cross_section提供了非常便捷的绘制剖面图的方法，具体可见网：

https://unidata./MetPy/latest/examples/cross_section.html#sphx-glr-examples-cross-section-py

如果你只需要简单地画个剖面，那么代码如下：

# %%import xarray as xrimport proplot as ppltfrom metpy.interpolate import cross_section# %%era5fnm = r'F:/era5/era5_2021071900_2021072023/era5_2021071900UTC.nc'era5f = xr.open_dataset(era5fnm)era5f = era5f.metpy.parse_cf().squeeze()# %%# create cross sectionstart = (12.0, 96.0)end = (27.0, 120.0)cross = cross_section(era5f, start, end)# %%fig = pplt.figure( refwidth = 5.0, )ax = fig.subplot( )# plot cross sectionax.contourf( cross['g0_lat_2'], cross['lv_ISBL1'][::-1], cross['V_GDS0_ISBL'][::-1, :], )ax.quiver( cross['g0_lat_2'][::4], cross['lv_ISBL1'][::-1], cross['V_GDS0_ISBL'][::-1, ::4], cross['W_GDS0_ISBL'][::-1, ::4]*10.,    )

如果你想要为剖面添加地形并且显示剖面的平面位置和一些场变量，那么代码如下：

# %%import xarray as xrimport proplot as ppltfrom metpy.interpolate import cross_sectionimport cartopy.crs as ccrsimport cartopy.feature as cfeatfrom cartopy.io.shapereader import Reader# %%era5fnm = r'F:/era5/era5_2021071900_2021072023/era5_2021071900UTC.nc'era5f = xr.open_dataset(era5fnm)era5f = era5f.metpy.parse_cf().squeeze()terfnm = r'F:/terrain/dixingdata.nc'terf = xr.open_dataset(terfnm)terf = terf.metpy.parse_cf().squeeze()c = terf['Y'].loc[30:10]# %%# create cross sectionstart = (12.0, 96.0)end = (27.0, 120.0)cross = cross_section(era5f, start, end)tercross = cross_section(terf, start, end)# %%fig = pplt.figure(    refwidth = 5.0,    )ax = fig.subplot(    )# plot cross sectionax.contourf(    cross['g0_lat_2'],    cross['lv_ISBL1'][::-1],    cross['V_GDS0_ISBL'][::-1, :],    )ax.quiver(    cross['g0_lat_2'][::4],    cross['lv_ISBL1'][::-1],    cross['V_GDS0_ISBL'][::-1, ::4],    cross['W_GDS0_ISBL'][::-1, ::4]*10.,    )# plot terrainoy = ax.alty()oy.area(    cross['g0_lat_2'],    tercross['elev'],    c = 'grey',    )oy.format(    ylim = (0, 10000),    yticks = 'none',    ylabel = '',    )# plot inset paneldata_crs = era5f['V_GDS0_ISBL'].metpy.cartopy_crsix = fig.add_axes(              # do note that using ax.inset() of proplot will cause a geo stale error    [0.105, 0.685, 0.4, 0.3],   # thus using fig.add_axes alternatively    projection = data_crs,    )ix.format(    lonlim = (95, 125),    latlim = (10, 30),    )ix.contour(    era5f['g0_lon_3'].loc[95:125],    era5f['g0_lat_2'].loc[10:30],    era5f['Z_GDS0_ISBL'].loc[850, 10:30, 95:125],    c = 'blue',    )ix.line(    cross['g0_lon_3'],    cross['g0_lat_2'],    c = 'red',    )# add geographic informationix.coastlines()provinces = cfeat.ShapelyFeature(    Reader(r'F:/ngcc/bou2_4m/bou2_4l.shp').geometries(),    ccrs.PlateCarree(),     edgecolor='black',    facecolor='none',    )river = cfeat.ShapelyFeature(    Reader(r'F:/Chinamap-master/cnmap/rivers.shp').geometries(),    ccrs.PlateCarree(),     edgecolor='lightblue',    facecolor='none',    )ix.add_feature(provinces, linewidth=0.5, zorder=2)ix.add_feature(river, linewidth=1.0, zorder=2)