pt_timvar程式說明

Table of contents

版本說明
- 粗細網格pt_timvar程式版本之比較
執行方式
TWN_3X3版本分段說明
EAsia_81K 版本差異
程式下載
Reference

有關pt_timvar的程式背景、策略構思、後處理及運作流程等等，詳見CMAQ點源變數檔案之準備。

版本說明

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

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

項目	EAsia_81K	TWN_3X3	說明
資料庫	REAS+TEDS	TEDS	前者範圍較大需要REAS資料庫
引數及主要檔	逐月CAMx高空點源nc檔案	逐月CAMx高空點源nc檔案	(same)
其他輸入檔	point_reas16’+mm+’.csv’, reas2cmaq.json, template.timvar.nc	template.timvar.nc	csv與json詳見https://sinotec2.github.io/FAQ/2022/07/06/REASPtRd.html
結果檔名	teds11.YYMM.timvar.nc	teds11.YYMM.timvar.nc	(same)
單月結果檔案容量	4.1G	3.9G

執行方式

目的：產生CMAQ所需的點源排放量中隨時間變化的部分，逐時數據
執行：
- 引數：逐月CAMx高空點源nc檔案(相對其他輸入檔，此檔為主要輸入檔)

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

MM=01~12
輸入檔(ptse檔案)：ptsEnHRBE.py結果
結果：teds11.YYMM.timvar.nc

TWN_3X3版本分段說明

模組及時間標籤轉換

副程式(dt2jul及jul2dt)

from pandas import *
import numpy as np
import netCDF4
import PseudoNetCDF as pnc
import datetime
import os,sys,json
import subprocess
     8
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)
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)
    19

讀入檔案名稱fname1

由其中讀取月份，計算月份的起迄日期
小時
- CAMx月份起始自前月最末日之8時(LST)
- CMAQ月份起始自0時(UTC)，亦即二者時間一致，因此不另調整。
日期
- CAMx月份起始自前月最末日
- CMAQ起始與結束日期則要看批次的定義，起始自前月15日後的16日，結束於第48日，共執行32天。

fname1=sys.argv[1]
mm=fname1[-5:-3]
mon=int(mm)
last=datetime.datetime(2019,mon,1)+datetime.timedelta(days=-1)
begP=last#+datetime.timedelta(days=0.5) #pt files begin last day 12:00(UTC)
begd=datetime.datetime(last.year,last.month,15)
start=begd+datetime.timedelta(days=4*(5-1))
idx=int((begP-start).seconds/3600)
enddt=begd+datetime.timedelta(days=4*12)
totalH=((enddt-start).days+1)*24+1

開啟CAMx點源處理結果

mpsp：部分CMAQ對照CAMx的名稱

#input the CAMx ptsource file
pt=netCDF4.Dataset(fname1, 'r')
v2=list(filter(lambda x:pt.variables[x].ndim==2, [i for i in pt.variables]))
nt,nopts=pt.variables[v2[0]].shape
mpsp={'PNA':'NA','POC':'POA','XYLMN':'XYL'}
    35
    36

定義輸出檔案名稱

使用ncks切割模版成為輸出檔並開啟之。
記錄全域屬性變數

#nc=/home/cmaqruns/2016_12SE1/emis/inln_point/ptnonipm/inln_mole_ptnonipm_20160701_12US1_cmaq_cb6_2016ff_16j.nc
#ncks -d TSTEP,0,0 -d ROW,1,25000 $nc template.timvar.nc
fname0='template.timvar.nc'
fname=fname1+'.timvar.nc'
path={'114-32-164-198.HINET-IP.hinet.net':'/opt/anaconda3/bin/', 'node03':'/usr/bin/','master':'/cluster/netcdf/bin/'}
path.update({'DEVP':'/usr/bin/','node01':'/cluster/netcdf/bin/','node02':'/cluster/netcdf/bin/'})
hname=subprocess.check_output('echo $HOSTNAME',shell=True).decode('utf8').strip('\n')
if hname not in path:
  sys.exit('wrong HOSTNAME')
if nopts>25000:sys.exit('nopts>25000')
os.system(path[hname]+'ncks -O -d ROW,1,'+str(nopts)+' '+fname0+' '+fname)
XCELL=pt.XCELL
XORIG=pt.XORIG
YORIG=pt.XORIG
V=[list(filter(lambda x:pt.variables[x].ndim==j, [i for i in pt.variables])) for j in [1,2,3,4]]
    52

開啟輸出檔

給定全域屬性變數

原點及點源座標會因網格系統的定義而改變。

nc = netCDF4.Dataset(fname,'r+')
jt=nt
v4=list(filter(lambda x:nc.variables[x].ndim==4, [i for i in nc.variables]))
bdate=jul2dt([pt.SDATE,pt.STIME])
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=dt2jul(start)[0]
nc.STIME=dt2jul(start)[1]
70

給定全月的時間標籤，藉此擴張檔案長度。

將所有變數矩陣全數歸0。

for t in range(totalH):
  dt=start+datetime.timedelta(hours=t)
  nc.variables['TFLAG'][t,:,0]=[dt2jul(dt)[0] for j in range(nc.NVARS)]
  nc.variables['TFLAG'][t,:,1]=[dt2jul(dt)[1] for j in range(nc.NVARS)]
    75
