pt_const程式說明

說明

整體工作流程詳見CMAQ點源常數檔案之準備([[2022-07-06-pt_const]])
目的：產生CMAQ所需的點源排放量中常數、不隨時間變化的部分
相關程式見[[2022-07-07-pt_timvar.py]]及[[2022-07-06-pt_timvar]]說明

執行

引數：CAMx分月高空點源，MM=01~12

for MM in {01..12};
do python pt_constLL.py fortBE.413_teds11.ptse$MM.nc
done

輸入檔(ptse檔案)：ptsEnHRBE.py港區計算結果（[[2022-07-15-ptsEnHRBE.py程式說明]]）
結果：teds11.YYMM.const.nc

版本

點源座標系統為LL
版本：
- TWN_3X3範圍、資料來源為TEDS資料庫，點源座標系統為LL
- EAsia_81K範圍，除了TEDS之外，還讀取REAS處理結果

粗細網格pt_const程式版本之比較

項目	EAsia_81K	TWN_3X3	說明
資料庫	REAS+TEDS	TEDS	前者範圍較大需要REAS資料庫
引數及主要檔	逐月CAMx高空點源nc檔案	逐月CAMx高空點源nc檔案	(same)
其他輸入檔	point_SPEC.csv, ^*nc檔模版	^*nc檔模版	csv詳見REASv3.1點源之讀取、格式轉換與合併
結果檔名	teds11.YYMM.const.nc, point_reasYYMM.csv(讓pt_timvar.py讀取)	teds11.YYMM.const.nc	(same)

^*nc檔模版來自於USEPA之範例檔。

分段說明

模組

時間標籤轉換副程式(dt2jul及jul2dt)

import numpy as np
import netCDF4
import PseudoNetCDF as pnc
import os,sys,twd97
from pyproj import Proj
import subprocess
from pandas import *
import datetime
     9
def dt2jul(dt):
  yr=dt.year
  deltaT=dt-datetime.datetime(yr,1,1)
  deltaH=int((deltaT.total_seconds()-deltaT.days*24*3600)/3600.)
  return (yr*1000+deltaT.days+1,deltaH*10000)
    15
def jul2dt(jultm):
  jul,tm=jultm[:]
  yr=int(jul/1000)
  ih=int(tm/10000.)
  return datetime.datetime(yr,1,1)+datetime.timedelta(days=int(jul-yr*1000-1))+datetime.timedelta(hours=ih)
    21
    22

讀入檔案名稱的月份
準備座標轉換之套件

mon=int(sys.argv[1][-5:-3])
#join the pollutants for this month
pth='/nas1/cmaqruns/2016base/data/ptse/'
#prepare interceptions
Latitude_Pole, Longitude_Pole = 23.61000, 120.9900
pnyc = Proj(proj='lcc', datum='NAD83', lat_1=10, lat_2=40, lat_0=Latitude_Pole, lon_0=Longitude_Pole, x_0=0, y_0=0.0)
    29

讀入CAMx點源檔案

讀取各維度變數名稱
讀取全域變數

fname1=sys.argv[1]
pt=netCDF4.Dataset(fname1, 'r')
v3=list(filter(lambda x:pt.variables[x].ndim==3, [i for i in pt.variables]))
v2=list(filter(lambda x:pt.variables[x].ndim==2, [i for i in pt.variables]))
v1=list(filter(lambda x:pt.variables[x].ndim==1, [i for i in pt.variables]))
nhr,nvar,dt=pt.variables[v3[0]].shape
nt,nopts=pt.variables[v2[0]].shape
XCELL=pt.XCELL
XORIG=pt.XORIG
YORIG=pt.XORIG
    40

讀取起始時間(0800 LST)
準備ncks程式

tb=pt.STIME-80000 #UTC
bdate=jul2dt([pt.SDATE,pt.STIME])
fname1=fname1.replace('fortBE.413_','').replace('ptse','19').replace('.nc','')
fname=fname1+'.const.nc'
fname0=pth+'stack_groups_ptnonipm_12US1_2016ff_16j.nc' #as template
    46
#ncks path
path={'114-32-164-198.HINET-IP.hinet.net':'/opt/anaconda3/bin/', 'node03':'/usr/bin/', \
        'master':'/cluster/netcdf/bin/','DEVP':'/usr/bin/'}
