本节提要：简要谈谈地形剖面图、纬度高度剖面图、时间纬度图的绘制方法。

提要中提到的这几种图形都是在气象上比较常用的，地形剖面主要研究地貌对降雨、气流的影响作用；纬度高度剖面图可以用来分析降雨的某些条件，如湿层深厚、上干下湿、风向风速等；时间纬度图研究某个固定经度上的值随时间的演变(这是和大气环流一般自西向东相匹配的，所以时间经度图比较少见)。

一、地形剖面图

绘制地形剖面图之前，需要了解自己使用的地形文件的格式与属性。我使用的是从气象家园巨佬Masterpiece处白嫖来的地形文件。文件为.nc格式，需要使用Python中的netCDF4或者xarray库包来读取。

首先我们先来读取一下文件，并print出来，看看其属性：

import numpy as npimport matplotlib.pyplot as pltimport matplotlib as mplimport cartopy.crs as ccrsimport cartopy.feature as cfeaturefrom cartopy.io.shapereader import Reader as shpreaderimport xarray as xrplt.rcParams['font.sans-serif']=['SimHei']#显示中文filename=r'E:\aaaa\world_geo.nc'#地形文件储存地址f=xr.open_dataset(filename)#读取文件print(f)#打印其属性

可以看出这个文件主要由x,y,z三个变量组成。

其中x表示经度，将全球东西360经度分为了10800刻度，相当于一个经度被分为30份；y表示纬度，将全球南北180纬度分为了5400份，也是将一个纬度分为30份。那么这个nc文件的精度就是0.0333°×0.0333°，由于是地形文件，显然要比常用的再分析资料精度更高；z表示高度，也就是地形。可以看出，z仅仅与y,x有关，且第一相关量为y而不是x，这与我们习惯不同，在取值时需注意。

因为是二维的数据，那么按照绘制平面填色图的ax.contourf命令是可以直接读取数据绘图的。接下来我们先绘制一个平面的地形图试试成色：

import numpy as npimport matplotlib.pyplot as pltimport matplotlib as mplimport cartopy.crs as ccrsimport cartopy.feature as cfeaturefrom cartopy.io.shapereader import Reader as shpreaderimport xarray as xrplt.rcParams['font.sans-serif']=['SimHei']#显示中文#######################################filename=r'E:\aaaa\world_geo.nc'#地形文件地址proj=ccrs.PlateCarree()#缩写投影extent=[70,140,5,75]#绘图范围def truncate_colormap(cmap, minval=0.0, maxval=1.0, n=128):new_cmap = mpl.colors.LinearSegmentedColormap.from_list('trunc({n},{a:.2f},{b:.2f})'.format(n=cmap.name,a=minval, b=maxval),cmap(np.linspace(minval, maxval, n)))return new_cmap#新的等高图色条，比较符合地理地形图的样子，来源气象家园巨佬############################################f=xr.open_dataset(filename)#读取文件lon=f['x'][:]#将文件中的x变量赋值为经度lat=f['y'][:]#赋值为纬度height=f['z'][:]#将z变量赋值为高度fig=plt.figure(figsize=(10,9),dpi=500)ax=fig.add_subplot(projection=proj)cmap_new = truncate_colormap(plt.cm.terrain, 0.23, 1.0) #截取colormap，要绿色以上的(>=0.23)cmap_new.set_under([198/255,234/255,250/255]) #低于0的填色为海蓝lev=np.arange(0,4000,200)norm3 = mpl.colors.BoundaryNorm(lev, cmap_new.N) #标准化level，映射色标cf=ax.contourf(lon[7500:9600],lat[2850:4950],height[2850:4950,7500:9600],levels=lev,norm=norm3,cmap=cmap_new,extend='both')ax.set_extent(extent)ax.set_title('一个框框地形图',fontsize=20)plt.show()

能画出地图，那就足够了，我们不从地理学和美学的角度再去探讨他。现在你能明白，这个图和你绘制的等温图其实是一个原理的，都是ax.contourf这个命令出来的。

接下来最重要的是绘图的这一句：