z=np.zeros(shape=nopts,dtype='float32')
for v in v4:
  for t in range(totalH):
    nc.variables[v][t,:,:,:]=z

給定CMAQ之PMOTHR排放量

給定其他特殊的粒狀物、VOC排放量(變數名稱對應不到的部分)

轉換單位由小時到秒

#CAMx_pt part
pmother=(pt.variables['CPRM'][:,:]+pt.variables['FPRM'][:,:])/ 3600.
nc.variables['PMOTHR'][idx:idx+jt,0,:nopts,0]=pmother
for v in mpsp: #mpsp={'PNA':'NA','POC':'POA','XYLMN':'XYL'}
  if mpsp[v] not in V[1]:continue
  nc.variables[v][idx:idx+jt,0,:nopts,0]=np.array(pt.variables[mpsp[v]][:,:nopts], dtype='float32')/ 3600.

將CAMx檔案內容倒入新檔

for v in v2: #ptse from CAMx
  if v not in v4:continue #outside of v4 (the template)
  print(v)
  var=pt.variables[v][:,:]/3600. #gmole/hr -> gmole/sec
  nc.variables[v][idx:idx+jt,0,:nopts,0]=var[:]

處理NOX的部分

-如果CMAQ要求的時間範圍較CAMx大，則以最靠近的日變化代入。

nox= nc.variables['NO2'][:,0,:,0]+nc.variables['NO'][:,0,:,0]
    92
nc.variables['NO2'][:,0,:,0]=nox *1./10.
nc.variables['NO'][:,0,:,0] =nox *9./10.
#np.zeros(shape=(jt,nopts),dtype='float32') #
    96
for v in v4:
  if np.sum(nc.variables[v][:])==0:continue
  if idx>0:
    for t in range(idx):
      nc.variables[v][t,0,:,0]=nc.variables[v][t+24,0,:,0]
  if totalH>idx+jt:
    for t in range(idx+jt,totalH):
      nc.variables[v][t,0,:,0]=nc.variables[v][t-24,0,:,0]
# print(v)
nc.close()

EAsia_81K 版本差異

讀取point_reas16MM.csv長度

$ diff ./REAS/pt_timvar.py ./twn/pt_timvarLL.py
...
28,30c29
< totalH=(enddt-start).days*24+1
< #input the results of pt_const.py fortBE.14.teds.baseMM
< len_df=int(subprocess.check_output('wc -l '+'point_reas16'+mm+".csv |awk '{print $1'}",shell=True).decode('utf8').strip('\n'))-1
---
> totalH=((enddt-start).days+1)*24+1

總點源個數與座標系統

50,53c46,52
< #if nopts+len(df)>10000:sys.exit('nopts+len(df)>10000')
< os.system(path[hname]+'ncks -O -d ROW,1,'+str(nopts+len_df)+' '+fname0+' '+fname)
< XCELL=81000.000
< XORIG=-2389500.
---
> if nopts>25000:sys.exit('nopts>25000')
> os.system(path[hname]+'ncks -O -d ROW,1,'+str(nopts)+' '+fname0+' '+fname)
> XCELL=pt.XCELL
> XORIG=pt.XORIG
> YORIG=pt.XORIG
> V=[list(filter(lambda x:pt.variables[x].ndim==j, [i for i in pt.variables])) for j in [1,2,3,4]]
57,58c56,58
< nc.NROWS=nopts+len_df
< nc.GDNAM='EAsia_81K'
---
> bdate=jul2dt([pt.SDATE,pt.STIME])
> nc.NROWS=nopts
> nc.GDNAM='TWN_3X3'

展開模版檔案

在nc檔案內使用迴圈非常費時，應儘量避免。
- 54～79是展開模版的歷程，無法避免。但是讓時間標籤先填入0，展開後再整批帶入，會比一一對照要快很多。
- 77～78開啟mask的過程，是最花時間的段落，因為nc基本上是個DICT、不適用整批帶入，沒辦法再更快。
- 關閉檔案再開啟：經測試並不會減省大量的記憶體。反而過多IO拖慢程式的速度
- 模版之展開（line 71～72）
- TFLAG填入日期時間（73～76）
  - 使用dt2jul函數，將datetime.timedelta計算結果轉成日期及時間陣列
  - 整理成nt2的陣列，轉到ntNVARS*2的大陣列
  - 一併填入TFLAG內
  - 如果有ETFLAG可以參考mod_tflag.py（適用uamiv檔案，如果沒有ETFLAG、aok等等很多程式會出錯）
- 防止陣列被mask鎖住（77~78）
- 儘量不要用2個nc檔案直接傳遞數值。

  nc = netCDF4.Dataset(fname,'r+')
...
  for t in range(totalH):
    nc.variables['TFLAG'][t,:,:]=0
  sdt=np.array([dt2jul(start+datetime.timedelta(days=t/24.)) for t in range(totalH)]).flatten().reshape(totalH,2)
  ARR=np.zeros(shape=(totalH,nc.NVARS,2))
  ARR[:,:,:]=sdt[:,None,:]
  nc.variables['TFLAG'][:,:,:]=ARR[:,:,:]
  for v in v4:
    nc.variables[v][:,0,:,0]=0. 
  print('save blank file')

