m3nc檔案轉GIF

Table of contents

背景

  • USEPA Model3 的nc檔案格式(IOAPI)雖然可以用VERDI或其他軟體打開、檢視,但還是不夠簡潔,尤其在公版模式、自動執行、遠端計算服務的過程中,如果能在模式計算完、後處理都結束後,自動產出GIF檔案,公布在網站上,會非常有幫助。
  • 畫等值線圖對matploglib不是問題,但因為座標投影的問題,很容易發生位置扭曲的情形,還是需要用經緯度來校正。這也是wrf-python的強項。麻煩的是必須從wrf_out來讀取座標系統的轉換參數。
  • matplotlib等值圖並不是raster圖檔,而是有內插效果的,不會像VERDI一樣鋸齒狀。
  • 圖面上的濃度可能差異好幾個數量級,用線性等級等值線會集中在高值,圖面上太稀疏無資訊。VERDI也有線性與對數2種作法因應。分開作業才是合理的方式。
  • 沿襲環保署提供公版後製工具
    • extend=’max’
    • 在footnote提供圖面上最大值之數字
    • (類)彩虹色階
  • 此處也同時解決了該工具的缺失:
    • 直角座標系統改成蘭伯投影系統
    • 256色階改成10~15階
    • 取消白色~淡藍色之間的漸層
    • 色階改成彈性判斷:線性或對數、階層數及有效位數自動調整(過去grads也是這樣設定的)。

輸入及輸出

引數

  • m3nc檔案名稱
    • 檔名會寫在Title上,雖會去掉路徑,但字數長度還是有限,
    • 最好控制在10個字元之內。

輸入檔案

  • m3nc檔案
  • wrfout_d04
    • 需存在同一工作目錄下
    • 注意水平網格之裁剪,一般wrfout會比MCIP/CCTM的網格範圍還大一些,需按照mcip的設定進行裁切,如
      • ncks -O -d Time,0 -d bottom_top,0 -d bottom_top_stag,0 -d west_east,2,201 -d south_north,2,201 -d west_east_stag,3,202 -d south_north_stag,3,202 $nc wrfout_d04
      • ncks -O -d Time,0 -d bottom_top,0 -d bottom_top_stag,0 -d west_east,7,101 -d west_east_stag,7,102 $nc /nas2/cmaqruns/2022fcst/grid03/cctm.fcst/daily/wrfout_d04

輸出檔案

  • v_NN.png、v為模擬空氣品質項目,NN為00~99的時間序列。
  • v.gif:最後整併結果

png2gif檔案轉換(convert)

  • imageMagics的convert指令是最佳的選擇,只需要將bash中的執行,放在python中直接進行就好了
  • 基本程序:去邊、統一增加5%的邊、整併成gif檔
  • png合併成gif需給定背景與像素尺寸的理由詳見格式轉換
  • 存檔格式png-or-jpg的選擇

程式外進行方式

for s in $(ls *_24.png);do 
  v=$(echo $s|cut -d'_' -f1)
  test $v=='PM25' && v=$(echo $s|cut -d'_' -f2)
  for png in $(ls ${v}*.png);do 
    convert -bordercolor white -trim $png tmp.png
    convert -bordercolor white -border 5%x5% tmp.png $png
  done
  convert  -dispose 2 -coalesce +repage -background none ${v}_*.png -size 607x774 ${v}.gif
done

size的問題

  • png的像素size會隨著圖形大小、邊界裁切等條件而異,應該不要設成定值比較保險
  • 詢問png檔案的size可以用convert '+v+'_00.png -format "%wx%h" info:指令(參考DETECT IMAGE SIZE)
  • python 可以用subprocess.checkout將結果讀成string,再寫成命令列即可確保gif的大小能夠一致。

m3nc2gif.py程式設計重點

座標系統

  • 雖然說一般d4範圍並不大、搭配了行政區的邊界線也沒有造成誤解的可能,經緯度系統似乎沒有必要,
    • 但有經緯度格線看起來確實還是比較專業一點。
    • 這也是必須套用wrf-python的理由。
  • 有座標值的內插可能是讓圖面看起來不再是raster圖,而是真正shaded圖的原因
  • VERDI圖有網格數的XY座標,但還是長度為單位比較直覺。如果可以用公里為單位就更好了。
