將CWB數據填入WRF客觀分析場

背景

在預備WRF的初始與邊界場時，REAL會讀取OBSGRID的客觀分析結果(metoa_em、auxinput1)，一般是以點狀觀測值，內插到3度空間的格點位置。
由於中央氣象局已經完成WRF模擬，不需要重新進行積分，反而有需要取得WRF執行結果中雲、邊界層數據，因此，需要將CWBWRF結果轉換成metoa_em檔案，以將其結果納入WRF模式中進行FDDA。
客觀分析結果檔(metoa_em)的內容
- 基本上的格式、變數定義等，與metgrid結果檔(met_em)完全一樣，只差後者可能沒有海溫(SST)。
- 垂直軸名稱為num_metgrid_levels，涉及6項變數：['PRES', 'GHT', 'RH', 'VV', 'UU', 'TT']
  - 壓力：除了地面之外，其餘33層皆為定壓層。範圍由1000mb~1mb。
  - RH相對濕度，單位為%
  - UU,VV為風速，單位為m/s
  - TT為氣溫(非位溫)，單位為K。
- 其他項目、或這些項目100mb以上高空，與met_em相同即可

CWBWRF與metoa_em壓力層之比較

項目	CWBWRF	metoa_em	說明
層數	11	34
第1層	1000mb	地面	風速可以使用uv10、溫度可以使用T2、相對濕度不另計算，令為1000mb值
最高層	100mb	1mb	CWBWRF僅到34層中的22層
相容性	-	-	CWBWRF每層都在metoa_em中

內插方法

水平網格
1. 讓metoa_em與CWBWRF(儘可能)完全一樣。不進行內插(CWBWRF_15k及45k版本)
2. 略有差異、但因網格較粗，選擇最近點。(mk_metoaD12.py)
3. 需精確對準座標值、採scipy.griddata內插([mk_metoaT.py][mk_metoaT.py]及mk_metoaTD3.py平行運作版本)
垂直內插
- 地面無法內插。另外處理
- CWBWRF定壓層：不必計算，直接引用
- 其他壓力層：使用wrf-python的內插函數 wrf.interplevel，不另外設計內插法。
時間內插
- 雖然CWBWRF為6小時輸出、預報84小時，但在ROOT/Data/cwb/WRF_3Km/get_M-A0064.cs中已使用Cubic Spline將其內插到逐時，並且轉成wrfout的格式。
- 為提供FDDA對模擬的影響，此處使用3小時為頻率、預報長度在GFS之180小時之內，採168小時(7天)。

mk_metoa.py程式說明

高度之對應

p：使用wrf-python的getvar函數讀取CWBWRF的壓力(定壓層)、詳參[[2022-08-11-wrf_pythonTAB]]
p_met：metoa_em的模版就是met_em檔案，也讀取其中34層壓力。除地面外，其餘也是定壓層，在metgrid.exe中所定義。
找到p在p_met的序號、並做成kdic備用

fname='CWB_wrfout_d01'
wrfin = Dataset(fname,'r')
p = getvar(wrfin, "pressure") #p is constant and 
...
fname='met_em.d01.'+fn_dt[0]+'.nc'
metin = Dataset(fname,'r')
p_met=getvar(metin, "pressure")
k_met=[list(p_met[:,0,0].values).index(k) for k in p[:,0,0].values]
kdic={k_met[k]:k for k in range(11)}

時間之對應

met_em及metoa_em都是單一筆時間在一個檔案內，其時間標籤被metgrid.exe設計為檔名的一部分。只需知道起迄時間(staend)，所有檔名都可以datetime計算而得(fn_dt)。
- 如果此時間在CWBWRF檔案中找得到，則讀取wrfout將其值內插到metoa_em檔中
- 如否：met_em也可以做為FDDA的依據，只是是GFS 1度內插的結果，解析度較差。
起迄時間是namelist.input的內容，以python讀取出來成為一個dict(staend)
以時間差計算nt，要記得前後都算所以必須+1。
fn_dt：以wrf系統特有的時間標籤。每3小時輸出備用。

#get the start and end time of simulation
with open('namelist.input','r') as f:
  lines=[i for i in f]
