Link Search Menu Expand Document
Focus on Air Quality

ISAM結果檔案之讀取(PM25_IONS)

Table of contents

背景

  • 雖然ISAM設計輸出PM25_IONS的成分與其來源追蹤，但結果並沒有一個變數項目稱為PM25_IONS，找不到類似combine的後處理程式，可以整合CMAQ-ISAM執行的成果。只能自行撰寫程式。
    • CMAQ-ISAM執行的成果CCTM_SA_ACONC變數的命名規則：以底線(_)區隔之複合變數spec_group，其中：
      • spec=CCTM_ACONC的污染項目，詳cb6r3_ae7_aq之物質名稱
      • group=isam_control.txt檔案裏定義的TAG_NAME，另外還包括ICON、BCON、OTHR等固定內設的標籤。
  • 綜合性污染物的定義：可以參考$REPO_HOME/POST/combine/scripts/spec_def_files/SpecDef_${MECH}.txt的內容
  • ISAM可能同時執行多個分區之排放分析，因此程式輸出結果須能有所區別
    • 分區的定義在isam_control.txt檔案裏的REGION(S)，指定到$ISAM_REGIONS檔案的分區名稱內容
    • 此處以中國大陸的空氣質量預報分區為例
空氣質量預報分區$ISAM_REGIONS(REGION(S))名稱ISAM結果檔名標籤
華北(河北山東)AQFZ1JJZ
汾渭平原、山西、陜西、河南AQFZ2FWS
東北區域AQFZ3NEC
甘肅塞外AQFZ4NWC
兩湖華南AQFZ5SCH
四川廣西青康藏AQFZ6SWC
華東江蘇浙江安徽AQFZ7YZD

程式說明

ISAM執行方式

連續執行批次

kuang@DEVP /home/cmaqruns/2018base
$ cat do_isam.csh
foreach BSN ('JJZ' 'SCH' 'YZD' 'FWS' 'NEC' 'NWC')
  source run_isamMM_RR_DM.csh 04 6 d01 $BSN
end

執行腳本

腳本說明

程式名稱

執行方式

引數

  • CCTM_SA_ACONC檔案名稱
  • eg. CCTM_SA_ACONC_v53_gcc_1804_run6_20180404_EAsia_81K_11_SCH.nc

I/O檔案

  • CCTM_SA_ACONC檔案：為執行CMAQ-ISAM的結果
  • PM25_IONS模版檔案：template_PM25_IONS.nc
    • 單一污染物、單一時間、地面層
  • 輸出結果檔案：'PM25_IONS'+path+'_'+ymdh+'_'+g+'.nc'
    • path：ISAM分區(詳前表ISAM結果檔名標籤及執行批次)
    • ymdh：年月日時
    • g：排放類別標籤(TAG_NAME)

執行腳本(proc.cs)

  • 將ISAM結果與模版連結到同一目錄
  • 每個結果檔案都執行一遍SA_PM25_IONS.py，
    • 會拆分成逐日、每個排放標籤、每個分區之分析濃度
    • 再逐層加總或合併
  • 加總使用到python程式addNC
    • 欲加總的內容項目太多了，以ncs變數累加之
    • 加總結果檔名：PM25_IONS${z}_2018040${d}.nc
  • 合併(照日期附加)
#isam job for 20180404~8, 6 AirQualityForecastZone(AQFS), only GR1~4 and PTA are taken into account
#sum up aerosol results using python program, to be PM25_IONS, see SA_PM25_IONS.py
in $(ls CCTM*nc);do python ../SA_PM25_IONS.py $nc;done >& /dev/null
for z in FWS JJZ NEC NWC SCH YZD;do 
  for d in {4..8};do
    #only GR13、GR24、PTA results(_[GP]*) are summed, see addNC, the python program
    ncs='';for nc in $(ls PM25_IONS${z}_2018040${d}_[GP]*.nc);do ncs=${ncs}" "$nc;done;
    python ~/bin/addNC $ncs PM25_IONS${z}_2018040${d}.nc
  done
done
# combine all days for each zone
for z in FWS JJZ NEC NWC SCH YZD;do ncrcat -O PM25_IONS${z}_2018040?.nc PM25_IONS${z}.nc;done

程式分段說明

  • 調用模組
import numpy as np
import netCDF4
import sys, os
  • 讀取引數(fname)
    • 由fname中讀取分區路徑(path)、小時(ymdh)
