Dr. Kuang's Utilities

預報系統圖檔之拮取

 
  • Oct 20, 2022
  • views
>

背景

  • [地區空品預報系統][fcst]擁有多污染物、多時間、多解析度的龐雜檔案系統，因此當在Earth展示畫面發現有特殊的污染事件，要從中拮取高清畫面，似乎不是那麼友善。需要一個快捷的工具。
  • 此處以類似拮取環保署數據的[specHrSlider][slider]輸入方式，給定起迄時間(YYYMMDDHH格式)、污染物名稱、以及模擬網格編號，即可使用[m3nc2gif][m3nc2gif]程式（see [[4.m3nc2gif]] )，將CMAQ結果轉成gif檔案。
    • 不同的是[specHrSlider][slider]是個python程式，此處是bash的腳本。
    • [m3nc2gif][m3nc2gif]也是個python程式，需在本地安裝[wrf-python][wrf-python]

程式腳本IO

命令列引數

使用範例：$fcst/aconc2gif.cs 2022102209 2022102409 PM10 1

  1. 開始日期時間：UTC、YYYMMDDHH格式，小時自0時~23時
  2. 結束日期時間：UTC、YYYMMDDHH格式，可與前述開始時間相同，但不能較早(test $b > $s && exit 0)
  3. 污染項目：目前以PMs、VOC等為選擇
  4. 模擬網格編號：1~3

輸入檔案

輸出gif檔案

  • 會存在root=${fcst}/${GRD[$d]}/cctm.fcst/daily目錄下
  • 檔名為${s}_${b}.gif，如PM10_2022102409.gif

呼叫程式

  • ~/bin/j2c：將julian day轉換成calendar day以便串連到模式輸出檔名，原理與使用可以參考[儒略日to月曆日][j2c]
  • NCO程式
    • ncks
      • 模擬結果檔案的切割
      • 按序號進行切割，序號(自0起算)，與時間(0~23)定義一致
    • ncrcat
      • 切割結果循序連接
      • 此處將需連接的每個檔名以序列變數明確指定，以避免錯誤。
    • ncatted
      • 正常檔案的STIME全域變數為0，如果起始小時($bh)非為0，則需要修改。
      • 否則以VERDI開啟時，會有錯誤的時間註記。
      • test $bh >0 && NCATTED -a STIME,global,o,f,${bh}0000. $s.nc
    • 個別程式之使用可參考NCO使用手冊、或教學範例
  • ~/bin/join_nc.py([[2022-09-15-join_nc]])，這個版本用6個小時進行日期轉換時間軸的漸變與均勻化，一樣不作用在原始模擬結果，只作用在ncrcat後之暫存檔，說明詳下。
  • ~/bin/[m3nc2gif.py][m3nc2gif]：進行[wrf-python][wrf-python]製圖與convert。

程式設計