hname=subprocess.check_output('echo $HOSTNAME',shell=True).decode('utf8').strip('\n')
if hname not in path:
  sys.exit('wrong HOSTNAME')
ncks=path[hname]+'ncks'

由既有模版使用ncks切出所要的長度，成為輸出檔案的底版。

x=list(pt.variables['xcoord'][:nopts])
y=list(pt.variables['ycoord'][:nopts])
lon, lat = pnyc(x, y, inverse=True)
    57
os.system(ncks+' -O -d ROW,1,'+str(nopts)+' '+fname0+' '+fname)

開啟輸出檔，寫入全域變數。

nc = netCDF4.Dataset(fname,'r+')
nc.NROWS=nopts
nc.GDNAM='TWN_3X3'
nc.P_ALP = np.array(10.)
nc.P_BET = np.array(40.)
nc.P_GAM = np.array(120.98999786377)
nc.XCENT = np.array(120.98999786377)
nc.YCENT = np.array(23.6100196838379)
nc.XCELL=XCELL #27000.000
nc.YCELL=XCELL #27000.000
nc.XORIG=XORIG #-877500.0
nc.YORIG=XORIG #-877500.0
nc.SDATE,nc.STIME=dt2jul(bdate+datetime.timedelta(hours=-8))
    72

參數名稱、單位、其他對照與轉換

點源變數名稱：CMAQ對照到CAMx
單位轉換：每小時->每秒
給定經緯度、網格位置、時間標籤等

mp={'STKDM':'stkdiam','STKHT':'stkheight','STKTK':'stktemp','STKVE':'stkspeed','XLOCA':'xcoord', 'YLOCA':'ycoord',}
    74
for v in mp:
  nc.variables[v][0,0,:,0]=pt.variables[mp[v]][:]
#velocity is m/hr in CAMx pt, but in m/s in CMAQ_pt
nc.variables['STKVE'][0,0,:,0]=pt.variables[mp[v]][:]/3600.
nc.variables['IFIP'][0,0,:,0]=[1000+i for i in range(nopts)]
nc.variables['ISTACK'][0,0,:,0]=[1+i for i in range(nopts)]
    81
nc.variables['COL'][0,0,:,0]=[int((i-nc.XORIG)/nc.XCELL) for i in x]
nc.variables['ROW'][0,0,:,0]=[int((i-nc.YORIG)/nc.YCELL) for i in y]
nc.variables['TFLAG'][0,:,0]=[nc.SDATE for i in range(nc.NVARS)]
nc.variables['TFLAG'][0,:,1]=[nc.STIME for i in range(nc.NVARS)]
nc.variables['LATITUDE'][0,0,:,0]=lat
nc.variables['LONGITUDE'][0,0,:,0]=lon
nc.close()

EAsia_81K 版本差異

輸入全年分項檔案

讀出其中的座標以進行轉換
讀出當月數據並寫成檔案以供pt_timvar.py使用

28a16,55
> XCELL=81000.000
> XORIG=-2389500.
>
> dfa=DataFrame({})
> for fn in fnames:
>   spec=fn.split('/')[-1].replace('point_','').replace('.csv','').replace('./','')
>   df=read_csv(fn)
>   #interceptions
>   lat,lon=(list(df.lat),list(df.lon))
>   x,y=pnyc(lon, lat, inverse=False)
>   df['x_m'],df['y_m']=(x,y)
>   boo=(df['x_m']>=XORIG) & (df['x_m']<=-XORIG) & (df['y_m']>=XORIG) & (df['y_m']<=-XORIG)
>   df=df.loc[boo].reset_index(drop=True)
>   #reduction
>   df.lat=[round(i,1) for i in df.lat]
>   df.lon=[round(i,1) for i in df.lon]
>   c=df.iloc[:,[0,1,mon+1]]
>   c.columns=['lon','lat',spec]
>   c=pivot_table(c,index=['lon','lat'],values=spec,aggfunc=sum).reset_index()
>   if fn==fnames[0]:
>     dfa=c
>   else:
>     dfa=merge(dfa,c,on=['lon','lat'],how='outer')
> df=dfa.fillna(0)
>
> fac=1000.*1000./monthrange(2015,mon)[1]/24/3600.
> mws={'BC':1,'BENZENE':78.11,'BUTANES':58.12,'CO':28, 'CO2':44, 'ETHANE':30.07,
>        'FORMALDEHYDE':30.031, 'INTERNAL_ALKENES':28.05, 'KETONES':86.13, 'NH3':17, 'NOX':46, 'OC':1,
>        'OTHER_ALKANES':58.12, 'OTHER_AROMATICS':106.16, 'PENTANES':72.15, 'PM10':1, 'PM25':1,
>        'PROPANE':44.1 , 'SO2':64, 'TERMINAL_ALKENES':28.05, 'TOLUENE':92.14, 'TOTAL_NMV':58.12,
>        'XYLENES':106.16}
>
> #units of gmole/s or gram/s
> for c in [ i for i in df.columns if i not in ['lon','lat']]:
>   if c not in mws:continue
>   df[c]=np.array(df[c])*fac/mws[c]
> lat,lon=(list(df.lat),list(df.lon))
> x,y=pnyc(lon, lat, inverse=False)
> df['x_m'],df['y_m']=(x,y)
> df.set_index('lon').to_csv('point_reas16'+sys.argv[1][-5:-3]+'.csv')