ax.contourf(lon[7500:9600],lat[2850:4960],height[2850:4960,7500:9600],levels=lev,norm=norm3,cmap=cmap_new,extend='both')ax.set_extent(extent)############################################ax.contourf(lon,lat,height,levels=lev,norm=norm3,cmap=cmap_new,extend='both')ax.set_extent(extent)

上下两种命令，出图应该是完全一样的，但是最好用上面这种，理由如下：第二种不对导入的数据做取舍，那么程序在绘图时，就会将全球都绘制出来，然后再裁剪边界，这样出图效率大概慢十倍。第一种本质上是将数据扣出一块，只绘制这一块，速度大大提高。

为什么要插这一句嘴，实际上有助于我们在接下来绘制剖面图时理解切片操作。

以经度为例，前面已经讲到将一个经度分为30份，那么我们要画东经70-140的图，那就需要对经度数据切片，原理如下(纬度同理)：

起始：(180+70)×30=7500(在前面属性可知，切片是需加上西经180)

终止：(180+140)×30=9600

接下来就是z的切取了，前面读取属性时我们已经知道，纬度为第一相关量，经度为第二相关量，所以应该先切纬度，后切经度：

height [ 2850:4960 , 7500:9600 ]

接下来，就是本节关键，怎么绘制地形剖面图。在绘制地形填色时，我们使用的是ax.contourf命令，他要求输入横坐标，纵坐标，与横纵坐标有关系的z值。这样z就必须是二维的，以与横纵坐标相关，所以切片时，我们必须使z切取范围与x,y完全一致，否则报错。

但是绘制剖面图，我们还需不需要contourf命令呢？显然是不需要的，我们只想知道沿某个经度(或纬度)的地形变化如何，用ax.plot命令结合fill_between命令即可。而这两个命令，只需要传入一个一维的横坐标，和一维的纵坐标即可。关键就在怎么把z从二维的变为一维的。

这就需要上面的切片方法了，比如我要画东经108.98°这个经线的剖面，那就直接在z取值时，将其x取值设置为固定的8669.

import numpy as npimport matplotlib.pyplot as pltimport matplotlib as mplimport cartopy.crs as ccrsimport cartopy.feature as cfeaturefrom cartopy.io.shapereader import Reader as shpreaderimport matplotlib.ticker as mtickerimport xarray as xrplt.rcParams['font.sans-serif']=['SimHei']plt.rcParams['axes.unicode_minus']=False#######################################filename=r'E:\aaaa\world_geo.nc'f=xr.open_dataset(filename)lat=f['y'][3591:3621]height=f['z'][3591:3621,8669]fig=plt.figure(figsize=(4,1.5),dpi=700)ax=fig.add_axes([0,0,1,1])ax.plot(lat,height,c='k',lw=1)ax.fill_between(lat,height,facecolor='white',hatch='///')#填充阴影ax.set_xlim(29.7,30.6)ax.set_xlabel('北纬(N)',fontsize=7)ax.set_ylim(700,1650)ax.set_ylabel('海拔高度(m)',fontsize=7)ax.tick_params(which='both',labelsize=5)plt.show()

可以进一步print一下这个切片操作：

lat=f['y'][3591:3621]height=f['z'][3591:3621,8669]print(lat)print(height)print(len(lat))print(len(height))

可以看出，两个都变为长度为30的一维数组了。理解这个，就为后面更多维度的切片打下基础。

二、利用NECP的再分析资料绘制纬度高度剖面图

由于笔者的电脑没办法安装pygrib，所以不能直接读取grib2格式的再分析资料，故委托基友杨雨蒙利用NCL转成nc格式的数据了。可能大家目前最需要的是解决在win上读grib2问题，笔者暂时还不能给出满意的解答，气象家园已有xarray配合eccodes和cfgrib或者李开元老师的方法wgrib转换方法，大家可以参考。

仍然和第一节中一样，我们先读取查看一下属性：(由于再分析资料的特性，其属性又臭又长)：

import xarray as xrds = xr.open_dataset(r'E:\aaaa\datanc\he\fnl_20200715_12_00.nc')print(ds)