inputs={}
for i in lines[1:]:
  itm=i.split()
  if len(itm)<3:continue
  try:
    inputs.update({itm[0]:float(itm[2].replace(',',''))})
  except:
    continue
staend={}
for var in ['sta','end']:
  tm=[i for i in inputs if var == i[:3]]
  for t in tm:
    inputs.update({t:int(inputs[t])})
  staend.update({var:datetime.datetime(inputs[tm[0]],inputs[tm[1]],inputs[tm[2]],inputs[tm[3]],)})
nt=int((staend['end']-staend['sta']).total_seconds()/3600./3)+1
fn_dt=[(staend['sta']+datetime.timedelta(hours=3*t)).strftime('%Y-%m-%d_%H:00:00') for t in range(nt)]

CWBWRF基本上是一個wrfout檔案，因此其時間標籤是個變數，將其讀出成序列以利比對。

fname='CWB_wrfout_d01'
wrfin = Dataset(fname,'r')
nt_cwb=wrfin.dimensions['Time'].size
strT=[''.join([i.decode('utf-8') for i in wrfin.variables['Times'][t,:]]) for t in range(nt_cwb)]

時間標籤匹配：tcwb為CWBWRF中的時間序位

...
  if fn_dt[t] in strT:
    tcwb=strT.index(fn_dt[0])
  else:
    continue

氣象要素的讀取

使用wrf-python的getvar有其便利性，且不必再進行複雜的計算、單位轉換等過程，也可以避免錯誤的發生。
- getvar項目內容可以參考官網，或分類列表
地面物理量
- uv10為U10及V10的合併
- T2單位為K、tc單位為度C

  p_met=getvar(metin, "pressure")
  uv10=getvar(wrfin,"uvmet10",timeidx=tcwb)
  t2=getvar(wrfin,"T2",timeidx=tcwb)
  for itm in ['ua','va','rh','tc']:
    cmd='cwb=getvar(wrfin, "'+itm+'",timeidx=tcwb)'
    exec(cmd)

高度內插

地面值：因每個項目的處理方式都不太一樣，只好以if…else來處理。
1~22層依序處理
- 如果是CWBWRF的11層，直接移轉矩陣。不要再內插。內插反而會出錯(1000mb以下無值)。
- 雖然官網手冊說可以整批一起內插，但經T/E發現只能一層層進行。
回存時要注意metoa_em的風速UU(在x方向)及VV(在y方向)的維度會多1格，與getvar讀的結果不同，因此需限制其最大維度。

    var=np.zeros(shape=p_met.shape)
    if itm[1]=='a':
      var[0,:,:]=uv10[uvi[itm],:,:]
    elif itm=='tc':
      cwb+=273.
      var[0,:,:]=t2[:,:]
    else:
      var[0,:,:]=cwb[0,:,:]
    for k in range(1,22): #34-12):
      if k in kdic:
        var[k,:,:]=cwb[kdic[k],:,:]
      else:
        var[k,:,:]=interplevel(cwb,p,p_met[k,:,:])
    metin.variables[vdic[itm]][0,:22,:ny,:nx]=var[:22,:,:]

CWBWRF_15k vs 45k(mk_metoa45.py)

15k版本裏新舊網格系統完全相同，只有高度時間的對應。
在45k版本，新網格是每3格採1次舊網格(CWBWRF_15k)數據，只要在水平xy平面上，以間格方式對照即可。

31a32
> nz0,ny0,nx0=p.shape
41c42
<
---
> N=nx0//nx
60c61
<       var[0,:,:]=uv10[uvi[itm],:,:]
---
>       var[0,:,:]=uv10[uvi[itm],:-1:N,:-1:N]
63c64
<       var[0,:,:]=t2[:,:]
---
>       var[0,:,:]=t2[:-1:N,:-1:N]
65c66
<       var[0,:,:]=cwb[0,:,:]
---
>       var[0,:,:]=cwb[0,:-1:N,:-1:N]
68c69
<         var[k,:,:]=cwb[kdic[k],:,:]
---
>         var[k,:,:]=cwb[kdic[k],:-1:N,:-1:N]
70c71
<         var[k,:,:]=interplevel(cwb,p,p_met[k,:,:])
---
>         var[k,:,:]=interplevel(cwb[:,:-1:N,:-1:N],p[:,:-1:N,:-1:N],p_met[k,:,:])

