Link Search Menu Expand Document
Focus on Air Quality

垂直向量及純量之分布

Table of contents

背景

NCL、wrf-python與CCTM純量場

  • NCL 雖然已經有很多範例，但是在垂直面的向量與純量分布圖這項仍然很弱。
    • 主要因為wrf或CCTM系統的垂直座標並非等間距，還需要垂直內插。
  • NCL-for-CMAQ雖然有些範例是CCTM純量的垂直分布，但X軸是時間，並非水平座標，無法套用地形。
  • 最困難的還是內插到統一的垂直網格系統、與地形地圖之間能沒有間隙。
    • 這點wrf-python可以做到最好。
    • 但是目前wrf-python的getvar還不能讀CMAQ系統的檔案，必須自行轉接。

垂直剖面內插程式

引用模版及副程式

  • 模版程式：cross-section-with-mountains
  • 副程式
    • None

      IO Files and Coordinate translation

  • Input Files
    • wrfout檔案：作為空間定位的模版、地形數據
    • CCTM_ACONC檔案：與wrfout解析度相同、只是座標平移、套入wrfout的網格系統中，以便應用wrf-python內插程式庫
...
wrf_file = Dataset("wrfout_d04_2019-01-29-31")
nc=Dataset('BASE3_K24.nc')
t=int(sys.argv[1])
pm=nc["PM25_TOT"][t,:,:,:]
nc=Dataset('METDOT3D_Taiwan.nc')
t0=35*24+t
UWIND=nc["UWIND"][t0,:,:-1,:-1]
VWIND=nc["VWIND"][t0,:,:-1,:-1]
...
# Get the WRF variables
ht = getvar(wrf_file, "z", timeidx=t)
ter = getvar(wrf_file, "ter", timeidx=-1)
u = getvar(wrf_file, "ua", timeidx=t)
v = getvar(wrf_file, "va", timeidx=t)
  • 平移CMAQ vs WRF
#u[:24,8:8+131,:92-12] = UWIND[:,:,12:]
#v[:24,8:8+131,:92-12] = VWIND[:,:,12:]
w = getvar(wrf_file, "wa", timeidx=t)
#U,W= 0.1**(u*cost+v*sint),0.1**(w*34)
U,W= u*cost+v*sint,w*34
dbz = getvar(wrf_file, "dbz", timeidx=-1)
dbz[:11,8:8+131,:92-12] = pm[:,:,12:]
  • Output Files
    • pm25_nn.png：nn=00~71
    • 為典型matplotlib存檔方式
pyplot.savefig('pm25_'+'{:02d}'.format(t)+'.png')

設定

  • 垂直剖線的起訖點：不限定是X或Y方向、可以是任意點
    # Define the cross section start and end points
cross_start = CoordPair(lat=24.453917871182558, lon=120.225062233815342,)
cross_end = CoordPair(lat=23.41214780981217, lon=121.79586963982419)
  • 高度：因wrfout有40層之多，大多是自由流範圍，因此必須限定繪圖的上邊界，凸顯邊界層現象。
dbz_contours = ax_cross.contourf(xs,
                                 ys[:-80],
                                 to_np(dbz_cross_filled)[:-80,:],
                                 levels=dbz_levels,
                                 cmap="rainbow",
                                 norm=dbz_norm,
                                 extend="max")
  • X軸標籤的有效位數
    • wrf-python使用經緯度tuple作為標籤，有效位數未經修剪長短不齊。此處統一修為2碼
    • 由於原程式是運用to_np之解讀程式，將coord_pairs經緯度值讀成字串，需將其拆解後方能改變有效位數。
x_labels = [pair.latlon_str() for pair in to_np(coord_pairs)]
lab=[[round(float(i),2) for i in j.split(',')] for j in x_labels[:]]
x_labels=[str(i[0])+','+str(i[1]) for i in lab]
  • 任意方向上的沿流方向風速分量
    • 計算分量後，U及W與濃度一樣進行垂直內插到等間距網格上
    • $ U = u \times cost + v \times sint $;
      • lent = np.sqrt((x1-x0)**2+(y1-y0)**2)
      • cost = (x1-x0)/lent
      • sint = (y1-y0)/lent
  • 貼上向量
    • 參考stackoverflow
    • 使用matplotlib.quiver指令來畫箭頭，中文可參考程式人生
    • 同樣使用[:-80]來限定高度範圍、[::2]來控制箭頭的密度