可以看出，数据由两部分组成。第一部分为经纬度与层次，第二部分为各种物理量。

前面这部分前缀为lv的表示层次，最后两个为经纬度，层次有各种划分方法，学气象的同学大概都知道。

这个文件中有很多基础物理量，举例子常用的：TMP温度 Pres气压 HGT位势高度 RH相对湿度 UGRD纬向风 VGRD经向风 CAPE 对流有效位能。

最前面的TMP表示温度，但是有9种，有的与海平面相关，有的与各层气压相关。

如第一个气温，后面的说明中表示这个只与(lat_0,lon_0)有关；第四个气温与(lv_ISBL0,lat_0,lon_0)有关。这样第一个就是二维的，可以直接绘制等值线填色图，第四个就是三维的，不能直接绘制等值线填色图，而只能在提取了某一层之后，变为二维的，才能绘制等值线填色图，如：

import xarray as xrds = xr.open_dataset(r'E:\aaaa\datanc\he\fnl_20200715_12_00.nc')single_temp=ds['TMP_P0_L1_GLL0'][:][:]#这是二维的many_temp=ds['TMP_P0_L100_GLL0'][:][:][:]#这是三维的single_many_temp=many_temp[0][:][:]#这就变为二维的，只取了一层次

根据第一节中提到的，我们现在要绘制一个某个经度的垂直剖面图，说明我们的横坐标应该是纬度，纵坐标应该是高度，但是在气象上一般不使用高度，而是气压层，如925hPa、850hPa、700hPa、500hPa、200hPa等，而经度就取一个固定值，这样也能变成二维数组。下面通过一个程序讲明，注释直接在程序中：

import numpy as npimport matplotlib.pyplot as pltimport xarray as xrimport cartopy.crs as ccrsimport cartopy.feature as cfimport cartopy.io.shapereader as shpreaderfrom cartopy.mpl.ticker import LongitudeFormatter,LatitudeFormatterfrom cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTERimport matplotlib.ticker as mtickerplt.rcParams['font.sans-serif']=['SimHei']filename=r'E:\aaaa\datanc\he\fnl_20200715_18_00.nc'f=xr.open_dataset(filename)lon=f['lon_0'][:]#读取经度，一维的lat=f['lat_0'][:]#读取纬度，一维的RH_P0_L100_GLL0=f['RH_P0_L100_GLL0'][:][:][:]#读取相对湿度，三维的UGRD_P0_L100_GLL0=f['UGRD_P0_L100_GLL0'][:][:][:]#读取纬向风，三维的VGRD_P0_L100_GLL0=f['VGRD_P0_L100_GLL0'][:][:][:]#读取经向风，三维的pres= f['lv_ISBL5'][:]*0.01#读取气压层数，一维的fig=plt.figure(figsize=(5,3),dpi=700)#添加画布ax=fig.add_axes([0,0,1,1])#添加子图ax.invert_yaxis()#反转纵轴，使1000hPa作为起点ax.set_yticks([1000,925,850,700,500,300])#突出我们感兴趣的气压层ax.set_xticks(np.arange(28,36,1))ax.set_xticklabels([r'28$^\degree$N',r'29$^\degree$N',r'30$^\degree$N',r'31$^\degree$N',r'32$^\degree$N',r'33$^\degree$N',r'34$^\degree$N',r'35$^\degree$N'])#转换为纬度格式ax.set(ylim=(1000,300))#设置气压层绘图范围，此处我们只显示到300hPaax.set_xlabel('纬度',fontsize=7)ax.set_ylabel('层次(hPa)',fontsize=7)ax.tick_params(axis='both',which='both',labelsize=7)##################################################################################ac=ax.contourf(lat[55:63],pres[:],RH_P0_L100_GLL0[:,55:63,109],cmap='Greens',levels=np.arange(60,101,10),extend='both',alpha=0.75)ax.barbs(lat[55:63],pres[:],UGRD_P0_L100_GLL0[:,55:63,109],VGRD_P0_L100_GLL0[:,55:63,109],barb_increments={'half':2,'full':4,'flag':20},length=5)cb=fig.colorbar(ac,extend='both',shrink=1,label='相对湿度(%)',pad=0.01)cb.ax.tick_params(axis='both',which='both',length=1,labelsize=6)ax.set_title('2020年7月15日20时相对湿度与风场剖面图',fontsize=15)