東海台海版本

因為座標系統需要精準校對，此處改以scipy.griddata內插取代前述直接引用。
CWBWRF_3Km格點多、檔案龐大，在get_M-A0064.cs中已經處理成逐日檔(約8.1G)，不是所有預報時間都在同一檔案內，需要一個目錄對照表，以讀取到正確的時間。

CWBWRF_3Km網格定義

使用getvar取得wrfout的經緯度(為YX 2維矩陣)
使用pyproj.Proj模組轉到目標網格系統(metin)

fname='CWB_wrfout_d03'
wrfin = Dataset(fname,'r')
p = getvar(wrfin, "pressure") #p is constant and homogeneous in x and y
latm,lonm=getvar(wrfin,'lat'),getvar(wrfin,'lon') #CWB net
...
pnyc = Proj(proj='lcc', datum='NAD83', lat_1=metin.TRUELAT1, lat_2=metin.TRUELAT2, lat_0=metin.CEN_LAT, lon_0=metin.CEN_LON, x_0=0, y_0=0)
x,y=pnyc(lonm,latm, inverse=False) #CWB net

目標網格系統

為要套用環保署公版模式的排放數據，需要將網格系統設計與其完全一致。依據其grid03/mcip/GRIDDESC內容如下：

$ cat /nas2/cmaq2019/download/input/201901/grid03/mcip/GRIDDESC
' '
'Taiwan'
  2        10.000        40.000       120.000       120.000        25.000
' '
'Taiwan'
'Taiwan'    -72000.000   -345000.000      3000.000      3000.000  92 131   1
' '

其中心位置(25,120)位在台海中線偏北方向，並非以臺灣為中心。
- 之所以會以此為模擬中心，應為風場有更大範圍之需求。
- 此處因主要援引中央氣象局的數值預報，做為4階同化的依據，並不重新進行動力模式模擬。範圍較小一些，可以取得mcip所需wrfout項目即可。
如要含概原點範圍，至少需有(201 × 301)點，
- 範圍如圖所示
- 取單數點，讓GRIDDESC原點可以保持是dx、dy的整數倍。


東海台海模擬範圍

fname='met_em.d01.'+fn_dt[0]+'.nc'
metin = Dataset(fname,'r')
nz,ny,nx=p_met.shape

