逐日排放檔案之切割程式

Table of contents

逐日排放檔案之切割程式

背景

雖然空品預報是模擬當天的空氣品質，但是大多數排放數據資料庫仍然是舊有數據，除了需要校正到正確的月份與星期，以符合氣候與星期週期的排放變化，也需要將舊的日期標籤改換成新的標籤，以符合CMAQ的要求。
- d01、d02仍然使用REAS v3.2數據，
- d03則使用TEDS11數據
- 即期排放(如電廠運轉率、CEMS數據、即時交通數據或預報數據)仍然有待發展

作業流程

在整體預報作業流程中的位置

因排放檔案與氣象處理沒有相關，因此可以在背景執行。
將網格範圍的選擇納入python程式中，讀取所在目錄而定，因d01、d02的方式很接近，並列在同一個程式內，用if block來切換

#execute ungrib and metgrid in background
...
#background executions of mk_emis and mk_ptse
for i in 0 1;do
  cd $fcst/${GRD[$i]}/smoke
  ~/bin/sub ../../mk_emis.py $BEGD
done
~/bin/sub $fcst/em3.cs

#CMAQ stream
...

d03有地面及高空汙染源之分，處理程序較為不同，使用下列腳本($fcst/em3.cs)。

kuang@dev2 /u01/cmaqruns/2022fcst/grid03/smoke
$ cat ../../em3.cs
cd $fcst/grid03/smoke
if [[ $BEGD == "" ]];then
  today=$(date -d -0day +%Y%m%d)
  export BEGD=$(date -d "$today -0days" +%Y-%m-%d)
  HR=10#$(date +%H)
  if [[ $HR -ge 7 && $HR -lt 16 ]];then export BEGD=$(date -d "$today -1days" +%Y-%m-%d);fi
fi
../../mk_emis.py $BEGD
/usr/bin/ncks -O -d LAY,0 TEDS.ncf TEDS0.ncf
/usr/bin/ncatted -a NLAYS,global,o,i,1 TEDS0.ncf
./mk_ptse.py $BEGD

程式說明

IO

引數：模擬開始日期$BEGD，格式為+%Y-%m-%d
程式也會讀取工作目錄，從中讀取網格設定，以開啟正確的排放資料母體檔案
fnames.txt這個檔案將TEDS11排放檔路徑集中，由檔名中讀取起始日期(bdate7s)與檔案天數(days)，經對齊星期日數後，篩選含有$BEGD之檔案。詳下說明。

網格	工作目錄	排放資料母體檔案	說明
d01	/nas2/cmaqruns/2022fcst/grid45/smoke	/nas1/TEDS/REAS3.2/origins/2015_D7.nc
d02	/nas2/cmaqruns/2022fcst/grid09/smoke	/nas1/TEDS/REAS3.2/origins/2015_D8.nc
d03	/nas2/cmaqruns/2022fcst/grid09/smoke	/nas2/cmaq2019/download-20220503/input/2019+mm+/grid03/smoke/	mm=01~12，檔名詳下說明

kuang@dev2 /u01/cmaqruns/2022fcst/grid03/smoke
$ cat fnames.txt
.../201901/grid03/smoke/b3gts_l.20181225.38.d4.ea2019_d4.ncf
.../201902/grid03/smoke/b3gts_l.20190125.35.d4.ea2019_d4.ncf
.../201903/grid03/smoke/b3gts_l.20190222.38.d4.ea2019_d4.ncf
.../201904/grid03/smoke/b3gts_l.20190325.37.d4.ea2019_d4.ncf
.../201905/grid03/smoke/b3gts_l.20190424.38.d4.ea2019_d4.ncf
.../201906/grid03/smoke/b3gts_l.20190525.37.d4.ea2019_d4.ncf
.../201907/grid03/smoke/b3gts_l.20190624.38.d4.ea2019_d4.ncf
.../201908/grid03/smoke/b3gts_l.20190725.38.d4.ea2019_d4.ncf
.../201909/grid03/smoke/b3gts_l.20190825.37.d4.ea2019_d4.ncf
.../201910/grid03/smoke/b3gts_l.20190924.38.d4.ea2019_d4.ncf
.../201911/grid03/smoke/b3gts_l.20191025.37.d4.ea2019_d4.ncf
.../201912/grid03/smoke/b3gts_l.20191124.38.d4.ea2019_d4.ncf

公版模式排放檔案命名規則