在TEDS之後再加上REAS數據

53,55c76,78
< ncks=path[hname]+'ncks'
< x=list(pt.variables['xcoord'][:nopts])
< y=list(pt.variables['ycoord'][:nopts])
---
>
> x=list(pt.variables['xcoord'][:nopts])+list(df.x_m)
> y=list(pt.variables['ycoord'][:nopts])+list(df.y_m)
58c81
< os.system(ncks+' -O -d ROW,1,'+str(nopts)+' '+fname0+' '+fname)
---
> os.system(path[hname]+'ncks -O -d ROW,1,'+str(nopts+len(df))+' '+fname0+' '+fname)
60,61c83,84
< nc.NROWS=nopts
< nc.GDNAM='TWN_3X3'
---
> nc.NROWS=nopts+len(df)
> nc.GDNAM='EAsia_81K'

填入REAS煙囪參數

73a97
> va={'STKDM':17.0  ,'STKHT':250   ,'STKTK':373.  ,'STKVE':19.0}
75,76d98
< for v in mp:
<   nc.variables[v][0,0,:,0]=pt.variables[mp[v]][:]
78,80c100,122
< nc.variables['STKVE'][0,0,:,0]=pt.variables[mp[v]][:]/3600.
< nc.variables['IFIP'][0,0,:,0]=[1000+i for i in range(nopts)]
< nc.variables['ISTACK'][0,0,:,0]=[1+i for i in range(nopts)]
---
> v='STKVE'
> STKVE=list(pt.variables[mp[v]][:]/3600.)
>
> for i in va:
>   val=va[i]
>   va.update({i:[val for i in range(len(df))]})
> va.update({'XLOCA':list(df.x_m),'YLOCA':list(df.y_m)})
> nc.variables['LMAJOR'][0,0,:,0]=[0 for i in range(nopts+len(df))]
> nc.variables['LPING'][0,0,:,0]=[0 for i in range(nopts+len(df))]
> for i in range(nopts+len(df)):
>   HSTK=250.
>   if i<nopts:
>     HSTK=pt.variables['stkheight'][i]
>   if HSTK>150.:
>     nc.variables['LPING'][0,0,i,0]=1
>     nc.variables['LMAJOR'][0,0,i,0]=1
> for v in mp:
>   if v=='STKVE' :continue
>   nc.variables[v][0,0,:,0]=np.array(list(pt.variables[mp[v]][:])+va[v])
> v='STKVE'
> nc.variables[v][0,0,:,0]=np.array(STKVE+va[v])
> nc.variables['IFIP'][0,0,:,0]=[1000+i for i in range(nopts+len(df))]
> nc.variables['ISTACK'][0,0,:,0]=[1+i for i in range(nopts+len(df))]

程式下載

Download: EAsia_81K版本：pt_const.py

Download: TWN_3X3版本：pt_constLL.py

說明

執行

版本

粗細網格pt_const程式版本之比較

分段說明

模組

讀入CAMx點源檔案

參數名稱、單位、其他對照與轉換

EAsia_81K 版本差異

輸入全年分項檔案

在TEDS之後再加上REAS數據

填入REAS煙囪參數

程式下載

Reference