CAMx高空點源排放檔案之轉寫

Table of contents

背景
程式說明
- 程式分段說明
檔案下載
Reference

背景

因程式記憶體耗用太大(多達150G、DEVP之50%)，nc檔案在輸出時遭到困難。最終決定先寫成.fth檔案，同時也有提供中間檢查、確認的方便性。
feather檔案可以說是目前壓縮比例最高、存取速度最快的資料庫檔案形式，詳情可以參考知乎。

程式說明

程式分段說明

調用pandas的版本必須含有pyarrow，會連動到numpy、xarray的版本
str2lst：netCDF檔案內的字串內容必須是bytes格式，一般的字串要經過轉換。

$ cat -n wrtE.py
'''
Purpose: Generate CAMx Elev. PtSe. NC file from dfMM.fth (MM=01~12)
Usage: python wrtE.py YYMM
see descriptions at https://sinotec2.github.io/Focus-on-Air-Quality/EmisProc/ptse/wrtE/
'''
#! crding = utf8
from pandas import *
import numpy as np
import os, sys, subprocess
import netCDF4
import datetime
from calendar import monthrange
    13
from ioapi_dates import dt2jul
if not sys.warnoptions:
    import warnings
    warnings.simplefilter("ignore")
    18
    19
def str2lst(A):
    return [bytes(i,encoding='utf-8') for i in A[1:-1].replace("b'","").replace("'","").replace(" ","").split(',')][:8]
    22
    23

外部程式之相依性

#Main
#locate the programs and root directory
ncks=subprocess.check_output('which ncks',shell=True).decode('utf8').strip('\n')
P='./'
    28

讀取引數(4碼年月)、定義起迄時間
- CAMx的時間是LST，因此自上月最末日的早上8時開始。以同時符合CMAQ的要求(自00UTC開始)。

#time and space initiates
ym=sys.argv[1]
mm=ym[2:4]
mo=int(mm)
yr=2000+int(ym[:2])
ntm=(monthrange(yr,mo)[1]+2)*24+1
bdate=datetime.datetime(yr,mo,1)+datetime.timedelta(days=-1+8./24)
edate=bdate+datetime.timedelta(days=ntm/24)#monthrange(yr,mo)[1]+3)

讀取模版並進行時間軸的延長
- 一般nc檔案的矩陣會自動以masked arraynumpy.ma.array型式儲存，模版內容如果被遮罩遮蔽了，延長放大之後會是個災難。有關模版的mask array的檢查與修正見另文

#prepare the uamiv template
print('template applied')
NCfname='fortBE.413_teds10.ptsE'+mm+'.nc'
try:
  nc = netCDF4.Dataset(NCfname, 'r+')
