水平向量及純量之分布

Table of contents

背景
程式IO與執行
程式設計
Results
Reference

背景

wrf-python官網提供了SLP與等壓面高度分布圖的基本作法。只需進行必要的修改即可使用。
除了壓力，此處也應用到垂直溫差(KO index)分布圖的繪製，以利空污事件過程氣象數據的分析。
- 垂直擴散能力的指標有幾，經嘗試大多並不理想
  - CAPE、CIN：似與雲雨有關，並非連續的場，解析能力有限
  - K值
    - K=(T₈₅₀-T₅₀₀)+T_d850-(T₇₀₀-T_d700)
    - T=temp
    - 沒有0值的概念，不如KO直觀。
  - KO index
    - KO = ((T₇₀₀ + T₅₀₀) - (T₈₅₀ + T₁₀₀₀))/2.
    - T=ETH (Equivalent Potential Temperature)
    - 事件會伴隨KO=0通過臺灣地區(GFS與WRF皆然)

程式IO與執行

wrf檔案之準備

因跨日處理，需先將逐日wrfout連成一個檔案備用

d=/nas1/WRF4.0/WRFv4.2/201903/wrfout/wrfout_d03_2019-03-
t=_00\:00\:00
ncrcat -O ${d}29$t ${d}30$t ${d}31$t ${d}29-31
ln -s ${d}29-31 .

原程式ncfile = Dataset("wrfout_d01_2016-10-07_00_00_00")
因應不同月、日、與domain的wrfout，將輸入檔變成可程式化檔名

d=sys.argv[1] #domain
t=int(sys.argv[2]) #iobaric value
n=2     #frequency of wind bars
if d=='4':n=4
mt='03'
dt='29'
dt2='31'
ncfile = Dataset("wrfout_d0"+d+"_2019-"+mt+"-"+dt+"-"+dt2)

NaturalEarthFeature檔案之連線

這項cartopy的底圖檔，並不需要特別準備。
- 第一次使用時，程式會自動從Natural Earth官網自行下載。只需保持網路暢通即可。
- 後續同一使用者就不會需要再下載。
除了行政區界如有別的底圖需要，可以參考
- Python feature.NaturalEarthFeature方法代碼示例，或
- 官網說明。
- zip檔下載點
原程式：load檔名是admin_1_states_provinces_shp

states = NaturalEarthFeature(category="cultural", scale="50m",
                             facecolor="none",
                             name="admin_1_states_provinces_shp")
ax.add_feature(states, linewidth=.5, edgecolor="black")
ax.coastlines('50m', linewidth=0.8)

因臺灣範圍較小，需有縣市界，因此增加解析度到10m，

從臺灣屬中國範圍，檔案名稱沒有_shp

修改後：

# Download and add the states and coastlines
states = NaturalEarthFeature(category="cultural", scale="10m",
                         facecolor="none",
                         name="admin_1_states_provinces")
ax.add_feature(states, linewidth=0.5, edgecolor="black")
ax.coastlines('50m', linewidth=0.8)

輸出檔案

原範例是show圖在x window上，因為需要作的圖頗多，如果能存檔將可有利檔案管理。
檔名變數
- iso：壓力值(hPa)
- t：wrfout檔案內的小時順序，0~nt-1
- d：domain,1~4

...
#plt.show()
plt.savefig(str(iso)+'_'+'{:02d}'.format(t)+'D'+d+'.png')

程式執行方式

fig = plt.figure(figsize=(12,9))指令似乎需要close，不能在時間的do loop中重複開啟。目前仍採一次一圖方式執行
好在不同時間作圖彼此不會干擾，可以平行運作。
sub(mit)的用法，可以詳unix系統小工具-sub

for d in 1 2 3 4;do for t in {4..40..6};do sub python iso.py $d $t;done;done

for d in 1 2 3 4;do for t in {4..40..6};do sub python KO_index.py $d $t;done;done

程式設計

等值區間的設定

因程式同時應用在不同的小時序，因此最大、最小值的範圍要先讀完所有時間的數據後，統一設定
讀取特定變數場：使用getvar，可以參考變數分類定義，或官網程式範例。
tlst：time frame list，共7筆
fac：因最大濃度會超出最後一個等候，wrf-python的contourf程式會將其空白，因此範圍需拉大一些，不能正好等於最大值。
iso.py
取exp再取log值：因為等值線偏向高值，在高值部分較為平緩，無法辨識。

...
shp = [7]+list(getvar(ncfile, "pressure").shape)
p0=np.zeros(shape=shp)
z0=np.zeros(shape=shp)
tlst=[i for i in range(4,41,6)]
for it in tlst:
  i=tlst.index(it)
  p0[i,:,:,:] = getvar(ncfile, "pressure",timeidx=it)
  z0[i,:,:,:] = getvar(ncfile, "z", units="dm",timeidx=it)
ht_iso0 = interplevel(z0, p0, iso)
fac=1.01
n=2     #frequency of wind bars
if d=='4':
  n=5
  fac=1.0007
elif d=='3':
  fac=1.002
mx,mn=10**(np.max(ht_iso0)*fac/100.),10**(np.min(ht_iso0)/100)

levels0= np.arange(mn, mx, (mx-mn)/10)
levels = [np.log10(i)*100 for i in levels0]
...

KO_index.py

