在氣象中，我們常常需要用到剖面圖。地形剖面主要研究地貌對(duì)降雨、氣流的影響作用；緯度高度剖面圖主要用來(lái)分析降雨的某些條件，如濕層深厚、上干下濕、風(fēng)向風(fēng)速等。

一、地形剖面圖
繪制地形剖面圖之前，需要了解自己使用的地形文件的格式與屬性。文件為.nc格式，需要使用Python中的netCDF4或者xarray庫(kù)包來(lái)讀取。
首先我們先來(lái)讀取一下文件，并print出來(lái)，看看其屬性：

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'#地形文件儲(chǔ)存地址
f=xr.open_dataset(filename)#讀取文件
print(f)#打印其屬性

可以看出這個(gè)文件主要由x,y,z三個(gè)變量組成。
其中x表示經(jīng)度，將全球東西360經(jīng)度分為了10800刻度，相當(dāng)于一個(gè)經(jīng)度被分為30份；y表示緯度，將全球南北180緯度分為了5400份，也是將一個(gè)緯度分為30份。那么這個(gè)nc文件的精度就是0.0333°×0.0333°；z表示高度，也就是地形。可以看出，z僅僅與y,x有關(guān)，且第一相關(guān)量為y而不是。
因?yàn)槭嵌S的數(shù)據(jù)，那么按照繪制平面填色圖的ax.contourf命令是可以直接讀取數(shù)據(jù)繪圖的。
接下來(lái)我們先繪制一個(gè)平面的地形圖：

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()#縮寫(xiě)投影
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變量賦值為經(jīng)度
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的填色為海藍(lán)
lev=np.arange(0,4000,200)
norm3 = mpl.colors.BoundaryNorm(lev, cmap_new.N) #標(biāo)準(zhǔn)化level，映射色標(biāo)
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)

上下兩種命令，出圖應(yīng)該完全一樣（幾句前extent語(yǔ)句已經(jīng)限制了繪圖范圍），但是最好用上面這種，理由如下：
第二種不對(duì)導(dǎo)入的數(shù)據(jù)做取舍，那么程序在繪圖時(shí)，就會(huì)將全球都繪制出來(lái)，然后再裁剪邊界，這樣出圖效率大概慢十倍。第一種本質(zhì)上是將數(shù)據(jù)扣出一塊，只繪制這一塊，速度大大提高。

提到這里，我們不得不提下剖面圖繪制中的切片操作：
以經(jīng)度為例，前面已經(jīng)講到將一個(gè)經(jīng)度分為30份，那么我們要畫(huà)東經(jīng)70-140的圖，那就需要對(duì)經(jīng)度數(shù)據(jù)切片，原理如下（緯度同理）：
起始：（180+70）×30=7500（在前面屬性可知，切片是需加上西經(jīng)180）
終止：（180+140）×30=9600

接下來(lái)就是z的切取了，前面讀取屬性時(shí)我們已經(jīng)知道，緯度為第一相關(guān)量，經(jīng)度為第二相關(guān)量，所以應(yīng)該先切緯度，后切經(jīng)度：
height [ 2850:4960 , 7500:9600 ]

接下來(lái)，就是本節(jié)關(guān)鍵，怎么繪制地形剖面圖。

在繪制地形填色時(shí)，我們使用的是ax.contourf命令，他要求輸入橫坐標(biāo)，縱坐標(biāo)，與橫縱坐標(biāo)有關(guān)系的z值。這樣z就必須是二維的，以與橫縱坐標(biāo)相關(guān)，所以切片時(shí)，我們必須使z切取范圍與x,y完全一致，否則報(bào)錯(cuò)。

但是繪制剖面圖，我們還需不需要contourf命令呢？
顯然是不需要的，我們只想知道沿某個(gè)經(jīng)度（或緯度）的地形變化如何，用ax.plot命令結(jié)合fill_between命令即可。而這兩個(gè)命令，只需要傳入一個(gè)一維的橫坐標(biāo)，和一維的縱坐標(biāo)即可。關(guān)鍵就在怎么把z從二維的變?yōu)橐痪S的。
這就需要上面的切片方法了，比如我要畫(huà)東經(jīng)108.98°這個(gè)經(jīng)線的剖面，那就直接在z取值時(shí)，將其x取值設(shè)置為固定的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為圖形長(zhǎng)，dpi為設(shè)置圖形每英寸的點(diǎn)數(shù)
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()

進(jìn)一步print一下這個(gè)切片操作：

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

可以看出，兩個(gè)都變?yōu)殚L(zhǎng)度為30的一維數(shù)組了。理解這個(gè)，就為后面更多維度的切片打下基礎(chǔ)。

二、利用NECP的再分析資料繪制緯度高度剖面圖
先讀取查看一下屬性：

import xarray as xr
ds = xr.open_dataset(r'C:\Users\86132\fnl_20200717_00_00.nc')
print(ds)