ax_cross.quiver(xs[::2], ys[:-80],
          to_np(u_cross_filled[:-80, ::2]), to_np(w_cross_filled[:-80, ::2]))

色標

  • 原程式為自主設定。不但顏色沒有連續、也不具辨識能力。
  • cmap選項有：”jet”、”rainbow”、
  • 適用所有matplotlib選項
    • plot_color_gradients(‘Perceptually Uniform Sequential’, [‘viridis’, ‘plasma’, ‘inferno’, ‘magma’, ‘cividis’])
    • plot_color_gradients(‘Sequential’, [‘Greys’, ‘Purples’, ‘Blues’, ‘Greens’, ‘Oranges’, ‘Reds’, ‘YlOrBr’, ‘YlOrRd’, ‘OrRd’, ‘PuRd’, ‘RdPu’, ‘BuPu’, ‘GnBu’, ‘PuBu’, ‘YlGnBu’, ‘PuBuGn’, ‘BuGn’, ‘YlGn’])
    • plot_color_gradients(‘Sequential (2)’, [‘binary’, ‘gist_yarg’, ‘gist_gray’, ‘gray’, ‘bone’, ‘pink’, ‘spring’, ‘summer’, ‘autumn’, ‘winter’, ‘cool’, ‘Wistia’, ‘hot’, ‘afmhot’, ‘gist_heat’, ‘copper’])
    • plot_color_gradients(‘Diverging’, [‘PiYG’, ‘PRGn’, ‘BrBG’, ‘PuOr’, ‘RdGy’, ‘RdBu’, ‘RdYlBu’, ‘RdYlGn’, ‘Spectral’, ‘coolwarm’, ‘bwr’, ‘seismic’])
    • plot_color_gradients(‘Cyclic’, [‘twilight’, ‘twilight_shifted’, ‘hsv’])
    • plot_color_gradients(‘Qualitative’, [‘Pastel1’, ‘Pastel2’, ‘Paired’, ‘Accent’, ‘Dark2’, ‘Set1’, ‘Set2’, ‘Set3’, ‘tab10’, ‘tab20’, ‘tab20b’, ‘tab20c’])
    • plot_color_gradients(‘Miscellaneous’, [‘flag’, ‘prism’, ‘ocean’, ‘gist_earth’, ‘terrain’, ‘gist_stern’, ‘gnuplot’, ‘gnuplot2’, ‘CMRmap’, ‘cubehelix’, ‘brg’, ‘gist_rainbow’, ‘rainbow’, ‘jet’, ‘turbo’, ‘nipy_spectral’, ‘gist_ncar’])
  • 反轉色標：XXX_r。
    • rainbow -> rainbow_r
cmap="rainbow",

$grep color $(findc "*.py")
./Air_Increment_tool/Lib/plot2D.py...
...
import matplotlib.colors as colors
...
      colorlist = ['white','deepskyblue','forestgreen','gold','red','purple']
      cmap = colors.LinearSegmentedColormap.from_list('AAA', colorlist)

Parallel Operation

  • 因程式沒有時間前後的交互作用，各個時間可以獨立運作。
  • 將時間(小時順序)作為引數
for i in {4..40..6};do sub ver_ZhonBu.py $i;done

GIF Producing

convert pm2.5*.png PMF.gif

程式碼

Results

gif player

PM2.5之垂直分布

wrf-python.png
圖 CCTM模擬d04範圍中部地區垂直PM2.5濃度分布(wrf-python繪製)，動態播放，單位μg/M3

位溫之垂直分布

  • 位溫垂直分布是地區穩定度的重要指標。
  • 受到500hPa高空風的影響，高空氣流也會在中央山脈的阻擋而形成重力波的現象，造成垂直位溫產生傾斜。
  • 右圖為一般天情況，500hPa為西南風，與山脈走勢相同，重力波現象較為輕微。
  • 左圖為事件日，山脈東側3 ~ 5Km下沉嚴重，穩定度高，限制850hPa以下東風過山，產生繞流渦漩現象。山脈西側1.5 ~ 3Km高度氣層嚴重壓縮。
臺灣中部西北至東南剖面位溫之垂直分布(動態播放，單位K)

Reference