Link Search Menu Expand Document
Focus on Air Quality

CAMS預報數據寫成CMAQ邊界檔

Table of contents

背景

  • 過去執行中期天氣預報中心再分析空品數據寫成CMAQ邊界檔有2個作法
    1. 先將全月數據寫成m3nc檔案，在用CMAQ系統的bcon程式進行切割。由於全月m3nc檔案非常龐大，此一作法不單浪費空間，同時也耗費計算時間。
    2. 直接將grib檔案寫bcon檔案。此法也需處理龐大的grib檔案，效率不彰。
  • 此處除了針對ecmwf數值預報數據檔與再分析數據檔的差異，改進EAC4檔案轉成5階m3.nc作業方式，也需要改進前述作業的作法，提升計算效率，以運用在日常作業中。
  • 改進構想
    1. 是否直接讀取grib檔案，不再經轉換成nc檔案？這題在CAMS_ic內詳細討論，此處就不贅述了。對邊界檔而言，需處理41個time step，如果是要循序讀取grib檔案、將每個2維陣列放在正確的5度空間矩陣，那將會是個災難。
    2. 由於是廣大範圍、遙遠位置的邊界濃度，對臺灣地區模擬結果的影響不是很大，過去詳細進行網格內插(scipy.griddata)似乎沒有必要。改成最近點直接代替，會比較有效率，同時也可以避免找不到適合進行內插的值而發生nan的結果。
    3. 時間內插：如果CMAQ可以接受解析度3小時的邊界條件就最方便了。但這項因素似乎是模式整體控制，不能個別差異。如果要改，連mcip、排放等都要系統性更換。
  • 除了此處之邊界檔案之外，CMAQ還需要初始檔，可以進行grb2icon.py，與邊界檔作業沒有前後關連性。

建立CMAQ邊界點位置與CAMS網格系統的對照關係(BconInGrb.py)

  • 目標是產生BconInGrb.csv檔案，將每一邊界點最靠近的CAMS經緯度位置記錄下來，供後續使用。
  • 因為對照關係不會因為時間日期而異，因此可以從日常的下載及處理程序中獨立出來，只做一次。待網格系統有異動時再執行即可。

CMAQ邊界位置的座標

  • CMAQ的邊界定義有點特別，在水平面上共有(nrow+ncol)*2+4點，自模式外圍西南角落開始向東、向北、向西、再向南完成一圈。詳細的程式碼：
nt1,nlay1,nbnd1=nc1.variables[V1[2][0]].shape
nrow0,ncol0=nc1.NROWS+5,nc1.NCOLS+5
x1d=[nc1.XORIG+nc1.XCELL*(i-2) for i in range(ncol0)]
y1d=[nc1.YORIG+nc1.YCELL*(i-2) for i in range(nrow0)]
x1,y1=np.meshgrid(x1d, y1d)
i0,j0=1,1
i1,j1=i0+ncol1+1,j0+nrow1+1
idx=[(j0,i) for i in range(i0+1,i1+1)]  +   [(j,i1) for j in range(j0+1,j1+1)] + \
    [(j1,i) for i in range(i1-1,i0-1,-1)] + [(j,i0) for j in range(j1-1,j0-1,-1)]
idxo=np.array(idx,dtype=int).flatten().reshape(nbnd1,2).T
x1,y1=x1[idxo[0],idxo[1]],y1[idxo[0],idxo[1]]
  • 邊界比正常模式模擬範圍還要大1圈，x1d,y1d則設成比模擬範圍多2圈、nrow0,ncol0多5，以擴大內插範圍(沿用grb2bc的作法，雖然此處不再進行內插)。
  • 有了目標網格系統的座標，還需進行CAMS經緯度座標的轉換

CAMS經緯度座標的轉換

  • 應用Proj進行座標轉換。參數則直接引用bcon模版，不另外內設以避免錯誤。
  • 經緯度在grib檔與nc檔案裏都有，此處版本使用grib檔，經緯度是任一grb的屬性，直接呼叫比較方便。
  • 注意：此處並不翻轉南北方向之y軸，未來應用時也不必翻轉。
...
pnyc = Proj(proj='lcc', datum='NAD83', lat_1=nc1.P_ALP, lat_2=nc1.P_BET, lat_0=nc1.YCENT, lon_0=nc1.XCENT, x_0=0, y_0=0.0)
...
#grib file
fname='allEA_1.grib'
grbs = pygrib.open(fname)
grb=grbs[1]
lats, lons = grb.latlons()
x,y=pnyc(lons,lats, inverse=False) #no flip

計算距離並選出最近距離

  • 先產生座標系統的JI值(seq)。seq, x, y都壓縮成1維備用。
ny,nx=x.shape
seq=np.zeros(shape=(ny,nx),dtype=int)
for j in range(ny):
    for i in range(nx):
        seq[j,i]=j*1000+i
