CWB_WRF grib數據檔的時間內插

背景

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

grib與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

• sinotec2(2022) 中央氣象局WRF_3Km數值預報產品 at Github/Focus-on-Air-Quality
• Jeff Whitaker(2021), pygrib官網 at pypi.org
• jiangleads(2018)， pygrib学习 @博客園，posted @ 2018-09-26 12:04.
• Jeff Whitaker(2021), Matplotlib/Basemap and pygrib example @nbviewer