REAS分區修正之程式說明

背景

REAS排放量含蓋中國大陸及其他國家，如何有效掌握其空間的特性，針對特定地區進行排放敏感性分析，有賴東亞(主要是中國大陸)[地理分區的網格遮罩(gridmask)檔案][withinD1]之應用。
mod11.py以大陸京津冀地區排放檔案的敏感性為範例，進行點源及面源排放之關閉與乘數應用。
雖然前述遮罩檔案是為CMAQ-ISAM使用，但同樣原理也可以應用在CAMx。此處應用對象主要是CAMx REAS的點源([[2022-07-06-REASPtRd]])及面源([[2022-07-06-REASgrnd]])檔案修正。

分區對照關係

分區的網格遮罩(gridmask)檔案
- fname='/'+hmp+'/cmaqruns/2016base/data/land/epic_festc1.4_20180516/gridmask/AQFZones_EAsia_81K.nc'
- 其產生方式可以參考[d01地理分區檔案之準備][withinD1]
分區代碼、乘數及名稱對照
- AQFZ:名稱與序號之對照
- R_dsgn：名稱與乘數之對照
- zone：所有受修正的空氣質量預報區代碼，其遮罩矩陣則讀成矩陣cnty_frac備用

#AQFZ0~7 0:sea_other_country/1:jinjinji/2:huabei/3:dongbei/4:xibei/5:huanan_taiwan/6:xinan/7:huadong
AQFZ={'SOC':0,'JJZ':1,'HB':2,'DB':3,'XB':4,'HN':5,'XN':6,'HD':7}
ic=0
dsgn=['JJZ']
R_dsgn={'JJZ':0.5}
zone=['AQFZ'+str(AQFZ[i]) for i in dsgn]
nv=len(zone)
izone=[V[3].index(v) for v in zone]
cnty_frac=np.zeros(shape=(nv,nrow,ncol))
cnty_frac[ic,:,:]=nc.variables[zone[ic]][0,0,:,:]

點源檔案之修正

CAMx點源與面源檔案最大的差異就是沒有4階矩陣(詳[點源nc檔案「排放量」數據之版本差異][CAMx67_timvar])
點源位置是在常數檔頭部分，CAMx6及CAMx7的座標變數名稱不同，此處為CAMx7版本，詳[點源nc檔案「煙囪參數」之版本差異][CAMx67_const]。
將網格化之點源位置之遮罩儲存在新的矩陣cnty_frac1備用
將所有排放量數據記存在矩陣var內

#ptse file
if len(V[3])==0 and len(V[1])>2:
  nt,nopts=(nc.variables[V[1][0]].shape[i] for i in range(2))
#determine which stacks remain in the domain
  IX=[int((x-XORIG)/XCELL) for x in nc.variables['xcoord'][:]]
  IY=[int((y-YORIG)/YCELL) for y in nc.variables['ycoord'][:]]
  cnty_frac1=np.zeros(shape=(nopts,nv))
  for i in range(nopts):
    if IY[i]*(IY[i]-nrow)<=0 and IX[i]*(IX[i]-ncol)<=0:
      cnty_frac1[i,:]=cnty_frac[:,IY[i],IX[i]]
  nv1=len(V[1])
  var=np.zeros(shape=(nv1,nt,nopts))
  v0=['CP_NO','plumerise']
#store the variable matrix
  for v in V[1]:
    if v in v0:continue
    iv=V[1].index(v)
    var[iv,:,:]=nc.variables[v][:,:]
    if np.sum(var[iv,:,:])==0:v0.append(v)
  nc.close()

開啟新檔案
將所有排放量乘上遮罩與乘數，可以得到指定地區減量後之排放量(var[iv,:,:]*(1-cnty_frac1[:,ic]*R_dsgn[dsgn[ic]]))

  nc = netCDF4.Dataset(fnameO, 'r+')
  for v in V[1]:
    if v in v0:continue
    iv=V[1].index(v)
    nc.variables[v][:,:]=var[iv,:,:]*(1-cnty_frac1[:,ic]*R_dsgn[dsgn[ic]])
  nc.close()

面源部分

將面源排放量存在var矩陣之內

#area files
else:
  nt,nlay,nrow,ncol=(nc.variables[V[3][0]].shape[i] for i in range(4))
  nv=len(V[3])
  var=np.zeros(shape=(nv,nt,nrow,ncol))
  v0=[]
  for v in V[3]:
    iv=V[3].index(v)
    var[iv,:,:,:]=nc.variables[v][:,0,:,:]
    if np.sum(var[iv,:,:,:])==0:v0.append(v)
  nc.close()

開啟新檔
cf為指定地區減量比例矩陣，相乘後存檔離開

  nc = netCDF4.Dataset(fnameO, 'r+')
  cf=1-cnty_frac[ic,:,:]*R_dsgn[dsgn[ic]]
  for v in V[3]:
    if v in v0:continue
    iv=V[3].index(v)
    nc.variables[v][:,0,:,:]=var[iv,:,:,:]*cf
  nc.close()

程式碼下載

Download: REAS排放分區修正之程式：mod11.py

系列其他程式

mod28-41.py程式將刪除位於(28,41)網格之點源排放量，以進行敏感性測試
- 以python語法此位置為cnty_frac[ic,40,27]

kuang@master /nas1/camxruns/2016_v7/ptse
$ diff mod11.py mod28_41.py
13,16c13,16
< ic=0
< dsgn=['JJZ']
< R_dsgn={'JJZ':0.5}
< zone=['AQFZ'+str(AQFZ[i]) for i in dsgn]
---
> ic=1
> dsgn=['JJZ','YP']
> R_dsgn={'JJZ':0.5,'YP':0.}
> zone=['AQFZ'+str(AQFZ[i]) for i in dsgn[:1]]+[0]
18c18
< izone=[V[3].index(v) for v in zone]
---
> #izone=[V[3].index(v) for v in zone]
20c20,21
< cnty_frac[ic,:,:]=nc.variables[zone[ic]][0,0,:,:]
---
> #nty_frac[ic,:,:]=nc.variables[zone[ic]][0,0,:,:]
> cnty_frac[ic,40,27]=1.

Download: 刪除特定網格位置REAS點源排放檔：mod28-41.py