fname=sys.argv[1]
intv=fname.split('_')[-1]
path=intv.strip('.nc')
ymdh=fname.split('_')[7]
  • 開啟ISAM執行成果檔案
    • 讀取變數名稱，從其中區分成分項目v4，與TAG_NAME標籤grp
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]]
nt,nlay,nrow,ncol=(nc.variables[V[3][0]].shape[i] for i in range(4))

v4=list(set([i.split('_')[0] for i in V[3]]));v4.sort()
grp=list(set([i.split('_')[1] for i in V[3]]));grp.sort()
  • PM的所有可能成分。由於ISAM只會輸出PM25_IONS的成分，包括在aer集合之內。
Sulfate='ASO4I ANO3I ANH4I ACLI ASO4J ANO3J ANH4J ASO4K ANO3K ANH4K'.split()
BC='AECI AECJ'.split()
OC='APOCI APNCOMI APOCJ AOTHRJ AXYL1J AXYL2J AXYL3J ATOL1J  ATOL2J  ATOL3J  ABNZ1J  ABNZ2J  ABNZ3J\
  AISO1J AISO2J  AISO3J  ATRP1J  ATRP2J  ASQTJ  AALK1J  AALK2J AORGCJ  AOLGBJ  AOLGAJ  APAH1J  APAH2J\
  APAH3J  APNCOMJ AOTHRI'.replace('_','').split()
Dust='AFEJ AALJ ASIJ ACAJ AMGJ AKJ AMNJ ACORS ASOIL ATIJ'.split()
SS='ANAI ANAJ ACLJ ACLK ASEACAT'.split()
Semi='ASVPO1J ASVPO2J  ASVPO3J  ASVPO2I  ASVPO1I ALVPO1J   ALVPO1I  AIVPO1J'.split()
aer=Sulfate+BC+OC+Dust+SS+Semi
  • 逐一拆解檔案
    • 定義輸結果檔案名稱：'PM25_IONS'+path+'_'+ymdh+'_'+g+'.nc'
for g in grp:
  fname = 'PM25_IONS'+path+'_'+ymdh+'_'+g+'.nc'
  os.system('cp template_PM25_IONS.nc '+fname)
  nco= netCDF4.Dataset(fname, 'r+')
  • 模版為單一時間檔案，須延長時間，以配合輸入檔長度。
  for t in range(nt):  
    nco['TFLAG'][t,0,:]=nc['TFLAG'][t,0,:]
  • 累加PM25_IONS
    • 因ISAM未處理氯離子等很多項目，矩陣無內容，因此須研判是否為nan，以避免被netCDF4遮蔽
  var=np.zeros(shape=(nt,nlay,nrow,ncol))
  for v in set(aer) & set(v4):
    vg=v+'_'+g
    if np.isnan(np.max(nc[vg][:])):continue
    var[:]+=nc[vg][:]
  nco['PM25_IONS'][:]=var[:]
  nco.close()

程式下載

成果檢視

20180404ISAM-wanli.PNG
2018/03/31-04/07沙塵暴期間萬里測站PM10各分區人為污染貢獻濃度之比較
CMAQ_WBDcomp.PNG
2018/03/31-04/07沙塵暴期間萬里測站PM10風吹揚砂濃度之比較
  • 時間序列以下列簡易程式進行解讀，轉成csv檔案，以excell進行繪圖。
$ cat rd_ncs.py
import netCDF4
from pandas import *
import datetime

IX,IY=28-1,29-1

zs=['NWC','JJZ','YZD','FWS']

for z in zs:
  fname='PM10'+z+'.nc'
  nc = netCDF4.Dataset(fname, 'r')
  exec(z+'=nc["PM10"][:,0,IY,IX]')
  if z==zs[0]:
    tflag=[datetime.datetime.strptime(str(i),'%Y%j').strftime('%Y%m%d')+'{:02d}'.format(j//10000) for i,j in zip(nc['TFLAG'][:,0,0],nc['TFLAG'][:,0,1])]
DD={}
for z in zs+['tflag']:
  exec('DD.update({"'+z+'":'+z+'})')
df=DataFrame(DD)
df.set_index('tflag').to_csv('ntw.csv')
  • note
    • IX,IY為VERDI畫面上北台灣的位置座標，是FORTRAN自1開始的標籤習慣，與python之標籤習慣從0開始相差1。

```

Reference