m3nc檔案轉GIF

背景

• 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_xlim及ax.set_ylim
  • 當模擬範圍較小時，為取整數刻度，matplotlib會將範圍向外多取一些，導致圖面出現空白框。可以使用此2指令強制縮小範圍。
  • cartopy_xlim(p)等為wrf-python指令，結果為下限與上限值之陣列。因此如需增減，在其後加上一陣列即可。
  • 上下限與增減單位為公尺
  • 如果內縮太多，可另行調整。

浮動的濃度等級

• 如前所述，濃度等級需要幾個項目是浮動的(至少想到5個)：
  • 1.levels極值的設定：全時間、全域、正值之
    • 99.99 percentile (等級之最大值mxv)
    • 0.01 percentile (等級之最小值mnv)
    • 不服務負值的理由是：
      • 正負值同框出現，不符合一般直覺與作圖用意。
      • 取log10會出錯
      • 取低值會差很多(0值可能在增量模擬結果會非常多)。
    • 如果負值一定要套用，可以先在nc檔案階段處理，將負值轉正，再應用繪圖程式。
    • 某時間的最大值低於mnv，contourf將會出錯、跳出程式，不會產生gif檔：
      • 發生在夏季午後降雨事件，濃度會在短時間內消失。
      • 如果是前面階段，從未產生過圖檔，只能捨棄這些小時。
      • 如果是中間、後面階段，已經產生過圖檔，還可以沿用上一小時之contourf結果
      • 檢討：因最大值低於mnv是非常極端個案，整體時間變化看起來是低濃度持續一陣子，還算順暢。
  • 2.線性或對數色階，改成彈性判斷
    • 最大/最小比值高於15：非線性、對數色階，適用在煙流增量分布，
    • 比值低於15：線性色階，適用在濃度變化不大的空品項目。
  • 3.濃度等級的階層數(nlev)
    • 對數色階的數字不是等間距，沒有階層數的問題，一律使用15階即可
    • 線性色階間距的數字需要是整數，以便能直接辨識濃度值
    • 此處分為2類，最大值第1碼為[3,6,9]者取15階，其餘(7視為6)取10階
  • 4.對數等級對應到顏色的分布(norm=nm)
    • 正常內設狀況下線性的等級對應到均勻的色階顏色，是沒有問題的。
    • 如果是對數濃度等級對應到內設色階分布，結果還是集中在高濃度，並沒有改變。需要改變顏色的分布。
    • matplotlib.colors至少有Normalize、LogNorm、SymLogNorm、PowerNorm、BoundaryNorm、OffsetNorm等分布型態。
    • 此處之nm值應用到前2者。顏色分布的極值(vmin,vmax)調整則以圖面上保持最多顏色而且均勻為原則。
  • 5.同時應用在小數及大數的有效位數 N
    • 濃度太低、太高的情形應考慮轉換單位、重新整理nc檔案，而不是在圖面上用有效位數調整。因此不考慮用科學記號。
    • N=int(3-np.log10(mxv))，大致上可以保持有3~4個有效數字，
    • 有效位數不足：低值部分因去尾造成某些區間的間格為0 (len(level)!=len(set(level)))，這會造成程式錯誤。再增加一位有效位數即可。
    • 測試結果如下

• 相關程式碼

#判斷線性或對數濃度色階
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-or-jpg的選擇
• 程序：去邊、統一增加5%的邊、取得png尺寸、整併成gif檔

...
    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模擬最大小時值之濃度分布

公版後製工具(線性等間距格柵)	VERDI(對數等間距格柵)	m3nc2gif(不等間距等濃度)

檢討

• 其他有待精進項目
  • 左右側範圍
  • XY軸單位
  • 中文Title/subTitle?

降雨的效果

• NullSchool 3HPA

更新色階至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/M³，讓10~15μg/M³落在淺綠色段落，可以用來與其他單位模擬結果進行比較。

結果

討論

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')