地形剖面图、纬度高度剖面图如何绘制 |
您所在的位置：网站首页 › 绘制地形剖面图的方法有哪些 › 地形剖面图、纬度高度剖面图如何绘制 |
地形剖面图、纬度高度剖面图如何绘制
|
在气象中，我们常常需要用到剖面图。地形剖面主要研究地貌对降雨、气流的影响作用；纬度高度剖面图主要用来分析降雨的某些条件，如湿层深厚、上干下湿、风向风速等。 一、地形剖面图 绘制地形剖面图之前，需要了解自己使用的地形文件的格式与属性。文件为.nc格式，需要使用Python中的netCDF4或者xarray库包来读取。 首先我们先来读取一下文件，并print出来，看看其属性： import numpy as np import matplotlib.pyplot as plt import matplotlib as mpl import cartopy.crs as ccrs import cartopy.feature as cfeature from cartopy.io.shapereader import Reader as shpreader import xarray as xr plt.rcParams['font.sans-serif']=['SimHei']#显示中文 filename=r'C:\Users\86132\world_geo.nc'#地形文件储存地址 f=xr.open_dataset(filename)#读取文件 print(f)#打印其属性
出图如下：
上下两种命令，出图应该完全一样（几句前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 np import matplotlib.pyplot as plt import matplotlib as mpl import cartopy.crs as ccrs import cartopy.feature as cfeature from cartopy.io.shapereader import Reader as shpreader import matplotlib.ticker as mticker import xarray as xr plt.rcParams['font.sans-serif']=['SimHei'] plt.rcParams['axes.unicode_minus']=False ####################################### filename=r'C:\Users\86132\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)#a为图形宽，b为图形长，dpi为设置图形每英寸的点数 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()
二、利用NECP的再分析资料绘制纬度高度剖面图 先读取查看一下属性： import xarray as xr ds = xr.open_dataset(r'C:\Users\86132\fnl_20200717_00_00.nc') print(ds)
根据之前提到的，我们现在要绘制一个某个经度的垂直剖面图，说明我们的横坐标应该是纬度，纵坐标应该是高度，但是在气象上一般不使用高度，而是气压层，如925hPa、850hPa、700hPa、500hPa、200hPa等，而经度就取一个固定值，这样也能变成二维数组。下面通过一个程序讲明，注释直接在程序中： import numpy as np import matplotlib.pyplot as plt import xarray as xr import cartopy.crs as ccrs import cartopy.feature as cf import cartopy.io.shapereader as shpreader from cartopy.mpl.ticker import LongitudeFormatter,LatitudeFormatter from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER import matplotlib.ticker as mticker plt.rcParams['font.sans-serif']=['SimHei'] filename=r'C:\Users\86132\fnl_20200717_00_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,4),dpi=200)#添加画布 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))#设置气压层绘图范围，此处我们只显示到300hPa ax.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)
我们使RH_P0_L100_GLL0取为[ : , 55:63 , 109 ],这表示什么呢？在前面读取阶段，我们知道了RH_P0_L100_GLL0这个物理量与三个参量有关，按顺序分别是：
风场这个也是同样的道理，经向风与纬向风按同样方法切片取值。 接下来是一个五维数据的剖面图绘制，可以帮助各位读者更好理解切片降维方法。
现在各位应该知道绘制剖面图技巧了，无论有多少维度，只保留感兴趣的两维，其他维度都做降维处理，处理完的数据变为二维，二维数据直接传入ax.contourf（）中画图。 本文摘自：https://my.oschina.net/u/4579644/blog/4518065 |
可以看出这个文件主要由x,y,z三个变量组成。 其中x表示经度，将全球东西360经度分为了10800刻度，相当于一个经度被分为30份；y表示纬度，将全球南北180纬度分为了5400份，也是将一个纬度分为30份。那么这个nc文件的精度就是0.0333°×0.0333°；z表示高度，也就是地形。可以看出，z仅仅与y,x有关，且第一相关量为y而不是。 因为是二维的数据，那么按照绘制平面填色图的ax.contourf命令是可以直接读取数据绘图的。 接下来我们先绘制一个平面的地形图：
其中，最重要的是绘图的这一句：
进一步print一下这个切片操作：
可以看出，两个都变为长度为30的一维数组了。理解这个，就为后面更多维度的切片打下基础。
可以看出，数据由两部分组成。第一部分为经纬度与层次，第二部分为各种物理量。 前面这部分前缀为lv的表示层次，最后两个为经纬度，层次有各种划分方法。 这个文件中有很多基础物理量，举例子常用的：TMP温度 Pres气压 HGT位势高度 RH相对湿度 UGRD纬向风 VGRD经向风 CAPE 对流有效位能。 最前面的TMP表示温度，但是有9种，有的与海平面相关，有的与各层气压相关。 如第一个气温，后面的说明中表示这个只与（lat_0,lon_0）有关；第四个气温与（lv_ISBL0,lat_0,lon_0）有关。这样第一个就是二维的，可以直接绘制等值线填色图，第四个就是三维的，不能直接绘制等值线填色图，而只能在提取了某一层之后，变为二维的，才能绘制等值线填色图，如：
最关键的仍然是这一句：
而在文件属性界面，我们可以知道，lv_ISBL5表示气压层，lat_0表示纬度，lon_0表示经度。 所以[ : ]表示取全部的气压层次高度，[ 55:63 ]表示取第55至63个纬度的值（不是北纬55-63，这 个是切片序号，不是其 存 放纬度值， 具体纬度值是多少需要你去算，我选的纬度是28-35），[ 109 ]表示取第109个经度值（也是切片序号，但是恰恰其存放值为109°E），经过切片后，经度因为只取了一个值，所以被降维， 由于经度被降维了，这个相对湿度物理量只剩纬度，气压层次两维了，而两维数据就可以直接绘图了。
可以看出，这个文件里的z与五个参数有关，所以可以称为五维变量，下面是绘制其剖面图的方法：
在z[ 1 , 1 , : , 90:181 , 100 ]里，按顺序分别表示years取第一个切片值；time取第一个切片值；层次level从上至下全部取完；纬度取第90到181个切片值；经度取第100个切片值。所以，years、time、经度这三个维度遭到降维打击，那么z变为仅与lat与level相关的二维数据，就可以画等值线填色图了。【本文地址】
今日新闻 |
推荐新闻 |