CWB_WRF grib數據檔的時間內插

 

背景

  • CWB_WRF數據偶會發生檔案喪失(data missing)的情形。原因不明,或許在下載時正好遇到電腦正在更新檔案,還是什麼環節出問題,總之,會有需要進行補遺的情況。
  • 幸好大多喪失情況是缺少某一個特定時間的預報檔案。經由前後時間的內插即可解決此一問題。
  • 雖然CWB_WRF數據並沒有統一的維度、雨量也是累加值,所幸前後時間的平均尚能符合定義。
  • 一般使用pygrib來進行讀取grib檔案,但不能用在寫出檔案,pygrib.open也沒有控制讀寫的選項,一律唯讀。必須另外寫出binary檔,這部分可以參考網友jiangleads的範例。
  • 由於每個grib檔案的內涵多有差異,本程式無法適用所有的grib格式檔,還是需一一檢視。
  • grib檔案的時間標籤可以使用[[EAC4_Times]]顯示。

grib與netCDF格式之比較

項目 grib netCDF 說明
版本 2 4  
字元 binary binary  
主要應用 氣象模式 氣象、大氣化學、</br>地球科學  
順序 循序讀寫 階層讀寫  
檔頭  
python模組 pygrib netCDF4  
讀寫 唯讀 可讀寫 前者另開binary檔案
官網 WMO Guide www.unidata.ucar.edu</br>/software/netcdf  

interp_grbs.py程式說明

IO

  • 引數:前後正常下載之檔案名稱(間隔12小時),喪失檔案為其中間時間
    • interp_grbs.py M-A0064-054.grb2 M-A0064-066.grb2,54與66為正常下載(間隔12小時) 之時間,中間的60則為喪失檔時間。
  • 結果:M-A0064-060.grb2

    重要參數

    1. grbs[0].messages:grib可以視為傳統的binary檔案,也是循序一批批寫入的檔案,每一批grib的習慣稱之為message。因此有總批次筆數,就是messages。(相對netCDF為階層式讀寫檔案,netCDF從來不需要、也不知道有rewind,grib就需要rewind)
    2. range(1,grbs[0].messages+1):每一批message的批號數是其內容之一,批號自1開始編號,一直到最後再加1(python習慣)。
    3. grbs[0].message(i):呼叫(抽取)批號i的message,不是序列也不是字典,而是函數關係。
    4. grb.analDate:grib檔案的分析時間(模式的起始時間)。這個時間在不同的預報輸出檔內都是一樣,在於區分預報執行的批次。因此不需要變動。
    5. grb.validDate:這個時間標籤不能修改,似乎是grb方法自行計算出來的,從grb.analDate以及以下要介紹的grb['forecastTime']計算而得。
    6. grb['dataDate']:grb是方法,也是個dict,有幾個索引(keys),然而與我們此次需要修改的項目只有時間。dataDate是8碼的整數,似乎沒有其他作用,後續程式也不會用到,還是參照網友範例將其修正。
    7. grb['forecastTime']:預報時間,單位為小時,非負整數。這項需要修改。
    8. grb.values[:]:此一方法的結果是個numpy.array。類似方法grb.data()的結果是個tuple。同樣是取值為什麼需要2種方法,不是很瞭解,只知道grb.data()目前還沒有看到實際的範例就是了。
    9. grb.tostring():此方法結果是個很長的字串,類似dump的作用,也是grib檔案message的意涵。將其按照批次順序寫進binary檔案就好了。

批次內插並寫進檔案之迴圈

grbout = open(fnames[0].replace(fcst_hrA[0],s),'wb')

for i in range(1,grbs[0].messages+1):
  grb=grbs[0].message(i)
  s=(grb.validDate+datetime.timedelta(hours=6)).strftime("%Y%m%d")
  grb['dataDate']=int(s)
  grb['forecastTime']=fcst_hr
  grb.values[:]=(grbs[0][i].values+grbs[1][i].values)/2.
  msg = grb.tostring()
  grbout.write(msg)

程式下載

Reference