x1d=[metin.DX*i for i in range(-nx//2,nx//2)]
y1d=[metin.DY*i for i in range(-ny//2,ny//2)]
x1,y1=np.meshgrid(x1d, y1d)
maxx,maxy=x1[-1,-1],y1[-1,-1]
minx,miny=x1[0,0],y1[0,0]
boo=(abs(x) <= (maxx - minx) /2+metin.DX*10) & (abs(y) <= (maxy - miny) /2+metin.DY*10)
idx = np.where(boo)
mp=len(idx[0])
xyc= [(x[idx[0][i],idx[1][i]],y[idx[0][i],idx[1][i]]) for i in range(mp)]

griddata內插副程式

因為需要內插的2維變數太多樣了，寫成簡單的副程式會比較方便

def interpXY(arr):
  c = [arr[idx[0][i], idx[1][i]] for i in range(mp)]
  return griddata(xyc, c, (x1, y1), method='linear')

內插後存到var矩陣內，再一併回存到metin。
地面風速(uv10)、地面溫度(t2)、相同定壓層變數、需垂直內插之變數(先執行垂直內插)等，都呼叫interpXY

  for itm in ['ua','va','rh','tc']:
    cmd='cwb=getvar(wrfin, "'+itm+'",timeidx=tcwb)'
    exec(cmd)
    var=np.zeros(shape=(nz,ny,nx))
    if itm[1]=='a':
      var[0,:,:]=interpXY(uv10[uvi[itm],:,:])
    elif itm=='tc':
      cwb+=273.
      var[0,:,:]=interpXY(t2[:,:])
    else:
      var[0,:,:]=interpXY(cwb[0,:,:])
    for k in range(1,22): #34-12):
      if k in kdic:
        var[k,:,:]=interpXY(cwb[kdic[k],:,:])
      else:
        var[k,:,:]=interpXY(interplevel(cwb[:,:,:],p[:,:,:],p_met[k,0,0]))
    metin.variables[vdic[itm]][0,:22,:ny,:nx]=var[:22,:,:]

時間標籤與CWBWRF檔案目錄的對照(dates)

fn_dt為檔名中的時間標籤，與前述版本都一樣，沒有更動。逐3小時一筆
fn_dt對照到WRF_3Km下存檔目錄(日期)，dates的關係如下所示

...
fn_dt=[(staend['sta']+datetime.timedelta(hours=3*t)).strftime('%Y-%m-%d_%H:00:00') for t in range(nt)]
#path dic of CWB_wrfouts, searched by fn_dt
dates={fn_dt[t]:(staend['sta']+datetime.timedelta(hours=t*3-6)).strftime('%Y%m%d') for t in range(nt)}
...

WRF_3Km自06Z開始存檔，因此時間需扣6小時
日期格式為YYYYMMDD，即’%Y%m%d’
對照表應用在時間的迴圈內，因為是唯讀檔，不需特別予以關閉。strT也是隨著檔案開啟後再產生，以找到正確的序位(tcwb，用法與前述版本都一樣)。

for t in range(nt):
  fname='met_em.d01.'+fn_dt[t]+'.nc'
  fnameO=fname.replace('met','metoa')
  os.system('cp '+fname+' '+fnameO)
  # only if fn_dt in strT(time overlaped)
  fname='/nas1/Data/cwb/WRF_3Km/'+dates[fn_dt[t]]+'/wrfout_d03_:00:00'
  wrfin = Dataset(fname,'r')
  nt_cwb=wrfin.dimensions['Time'].size
  strT=[''.join([i.decode('utf-8') for i in wrfin.variables['Times'][t,:]]) for t in range(nt_cwb)]
  if fn_dt[t] in strT:
    tcwb=strT.index(fn_dt[t])
  else:
    continue

平行計算版本(mk_metoaT.py)

因griddata計算需時較久(處理一個metoa檔約需20分鐘)，如循序處理，會需要近19小時，因此必須同步處理不同時間的metoa檔。
設計將序號當成引數，前述時間迴圈只進行該小時之內插。

$ diff mk_metoa45.py mk_metoaT.py
0a1
> #argument with time t=0~nt-1
5c6
< import os
---
> import os, sys
63c64,66
< for t in range(nt):
---
> tin=int(sys.argv[1])
> if tin>=nt:sys.exit('wrong T')
> for t in range(tin,tin+1):

因同時進行開啟、複製檔案，會造成記憶體嚴重耗盡，讓每批次工作之間相隔10秒鐘：

for i in {0..56};do sub python mk_metoaT.py $i;sleep 10s;done

程式碼精進：去除if block

雖然不會浪費判斷的時間，但是python語法有dictionary可以用，還用if block的寫法似乎就沒有必要，也缺乏簡捷性與系統性。
需要整併的項目有2
1. 地面變數的內容指定到t2、uv10，這會用到exec指令
2. 溫度還需要加273、其他項則不需要

修改前

‘ua’及’va’地面風指向uv10、地面溫度指向t2
其餘：地面則為cwb[0,:,:]

    if itm[1]=='a':
      var[0,:,:]=interpXY(uv10[uvi[itm],:,:])
    elif itm=='tc':
      cwb+=273.
      var[0,:,:]=interpXY(t2[:,:])
    else:
      var[0,:,:]=interpXY(cwb[0,:,:])

修改後

用dict的update來替代if block的特例狀況
變數名稱與維度全改成字串(surf_var)，用exec來執行
cmd的內容有2個itm，意義分別是字串itm與變數itm。

...
base={i:0 for i in vdic} #used for temp, degC to degK
base.update({'tc':273})

surf_var={i+'a':'uv10[uvi[itm],' for i in 'uv'} #surface wind
surf_var.update({'tc':'t2[','rh':'cwb[0,' for i in vdic}) #temp and rh
...
    cmd='cwb=getvar(wrfin, "'+itm+'",timeidx=tcwb)+base[itm]'
    ...
    exec('var[0,:,:]=interpXY('+surf_var[itm]+':,:])') #surface mapping

三層套疊版本

個別網格系統可以執行後，遭遇邊界條件解析度不足的問題、45公里解析度結果直接引用到3公里解析度的模擬，很難模擬跨境傳輸的現象。
中間增加一層9公里解析度網格，模擬CWBWRF_3Km範圍(南中國與台灣)，再將結果做為台灣範圍的邊界條件，似乎為合理可行的方案。各層範圍詳見下表：

網格層次	中央氣象局WRF之引用	引用方式	CMAQ解析度	用途說明
d01	CWBWRF_15Km範圍(d01)	取最近值(mk_metoaD12.py)	45Km	中亞到日本、印尼到蒙古範圍，與CAMS預報解析度接近
d02	CWBWRF_3Km範圍(d03)	取最近值(mk_metoaD12.py)	9Km	南中國～台灣之中間層
d03	CWBWRF_3Km範圍	線性內插值(mk_metoaTD3.py)	3Km	環保署公版模式範圍。

策略考量
1. d01主要是引用CAMS的預報結果，d02是中間過渡，d03才是真正需要仔細預報空品的範圍。
2. 因此，CWBWRF結果的引用，前二者使用較粗略但快速的方式，直接取最近點的模擬結果進行FDDA，
3. 第三層則使用內插值，一方面有較高的精度，一方面網格數相對(前述東海台海方案)較少，計算時間還可以省一些(還是必須平行計算，作業化的wait程式還需進一步撰寫)。

WPS/geogrid 之網格設定

$ cat /nas1/WRF4.0/WPS/202208_CWB3nests/namelist.wps
&share
 wrf_core = 'ARW',
 max_dom = 3,
 start_date = '2022-08-10_00:00:00','2022-08-10_00:00:00','2022-08-10_00:00:00','2020-06-15_00:00:00'
 end_date   = '2022-08-17_00:00:00','2022-08-17_00:00:00','2022-08-17_00:00:00','2020-08-05_00:00:00'
 interval_seconds = 10800
 io_form_geogrid = 2,
/

&geogrid
 parent_id         =   1,    1,    2,  3,
 Parent_grid_ratio =   1,    5,   3,
 i_parent_start    =   1,   90,  98,
 j_parent_start    =   1,   43,  71,
 e_we              = 221,  206, 103,
 e_sn              = 129,  206, 142,

相對而言d02的格點數有點多，
1. 這是因為原本CWBWRF_3Km的網格點數就比CWBWRF_15Km多，此處僅將其解析度減少3分之1，但仍能保持原CWB模擬範圍，儘量使用所有資訊。
2. 中間過渡何不考慮前述東海台海方案：原因是網格數太少(<40)，WRF及CMAQ都拒絕執行。
而d03原本就不是一個wrf模擬的合理範圍(將大陸的地形排除在模擬範圍之外)，僅是為CMAQ模擬所需的設定。
i(j)_parent_start：d02可以約略估算、d03則需試誤確定。
1. geogrid.exe結果之VERDI確認
2. 執行wrf ~ mcip以致確認GRIDDESC是正確的

最近點對照表之建立mk_nearst_csv.py

參考CAMS預報數據寫成CMAQ邊界檔的作法，找到CMAQ網格點最近的CWBWRF點，並在時間迴圈之外將其紀錄下來，寫成csv檔案備用，不必每次計算。
2個水平面網格點的總數雖然很多、彼此間的距離(點數之乘積)矩陣更大，然而使用python的numpy模組可以輕鬆解決。未來進一步使用pandas資料表代入矩陣的線性化，速度更快，因此不必擔心龐大矩陣的問題。
輸入檔
1. 中央氣象局WRF模版：CWB_wrfout_d0'+d，d=1,2，需連結到日常get_M-A0064.cs作業之成果。注意層數定義的對照(前表)。
2. 目標模版：fname='geo_em.d0'+d+'.nc'，d=1,2。需連結到前述WPS/geogrid.exe結果。
輸出檔：'nearstD'+d+'.csv'，d=1,2。
形成中央氣象局WRF及此處metoa_em檔等2個水平面網格系統：
- 因同時需2個系統的參數，程式設計成交錯讀取處理
- 0為CWB舊網格系統、1為目標網格系統，二者同時以後者的定義為基準
- 產生2維直角座標值後，將其壓縮成為1維方便後續使用

  fname='CWB_wrfout_d0'+d
  wrfin = Dataset(fname,'r')
  latm,lonm=getvar(wrfin,'lat'),getvar(wrfin,'lon')
  ny0,nx0=latm.shape
  fname='geo_em.d0'+d+'.nc'
  metin = Dataset(fname,'r')
  pnyc = Proj(proj='lcc', datum='NAD83', lat_1=metin.TRUELAT1, lat_2=metin.TRUELAT2, lat_0=metin.CEN_LAT, lon_0=metin.CEN_LON, x_0=0, y_0=0)
  x0,y0=pnyc(lonm,latm, inverse=False) #CWB net
  nz,ny,nx=metin['HGT_M'].shape
  x1d=[metin.DX*(i+0.5) for i in range(-nx//2,nx//2)]
  y1d=[metin.DY*(i+0.5) for i in range(-ny//2,ny//2)]
  x1,y1=np.meshgrid(x1d, y1d)
  for i in 'xy':
    for j in '01':
      exec(i+j+'='+i+j+'.flatten()')

計算最近點、並記錄其標籤(1維)，這段會需要一些時間。

  n=[]
  for i in range(ny*nx):
    dist=(x0-x1[i])**2+(y0-y1[i])**2
    idx=np.where(dist==np.min(dist))[0][0]
    n.append(idx)

轉成2維引數、寫成資料表並記錄下來備用。

df=DataFrame({'num':[i for i in range(ny*nx)],'J0':[n[i]//nx0 for i in range(ny*nx)], 'I0':[n[i]%nx0 for i in range(ny*nx)]})
  df['J1']=df.num//nx
  df['I1']=df.num%nx
  df.set_index('num').to_csv('nearstD'+d+'.csv')

mk_metoaD12.py

這個版本就是將前述interpXY副程式，改寫成最近位置標籤的引用(副程式pick_nearst)。

def pick_nearst(oldarr):
  n=oldarr.ndim
  shp=oldarr.shape
  if n==2:
    return oldarr[df.J0,df.I0].reshape(ny,nx)
  if n==3:
    k,j,i=oldarr.shape
    return oldarr[:,df.J0,df.I0].reshape(k,ny,nx)

因getvar結果不是2維就是3維矩陣(給定時間項)，因此只需針對此2種rank形態進行套用。
3維狀況的垂直範圍可能不同，需針對個別矩陣解析(k，不能用nz)、妥為因應。
np.array 1維之選用非常快速，但要確認進來的確是矩陣而不只是getvar的結果(型態為一方法，需以getvar().data形態才是矩陣)

mk_metoaTD3.py

這個版本與前述mk_metoaT.py一樣，只將met_em.d01改成了met_em.d01：

kuang@DEVP /nas1/WRF4.0/WRFv4.2/202208
$ diff mk_metoaTD3.py ../202208TWEPA3k/mk_metoaT.py
44c44
< fname='met_em.d03.'+fn_dt[0]+'.nc'
---
> fname='met_em.d01.'+fn_dt[0]+'.nc'
73c73
<   fname='met_em.d03.'+fn_dt[t]+'.nc'
---
>   fname='met_em.d01.'+fn_dt[t]+'.nc'

結果檢討

東海台海模擬範圍之檢討

含蓋了CMAQ模擬的範圍
風速性質可以反映CWBWRF_3Km的特性


東海台海模擬範圍之地面風氣流線與風速

執行mcip結果檢討

雖然前述網格系統定義並不是以台灣島為中心，mcip結果卻不會因此而有偏移。
邊界線確實非常靠近台灣島(或澎湖)，未來可以預期邊界濃度將大幅影響模式模擬結果。


CMAQ TWEPA公版模式模擬範圍的地形高程

mk_metoa.py程式下載

Download: 將CWB數據填入WRF客觀分析場之程式：mk_metoa.py

[2022-08-11-wrf_pythonTAB]: https://sinotec2.github.io/FAQ/2022/08/11/wrf_pythonTAB.html “wrf-python [//end]: # “Autogenerated link references”