最关键的仍然是这一句：

ax.contourf(lat[55:63],pres[:],RH_P0_L100_GLL0[:,55:63,109],cmap='Greens',levels=np.arange(60,101,10),extend='both',alpha=0.75)

我们使RH_P0_L100_GLL0取为[ : , 55:63 , 109 ],这表示什么呢？在前面读取阶段，我们知道了RH_P0_L100_GLL0这个物理量与三个参量有关，按顺序分别是

而在文件属性界面，我们可以知道，lv_ISBL5表示气压层，lat_0表示纬度，lon_0表示经度。

所以[ : ]表示取全部的气压层次高度，[ 55:63 ]表示取第55至63个纬度的值(不是北纬55-63，这

个是切片序号，不是其

存

放纬度值，

具体纬度值是多少需要你去算，我选的纬度是28-35)，[ 109 ]表示取第109个经度值(也是切片序号，但是恰恰其存放值为109°E)，经过切片后，经度因为只取了一个值，所以被降维，

由于经度被降维了，这个相对湿度物理量只剩纬度，气压层次两维了，而两维数据就可以直接绘图了。

ax.barbs(lat[55:63],pres[:],UGRD_P0_L100_GLL0[:,55:63,109],VGRD_P0_L100_GLL0[:,55:63,109],barb_increments={'half':2,'full':4,'flag':20},length=5)

风场这个也是同样的道理，经向风与纬向风按同样方法切片取值。

接下来是一个五维数据的剖面图绘制，当然其实没什么意义了，但是可以帮助各位读者更好理解切片降维方法

可以看出，这个文件里的z与五个参数有关，所以可以称为五维变量，听起来很可怕实际上很简单。下面是绘制其剖面图的方法：

import xarray as xrimport numpy as npimport matplotlib.pyplot as pltplt.rcParams['font.sans-serif']=['SimHei']filename=r'E:\aaaa\z1.nc'f=xr.open_dataset(filename)lat=f['lat'][:]level=f['level'][:]z=f['z'][:][:][:][:][:]fig=plt.figure(figsize=(5,4),dpi=500)ax=fig.add_axes([0,0,1,1])ax.invert_yaxis()ca=ax.contourf(lat[90:181],level[:],z[1,1,:,90:181,100],levels=np.arange(0,320000,10000))fig.colorbar(ca)ax.set(xlabel='北纬',ylabel='气压层(百帕)',title='多维数据的剖面图')plt.show()

在z[ 1 , 1 , : , 90:181 , 100 ]里，按顺序分别表示years取第一个切片值；time取第一个切片值；层次level从上至下全部取完；纬度取第90到181个切片值；经度取第100个切片值。所以，years、time、经度这三个维度遭到降维打击，那么z变为仅与lat与level相关的二维数据，就可以画等值线填色图了。

现在各位应该知道绘制剖面图技巧了，无论有多少维度，只保留感兴趣的两维，其他维度都做降维处理，处理完的数据变为二维，二维数据直接传入ax.contourf()中画图。

三、时间纬度图(或时间高度图)

时间纬度图已经有摸鱼咯大佬给出了一种方法，大家可以左转去参考。

时间高度图就是像下面这种：

我还没有画过，但是猜测应当是这个数据为四维数据，将经度、纬度做降维处理，从图上可以看出，这张图代表(30.28°E，108.93°N)这一个点的整层数据随时间的变化。

可能是这样取的[ 1414:1720 , all level , 30.28 , 108.93 ]。等后面我画出来了再给大家补充吧。

实验文件提取(包括一个地形文件、一个气温再分析资料文件、一个转换为nc格式的grib2文件)

链接：https://pan.baidu.com/s/1ZK47zL2XJjKn0e2ubDYNuw

提取码：oci0

感谢阅读！！！

欢迎关注云台书使公众号获取更多资讯