地形剖面图、纬度高度剖面图如何绘制

在气象中，我们常常需要用到剖面图。地形剖面主要研究地貌对降雨、气流的影响作用；纬度高度剖面图主要用来分析降雨的某些条件，如湿层深厚、上干下湿、风向风速等。

一、地形剖面图
绘制地形剖面图之前，需要了解自己使用的地形文件的格式与属性。文件为.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)#打印其属性

可以看出这个文件主要由x,y,z三个变量组成。
其中x表示经度，将全球东西360经度分为了10800刻度，相当于一个经度被分为30份；y表示纬度，将全球南北180纬度分为了5400份，也是将一个纬度分为30份。那么这个nc文件的精度就是0.0333°×0.0333°；z表示高度，也就是地形。可以看出，z仅仅与y,x有关，且第一相关量为y而不是。
因为是二维的数据，那么按照绘制平面填色图的ax.contourf命令是可以直接读取数据绘图的。
接下来我们先绘制一个平面的地形图：

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'#地形文件地址
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=(14,9),dpi=100)
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=8)
plt.savefig('demo.png',bbox_inches='tight') #存图
plt.show()

出图如下：

其中，最重要的是绘图的这一句：

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)

上下两种命令，出图应该完全一样（几句前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()

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

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

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

二、利用NECP的再分析资料绘制纬度高度剖面图
先读取查看一下属性：

import xarray as xr
ds = xr.open_dataset(r'C:\Users\86132\fnl_20200717_00_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 xr
ds = xr.open_dataset(r'C:\Users\86132\fnl_20200717_00_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 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)

最关键的仍然是这一句：

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 xr
import numpy as np
import matplotlib.pyplot as plt
plt.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（）中画图。

本文摘自：https://my.oschina.net/u/4579644/blog/4518065

地形剖面图、纬度高度剖面图如何绘制相关推荐

unity地形之splatalpha研究地形贴图导出更换与绘制
unity中的地图贴图的绘制常常使用的是paint texture里面的但是这个方式往往费时很多,却只能做出很少的效果,这里要介绍的就是通过外部绘制splatalpha 来替换,达到unity中地形 ...
matlab 什么是剖面图,什么是剖面图剖面图作用及注意事项
剖面图一般用于工程的施工图和机械零部件的设计中,补充和完善设计文件,是工程施工图和机械零部件设计中的详细设计,用于指导工程施工作业和机械加工.剖面图的作用是什么呢?在制作剖面图时有哪些需要注意的地方, ...
飞机3D轨迹绘制(经度-纬度-高度）
使用Python绘制 #绘制三维直线图,将飞机飞行的航迹用(经度,纬度和高度)来描述 #******************************************************** ...
python opengl 截图_初试PyOpenGL二 (Python+OpenGL)基本地形生成与高度检测
在上文中,讲述了PyOpenGL的基本配置,以及网格,球形的生成,以及基本的漫游.现在利用上一篇的内容,来利用高程图实现一个基本的地形,并且,利用上文中的第三人称漫游,以小球为视角,来在地形上前后左右 ...
Ogre 天龙八部地形 Heightmap（高度图）+GridInfo（地表信息）初步结果
分享一下我老师大神的人工智能教程!零基础,通俗易懂!http://blog.csdn.net/jiangjunshow 也欢迎大家转载本篇文章.分享知识,造福人民,实现我们中华民族伟大复兴! 刚研究出 ...
geotif 添加坐标_tiff和geotiff经度纬度高度值读取
使用 tiff3.8.2 和 geotiff1.2.5 读取一个 geotiff 格式的图像,获取图片的经度.纬度和高度值. 1. 经度和纬度可以通过 geotiff 读到栅格坐标和地理坐标的换算关系 ...
泰山OFFICE技术讲座：由行的布局高度，谈绘制高度的高度溢出、高度缩水（全网首发）
在之前的博文中,吾提出除,行视图除了字体高度,还有一个布局高度.现在又要提出绘制高度的概念.什么叫绘制高度?比如说,布局高度是18.1,绘制时的高度一般是18,也有可能不是.这个时候就需要绘制高度. ...
python绘制剖面图_Python气象绘图教程—（十九）剖面图
本节提要:简要谈谈地形剖面图.纬度高度剖面图.时间纬度图的绘制方法. 提要中提到的这几种图形都是在气象上比较常用的,地形剖面主要研究地貌对降雨.气流的影响作用:纬度高度剖面图可以用来分析降雨的某些条件 ...
【ArcGIS】绘制地形剖面图
数据准备 DEM高程(tiff格式) 1.使用ArcGIS打开 DEM 数据 2.右击工具条空白处打开3D Analyst工具条 3.使用插入线工具绘制剖面线 (双击结束编辑) 4.使用创建剖 ...

地形剖面图、纬度高度剖面图如何绘制

地形剖面图、纬度高度剖面图如何绘制相关推荐

最新文章

热门文章