grb2wrfout_d04轉檔

Table of contents

grb2wrfout_d04轉檔

背景

續樓上、以及下載程序之說明，此處詳述轉檔歷程。
由grb2轉到wrfout_d04的空間內插問題，詳見WRF_3Km空間內插說明。
雖然grb2格式也有其解讀、應用的軟體，然而在空氣污染領域還並不多。因此還是需要轉成wrfout的nc格式。
此處應用pygrib模組進行grb2檔案的解析

rd_grbCubicA.py分段說明

宣告與基本設定

引用模組

pygrib的建置比較特別，可以參考pygrib的安裝、重要語法

空間內插之標籤與權重檔案是以FortranFile儲存的，時間內插使用到CubicSpline，垂直空間內插使用到interp1d。

kuang@MiniWei /Users/Data/cwb/WRF_3Km
$ cat -n rd_grbCubicA.py
#!/opt/anaconda3/envs/gribby/bin/python
import pygrib
import numpy as np
import netCDF4
from scipy.io import FortranFile
from scipy.interpolate import CubicSpline
from scipy.interpolate import interp1d
import datetime
import os
import subprocess
11

溫度轉飽和蒸汽壓Arden Buck equation
- grb2及wrfout的溫度單位都是K，但Arden Buck equation溫度單位為C。
```
12  def buck(K):
13      C=K-273.
14      return 611.21*np.exp((18.678-C/234.5)*C/(257.14+C))
15
16
```

開啟模版、讀取d04的網格數

由於wrfout的變數較grb2檔案還多，需事先處理好(ncks -O -x -v VAR1,VAR2... $oldNC $newNC)。

當然做為模版，也要有正確的時間、空間結構，可以使用ncks -d指令裁切既有檔案。

fname='wrfout_d04'
tmax=36+1
nc = netCDF4.Dataset(fname,'r+')
Vs=[list(filter(lambda x:nc.variables[x].ndim==j, [i for i in nc.variables])) for j in [1,2,3,4]]
nt1,nrow1,ncol1=(nc.variables['U10'].shape[i] for i in range(3))
22

讀取空間內插之標籤檔案idxD4.bin、權重矩陣檔案wtsD4.bin

因2個網格系統的位置不會隨時間改變，可以事先處理好備用。

#path='/Users/kuang/MyPrograms/UNRESPForecastingSystem/CWB_data/raw/'
#path='/home/cpuff/UNRESPForecastingSystem/CWB_data/raw/'
#path='/nas1/backup/data/CWB_data/raw/'
path='/Users/Data/cwb/WRF_3Km/'
f=FortranFile(path+'idxD4.bin', 'r')
idx=f.read_record(dtype=np.int64)
idx=idx.reshape(nrow1,ncol1,2)
30
f=FortranFile(path+'wtsD4.bin', 'r')
wts=f.read_record(dtype=np.float64)
wts=wts.reshape(nrow1,ncol1,4)
34

時間迴圈前之準備

atbs:變數名稱對照表，短名稱為wrfout內的名稱，長名稱為grb2內的變數屬性

one=np.ones(shape=(nrow1,ncol1),dtype=np.int64)
atbs={'U10': '10 metre U wind component',
 'V10': '10 metre V wind component',
 'T': 'Temperature',
 'V': 'V component of wind',
 'U': 'U component of wind',
 'W': 'Geometric vertical velocity',
 'QVAPOR': 'Relative humidity',
 'T2': '2 metre temperature',
 'SST': 'Temperature',
 'TSK': 'Skin temperature',
 'SWNORM': 'Net short-wave radiation flux (surface)',
 'SWDOWN': 'Net short-wave radiation flux (surface)',
 'PHB': 'Geopotential Height',
 'PSFC': 'Surface pressure',
 'RAINC':'Total precipitation'}

atbs2:atbs中之3維變數

atbs2=[]
for v in atbs:
  if v in Vs[2]:atbs2.append(v)
#'U10','V10','T2','SWNORM','SWDOWN','PSFC','SST','TSK']
#atbs2.append('TSLB') #3d in grb but 4d in wrfout

開啟變數陣列

b為時間標籤
每個atbs項目都有自己的陣列(s開頭的全時間變數、有別於沒有s為某時間的變數)

grb2的高度順序是相反的，lv為其標籤。

b=[t for t in range(0,tmax)]
for a in atbs:
  if a in atbs2:
    exec('s'+a+'=np.zeros(shape=(tmax,nrow1,ncol1),dtype=np.float64)')
  else:
    exec('s'+a+'=np.zeros(shape=(tmax,11,nrow1,ncol1),dtype=np.float64)')
 #'GPH': 'Geopotential Height',
63
lv=[i for i in range(10,-1,-1)]
PB,level=[],[]
66

開啟檔案、讀取變量

開始時間迴圈、開啟grb2檔案

for t in range(0,tmax,6):
  fname='M-A0064-0'+'{:02d}'.format(t)+'.grb2'
  grbs = pygrib.open(fname)
  70

除累積雨量外、讀取每一grb2變數

