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

一、地形剖面圖
繪制地形剖面圖之前，需要了解自己使用的地形文件的格式與屬性。文件為.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

總結

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

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。