except:
  os.system('cp '+P+'template_v7.nc '+NCfname)
  nc = netCDF4.Dataset(NCfname, '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,nv,dt=nc.variables[V[2][0]].shape
nv=len([i for i in V[1] if i !='CP_NO'])
nc.SDATE,nc.STIME=dt2jul(bdate)
nc.EDATE,nc.ETIME=dt2jul(edate)
nc.NOTE='Point Emission'
nc.NOTE=nc.NOTE+(60-len(nc.NOTE))*' '
nc.NVARS=nv

注意name 屬性之調用，可能會遭遇問題。

#Name-names may encounter conflicts with newer versions of NCFs and PseudoNetCDFs.
nc.name='PTSOURCE  '
nc.NSTEPS=ntm
if 'ETFLAG' not in V[2]:
  zz=nc.createVariable('ETFLAG',"i4",("TSTEP","VAR","DATE-TIME"))
if nt!=ntm or (nc.variables['TFLAG'][0,0,0]!=nc.SDATE and nc.variables['TFLAG'][0,0,1]!=nc.STIME):
  for t in range(ntm):
    sdate,stime=dt2jul(bdate+datetime.timedelta(days=t/24.))
    nc.variables['TFLAG'][t,:,0]=[sdate for i in range(nv)]
    nc.variables['TFLAG'][t,:,1]=[stime for i in range(nv)]
    ndate,ntime=dt2jul(bdate+datetime.timedelta(days=(t+1)/24.))
    nc.variables['ETFLAG'][t,:,0]=[ndate for i in range(nv)]
    nc.variables['ETFLAG'][t,:,1]=[ntime for i in range(nv)]
for c in V[1]:
  nc.variables[c].set_auto_mask(False)
  if c=='CP_NO': continue
  nc.variables[c][:]=0.
nc.close()
#template OK
    72

讀取逐時、逐煙道之排放量.fth檔案
讀取煙囪條件

df=read_feather('df'+mm+'.fth')
pv=read_csv('pv'+mm+'.csv')
pv.CP_NOb=[str2lst(i) for i in pv.CP_NOb]
nopts=len(pv)
#item sets definitions
XYHDTV=['UTM_E','UTM_N','HEI','DIA','TEMP','VEL','ORI_QU1']
    79

CAMx版本在此設定

#determination of camx version
ver=7
if 'XSTK' in V[0]:ver=6
col_mn=['ORI_QU1','TEMP','VEL','UTM_E', 'UTM_N','HY1','HD1','DY1']
lspec=[i for i in df.columns if i not in col_mn+['index','C_NO']]
    85
dimn={6:'NSTK',7:'COL'}

為延長煙道數(COL軸)，須先將模版的COL軸設為可增加之記錄軸(--mk_rec_dmn)

print(dimn[ver]+' expanding and reopening')
res=os.system(ncks+' -O --mk_rec_dmn '+dimn[ver]+' '+NCfname+' tmp'+mm)
if res!=0: sys.exit(ncks+' fail')
res=os.system('mv tmp'+mm+' '+NCfname)
if res!=0: sys.exit('mv fail')

準備煙囪參數的對照表
- CAMx 6與CAMx 7的煙囪參數變數名稱有很大的改變
- CAMx設有開啟pig的機制。對CMAQ沒有效果。

#CP_NO in S1(byte) format
print('ncfile Enlargement')
#prepare the parameter dicts
PRM='XYHDTV'
v2n={PRM[i]:XYHDTV[i] for i in range(6)}
names={7:['xcoord','ycoord','stkheight','stkdiam','stktemp','stkspeed'],
       6:[v+'STK' for v in PRM]}
v2c={PRM[i]:names[ver][i] for i in range(6)}
a=DataFrame({'SN':df.SO2+df.NO2})
a=a.sort_values('SN',ascending=False)
pig=[]#a.index[:100]
#filling the stack parameters for camx700nc
   104

再次開啟模版。關閉再開的理由是節省不必要的暫存記憶體。
- set_auto_mask(False)或set_always_mask(False)的原因，是因為模版中有被遮蔽的矩陣內容。須先在模版階段解決。

nc = netCDF4.Dataset(NCfname, 'r+')
#enlarge the record dimension (COL)
z=np.zeros(shape=ntm)
for c in V[1]:
  nc.variables[c].set_auto_mask(False)
  if c in ['CP_NO']:continue
  for i in range(nopts):
    nc.variables[c][:ntm,i]=z
if ver==7:nc.variables['pigflag'][:nopts]=0
nc.close()

如不再增加時間的長度，不必再將記錄軸設回來。

#res=os.system(ncks+' -O --mk_rec_dmn TSTEP '+NCfname+' tmp'+mm)
#if res!=0: sys.exit(ncks+' fail')
#res=os.system('mv tmp'+mm+' '+NCfname)
#if res!=0: sys.exit('mv fail')
   119

寫入煙囪參數，CAMx的單位是每小時

nc = netCDF4.Dataset(NCfname, 'r+')
for v in PRM:
  var=v2c[v]
  nc.variables[var].set_auto_mask(False)
  nc.variables[var][:nopts]=np.array(pv[v2n[v]])
nc.variables[v2c['V']][:nopts]=nc.variables[v2c['V']][:]*3600.
nc.variables[v2c['T']][:nopts]=nc.variables[v2c['T']][:]+273.

開啟pig設定

#first 100 for PiG
if len(pig)>0:
  if ver==7:
    nc.variables['pigflag'][pig]=1
  else:
    nc.variables[v2c['D']][pig]=nc.variables[v2c['D']][pig]*-1.

寫入每污染物之排放量

for c in V[1]:
  if c not in lspec:continue
  if c not in df.columns:continue
  if c in ['CP_NO']: continue
  ic=lspec.index(c)
  nc.variables[c][:,:nopts]=np.array(df[c]).reshape(ntm,nopts)
  print(c)

將管煙編號(=管編+煙道編號)寫進檔案中以備後續增、減排放源使用。
- CMAQ雖然有FIP編號可供追蹤，但因FIP編號是整數，與國內的編碼系統不符合，且為管編，不能指定到個別煙道，因此建立自己的編號系統是必須的。
- 且因非CAMx或CMAQ控制之變數名稱，程式會跳開不讀，並不會報錯。

nc.variables['CP_NO'][:nopts,:8]=np.array(list(pv.CP_NOb)).flatten().reshape(nopts,8)
nox=nc.variables['NO2'][:,:nopts]
nc.variables['NO'][:,:nopts]=nox[:,:nopts]*0.9
nc.variables['NO2'][:,:nopts]=nox-nc.variables['NO'][:,:nopts]
nc.NOPTS=nopts
nc.close()
   146

檔案下載

Download: python程式：wrtE.py

Reference

The NumPy community, numpy.ma.array, Jun 22, 2021
数据如琥珀, 轻如“鸿毛（Feather）”的文件格式却重于泰山, 知乎, 2020-09-16