Link Search Menu Expand Document

水平向量及純量之分布

Table of contents

背景

  • wrf-python官網提供了SLP與等壓面高度分布圖的基本作法。只需進行必要的修改即可使用。
  • 除了壓力,此處也應用到垂直溫差(KO index)分布圖的繪製,以利空污事件過程氣象數據的分析。
    • 垂直擴散能力的指標有幾,經嘗試大多並不理想
      • CAPE、CIN:似與雲雨有關,並非連續的場,解析能力有限
      • K值
        • K=(T850-T500)+Td850-(T700-Td700)
        • T=temp
        • 沒有0值的概念,不如KO直觀。
      • KO index
        • KO = ((T700 + T500) - (T850 + T1000))/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檔案之連線

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)

空間座標位置與投影

# 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")

向量圖

  • 因高層風速較大,因此選擇用風標(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結果比較
wrf-python.png
Domain3範圍事件日(左)及非事件日(右)500hPa等壓面高度
  • KO_index.py結果比較
wrf-python.png
Domain3範圍事件日及非事件日800hPa等壓面風標、500hPa垂直運動Omega及KO等值線

Reference