for i in ['x','y','seq']:
    exec(i+'='+i+'.flatten()')
  • 計算舊網格到邊界上的距離，找到最近點，將其記錄下來，以對照前述seq及x,y。
n=[-1 for i in range(nbnd1)]
for i in range(nbnd1):
    dist=(x-x1[i])**2+(y-y1[i])**2      #nearest grib data for bcon
    idx=np.where(dist==np.min(dist))[0]
    if type(idx)==list and len(idx)>1: idx=idx[0]
    n[i]=idx

儲存結果

  • 以DataFrame形式儲存結果，方便檢查與日後讀取。
df=DataFrame({'NumOfBnd':[i for i in range(nbnd1)],'x1':x1,'y1':y1,'JinBCON':idxo[0],'IinBcon':idxo[1],'JIseqInGrb':[seq[n[i]][0] for i in range(nbnd1)]})
df.set_index('NumOfBnd').to_csv('BconInGrb.csv')
  • 由於df是正確順序的資料表，未來引用JIseqInGrb即可連到grib的數據。

grb2bcon.py

  • 首要任務是將3批下載的grib檔案進行整併，存成5維度的矩陣備用。
    • 不在作業環境使用ncks進行橫向污染物質整併的理由
      1. 檔案會非常大，沒有必要儲存。
      2. 累加的作業方式也不利平行化處理。
  • 讀進前述位置對照關係(BconInGrb.csv)之濾鏡，經此濾鏡可以快速將5維度矩陣降階為4度(水平x,y消失，簡化為邊界線形維度)。
    • 前述5維度矩陣之記憶體可以覆蓋、不再使用
  • 空間問題解決後，即可在污染物維度進行展開，由27項展開為39項。
  • 最後進行時間的內插，將CAMS 3小時之架構線性內插展開為逐時檔。時間由41項增加為121。存入舊檔模版。

grb2icon.py程式IO

  • AllEA_1.nc ~ AllEA_3.nc：3批次下載之grib檔案，經ncl_convert2nc轉換之結果
  • 與grb2icon.py相同
    1. BCON_v53_1912_run5_regrid_20191201_TWN_3X3：取粒狀物IJK比例(rate)之模版
    2. METCRO3D_2208_run8.nc：讀取空氣密度進行粒狀物的單位轉換
  • BconInGrb.csv：座標點對照關係(前述)
  • BCON_v53_2208_run8_regrid_20220811_CWBWRF_45k：模版及結果檔

grb2bcon.py流程

graph TD
    A((AllEA_?.nc)) --> B(bdate, 5d array var)
    C((BCON_old)) --> D(rate)
    J((BconInGrb.csv)) --> K{horizontal mirroring}
    E((METCRO3D)) --> K
    G((BCON_template)) --> I{time flags}
    B --> I
    B --> K
    K --> L(4d new array varbc, dens1d)
    L --> M{spec expansion}
    D --> M
    M --> Q(varv)
    Q --> O{time interpolation}
    I --> O
    O --> H(BCON_new)
    H --> P((save and quit))

污染項目對照與展開

  • 此部分之處理方式與CAMS_ic相同，僅輸入及輸出的矩陣因應調整
#iv in (gas+par) transfer to V1[2]
for v in list(nms_gas)+list(nms_part):
  print(v)
  iv=(gas+par).index(dic[v])
  if v in nms_gas:
    nm=nms_gas[v]
    if nm not in V1[2]:continue
    varv[V1[2].index(nm),:,:,:]=varbc[iv,:,:,:]*28.E6/mws[dic[v]] #mixing ratio to ppm
    continue
  skip=0
  nms=nms_part[v]
  for nm in nms:
    if nm not in V1[2]:skip=1
  if skip==1:continue
#    unit=dens[:] (kg/kg) to microgram/M3
  for nm in nms:
    im=nms.index(nm)
    varv[V1[2].index(nm),:,:,:]+=varbc[iv,:,:,:]*rate[v][im]*dens1d[None,:,:]

結果檢討

2023/06更新

程式更新

  • 因應grib_to_netcdf轉檔結果，時間標籤改在nc['time']，單位是1900-01-01 00:00:00.0之後的小時數。其餘沒有更動。
  • 對照表的內容更新了，但檔名沒有更動。
kuang@dev2 /u01/ecmwf/CAMS/CAMS_global_atmospheric_composition_forecasts/2022
$ diff grb2bcon.py grb2bconV2.py
44,45c44,47
<     fcst_hr=float(np.array(nc['forecast_time0'])[0])
<     bdatef=datetime.datetime.strptime(nc.variables[V[3][0]].initial_time,'%m/%d/%Y (%H:%M)')+datetime.timedelta(hours=fcst_hr)
---
>     v='time'
>     udate=datetime.datetime.strptime(' '.join(nc[v].units.split()[2:4]),'%Y-%m-%d %H:%M:%S.0')
>     fcst_hr=float(np.array(nc[v])[0])
>     bdatef=udate+datetime.timedelta(hours=fcst_hr)