Link Search Menu Expand Document
Focus on Air Quality

CMAQ/CAMx排放量檔案之轉換

Table of contents

背景

  • 網格模式的排放量檔案自然也是5d規格檔案，但是各個模式略有不同，差異比較如下表：

CMAQ/CAMx排放量檔案之差異比較

項目CAMx模式CMAQ模式
面源格式uamiv(可以用pncdump指令檢視，CAMx6以後可以接受nc格式，但6與7版本也有差異IOAPI-ncf(可以用ncdump指令檢視)
高空面源(not try)可接受
點源格式ptsource(可以用pncdump指令檢視)IOAPI-ncf(可以用ncdump指令檢視)
時間範圍個案～全月(LST)逐日分開(UTC)、v5.2+可以接受全月同一檔
污染物項目可與主程式略有不同必須與主程式完全一致
點源檔案個數1（標頭及時序同一檔），CAMx6以後可以接受nc格式，但6與7版本也有差異const+timvar = 2
排放量單位(氣狀物/粒狀物)gmole/hr, g/hrgmole/s, g/s
nc檔時間標籤TFLAG與ETFLAGTFLAG
nc檔全域屬性除ioapi標準外有另外的需求。< v7.10會嚴格檢查NAME內容ioapi標準項目
NAME, CAMx_NAME 指定內容EMISSION, PTSOURCE(無)
  • 即使是 >CAMx6 的nc檔案，與CMAQ的格式還是有很大的差異，可以參考進一步點源的煙囪參數及排放量比較。
  • 除了檔案格式之外，在時間定義、污染項目、排放量單位、檔案切割(合格)等，都有所差異。
  • 所以轉換作業的重點包括
    1. 面源檔案格式、污染項目的改變
    2. 點源檔案格式、污染項目、表頭表身的分割（合併）

CMAQ2CAMx程式之準備

  1. Environ RAMBOLL公司提供了轉換程式，可以轉換排放量(面源及點源)、IC、BC。
  2. 程式除了連結netCDF程式庫之外，還需要IOAPI程式庫。下載與編譯可以參考CMAQ2CAMx之單向轉換。
  3. 使用化學物質對照表來對應2個模式的物質種類

camx2ioapi程式之準備

  • 下載、編譯及執行範例可以詳見CAMx2ioapi說明。

Ramboll轉換工具的差異

  • 前述2個程式雖然是2016年版本，迄今有些歷史，然而Ramboll公司並沒有更新的計畫，主要是2個模式版本更新太快，這些converter似乎永遠都跟不上。
項目cmaq2camxcamx2ioapi說明
年代20162016 
版本v2v6 
檔案類型icon, bcon, emis, ipntuamiv格式檔包括AVRG, EMIS, DDEP, WDEPinline point source 現已廢棄
污染物對照有無 
垂直層對照無有 

面源檔案之轉換

  • 早期Environ公司曾經提供此一轉換程式(camx2ioapi)，雖然當時連結IOAPI3.1，但似乎仍可以與IOAPI3.2混合運用。
  • 下載、編譯、詳見CAMx2IOAPI檔案之轉換

執行CAMx(uamiv) to CMAQ(ioapi)

  1. 使用csh環境設定，如下表。source camx2cmaq.job 01 (01~12為月份)
  2. MAP_PROJ_OVERRIDE ：如投影方式改變，可以在此設定。
  3. TIMEZONE_OVERRIDE ：時區設定，CMAQ系統為UTC，CAMx可以接受LST，必須在此改變。
  4. UNITS3D_OVERRIDE ：單位轉換，CMAQ為mole/s，CAMx為mole/hour
kuang@master /nas1/cmaqruns/2016base/data/emis
$ cat camx2cmaq.job
setenv MAP_PROJ_OVERRIDE "2 10.000 40.000 120.990 120.990 23.610"
setenv VGTOP3D_OVERRIDE 5000.0
setenv TIMEZONE_OVERRIDE -8 # PST
setenv OUT_LAYER_INDEX 1
setenv GDNAM_OVERRIDE "sChina_27k"
setenv MON $argv[1]
setenv CAMx_IN  /home/camxruns/2016/emis/fortBE.213_teds10.base$MON
setenv IOAPI_OUT /nas1/cmaqruns/2016base/data/emis/fortBE.213.teds10.base$MON.nc
rm -f $IOAPI_OUT

/cluster/src/CAMx/camx2ioapi_v6/src/camx2ioapi << EOF
Input CAMx filename|$CAMx_IN
Data Type          |EMIS
Sigma Levels      |
EOF

CAMx(nc) to CMAQ(IOAPI)

  • 前述是CAMx(uamiv) to CMAQ(IOAPI)的轉換，近年來CAMx模式逐漸在IO格式上向CMAQ系統靠攏，除了傳統的uamiv格式之外，也容許nc格式的IO，在camx7版本中更像是IOAPI格式的nc檔案。
  • 然而畢竟是2個不同的模式，檔案內容有些微的差異，因為差異不大，2個模式的官網似乎也沒有強調，但使用者還是必須注意、處理之。
  • camx2cmaq.py程式的說明可以詳見D4範圍地面排放檔案之轉換。

點源檔案之轉換

CMAQ轉CAMx

  1. 使用csh環境設定，如下表點源範例。source conv_ipnt.job
  2. INPUT_CMAQ_EMIS  ：為點源排放量時間變化部分檔案(timvar檔)。
  3. INPUT_STACK_GRP ：為點源固定參數部分檔案(const檔)
  4. SPECIES_MAPPING ：為污染物名稱的對照表，按照模式機制選定
  5. OUTPUT_TIMEZONE ：時區設定，CMAQ系統為UTC，CAMx可以接受LST，必須在此改變。
  6. 單位轉換，CMAQ為mole/s，CAMx為mole/hour，CAMx不會檢查單位，因此轉接程式直接乘以3600(emsfac, vstk, flow)。
kuang@114-32-164-198 /Users/camxruns/src/cmaq2camx
$ cat conv_ipnt.job
#!/bin/csh -f

set INPUT_CMAQ_EMIS  = ./test_sample/inln_mole_ptipm_20020101_12US1_cmaq_cb05.ncf
set INPUT_STACK_GRP  = ./test_sample/stack_groups_ptipm_12US1_cmaq.ncf
set OUTPUT_CAMx_PNT  = ./test_sample/point.camx.20020101.bin
set SPECIES_MAPPING  = ./Species_Mapping_Tables/MAPTBL.CAMx5.2_CB05_CF.CMAQ_CB05_AE5_EMIS
set OUTPUT_TIMEZONE  = 0

setenv INFILE  $INPUT_CMAQ_EMIS
setenv OUTFILE tmp_emis
setenv MAPTBL  $SPECIES_MAPPING
rm -f $OUTFILE

./src/spcmap

setenv CMAQEMIS tmp_emis
setenv STKGROUP $INPUT_STACK_GRP
rm -f $OUTPUT_CAMx_PNT

./src/cmaq2uam << EOF
File Type          |IPNT
OutFile1 (IC,BC,AR)|
OutFile2 (TopC,PT) |$OUTPUT_CAMx_PNT
Output Timezone    |$OUTPUT_TIMEZONE
EOF

rm -f tmp_emis

CAMx2CMAQ檔案之轉換

由於python同時有ptsource及netcdf的模組，因此只要找到適合的IOAPI模版，將其適當減裁，填入ptsource的數據即可。

pt_constTWN.py常數部分(檔頭)

  1. 引數：欲轉換格式的ptsource檔案名稱(line 9)
  2. fname0 ：煙囪數夠多的模版檔案(line18)
  3. 煙囪的維度被安排在ROW，因此可以用ncks切出一樣大小的檔案(line 19)
  4. 填入global constants(line 21~33)
  5. mp：新舊檔案中變數名稱的對照字典(line 34)
  6. twd97.towgs84 ：批次進行座標轉換的做法(line51~56)
#kuang@114-32-164-198 /Users/cmaqruns/2016base/data/ptse
#$ cat pt_constTWN.py
  1 import numpy as np
  2 import netCDF4
  3 import PseudoNetCDF as pnc
  4 import os,sys
  5 import twd97
  6
  7
  8
  9 fname1=sys.argv[1].replace('fortBE.14.','').replace('base','16')
10 pt=pnc.pncopen(fname1,format='point_source')
11 v3=list(filter(lambda x:pt.variables[x].ndim==3, [i for i in pt.variables]))
12 v2=list(filter(lambda x:pt.variables[x].ndim==2, [i for i in pt.variables]))
13 v1=list(filter(lambda x:pt.variables[x].ndim==1, [i for i in pt.variables]))
14 nhr,nvar,dt=pt.variables[v3[0]].shape
15 nt,nopts=pt.variables[v2[0]].shape
16 tb=pt.STIME[0]-8 #UTC
17 fname=fname1+'.const.nc'
18 fname0='stack_groups_ptnonipm_12US1_2016ff_16j.nc'
19 os.system('/cluster/netcdf/bin/ncks -O -d ROW,1,'+str(nopts)+' '+fname0+' '+fname)
20 nc = netCDF4.Dataset(fname,'r+')
21 nc.NROWS=nopts
22 nc.GDNAM='sChina_27k'
23 nc.P_ALP = np.array(10.)
24 nc.P_BET = np.array(40.)
25 nc.P_GAM = np.array(120.98999786377)
26 nc.XCENT = np.array(120.98999786377)
27 nc.YCENT = np.array(23.6100196838379)
28 nc.XCELL=27000.000
29 nc.YCELL=27000.000
30 nc.XORIG=-877500.0-nc.XCELL/2.
31 nc.YORIG=-877500.0-nc.YCELL/2.
32 nc.SDATE=2000000+pt.SDATE[0]
33 nc.STIME=tb*10000
34 mp={'STKDM':'DSTK','STKHT':'HSTK','STKTK':'TSTK','STKVE':'VSTK','XLOCA':'XSTK', 'YLOCA':'YSTK',}
35 nc.variables['LMAJOR'][0,0,:,0]=[0 for i in range(nopts)]
36 nc.variables['LPING'][0,0,:,0]=[0 for i in range(nopts)]
37 for i in range(nopts):
38  if pt.variables['HSTK'][i]>150.:
39    nc.variables['LPING'][0,0,i,0]=1
40    nc.variables['LMAJOR'][0,0,i,0]=1
41 for v in mp:
42  nc.variables[v][0,0,:,0]=np.array(pt.variables[mp[v]][:],dtype='float32')
43 nc.variables['IFIP'][0,0,:,0]=[1000+i for i in range(nopts)]
44 nc.variables['ISTACK'][0,0,:,0]=[1+i for i in range(nopts)]
45
46 x=pt.variables['XSTK'][:]
47 y=pt.variables['YSTK'][:]
48 nc.variables['COL'][0,0,:,0]=[int((i-nc.XORIG)/nc.XCELL) for i in x]
49 nc.variables['ROW'][0,0,:,0]=[int((i-nc.YORIG)/nc.YCELL) for i in y]
50
51 Latitude_Pole, Longitude_Pole = 23.61000, 120.9900
52 Xcent, Ycent = twd97.fromwgs84(Latitude_Pole, Longitude_Pole)
53 x=x+Xcent
54 y=y+Ycent
55 ll = np.array([twd97.towgs84(i, j) for i, j in zip(x, y)])
56 lat, lon = (ll[:, i] for i in [0, 1])
57 nc.variables['TFLAG'][0,:,0]=[nc.SDATE for i in range(nc.NVARS)]
58 nc.variables['TFLAG'][0,:,1]=[nc.STIME for i in range(nc.NVARS)]
59 nc.variables['LATITUDE'][0,0,:,0]=lat
60 nc.variables['LONGITUDE'][0,0,:,0]=lon
61 nc.close()

pt_timvarTWN.py隨時間改變部分

  1. 引數：欲轉換格式的ptsource檔案名稱(line 12)
  2. fname0 ：煙囪數夠多的模版檔案(line13)
  3. 煙囪的維度被安排在ROW，因此可以用ncks切出一樣大小的檔案(line 15)
  4. 填入global constants(line 21~33)
  5. mpsp ：新舊檔案中變數名稱的對照字典(line 9, 53~54)，名字相同的直接倒入。
  6. 3600  ：ptsource排放量是逐時，CMAQ是每秒
#kuang@114-32-164-198 /Users/cmaqruns/2016base/data/ptse
#$ cat pt_timvarTWN.py
  1 import numpy as np
  2 import netCDF4
  3 import PseudoNetCDF as pnc
  4 import os,twd97,sys
  5 import datetime
  6 pt=pnc.pncopen(sys.argv[1],format='point_source')
  7 v2=list(filter(lambda x:pt.variables[x].ndim==2, [i for i in pt.variables]))
  8 nt,nopts=pt.variables[v2[0]].shape
  9 mpsp={'PNA':'NA','POC':'POA','XYLMN':'XYL'}
10
11 id=20
12 fname1=sys.argv[1].replace('fortBE.14.','').replace('base','16')
13 fname0='inln_mole_ptnonipm_20160701_12US1_cmaq_cb6_2016ff_16j.nc'
14 fname=fname1+'.timvar.nc'
15 os.system('/cluster/netcdf/bin/ncks -O -d ROW,1,'+str(nopts)+' '+fname0+' '+fname)
16
17 nc = netCDF4.Dataset(fname,'r+')
18 jt=nt
19 v4=list(filter(lambda x:nc.variables[x].ndim==4, [i for i in nc.variables]))
20 tb= int(pt.STIME[0]-8) #UTC
21 nc.NROWS=nopts
22 nc.P_ALP = 10.
23 nc.P_BET = 40.
24 nc.P_GAM = 120.98999786377
25 nc.XCENT = 120.98999786377
26 nc.YCENT = 23.6100196838379
27 nc.GDNAM='sChina_27k'
28 nc.XORIG=-877500.000-nc.XCELL/2.
29 nc.YORIG=-877500.000-nc.YCELL/2.
30 nc.XCELL=27000.000
31 nc.YCELL=27000.000
32 nc.SDATE=2000000+pt.SDATE[0]
33 nc.STIME=tb*10000
34 yyyyjjj=[nc.SDATE for i in range(tb,24)]
35 begD=datetime.datetime(int(nc.SDATE/1000),1,1)+datetime.timedelta(days=int(nc.SDATE%1000)-1)
36 for i in range((24-tb),jt,24):
37  nowD=begD+datetime.timedelta(days=int(i/24+1))
38  yr=nowD.year
39  jj=yr*1000+(nowD-datetime.datetime(yr,1,1)).days+1
40  for j in range(i,min(jt,i+24)):
41    yyyyjjj.append(jj)
42
43 for j in range(nc.NVARS):
44  nc.variables['TFLAG'][:,j,0]=yyyyjjj
45 nc.variables['TFLAG'][:,:,1]=[[(i+tb)%24*10000 for j in range(nc.NVARS)] for i in range(jt)]
46
47 for v in v4:
48  nc.variables[v][:,0,:,0]=np.zeros(shape=(jt,nopts),dtype='float32')
49
50 idx=0 #list(pt.variables['TFLAG'][:,0,0]).index(nc.SDATE)+(tb+8) #pt is in LST
51 pmothr=pt.variables['FCRS'][idx:idx+jt,:]+pt.variables['FPRM'][idx:idx+jt,:]
52 nc.variables['PMOTHR'][:,0,:,0]=np.array(pmothr,dtype='float32') / 3600.
53 for v in mpsp:
54  nc.variables[v][:,0,:,0]=np.array(pt.variables[mpsp[v]][idx:idx+jt,:], dtype='float32')/ 3600.
55 for v in v2:
56  if v not in v4:continue
57  nc.variables[v][:,0,:,0]=np.array(pt.variables[v][idx:idx+jt,:],dtype='float32') / 3600.
58 #nox= nc.variables['NO2'][:,0,:,0]+nc.variables['NO'][:,0,:,0]
59 #nc.variables['NO2'][:,0,:,0]= np.zeros(shape=(jt,nopts),dtype='float32') #nox *1./10.
60 #nc.variables['NO'][:,0,:,0] = np.zeros(shape=(jt,nopts),dtype='float32') #nox *9./10.
61 nc.close()
62

逐日拆分作業方式（brk_day.cs）

除了點源的const檔之外，其餘面源及點源timvar檔皆須拆分到逐日，可以用批次檔brk_day.cs運用ncks進行拆解。

  • 針對日期的確認，brk_day2.cs有更新作法，主要是起始時間的計算方式，前版是用前月的最後日，在run批次的計算方式中，如果前月是小月就會發生誤差。
  • brk_day2.cs(*)用最原始的定義，前月15日為起始日，按照run5的實際日期起算，逐日產生分日檔案。
$ diff brk_day.cs brk_day2.cs
12c12,15
< begd=$(date -d "20${yrmn}01 -1days" +%Y%m%d)
---
> last=$(date -d "20${yrmn}01 -1days" +%Y%m%d)
> y=$(date -d $last +%Y)
> m=$(date -d $last +%m)
> begd=$(date -d "${y}-${m}-15 +16days" +%Y%m%d)
15c18,22
< if [ $begj != $SDATE ]; then echo 'not ok in SDATE'; exit;fi
---
> if [ $begj != $SDATE ]; then
> echo $begj $SDATE 'not ok in SDATE';
> jj=$(( $SDATE - 2016001 ))
> begd=$(date -d "2016-01-01 +${jj}days" +%Y%m%d)
> fi

點源轉檔結果之確認

  1. uamiv及IOAPI ncf均可以用verdi進行展示，包括空間分布、時間序列、以及局部數字之探索(詳VERDI使用說明 )
  2. ptsource可以用pncdump或者執行pt2emNest轉成uamiv格式用verdi檢視
  3. pt2em也有python版本

pt2em.f

  1. camx600以後的版本將很多輸入檔都轉成uamiv之通用格式(濃度、沉降、氣象、地面排放、瞬間再啟動濃度檔等等)，但是點源檔案仍然有許多煙囪參數是在表頭，不能適用VERDI，因此在圖形展示上造成不小困擾。
  2. 此系列程式就是將ptsource檔案轉成uamiv檔案格式，以利VERDI可以檢視。系列程式包括2個，
  3. 簡單版本pt2em直接將輸入檔(forBE.14)以D1座標系統進行轉換(fortBE14.13)
    • pt2em File14 File13
  4. pt2emNest則可以選擇輸入檔d?.in做為條件，以進行不同範圍與解析度之圖形，網格層數將會出現在輸出檔案名稱(如fortBE.14.base.d44.in)
    • pt2emNest File14 Nest.in
  5. 注意：
    1. 由於程式將網格內所有點源直接累加，因此即使同一位置，若是解析度較低的網格，排放量會較大。
    2. VOCs無法顯示排放總量。
  6. Nest.in定義，其實是用mm5camx的控制檔，如下面d4.in的內容，只會使用檔案中格點數、間距、原點等數據。
  7. 程式設計詳FAQ
KV Method          |OB70
Minimum Kv        |2.
Projection        |LCP
Diag subgrid clouds|.true.
Start/end date    |10103120 10113023
MM5 output freq    |60
Grid time zone    |-8
CAMx grid size    |83,137,15
CAMx Dx, Dy        |3.,3.
CAMx orig & params |-124.500, -205.500, 120.99, 23.61, 10., 40.
Layer mapping      |-1,0, 1, 1, 2, 2, 3, 3, 4, 5, 6, 7, 8,10,23
CAMx lu file      |1011d43.lu
CAMx 3d file      |1011d43.3d
CAMx 2d file      |1011d43.2d
CAMx Kv file      |1011d43.kv
CAMx Cld/rain file |1011d43.cr
# MM5 files to proc| 1
MM5 filename      |./01/MMOUT_DOMAIN3_61

pt2em_d04.py（CAMx nc format）

  • 7版以後的CAMx 可以接受nc檔案，但其實點源nc檔案有2個版本：
    • 一個是用pncgen產生，維度和變數名稱仍然維持camx6版以前的協定，這個版本並不能用來執行任何模式，只是用來修改，此處定義為ver=6。
    • 一者是camx700新的維度和變數名稱，這個版本的nc檔案是繼續給camx700nc執行使用，有較高的穩定性，因此本程式以此版本為主。(ver=7)
    • 事實上本程式只有使用到點源的座標，33~37會先判斷nc檔案內的變數名稱為何，版本為何，再使用正確的變數名稱來引用數據。
    32  ver=7
    33  if 'XSTK' in Vt[0]:ver=6
    34  X={6:'XSTK',7:'xcoord'}
    35  Y={6:'YSTK',7:'ycoord'}
    36  IX=np.array([(i-nc.XORIG)/nc.XCELL for i in nct.variables[X[ver]][:]],dtype=int)
    37  IY=np.array([(i-nc.YORIG)/nc.XCELL for i in nct.variables[Y[ver]][:]],dtype=int)

pt2em_d01.py（CMAQ IOAPI nc format）

  • 雖然點源檔案與哪一層座標系統的解析度並沒有關聯，但是為了展示方便，還是以最完整的d1來承載排放量。如果要改變其他的網格系統，只要修改模版即可。
  • CMAQ的點源檔案有常數部分、有時間變化部分，因此兩個檔案都必須讀進來，前者是座標計算的依據，後者則是排放量部分。
  • 由於點源檔案是一維系統，整併進二維系統時切忌用多層次的if判別，會增加很多不必要的時間。
  • 此處先將J、I兩個值合併成一個數字，在一維情況下先做網格解析度的整併（line34~43)，將各組獨立位置的足標紀錄下來。
    • 先用d1的原點、網格間距，計算出點源在網格系統內的JI值配對。
    • pwrt 值是為了要讓J值有足夠的乘數可以同時紀錄這兩個整數。ncol如果是2位數，就乘上100。
    • 使用set技巧整併重複值
    • 用np.where找到重複位置的足標
    34    pwr=10**(int(np.log10(ncol))+1)
    35    ji=[int((y-nc.YORIG)/nc.YCELL)*pwr+int((x-nc.YORIG)/nc.XCELL) for x,y in zip(X,Y)]
    36    sji=set(ji)
    37    ji=np.array(ji)
    38    idx,idx0=[],[]
    39    for a in sji:
    40        j,i=int(a/pwr),int(a%pwr)
    41        if i>=ncol or j>=nrow or i<0 or j<0:continue
    42        idx.append(np.where(ji==a)[0])
    43        idx0.append(a)
  • 其次污染物種類的對照也是很大的問題，模版和點源檔案之間除了對照關係之外，還有很多污染物其實並沒有數值，可以規避。
  • 最後才做物種與時間的迴圈，針對前面紀錄的足標(idx)直接用np.sum進行加總，而不是累加，如此可以大大減省加總的時間。
    64    for ii in range(len(idx)):
    65        var1[:,:,ii]=np.sum(var0[:,:,idx[ii]],axis=2)
    66   
    67    for iv in range(nv):
    68        v=sint[iv]
    69        print(v)
    70        for ii in range(len(idx0)):
    71            a=idx0[ii]
    72            j,i=int(a/pwr),int(a%pwr)
    73            nc.variables[v][:,0,j,i]=var1[iv,:,ii]

Resources

  • Barron  Henderson, pseudonetcdf tutorial, http://www.barronh.com/software/tutorials/pseudonetcdf-tutorial
  • verdi
    • usage https://www.airqualitymodeling.org/index.php/VERDI_1.5_User_Manual#3.1_Installation_Instructions_for_Linux_and_Mac
    • VERDI使用說明 : http://www.evernote.com/l/AH3leuVQTuBEF7Vrs0D1C8Q-Iff5CpHl7eU
  • pt2emNest:https://github.com/sinotec2/CAMx_utility/blob/master/pt2emNest.f
  • Example https://github.com/sinotec2/cmaq_relatives/blob/master/run_cctmMM_RR_DM.csh
  • Notes
    • CAMx(UAM)的檔案格式, Yungchuan Kuang edited this page on 12 Jul 2016 · 2 revision, shttps://github.com/sinotec2/camxruns/wiki/CAMx(UAM)的檔案格式
    • Here: CMAQ/CAMx排放量檔案之轉換
    • Relatives:
      • D4範圍地面排放檔案之轉換
      • CMAQ compilations
      • CMAQ初始及邊界條件設定
      • brk_day2.cs(*)