from pyproj import Proj
from wrf import (getvar, to_np, get_cartopy)
...
pnyc = Proj(proj='lcc', datum='NAD83', lat_1=nc.P_ALP, lat_2=nc.P_BET,lat_0=nc.YCENT, lon_0=nc.XCENT, x_0=0, y_0=0.0)
...
x1d=[nc.XORIG+nc.XCELL*i for i in range(ncol)]
y1d=[nc.YORIG+nc.YCELL*i for i in range(nrow)]
X,Y=np.meshgrid(x1d,y1d)
lons, lats= pnyc(X,Y, inverse=True)
ncfile = Dataset('wrfout_d04')
p = getvar(ncfile, "pressure",timeidx=0)
cart_proj = get_cartopy(p)

...
    ax = plt.axes(projection=cart_proj)
    contours = plt.contourf(lons, lats, nc[v][t,0,:,:] ,
                       levels=level, #colors="rainbow",#"black",
                       norm=nm,
                       cmap=get_cmap("rainbow"),
                       transform=crs.PlateCarree(),
                       extend='max')
...
    ax.set_xlim(cartopy_xlim(p)+np.array([+0,-30000]))
    ax.set_ylim(cartopy_ylim(p)+np.array([+30000,-20000]))
    ax.gridlines()
    ax.set_xticks(xtics,cart_proj) # Grid
    ax.set_yticks(ytics,cart_proj) # Grid

    # Set the x-axis and  y-axis labels
    ax.set_xlabel("meter", fontsize=12)
    ax.set_ylabel("meter", fontsize=12)
  • 圖面範圍:ax.set_xlimax.set_ylim
    • 當模擬範圍較小時,為取整數刻度,matplotlib會將範圍向外多取一些,導致圖面出現空白框。可以使用此2指令強制縮小範圍。
    • cartopy_xlim(p)等為wrf-python指令,結果為下限與上限值之陣列。因此如需增減,在其後加上一陣列即可。
    • 上下限與增減單位為公尺
    • 如果內縮太多,可另行調整。

浮動的濃度等級

  • 如前所述,濃度等級需要幾個項目是浮動的(至少想到5個):
    • 1.levels極值的設定:全時間、全域、正值
      • 99.99 percentile (等級之最大值mxv)
      • 0.01 percentile (等級之最小值mnv)
      • 不服務負值的理由是:
        • 正負值同框出現,不符合一般直覺與作圖用意。
        • log10會出錯
        • 取低值會差很多(0值可能在增量模擬結果會非常多)。
      • 如果負值一定要套用,可以先在nc檔案階段處理,將負值轉正,再應用繪圖程式。
      • 某時間的最大值低於mnvcontourf將會出錯、跳出程式,不會產生gif檔:
        • 發生在夏季午後降雨事件,濃度會在短時間內消失。
        • 如果是前面階段,從未產生過圖檔,只能捨棄這些小時。
        • 如果是中間、後面階段,已經產生過圖檔,還可以沿用上一小時之contourf結果
        • 檢討:因最大值低於mnv是非常極端個案,整體時間變化看起來是低濃度持續一陣子,還算順暢。
    • 2.線性或對數色階,改成彈性判斷
      • 最大/最小比值高於15:非線性、對數色階,適用在煙流增量分布,
      • 比值低於15:線性色階,適用在濃度變化不大的空品項目。
    • 3.濃度等級的階層數(nlev)
      • 對數色階的數字不是等間距,沒有階層數的問題,一律使用15階即可
      • 線性色階間距的數字需要是整數,以便能直接辨識濃度值
      • 此處分為2類,最大值第1碼為[3,6,9]者取15階,其餘(7視為6)取10
    • 4.對數等級對應到顏色的分布(norm=nm)
      • 正常內設狀況下線性的等級對應到均勻的色階顏色,是沒有問題的。
      • 如果是對數濃度等級對應到內設色階分布,結果還是集中在高濃度,並沒有改變。需要改變顏色的分布
      • matplotlib.colors至少有NormalizeLogNormSymLogNormPowerNormBoundaryNormOffsetNorm等分布型態。
      • 此處之nm值應用到前2者。顏色分布的極值(vmin,vmax)調整則以圖面上保持最多顏色而且均勻為原則。
    • 5.同時應用在小數及大數的有效位數 N
      • 濃度太低、太高的情形應考慮轉換單位、重新整理nc檔案,而不是在圖面上用有效位數調整。因此不考慮用科學記號。
      • N=int(3-np.log10(mxv)),大致上可以保持有3~4個有效數字,
      • 有效位數不足:低值部分因去尾造成某些區間的間格為0 (len(level)!=len(set(level))),這會造成程式錯誤。再增加一位有效位數即可。
      • 測試結果如下