len(grb)==1：只有一層數據
有11層數據者：要記得高度須反向儲存

海溫SST設定為地面溫度

  for a in set(atbs)-set(['RAINC']):
    grb = grbs.select(name=atbs[a])
    if len(grb)==1:
      cmd=a+'=grb[0].values'
      exec(cmd)
    else:
      arr=[]
      for k in range(11):
        arr.append(grb[lv[k]].values)
      arr=np.array(arr)
      exec(a+'=arr')
    if len(grb)==11 and len(level)==0:
      level=[grb[lv[i]].level for i in range(11)]
    if len(grb)==12 and a=='SST':SST=grb[11].values

由gr2的相對濕度計算wrfout中的水汽量

  nlay,nrow,ncol=QVAPOR.shape
  if len(PB)==0:
    PB=np.ones(shape=(nlay,nrow,ncol),dtype=np.int64)
    for k in range(nlay):
      PB[k,:,:]=level[k]*100
  ps=buck(T)
  QVAPOR=QVAPOR/100.*(ps*18.)/(PB*28.)
  92

累積雨量之計算

grb2是逐6小時累積，wrfout是逐時累積

  a='RAINC'
  if t==0:rain_old=np.zeros(shape=(nrow,ncol))
  t6=min([t+6,84])
  fname='M-A0064-0'+'{:02d}'.format(t6)+'.grb2'
  try:
    rain_acc=pygrib.open(fname)[62].values
  except:
    rain_acc=rain_old
  arr=(rain_acc-rain_old)/6. #average along 6 hrs
  exec(a+'=arr')
  rain_old[:]=rain_acc[:]
 104

空間內插

空間內插，並儲存在全時間變數陣列中

因2個檔案的解析度一樣，簡單的線性內插即可，然因原點與投影略有差異，還是有內插模式，詳見WRF_3Km空間內插說明。

2度空間變數與3度空間變數分開處理，後者逐層處理

  for a in atbs:
    exec('var='+a)
    n=var.ndim
    if n==2:
      kk=0
      for jj in [0,1]:
        for ii in [0,1]:
          vr=var[idx[:,:,0]+one*jj,idx[:,:,1]+one*ii]
          exec('s'+a+'[t,:,:]+=vr[:,:]*wts[:,:,kk]')
          kk+=1
    elif n==3:
      for k in range(11):
        kk=0
        for jj in [0,1]:
          for ii in [0,1]:
            vr=var[k,idx[:,:,0]+one*jj,idx[:,:,1]+one*ii]
            exec('s'+a+'[t,k,:,:]+=vr[:,:]*wts[:,:,kk]')
            kk+=1

時間內插

取得分析時間的時間標籤

校正到UTC 0600

#time stamps
  if t==0:
    V=grbs[1]
    beg_time=V.analDate
    if beg_time.hour != 6:
      beg_time=beg_time-datetime.timedelta(days=beg_time.hour/24)+datetime.timedelta(days=6/24)
 129

結束時間迴圈、進行時間內插

產生wrfout逐時的時間標籤

X=[i for i in range(0,tmax,6)]
Hstart=14 #begin local time(20 for CAMx, 14 for fitting WRF GMT 0600)
for t in range(0,tmax):
  time=beg_time+datetime.timedelta(days=t/24.)
  b[t]=np.array([bytes(i,encoding='utf-8') for i in time.strftime("%Y-%m-%d_%H:%M:%S")])
wname=''
for i in b[0]:
  wname+=i.decode('utf-8')
 138

3維變數的時間內插

TSLB在grb2是3度，在wrfout是4度

for a in atbs2:
  exec('ss=np.array([s'+a+'[t,:,:]   for t in range(0,tmax,6)])')
  cs = CubicSpline(X,ss)
  for t in range(0,tmax):
    if t%6==0:continue
    exec('s'+a+'[t,:,:]=cs(t)')
  if a == 'TSLB':
    exec('nc.variables["'+a+'"][:,0,:nrow1,:ncol1]=s'+a+'[Hstart-14:Hstart+10,:,:]')
  else:
    exec('nc.variables["'+a+'"][:,:,:]=s'+a+'[Hstart-14:Hstart+10,:,:]')

4維變數的時間內插

同前TSLB之計算

for a in set(atbs)-set(atbs2):
  exec('ss=np.array([s'+a+'[t,:,:,:] for t in range(0,tmax,6)])')
  cs = CubicSpline(X,ss)
  for t in range(0,tmax):
    if t%6==0:continue
    exec('s'+a+'[t,:,:,:]=cs(t)')
  exec('nc.variables["'+a+'"][:,:nlay,:nrow1,:ncol1]=s'+a+'[Hstart-14:Hstart+10,:,:,:]')
 157

個別變數再加工

壓力的計算

wrfout的壓力有基準量、擾動量之分，grb2只有固定層壓力

#pressure
v='PB'
nt,nlay,nrow,ncol=nc.variables[v].shape
for k in range(nlay):
  nc.variables[v][:,k,:,:]=level[k]*100.
nc.variables['P'][:]=0.
 164

重力位高