日期的轉換

  • 用cut指令將輸入的起(b)迄(e)時間抽出日期及時間
  • 月曆日(calendar day)轉換成儒略日(julian day)：(date -d "$?d" +%Y%j)
    • 避免使用cal2jul.f([[2022-10-21-cal2jul]])，降低程式的複雜度、跨機器需重新編譯的困擾。
    • 在MacOS平台date指令略有不同，需注意修改。
  • 小時數也從字串(01~24)轉換成10進位之整數(#10)，以利四則計算。
bd=$(echo $b|cut -c 1-8)
bj=$(date -d "$bd" +%Y%j)
bh=$(echo $b|cut -c 9-10)
bh=$(( 10#$bh ))

ed=$(echo $e|cut -c 1-8)
ej=$(date -d "$ed" +%Y%j)
eh=$(echo $e|cut -c 9-10)
eh=$(( 10#$eh ))
  • 每天的起迄小時：切割時需要精確到逐時
    • 如果是端點日期，需要改成起迄時
    • 如果不是，則為0~23
  jbh=0; test $bd == ${dates[$jd1]} && jbh=$bh
  jeh=23; test $ed == ${dates[$jd1]} && jeh=$eh

可變動的迴圈範圍

  • bash不接受可變動的迴圈(for i in {${bj}..${ej});)，主要因為{…}是優先保留給環境變數使用。
  • bash可以接受for i in $(seq $bj $ej);do的寫法。$()只有執行一次，還算單純。
  • 本腳本一共用了4次
...
dates=();for i in $(seq $bj $ej);do dates=( ${dates[@]} $(~/bin/j2c $i) );done
...
for jd in $(seq 1 $nd);do
...
  for jh in $(seq $jbh $jeh);do
...
rm ${s}_*.png;for i in $( seq 0 $len);do rm ${fnames[$i]};done

可變動的序列內容

  • bash的序列可以是固定長度內容，明確寫在程式內，也可以在程式內逐一增加擴充。本腳本一共用了5個序列，3個是固定，2個是浮動(每次執行臨時產生)。
  • 產生所有日期之序列：
    • 因應檔案的命名系統都是使用日期，
    • 因此需要起、迄之間每一天的日期序列、與序列長度
dates=();for i in $(seq $bj $ej);do dates=( ${dates[@]} $(~/bin/j2c $i) );done
nd=${#dates[@]}
  • ncrcat需要的引數
    • 過去ncrcat會以符合條件的所有檔(通配符 wild card *)來指定需連結的檔案，不但不明確也容易出錯。以序列定義絕對不會出錯。
    • 切割產生的結果檔，需全部、依序排列進行ncrcat將其連結成新的檔案，此時的結果檔名也以序列將其儲存、
    • 序列內容再以echo及eval傳到ncrcat程式內，避免使用通配符。
fnames=()
for jd in $(seq 1 $nd);do
...
  for jh in $(seq $jbh $jeh);do
    fnameo=${root}/${s}${dates[$jd1]}_$jh
    $NCKS -O -v TFLAG,$s -d TSTEP,$jh $fname $fnameo
    fnames=( ${fnames[@]} $fnameo )
  done
done
a='$NCRCAT -O '$(echo ${fnames[@]})' $s.nc'
eval $a

join_nc.py

  • 因逐日重新進行濃度場的初始化，在範圍較小的(如d02東南中國及d03臺灣島範圍)的跨日串連結果中會出現突兀的變化，需要進行漸變或均勻化處理，以使gif不會產生嚴重的跳動。
  • 這個版本一樣不作用在原始模擬結果，只作用在ncrcat後之暫存檔。

程式之執行

  • 引數：連結好、需修改的nc檔案
  • 結果：以同一檔名覆蓋原來檔案
d=$4 #domain 1~3
...
d=$(( $d - 1 ))
...
if [[ $d == 1 || $d == 2 ]];then ~/bin/join_nc.py $s.nc;fi

程式執行結果檢討

  • 雖然漸變與均勻化確實對跨日連結的結果造成一些正面的效果，但也有不合理之處
    1. 如果次日初始化濃度較高，煙陣移動會有逆向的奇怪效果
    2. 相反如果次日濃度低，煙陣移動會有突然加速、無端消失等不合理情形。
  • 釜底抽薪之計，仍應該針對事件重新進行模擬，才會有順暢的跨日變化。

join_nc.py程式碼

  • 以0時為中心，前後tsm小時值的濃度將會被修改成線性漸變過程
  • tsm此處(經試誤)定為3
$ cat ~/bin/join_nc.py
#!/opt/anaconda3/envs/pyn_env/bin/python
import subprocess, os, sys
import netCDF4
import numpy as np

fname=sys.argv[1]
nc = netCDF4.Dataset(fname, 'r+')
V=[list(filter(lambda x:nc.variables[x].ndim==j, [i for i in nc.variables])) for j in [1,2,3,4]]
times=np.array(nc['TFLAG'][:,0,1]/10000)
nt=len(times)
if 0 not in times:sys.exit()
idx=np.where(times==0)[0]
tsm=3
ts2=tsm*2
for v in V[3]:
  for i in idx:
    if i-tsm<0 or i+tsm+1>nt :continue
    var0=nc[v][i-tsm,:,:,:]
    var1=nc[v][i+tsm,:,:,:]
    for t in range(i-tsm,i+tsm+1):
      tt=t-(i-tsm)
      nc[v][t,:,:,:]=var0*(ts2-tt)/ts2+var1*tt/ts2
nc.close()

腳本內容

kuang@master /nas2/cmaqruns/2022fcst
$ cat aconc2gif.cs
NCKS=/usr/bin/ncks
NCRCAT=/usr/bin/ncrcat
NCATTED=/usr/bin/ncatted
fcst=/nas2/cmaqruns/2022fcst
b=$1 #begin YYYYMMDDHH
e=$2 #end YYYYMMDDHH
s=$3 #species name
d=$4 #domain 1~3
PM_s=( "PM10" "PM1_TOT" "PM25_TOT" "PMC_TOT" "VOC" )
GRD=( 'grid45'  'grid09'  'grid03' )
DOM=( 'CWBWRF_45k' 'SECN_9k' 'TWEPA_3k')


bd=$(echo $b|cut -c 1-8)
bj=$(date -d "$bd" +%Y%j)
bh=$(echo $b|cut -c 9-10)
bh=$(( 10#$bh ))

ed=$(echo $e|cut -c 1-8)
ej=$(date -d "$ed" +%Y%j)
eh=$(echo $e|cut -c 9-10)
eh=$(( 10#$eh ))

dates=();for i in $(seq $bj $ej);do dates=( ${dates[@]} $(~/bin/j2c $i) );done
nd=${#dates[@]}

len=$(( ( $ej - $bj ) * 24 + $eh - $bh ))

js=-1
for i in {0..4};do
  if [ $s == ${PM_s[$i]} ];then
    js=$i
  fi
done
d=$(( $d - 1 ))
if [ $js -ge 0 ];then
  FN=PMs
else
  FN=CCTM_ACONC_v532_intel_${DOM[$d]}_
fi

root=${fcst}/${GRD[$d]}/cctm.fcst/daily
fnames=()
for jd in $(seq 1 $nd);do
  jd1=$(( $jd - 1))
  fname=${root}/PMs${dates[$jd1]}.nc
  jbh=0; test $bd == ${dates[$jd1]} && jbh=$bh
  jeh=23; test $ed == ${dates[$jd1]} && jeh=$eh
  for jh in $(seq $jbh $jeh);do
    fnameo=${root}/${s}${dates[$jd1]}_$jh
    $NCKS -O -v TFLAG,$s -d TSTEP,$jh $fname $fnameo
    fnames=( ${fnames[@]} $fnameo )
  done
done
pwd=$PWD
cd $root
a='$NCRCAT -O '$(echo ${fnames[@]})' $s.nc'
eval $a
test $bh >0 && NCATTED -a STIME,global,o,f,${bh}0000. $s.nc
rm ${s}_*.png;for i in $( seq 0 $len);do rm ${fnames[$i]};done
if [[ $d == 1 || $d == 2 ]];then ~/bin/join_nc.py $s.nc;fi
~/bin/m3nc2gif.py $s.nc
mv $s.gif ${s}_${b}.gif
cd $pwd
echo $root/${s}_${b}.gif has been generated
Search