mxv範例有效位數N有效數字
12345.12345-112350.
1234.123401234.
123.1230123.
12.12112.1
1.121.10
0.1230.120
0.12330.123
0.0123440.0123
  • 相關程式碼
#判斷線性或對數濃度色階
def get_lev(N):
  if mxv/mnv>15:    
    dc=(np.log10(mxv)-np.log10(mnv))/15    
    level=[round(10**(dc*i+np.log10(mnv)),N) for i in range(15)]
    nm=colors.LogNorm(vmin=level[0], vmax=level[-1])
  else:
    i=int('{:e}'.format(mxv)[0])
    if i==7:i=6
    dc=i*10**int(np.log10(mxv))/nlev[i]
    level=[round(dc*i,N) for i in range(nlev[i])]
    nm=colors.Normalize(vmin=mnv, vmax=mxv)
  return level,nm
#2種線性間隔數
nlev={i:10 for i in [1,2,4,8]}
nlev.update({i:15 for i in [3,6,9]})
...
  #極值只考慮正值範圍
  a=np.where(nc[v][:,0,:,:]>0,nc[v][:,0,:,:],0)  
  mxv=np.percentile(a,99.99)
  mnv=np.max([np.percentile(a,0.01),mxv/100])
  N=int(3-np.log10(mxv))
  level,nm=get_lev(N)    
  if len(level)!=len(set(level)):level,nm=get_lev(N+1)
  #格式必須在時間迴圈外設定好,避免有偏差,GIF會跳動
  fmt='%.'+str(N)+'f'     
...
    plt.colorbar(contours, ax=ax, orientation="vertical",pad=.05,format=fmt)

浮動大小比例的色標(colorbar)

  • 色標大小的設定如果要跟著主要圖面更動,必須使用fraction指令(參考Geek)。
  • 0.047 is a magic number
...
im_ratio=nrow/ncol
...
plt.colorbar(contours, ax=ax, orientation="vertical",pad=.05,format='%.3f',fraction=0.047*im_ratio)

程式內進行png2gif之convert

...
    png=v+'_'+'{:02d}'.format(t)+'.png'
    plt.savefig(png)
    plt.close()
    os.system('convert -bordercolor white -trim '+png+' tmp.png')
    os.system('convert -bordercolor white -border 5%x5% tmp.png '+png)
  if nt<6:continue #too short for GIF    
  size=subprocess.check_output('convert '+v+'_00.png -format "%wx%h" info:',shell=True).decode('utf8').strip('\n')
  os.system('convert  -dispose 2 -coalesce +repage -background none '+v+'_*.png -size '+size+' '+v+'.gif')

結果比較檢討

比較

  • 公版後製工具VERDI、m3nc2gif
  • 2019/01全月NO2模擬最大小時值之濃度分布
p1p2p3
公版後製工具(線性等間距格柵)VERDI(對數等間距格柵)m3nc2gif(不等間距等濃度)

檢討

項目公版後製工具VERDIm3nc2gif說明
顏色分布與解析重點低濃度段中段高濃度段m3nc2gif結果低濃度段解析較差、有顏色之總面積在三者中為最小,雖然模擬是為瞭解高值,但低值卻反映傳輸現象。似仍有檢討空間。
格柵可以放大圖面解析度
等值線平緩曲線
色階解析度低、無法精確判讀可、但需轉換計算高、直接判讀浮動濃度等級奏效、濃度值也不必再轉10的次方
底圖解析度Natural Earth scale.無法再更細, Valid scales are “110m”, “50m”, and “10m”,是否改內政部shape檔?(無大陸部分)
標示極值陸上所有範圍 
  • 其他有待精進項目
    • 左右側範圍
    • XY軸單位
    • 中文Title/subTitle?