mod_tflag.py之應用
- 本程式是為製作2019年各月CAMx模擬結果的模版，藉此執行aok以匯總EPA測站數據之ovm.csv/ovm.dat_camx，藉以快速繪製時序圖。
- d_hrs是TFLAG/ETFLA的差異，前者是模擬開始時間、後者是TFLAG+TSTEP，模擬結束時間
  - 此處TSTEP是1小時
- 應用前述批次倒入TFLAG之作法

REAS 排放對照及轉檔

讀取csv及json檔案
先跳過NOx隨後再處理
重複使用ARR：並使用空維度None，對速度提升與降低記憶體使用都很有效
執行完畢將變數歸0，對減少記憶體也不無幫助（98～99）

  #CAMx_pt part
  mpsp={'PNA':'NA','POC':'POA','XYLMN':'XYL'}
  ARR=(pt.variables['CPRM'][:,:]+pt.variables['FPRM'][:,:])/ 3600.
  nc.variables['PMOTHR'][idx:idx+jt,0,:nopts,0]=ARR[:,:]
  with_data=['PMOTHR']
  
  for v in mpsp: #mpsp={'PNA':'NA','POC':'POA','XYLMN':'XYL'}
    if v not in mpsp or mpsp[v] not in v2:continue
    ARR[:,:]=np.array(pt.variables[mpsp[v]][:,:nopts], dtype='float32')/ 3600.
    nc.variables[v][idx:idx+jt,0,:nopts,0]=ARR[:,:]
    with_data.append(v)
  
  for v in v2: #ptse from CAMx
    if v not in v4:continue #outside of v4 (the template)
    ARR[:,:]=np.array(pt.variables[v][:,:nopts],dtype='float32') / 3600. #gmole/hr -> gmole/sec
    nc.variables[v][idx:idx+jt,0,:nopts,0]=ARR[:,:]
    with_data.append(v)
  pt=0.
  ARR=0.

這部分容量並不增加、計算速度也不低，

  #REAS part
  print('REAS part')
  df=read_csv('point_reas16'+mm+'.csv')
  specs=[ i for i in df.columns if i not in ['lon','lat','x_m','y_m']]
  ARR=np.zeros(shape=(jt,len_df))
  with open('reas2cmaq.json','r') as f:
    r2c=json.load(f)
  for v in specs:
    if v == 'NOX':continue
    if sum(np.array(df[v]))==0.:continue
    if v not in r2c :continue
    if type(r2c[v])==str: #simple mapping of specs
      if r2c[v] not in v4:continue
      a=np.array(df[v])
      ARR[:,:]=a[None,:]
      nc.variables[r2c[v]][idx:idx+jt,0,nopts:,0]=ARR[:,:]
      with_data.append(r2c[v])

CBM成分累加部分（117~133）與nc檔案脫鉤，在時間迴圈之外執行，這樣才有效率

  cbm=set()
  sv=[]
  for v in specs:
    if v not in r2c :continue
    if type(r2c[v])!=str:
      sv.append(v)
      cbm=cbm.union(set(r2c[v]))
  cbm=list(cbm)
  ARR=np.zeros(shape=(len(cbm),len_df))
  for v in sv:
    for i in r2c[v]:
      ARR[cbm.index(i),:]+=np.array(df[v])
  Acbm=np.zeros(shape=(jt,len(cbm),len_df))
  Acbm[:,:,:]=ARR[None,:,:]
  for v in cbm:
    nc.variables[v][idx:idx+jt,0,nopts:,0]=Acbm[:,cbm.index(v),:]
    with_data.append(v)

REAS氮氧化物排放量之讀取

144,147d92
< ARR=np.zeros(shape=(jt,len_df))
< a=np.array(df['NOX'])
< ARR[:,:]=a[None,:]
< nox[idx:idx+jt,nopts:]=ARR[:,:]

前一天排放補值

CAMx 是逐月檔，但在CMAQ是跟著WRF的批次，因此可能在每月的前後會需要多一點時間（小月問題），過去以手動方式加值，現則直接以程式添加。

  for v in v4:
    if v not in with_data:continue
    a=nc.variables[v][24:idx+24,0,:,0]
    nc.variables[v][:idx,0,:,0]=a[:,:]
    a=nc.variables[v][idx-jt-24:totalH-24,0,:,0]
    nc.variables[v][idx-jt:totalH,0,:,0]=a[:,:]
  
  nox= nc.variables['NO2'][:,0,:,0]+nc.variables['NO'][:,0,:,0]
  
  ARR=np.zeros(shape=(jt,len_df))
  a=np.array(df['NOX'])
  ARR[:,:]=a[None,:]
  nox[idx:idx+jt,nopts:]=ARR[:,:]
  nc.variables['NO2'][:,0,:,0]=nox *1./10.
  nc.variables['NO'][:,0,:,0] =nox *9./10.
  # print(v)
  nc.close()