日期及日數(date_day)：公版模式自每月前7天起始(bdate7)，同時將模擬總日數寫在檔案名稱之中，而每月的日數不同，因此需事先讀取所有的序列備查。
月份：因起始日的月份與實體月份不同，為方便程式撰寫，也改成檔名序列備查(smk=fns[m].split('b3g')[0])。
該模式共讀取3個網格排放檔案：生物源('begts.ncf')、TEDS排放源(‘TEDS.ncf’)、以及粗網格貢獻量('egts.ncf')。依序將檔名個元件組裝後、寫成fnames序列。

  yr=tdy.split('-')[0]
  mm=tdy.split('-')[1]
  with open('fnames.txt','r') as f:
    fns=[i.strip('\n') for i in f]
  sbdates=[i.split('.')[1] for i in fns]
  days=[i.split('.')[2] for i in fns]
  bdate7s=[datetime.datetime.strptime(i,"%Y%m%d") for i in sbdates]
  deld=bdate7s[-1].weekday()-datetime.datetime(int(yr), bdate7s[-1].month, bdate7s[-1].day).weekday()
  if deld<0:deld+=7
  bdate7m=[i-datetime.timedelta(days=deld) for i in bdate7s]
  bdate7m=[datetime.datetime(int(yr), i.month, i.day) for i in bdate7m]
  bdate7m[0]=[datetime.datetime(int(yr)-1, i.month, i.day) for i in bdate7m[:1]][0]
  m=-1
  for i in range(12):
    if bdate>bdate7m[i]: m=i
  if m<0:sys.exit('date before Christmas too much')
  mm='{:02d}'.format(bdate7m[m].month)
  date_day=sbdates[m]+'.'+days[m]
  smk=fns[m].split('b3g')[0]
  doms=['d4.ea2019_d4','TW3-d4.BaseEms','d4.ea2019_d4']
  kind=['b3gts_l.','cmaq_cb06r3_ae7_aq.'+mm+'-','egts_l.']
  exts=['tar.xz',  'tar.gz',                    'tar.xz' ]
  opts=['xvfJ',    'xvfz',                      'xvfJ'   ]
  fnameO=['begts.ncf','TEDS.ncf','egts.ncf']
  fnames=[smk+kind[i]+date_day+'.'+doms[i]+'.ncf' for i in range(nf)]

使用公版原名稱序列、以查詢方式找到正確月份的必要性
1. 每次計算年代間星期差異，適用在所有任意年代。
2. 預報時間範圍是否跨越不同的月份，與預報時間長度有關。公版模式每月只有多出1~2天，因此跨越計算的邏輯並不容易撰寫。
3. 使用序列查詢，只須確保公版模式檔案有跨越到下月初前2天即可(倒數第8天~下月初前2天=10天、倒數第7天以後則會開啟下個月檔案)。

TFLAG內容之填入

雖然是ncks切割出來的檔案，已經具有SDATE、TFLAG等屬性或變數內容，但因為原檔案是2019年，因此所有的時間標籤都必須修正。
先產生一個datetime的序列，儲存所有模擬時間範圍內逐時的datetime內容。
datetime與整數的轉換，使用dt2jul小程式。
將第0個汙染物的TFLAG轉換成新的時間標籤。此處要避免對污染物維度進行迴圈、nc檔案內容如果一一置換，其工作效率很低。
1. 先將時間標籤另存成2維的矩陣(var=np.array(nc1.variables['TFLAG'][:,0,:]))備用。np.array的展開效率遠大於一一填入nc檔案。
2. 使用None將2維矩陣複製到3維矩陣。(var3[:,:,:]=var[:,None,:])
3. 3維矩陣一次整批倒入nc檔案內。

for jf in range(nf):
  nc1 = netCDF4.Dataset(fnameO[jf],'r+')
  V1=[list(filter(lambda x:nc1.variables[x].ndim==j, [i for i in nc1.variables])) for j in [1,2,3,4]]
#bdate=datetime.datetime.combine(datetime.date.today(),datetime.time(0,0,0,0))
  nc1.SDATE,nc1.STIME=dt2jul(bdate)
  nt1=nhr
  SDATE=[bdate+datetime.timedelta(hours=int(i)) for i in range(nt1)]
  for t in range(nt1):
    nc1.variables['TFLAG'][t,0,:]=dt2jul(SDATE[t])
  var=np.array(nc1.variables['TFLAG'][:,0,:])
  var3=np.zeros(shape=nc1.variables['TFLAG'].shape)
  var3[:,:,:]=var[:,None,:]
  nc1.variables['TFLAG'][:]=var3[:]
  nc1.TSTEP=10000

d03點源的處理

由於公版模式將點源併入網格之中，大幅降低其解析度，此處只取其第一層結果。
高空部分
- 以CMAQ點源檔案型式輸入，檔案引用Focus-on-Air-Quality/EmisProc/ptse之處理結果。
- 在mk_ptse.py中對時間標籤加以處理。

程式下載點

Download: 逐日排放檔案之切割程式mk_emis.py

Download: 逐日點源排放檔案之切割程式mk_ptse.py