降雨的效果

更新色階至RdYlGn(紅綠燈)

緣起

  • 為比較外單位(中央大學大氣系、中研院變遷中心)的模擬成果圖檔,將原來彩虹色階改成紅綠燈色階。
  • 除顏色之外,色階的分層也採一致化,而不是程式給定。

變更項目

$ diff ~/bin/m3nc2gif.py ~/bin/m3nc2gif03.py
22,23c22,24
<     dc=(np.log10(mxv)-np.log10(mnv))/15
<     level=[round(10**(dc*i+np.log10(mnv)),N) for i in range(15)]
---
>     nlev=25
>     dc=(np.log10(mxv)-np.log10(mnv))/nlev
>     level=[round(10**(dc*i+np.log10(mnv)),N) for i in range(nlev)]
70a72
>   if v=='PM25_TOT':mxv=204.
74c76,80
<   level,nm=get_lev(N)
---
>   level,nm=get_lev(N)
>   clrs=get_cmap("rainbow")
>   if v=='PM25_TOT':
>     level=[2,3,5,7]+[i for i in range(10,70,5)]+[75,90,100,125,150,200]
>     clrs=get_cmap("RdYlGn").reversed()
95c101
<                        cmap=get_cmap("rainbow"),
---
>                        cmap=clrs,
  • 說明
    1. 色階層數原來設定15層,此次改為25層,以因應21層之指定需要。
    2. 因程式的引數是nc檔名。而同一檔案內可能有其他汙染物存在,會在程式內依序產生,因此需按照汙染項目來決定是否啟動紅綠燈色階。
    3. 最大值取204,此值為試誤結果,因取log scale,無法確定最大值。
    4. 一般物質仍取彩虹色階,只PM2.5為比較而取紅綠燈色階
    5. 各階層濃度採取一樣的序列
    6. 紅綠燈色階的低值為紅色,因此需要做.reversed,讓低值是綠色段。
    7. 低濃度段的色階解析度太差,是各家的共同的缺點之一。這讓圖面的大部分都是空值、或單一值。此處除了[2,7]之間再加上3,5等2各階層,以便讓黃色部分代表15~20μg/M3,讓10~15μg/M3落在淺綠色段落,可以用來與其他單位模擬結果進行比較。

結果

討論

  1. 彩虹色階中雖然確實含有綠色,但因純色較多,真正分配給綠色階段的濃度等級很少或不存在。雖然能提供最多的訊息,但也只能放棄。
  2. “seismic”雖然只有藍紅二色,但沒有綠色,而且中間濃度皆為淺的藍或淺紅色,辨識非常困難。
  3. 減少顏色如jet色階(與rainbow差在靛紫色),低濃度段仍然沒有太高解析度。
  4. “afmhot”只有紅色~黃色,因沒有綠色不適合。
  5. 2各色階相連,可以得到較不連續的色階譜。還蠻突兀的。(參考ChatGPT建議)
  6. 有關3各單位的模擬結果學理討論詳見公版vs學版模式模擬結果比較
...
#cmap3 = plt.get_cmap("seismic").reversed()
cmap3 = plt.get_cmap("RdYlGn").reversed()
#cmap2 = plt.get_cmap("afmhot").reversed()
cmap1 = plt.get_cmap("winter")
#cmap1 = plt.get_cmap("jet")
#cmap1 = plt.get_cmap("rainbow")
n = 256  # ▒~O▒~I▒▒~U▒▒~G~O
colors1 = cmap1(np.linspace(0., 1., n))
colors2 = cmap2(np.linspace(0., 1., n))
colors3 = cmap3(np.linspace(0., 1., n))
combined_colors = np.concatenate((colors1, colors3), axis=0)
combined_cmap = ListedColormap(combined_colors)
...
    contours = plt.contourf(lons, lats, nc[v][t,0,:,:] ,
                        ...
                       cmap=combined_cmap,
                        ...
                       extend='max')