同樣wrfout的的高度也有基準量、擾動量之分。前者令為所有時間的平均值，後者即為差值。

地面設為地形高度

#geopotential heights
v='PHB'
nc.variables['PH'][:]=0.
PHB=nc.variables[v][:]*9.8
PHBm=np.mean(PHB,axis=0) #const in time
PH=PHB[:,:,:,:]-PHBm[None,:,:,:]
PHBm=PHB-PH
#surface geopotential height=HGT*9.8
HGT=nc.variables['HGT'][:,:,:]*9.8
for k in range(nlay):
  nc.variables['PH'][:,k+1,:,:]=PH[:,k,:,:]
nc.variables[v][:,0,:,:]=HGT
for k in range(nlay):
  nc.variables[v][:,k+1,:,:]=PHBm[:,k,:,:]+HGT[0,:,:]
 179

位溫

溫度基準為定值

#potential temperature
v='T' #WRF is temp "pertubation"
PB,TK=nc.variables['PB'][:],nc.variables[v][:]
nc.variables[v][:]=TK*(100000./PB)**0.286
nc.variables[v][:]-=nc.variables['T00'][0]
 185

表面溫度

令為T2

 186  #T2->TSK
 187  nc.variables['TSK'][:]=nc.variables['T2'][:]
 188

逐時雨量累加成為累積雨量

v='RAINC'
for t in range(1,24):
  nc.variables[v][t,:,:]+=nc.variables[v][t-1,:,:]
 192

行星邊界層高度

求解行星邊界層高度。由於grb2沒有相關任何訊息，此處以混合層高度代之。

先求出dt/dz值(使用np.diff指令)

v='PBLH'
#PBLH=nc.variables[v][:]
#for t in range(24):
#  t_wrf=(t+6)%24
#  nc.variables[v][t,:,:]=PBLH[t_wrf,:,:]
T0=nc.variables['T2'][:]*(100000./nc.variables['PSFC'][:])**0.286
T=nc.variables['T'][:]+290
H=(nc.variables['PH'][:]+nc.variables['PHB'][:])/9.8
nt,nlay,nrow,ncol=(H.shape[i] for i in range(4))
nk=nlay-1
HT=np.zeros(shape=T.shape)
for k in range(nk):
  HT[:,k,:,:]=(H[:,k,:,:]+H[:,k+1,:,:])/2.-H[:,0,:,:]
 206

地面溫度少於其他高度者：混合層限制在最低高度

mixH=np.zeros(shape=T0.shape)
minT=np.min(T,axis=1)
idx=np.where(T0<=minT)
HT0=np.array([max(35,i) for i in HT[:,0,:,:].flatten()]).reshape(HT[:,0,:,:].shape)
mixH[idx[:]]=HT0[idx[:]]
 212

無逆溫者使用interp1d指令內插到符合地面溫度的高度

#potential temperature vertical profile is increasingly sorted
dT=np.diff(T,axis=1)
dTmin=np.min(dT,axis=1)
idx=np.where((dTmin>=0)&(T0>minT))
for n in range(len(idx[0])):
  t,j,i=(idx[k][n] for k in range(3))
  f=interp1d(T[t,:,j,i],HT[t,:,j,i])
  mixH[t,j,i]=f(T0[t,j,i])
 221

有逆溫者找到逆溫層高度

#with inversions
idx=np.where((dTmin<0)&(T0>minT))
for n in range(len(idx[0])):
  t,j,i=(idx[k][n] for k in range(3))
  k=list(dT[t,:,j,i]).index(dTmin[t,j,i])+1
  mixH[t,j,i]=HT[t,k,j,i]
nc.variables[v][:]=mixH[:]
 229

時間標籤與屬性資料之寫入

寫入時間標籤

修改wrfout檔案屬性訊息

修改開始時間
修改垂直層數：因變數BOTTOM-TOP_GRID_DIMENSION含減號-，不能在python中處理，必須用ncatted指令在OS處理。

將wrfout_d04更名

v='Times'
nc.variables[v][:,:]=b[Hstart-14:Hstart+10][:]
btime= beg_time.strftime("%Y-%m-%d_%H:%M:%S")
nc.SIMULATION_START_DATE=btime
nc.START_DATE           =btime
yr=beg_time.year
nc.JULYR                =yr
nc.JULDAY               =int((beg_time-datetime.datetime(yr,1,1)).total_seconds()/24/3600.)+1
nc.close()
ncatted='/usr/local/bin/ncatted'#subprocess.check_output('which ncatted',shell=True).decode('utf8').strip('\n')
os.system(ncatted+' -a BOTTOM-TOP_GRID_DIMENSION,global,o,i,12 wrfout_d04')
os.system('mv wrfout_d04 wrfout_d04_'+wname)

下載程式碼

可以由github找到原始碼。

Download: rd_grbCubicA.py

檢核

可以使用MeteoInfo或CWB網站

Reference

sinotec2, pygrib的安裝、重要語法, evernote, 2021年4月1日
Yaqiang Wang, MeteoInfo Introduction, meteothink, 2021,10,16