可以看出，數(shù)據(jù)由兩部分組成。第一部分為經(jīng)緯度與層次，第二部分為各種物理量。 前面這部分前綴為lv的表示層次，最后兩個(gè)為經(jīng)緯度，層次有各種劃分方法。
這個(gè)文件中有很多基礎(chǔ)物理量，舉例子常用的：TMP溫度 Pres氣壓 HGT位勢(shì)高度 RH相對(duì)濕度 UGRD緯向風(fēng) VGRD經(jīng)向風(fēng) CAPE 對(duì)流有效位能。
最前面的TMP表示溫度，但是有9種，有的與海平面相關(guān)，有的與各層氣壓相關(guān)。
如第一個(gè)氣溫，后面的說(shuō)明中表示這個(gè)只與（lat_0,lon_0）有關(guān)；第四個(gè)氣溫與（lv_ISBL0,lat_0,lon_0）有關(guān)。這樣第一個(gè)就是二維的，可以直接繪制等值線填色圖，第四個(gè)就是三維的，不能直接繪制等值線填色圖，而只能在提取了某一層之后，變?yōu)槎S的，才能繪制等值線填色圖，如：

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][:][:]#這就變?yōu)槎S的，只取了一層次

根據(jù)之前提到的，我們現(xiàn)在要繪制一個(gè)某個(gè)經(jīng)度的垂直剖面圖，說(shuō)明我們的橫坐標(biāo)應(yīng)該是緯度，縱坐標(biāo)應(yīng)該是高度，但是在氣象上一般不使用高度，而是氣壓層，如925hPa、850hPa、700hPa、500hPa、200hPa等，而經(jīng)度就取一個(gè)固定值，這樣也能變成二維數(shù)組。下面通過(guò)一個(gè)程序講明，注釋直接在程序中：

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'][:]#讀取經(jīng)度，一維的
lat=f['lat_0'][:]#讀取緯度，一維的
RH_P0_L100_GLL0=f['RH_P0_L100_GLL0'][:][:][:]#讀取相對(duì)濕度，三維的
UGRD_P0_L100_GLL0=f['UGRD_P0_L100_GLL0'][:][:][:]#讀取緯向風(fēng)，三維的
VGRD_P0_L100_GLL0=f['VGRD_P0_L100_GLL0'][:][:][:]#讀取經(jīng)向風(fēng)，三維的
pres= f['lv_ISBL5'][:]*0.01#讀取氣壓層數(shù)，一維的
fig=plt.figure(figsize=(5,4),dpi=200)#添加畫(huà)布
ax=fig.add_axes([0,0,1,1])#添加子圖
ax.invert_yaxis()#反轉(zhuǎn)縱軸，使1000hPa作為起點(diǎn)
ax.set_yticks([1000,925,850,700,500,300])#突出我們感興趣的氣壓層
ax.set_xticks(np.arange(28,36,1))#橫坐標(biāo)
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'])#轉(zhuǎn)換為緯度格式
ax.set(ylim=(1000,300))#設(shè)置氣壓層繪圖范圍，此處我們只顯示到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='相對(duì)濕度（%）',pad=0.01)
cb.ax.tick_params(axis='both',which='both',length=1,labelsize=6)
ax.set_title('2020年7月15日20時(shí)相對(duì)濕度與風(fēng)場(chǎng)剖面圖',fontsize=15)

最關(guān)鍵的仍然是這一句：

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這個(gè)物理量與三個(gè)參量有關(guān)，按順序分別是：

而在文件屬性界面，我們可以知道，lv_ISBL5表示氣壓層，lat_0表示緯度，lon_0表示經(jīng)度。
所以[ : ]表示取全部的氣壓層次高度，[ 55:63 ]表示取第55至63個(gè)緯度的值（不是北緯55-63，這 個(gè)是切片序號(hào)，不是其 存 放緯度值， 具體緯度值是多少需要你去算，我選的緯度是28-35），[ 109 ]表示取第109個(gè)經(jīng)度值（也是切片序號(hào)，但是恰恰其存放值為109°E），經(jīng)過(guò)切片后，經(jīng)度因?yàn)橹蝗×艘粋€(gè)值，所以被降維， 由于經(jīng)度被降維了，這個(gè)相對(duì)濕度物理量只剩緯度，氣壓層次兩維了，而兩維數(shù)據(jù)就可以直接繪圖了。

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)

風(fēng)場(chǎng)這個(gè)也是同樣的道理，經(jīng)向風(fēng)與緯向風(fēng)按同樣方法切片取值。
接下來(lái)是一個(gè)五維數(shù)據(jù)的剖面圖繪制，可以幫助各位讀者更好理解切片降維方法。

可以看出，這個(gè)文件里的z與五個(gè)參數(shù)有關(guān)，所以可以稱(chēng)為五維變量，下面是繪制其剖面圖的方法：

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='多維數(shù)據(jù)的剖面圖')
plt.show()

在z[ 1 , 1 , : , 90:181 , 100 ]里，按順序分別表示years取第一個(gè)切片值；time取第一個(gè)切片值；層次level從上至下全部取完；緯度取第90到181個(gè)切片值；經(jīng)度取第100個(gè)切片值。所以，years、time、經(jīng)度這三個(gè)維度遭到降維打擊，那么z變?yōu)閮H與lat與level相關(guān)的二維數(shù)據(jù)，就可以畫(huà)等值線填色圖了。

現(xiàn)在各位應(yīng)該知道繪制剖面圖技巧了，無(wú)論有多少維度，只保留感興趣的兩維，其他維度都做降維處理，處理完的數(shù)據(jù)變?yōu)槎S，二維數(shù)據(jù)直接傳入ax.contourf（）中畫(huà)圖。

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

總結(jié)

以上是生活随笔為你收集整理的地形剖面图、纬度高度剖面图如何绘制的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問(wèn)題。

如果覺(jué)得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。