KO圖同時也有omega(omg)的等值圖，因此需要3個純量場
omg內設單位Pa/s，需轉成hPa/hr，以與國際接軌。
等位溫在低層有可能不存在，繪圖時會發生斷線。因此需將山區部分補滿。
- 此處假設山區範圍低層的等位溫與高空(上一層)相同。
- 在700hPa面(約3000M)少部分山頂，KO總是等於0。
KO以等值線表示，沒有空白的問題
omg以shade表示，設成固定的區間，以與wetter3.de圖形比較。不另計算等值區間。

p0=np.zeros(shape=shp)
T=np.zeros(shape=shp)
W=np.zeros(shape=shp)
tlst=[i for i in range(4,41,6)]
for it in tlst:
  i=tlst.index(it)
  p0[i,:,:,:] = getvar(ncfile, "pressure",timeidx=it)
  T[i,:,:,:] = getvar(ncfile, "eth",units="K",timeidx=it)
  W[i,:,:,:] = getvar(ncfile, "omg",timeidx=it)*3600./100.
lvls=[500,700,850,1000]
for lvl in lvls:
  ilvl=lvls.index(lvl)
  slvl=str(lvl)
  exec("T"+slvl+" = interplevel(T, p0, "+slvl+")")
#  check eth at mountain region, if nan, let equal to upper layer THE(adiabatic condition)
  if lvl>500:
    uplvl=str(lvls[ilvl-1])
    exec("T"+slvl+" = np.where(np.isnan(T"+slvl+"),T"+uplvl+",T"+slvl+")")
KO0 = ((T700 + T500) - (T850 + T1000))/2.
mx,mn=np.round(np.max(W)*1.2,0),np.round(np.min(W),0)
mn,mx=-54,1
levels = list(range(mn,-22,8))+list(range(-22,mx,2))

等壓面垂直內插

使用interplevel程式，用法也可以參考官網其他範例。
以下以iso.py為例，KO_index.py也類似作法。
wspd：風速，因getvar是將風速、風向放在同一個變數內，取出時要指定第一個rank([0,:])。(最後因500hPa面在中小尺度沒有太大的風速特性，因此沒有繪出等值圖)。

p = getvar(ncfile, "pressure",timeidx=t)
z = getvar(ncfile, "z", units="dm",timeidx=t)
ua = getvar(ncfile, "ua", units="kt",timeidx=t)
va = getvar(ncfile, "va", units="kt",timeidx=t)
wspd = getvar(ncfile, "wspd_wdir", units="kts",timeidx=t)[0,:]

# Interpolate geopotential height, u, and v winds to 500 hPa
ht_iso = interplevel(z, p, iso)
u_iso = interplevel(ua, p, iso)
v_iso = interplevel(va, p, iso)
wspd_iso = interplevel(wspd, p, iso)

空間座標位置與投影

藏在interplevel method結果內。

# Get the lat/lon coordinates
lats, lons = latlon_coords(u_iso)

# Get the map projection information
cart_proj = get_cartopy(u_iso)

Create the figure

因與投影方式有關，因此放在座標率定之後

...
fig = plt.figure(figsize=(12,9))
ax = plt.axes(projection=cart_proj)
...

等值圖

shade填滿圖：contourf為matplotlib程式。用法也可以參考wrf-python官網範例。
重要引數包括
- x, y, value, levels
- cmap或colors
- transform：座標轉換
- norm：均化指標。會將所有值填滿在色階範圍。
- extend：{‘neither’, ‘both’, ‘min’, ‘max’}, default: ‘neither’，超出範圍要怎麼填滿。內設是留空白。
色標legend：colorbar
- orientation：垂直(vertical)或水平(horizontal)

# Add the 500 hPa geopotential height contours
contours = plt.contourf(to_np(lons), to_np(lats), field,
                       levels=levels, #colors="rainbow",#"black",
                       cmap=get_cmap("plasma"),
                       transform=crs.PlateCarree())
plt.colorbar(contours, ax=ax, orientation="vertical",pad=.05) 

等值線圖：contour
- extend：等值線是在數據點位置或延伸到圖框邊界上
- linewidths：線寬，內設1.5
clabel：線標籤

contours = plt.contour(to_np(lons), to_np(lats), to_np(KO),
                             levels=levels,
                             linewidths=1.,
                             cmap=get_cmap("winter_r"), #colors="black",
                             transform=crs.PlateCarree())
plt.clabel(contours, inline=1, fontsize=10, fmt="%i")

此範例是分開設定，也可以用matplotlib.contour.ContourSet整個包在一起設定
色標
- 適用所有matplotlib選項，
- 選擇考量見wrf-python 安裝與基本指令色標

向量圖

因高層風速較大，因此選擇用風標(Barbs)來表示。
因為是等壓面，假設垂直方向可以忽略，直接以ua、va來繪製
貼上向量
各domain風標密度：需嘗試錯誤得到較佳效果。

...
plt.barbs(to_np(lons[::n,::n]), to_np(lats[::n,::n]),
          to_np(u_iso[::n, ::n]), to_np(v_iso[::n, ::n]),
          transform=crs.PlateCarree(), length=5)
...          

程式碼

Results

iso.py結果比較


Domain3範圍事件日(左)及非事件日(右)500hPa等壓面高度

KO_index.py結果比較


Domain3範圍事件日及非事件日800hPa等壓面風標、500hPa垂直運動Omega及KO等值線

水平向量及純量之分布

背景

程式IO與執行

wrf檔案之準備

NaturalEarthFeature檔案之連線

輸出檔案

程式執行方式

程式設計

等值區間的設定

iso.py

KO_index.py

等壓面垂直內插

空間座標位置與投影

Create the figure

等值圖

向量圖

程式碼

Results

Reference