m3nc檔案轉GIF

Table of contents

背景
輸入及輸出
png2gif檔案轉換(convert)
- 程式外進行方式
- size的問題
m3nc2gif.py程式設計重點
結果比較檢討
更新色階至RdYlGn(紅綠燈)

背景

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)))，這會造成程式錯誤。再增加一位有效位數即可。
  - 測試結果如下

mxv範例	有效位數N	有效數字
12345.12345	-1	12350.
1234.1234	0	1234.
123.123	0	123.
12.12	1	12.1
1.1	2	1.10
0.12	3	0.120
0.123	3	0.123
0.01234	4	0.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-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(不等間距等濃度)

檢討

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

其他有待精進項目
- 左右側範圍
- 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³落在淺綠色段落，可以用來與其他單位模擬結果進行比較。

結果

討論

彩虹色階中雖然確實含有綠色，但因純色較多，真正分配給綠色階段的濃度等級很少或不存在。雖然能提供最多的訊息，但也只能放棄。
“seismic”雖然只有藍紅二色，但沒有綠色，而且中間濃度皆為淺的藍或淺紅色，辨識非常困難。
減少顏色如jet色階(與rainbow差在靛紫色)，低濃度段仍然沒有太高解析度。
“afmhot”只有紅色~黃色，因沒有綠色不適合。
2各色階相連，可以得到較不連續的色階譜。還蠻突兀的。(參考ChatGPT建議)
有關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')

m3nc檔案轉GIF

背景

輸入及輸出

引數

輸入檔案

輸出檔案

png2gif檔案轉換(convert)

程式外進行方式

size的問題

m3nc2gif.py程式設計重點

座標系統

浮動的濃度等級

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

程式內進行png2gif之convert

結果比較檢討

比較

檢討

降雨的效果

更新色階至RdYlGn(紅綠燈)

緣起

變更項目

結果

討論

png2gif檔案轉換(`convert`)

程式內進行png2gif之`convert`