EAC4檔案轉成5階m3.nc

Table of contents

背景
逐日密度檔案之準備
grb2D1m3.py程式說明
結果檢視
程式下載
Reference

背景

歐洲中期天氣預報中心(ECMWF)之EAC4 (ECMWF Atmospheric Composition Reanalysis 4)數據下載整併後，此處將其轉成m3.nc檔案，以供VERDI等顯示軟體、以及後續光化模式所需。
由於EAC4粒狀物單位(重量混合比)轉換過程需要大氣的密度，不單是高度的函數，也隨著天氣系統而有時空的變化。可以由mcip計算結果（METCRO3D）中讀取，需在執行轉換前預備好。

逐日密度檔案之準備

由於此處EAC4檔案的時間範圍為全月，而mcip結果是彼此會有重疊的批次作業，因此需以brk_day.cs拆解、讓程式可以逐日讀取，以降低複雜度。程序如下：
- 先以ncks讀取METCRO3D檔案中的密度(DENS)及時間標籤
- 再以brk_day.cs拆解成逐日檔案備用
- 密度單位kg/M³，此處未更動，在主程式中進行轉換。

for r in {5..12};do 
  nc=/nas1/cmaqruns/2018base/data/mcip/1804_run$r/sChina_81ki/METCRO3D_1804_run$r.nc
  ncks -O -v TFLAG,DENS $nc RHO.1804.nc;
  brk_day2.cs RHO.1804.nc
done
ls 1804
RHO.20180331.nc  RHO.20180405.nc  RHO.20180410.nc  RHO.20180415.nc  RHO.20180420.nc  RHO.20180425.nc  RHO.20180430.nc  RHO.20211222.nc
...
RHO.20180404.nc  RHO.20180409.nc  RHO.20180414.nc  RHO.20180419.nc  RHO.20180424.nc  RHO.20180429.nc  RHO.20211221.nc

Note:brk_day2.cs的引數必須以YYMM做為主檔名的最後標籤(如範例中的RHO.1804.nc)。

grb2D1m3.py程式說明

引數

grb2D1m3.py YYMM.nc
- 為EAC4下載、轉檔、合併之結果

I/O檔案

Inputs
- YYMM.nc
- BCON模版
- m3.nc模版
- 逐日密度檔案
JSON’s
- dic.json={eac4污染物代碼：污染物名稱}
- mws.json={污染物名稱：分子量}
- nms_gas.json：{eac4污染物代碼：CMAQ污染物名稱}
- nms_part.json：{eac4污染物代碼：CMAQ污染物名稱序列}
Output
- YYMMD1.m3.nc

dic.json

{ 
  "VAR_192_210_123_P0_L105_GLL0":"carbon_monoxide", 
  "VAR_192_217_45_P0_L105_GLL0":"ethane", 
  "VAR_192_210_124_P0_L105_GLL0":"formaldehyde", 
  "VAR_192_217_16_P0_L105_GLL0":"isoprene", 
  "VAR_192_210_121_P0_L105_GLL0":"nitrogen_dioxide", 
  "VAR_192_217_27_P0_L105_GLL0":"nitrogen_monoxide", 
  "VAR_192_217_11_P0_L105_GLL0":"olefins", 
  "VAR_192_217_6_P0_L105_GLL0":"nitric_acid", 
  "VAR_192_217_15_P0_L105_GLL0":"organic_nitrates", 
  "VAR_192_210_203_P0_L105_GLL0":"ozone", 
  "VAR_192_217_9_P0_L105_GLL0":"paraffins", 
  "VAR_192_217_13_P0_L105_GLL0":"peroxyacetyl_nitrate", 
  "VAR_192_217_47_P0_L105_GLL0":"propane", 
  "VAR_192_210_122_P0_L105_GLL0":"sulphur_dioxide"
  "MASSMR_P40_L105_GLL0_A0":"ozone",
  "MASSMR_P40_L105_GLL0_A4":"carbon_monoxide",
  "MASSMR_P40_L105_GLL0_A5":"nitrogen_dioxide",
  "MASSMR_P40_L105_GLL0_A8":"sulphur_dioxide",
  "MASSMR_P40_L105_GLL0_A11":"nitrogen_monoxide",

  "VAR_192_217_21_P0_L105_GLL0":"ammonium", 
  "VAR_192_210_4_P0_L105_GLL0":"dust_aerosol_0.03-0.55um_mixing_ratio", 
  "VAR_192_210_5_P0_L105_GLL0":"dust_aerosol_0.55-0.9um_mixing_ratio", 
  "VAR_192_210_6_P0_L105_GLL0":"dust_aerosol_0.9-20um_mixing_ratio", 
  "VAR_192_210_9_P0_L105_GLL0":"hydrophilic_black_carbon_aerosol_mixing_ratio", 
  "VAR_192_210_7_P0_L105_GLL0":"hydrophilic_organic_matter_aerosol_mixing_ratio", 
  "VAR_192_210_10_P0_L105_GLL0":"hydrophobic_black_carbon_aerosol_mixing_ratio", 
  "VAR_192_210_8_P0_L105_GLL0":"hydrophobic_organic_matter_aerosol_mixing_ratio", 
  "VAR_192_217_51_P0_L105_GLL0":"nitrate", 
  "VAR_192_210_1_P0_L105_GLL0":"sea_salt_aerosol_0.03-0.5um_mixing_ratio", 
  "VAR_192_210_2_P0_L105_GLL0":"sea_salt_aerosol_0.5-5um_mixing_ratio", 
  "VAR_192_210_3_P0_L105_GLL0":"sea_salt_aerosol_5-20um_mixing_ratio", 
  "VAR_192_210_11_P0_L105_GLL0":"sulphate_aerosol_mixing_ratio", 
  "MASSMR_P48_L105_GLL0_A62003":"ammonium_aerosol_mass_mixing_ratio",
  "MASSMR_P48_L105_GLL0_A65533":"nitrate_coarse_mode_aerosol_mass_mixing_ratio",
  "MASSMR_P48_L105_GLL0_A65534":"nitrate_fine_mode_aerosol_mass_mixing_ratio",
  }

nms_part.json

dust_aerosol_0.03-0.55um為塵土中的PM₁，但是CMAQ並沒有處理這一塊，為求質量守恆，將其併入PM_2.5之中(J類別)

{
"VAR_192_217_21_P0_L105_GLL0": ["ANH4I", "ANH4J", "ANH4K"], 
"VAR_192_210_4_P0_L105_GLL0": ["AFEJ", "AALJ", "ASIJ", "ACAJ", "AMGJ", "AKJ", "AMNJ"], 
"VAR_192_210_5_P0_L105_GLL0": ["AFEJ", "AALJ", "ASIJ", "ACAJ", "AMGJ", "AKJ", "AMNJ"], 
"VAR_192_210_6_P0_L105_GLL0": ["ACORS", "ASOIL"], 
"VAR_192_210_9_P0_L105_GLL0": ["AECI", "AECJ"], 
"VAR_192_210_7_P0_L105_GLL0": ["APOCI", "APNCOMI", "APOCJ", "AOTHRJ", "AISO3J", "ASQTJ", "AORGCJ", "AOLGBJ", "AOLGAJ"], 
"VAR_192_210_10_P0_L105_GLL0": ["AECI", "AECJ"], 
"VAR_192_210_8_P0_L105_GLL0": ["APOCI", "APNCOMI", "APOCJ", "AOTHRJ", "AISO3J", "ASQTJ", "AORGCJ", "AOLGBJ", "AOLGAJ"], 
"VAR_192_217_51_P0_L105_GLL0": ["ANO3I", "ANO3J", "ANO3K"], 
"VAR_192_210_1_P0_L105_GLL0": ["ANAI", "ACLI", "ANAJ", "ACLJ"], 
"VAR_192_210_2_P0_L105_GLL0": ["ACLK", "ASEACAT"], 
"VAR_192_210_3_P0_L105_GLL0": [], 
"VAR_192_210_11_P0_L105_GLL0": ["ASO4I", "ASO4J", "ASO4K"],
"MASSMR_P48_L105_GLL0_A62003": ["ANH4I", "ANH4J", "ANH4K"],
"MASSMR_P48_L105_GLL0_A65533": ["ANO3K"],
"MASSMR_P48_L105_GLL0_A65534": ["ANO3I", "ANO3J"],
}

分段說明

使用scipy的griddata進行水平的內插
將dt2jul、jul2dt寫成副程式dtconvertor，簡化程式版面
對照表也寫成json檔案

#argumens numNest gribFname
import netCDF4
import numpy as np
import datetime
from scipy.interpolate import griddata
import json

from pyproj import Proj
import sys,os,subprocess
from dtconvertor import dt2jul, jul2dt



分別讀取分子量(mws)、EAC4物質代碼(dic)、與對照到CMAQ氣狀物(nms_gas)及粒狀物(nms_part)之名稱
- 層數的對照表(deprecated，改在wsite進行對照)

for v in ['mws','dic','nms_gas','nms_part']:
 with open(v+'.json', 'r') as jsonfile:
   exec(v+'=json.load(jsonfile)')
uts=['PPM',"ug m-3          "]
l34=['21', '22', '23', '24', '25', '26', '27', '28', '29', '30', '31', '32', '33', '34', '35', '36', '37', 
    '38', '39', '40', '42', '43', '44', '46', '47', '48', '49', '50', '51', '53', '54', '56', '57', '59']
l40=['21', '21', '22', '22', '23', '24', '24', '25', '25', '26', '27', '28', '28', '29', '30', '31', '32', '32', '33', '34', 
    '35', '36', '37', '38', '39', '40', '42', '43', '44', '46', '47', '48', '49', '50', '51', '53', '54', '56', '57', '59']
d40_23={39-k:l34.index(l40[k]) for k in range(40)}

讀取一個BCON的樣版。從中取得粒狀物之間的比例關係(rate)
- 因為EAC4的粒狀物沒有粒徑I、J之分，需要將其切分。

byr=subprocess.check_output('pwd',shell=True).decode('utf8').strip('\n')[-2:]
#read a BC file as rate base
fname='/nas1/cmaqruns/2019base/data/bcon/BCON_v53_1912_run5_regrid_20191201_TWN_3X3'
nc = netCDF4.Dataset(fname,'r')
Vb=[list(filter(lambda x:nc.variables[x].ndim==j, [i for i in nc.variables])) for j in [1,2,3,4]]
rate={}
for v in nms_part:
 nms=nms_part[v]
 for nm in nms:
   if nm not in Vb[2]:sys.exit(v+' not in BCON file')
 avg=[np.mean(nc.variables[nm][:]) for nm in nms]
 sum_avg=sum(avg)
 if sum_avg==0:sys.exit('sum_avg==0')
 ratev=[avg[i]/sum_avg for i in range(len(avg))]
 rate.update({v:ratev})
for v in nms_gas:
 rate.update({v:[1.]})


讀取合併後之EAC4檔案

#read the merged grib files (ncl_convert2nc)
#lay,row in reversed directions
fname=sys.argv[1]
nc = netCDF4.Dataset(fname,'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,nlay,nrow,ncol=(nc.variables[V[3][0]].shape[i] for i in range(4))
#read the timestamp in nc and store at /expand the nc1
SDATE=[datetime.datetime.strptime(''.join([str(i, encoding='utf-8') for i in list(nc.variables[V[1][0]][t, :])]),\
'%m/%d/%Y (%H:%M)') for t in range(nt)]
bdate,edate=SDATE[0],SDATE[-1]
delt=edate-bdate
ntA=int(delt.total_seconds()/3600.)
JuliHr=[int((bdate+datetime.timedelta(hours=t)).strftime("%Y%j%H")) for t in range(ntA)]


讀取m3.nc模版
- 污染項目需事先按照nms_gas和nms_part的對照表挑選污染項目。
- 先對時間進行展開。由於模版污染項目總數即為nms_gas和nms_part值的總數，因此不會有未給定而被遮蔽的情形，不必花時間將矩陣全部清空。

N='D1'
tmps={'D1':'templateD1.ncV49K34','D2':'templateD2.ncV49K34'}
path='./'
fnameO=fname.replace('.nc',N+'.m3.nc')
if not os.path.exists(fnameO):os.system('cp '+path+tmps[N]+' '+fnameO)
nc1= netCDF4.Dataset(fnameO,'r+')
V1=[list(filter(lambda x:nc1.variables[x].ndim==j, [i for i in nc1.variables])) for j in [1,2,3,4]]
nv1=len(V1[3])
nt1,nlay1,nrow1,ncol1=nc1.variables[V1[3][0]].shape
if nt1<ntA:
 print('expand the matrix')
 for t in range(ntA):
   nc1.variables['TFLAG'][t,0,0]=0
if nt1>ntA:
 nc1.close()
 ncks=subprocess.check_output('which ncks',shell=True).decode('utf8').strip('\n')
 os.system(ncks+' -O -d TSTEP,0,'+str(ntA-1)+' '+fnameO+' '+fnameO+'_tmp;mv '+fnameO+'_tmp '+fnameO) 
 nc1= netCDF4.Dataset(fnameO,'r+')
nc1.SDATE=JuliHr[0]//100
nc1.STIME=JuliHr[0]%100*10000
var=np.zeros(shape=(ntA,nc1.NVARS,2))
var[:,:,0]=np.array([i//100 for i in JuliHr])[:,None]
var[:,:,1]=np.array([i%100  for i in JuliHr])[:,None]*10000
nc1.variables['TFLAG'][:,:,:]=var[:,:,:]

空間內插之準備
- 因EAC4檔案為等經緯度、非直角座標系統，因此選擇以scipy.griddata方式內插

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)

xlon, xlat = nc.variables['lon_0'][:].flatten(), np.flip(nc.variables['lat_0'][:].flatten())
lonm, latm = np.meshgrid(xlon, xlat)
x,y=pnyc(lonm,latm, inverse=False)

#interpolation indexing
x1d=[nc1.XORIG+nc1.XCELL*i for i in range(ncol1)]
y1d=[nc1.YORIG+nc1.YCELL*i for i in range(nrow1)]
x1,y1=np.meshgrid(x1d, y1d)
maxx,maxy=x1[-1,-1],y1[-1,-1]
minx,miny=x1[0,0],y1[0,0]
boo=(abs(x) <= (maxx - minx) /2+nc1.XCELL*10) & (abs(y) <= (maxy - miny) /2+nc1.YCELL*10)
idx = np.where(boo)
mp=len(idx[0])
xyc= [(x[idx[0][i],idx[1][i]],y[idx[0][i],idx[1][i]]) for i in range(mp)]


讀取[空氣密度]數據，以備粒狀物濃度的計算

print('read the density of air')
dlay=np.array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16,
      17, 18, 19, 20, 21, 23, 24, 25, 26, 28, 29, 30, 32, 34, 35, 37, 39])
dens=np.zeros(shape=(ntA,40, nrow1, ncol1))
caldat=list(set([int((bdate+datetime.timedelta(hours=t)).strftime("%Y%m%d")) for t in range(ntA)]))
caldat.sort()
for c in caldat:
 iday=caldat.index(c)
 fname='/nas1/cmaqruns/20'+byr+'base/data/mcip/RHO/RHO.'+str(c)+'.nc'
 ncr = netCDF4.Dataset(fname,'r')
 ntr=min(24,ncr.dimensions['TSTEP'].size)
 t1=iday*24
 t2=min(ntA,t1+min(24,ntr))
 hrs=t2-t1
 dens[t1:t2,:,:,:]=ncr.variables['DENS'][:hrs,:,:,:] *1E9 #(kg to microgram)
dens2=np.zeros(shape=(ntA,nlay1, nrow1, ncol1))
for k in range(nlay1):
 dens2[:,k,:,:]=dens[:,dlay[k],:,:]

讀取EAC4濃度內容，進行空間內插
- EAC4的高度方向自頂到底、緯度自北向南，與m3.nc定義相反，需要進行翻轉(np.flip)，4個軸中高度、緯度分別是第1、2軸。
- 翻轉後矩陣、併同前述d01網格座標系統，一起輸入griddata模組進行(線性)內插。

var=np.zeros(shape=(nt, nlay, nrow, ncol))
zz=np.zeros(shape=(nt, nlay1, nrow1, ncol1))
var2=np.zeros(shape=(ntA,nlay1, nrow1, ncol1))
print('fill the matrix')
for v in list(nms_gas)+list(nms_part) :
 var[:,:,:,:]=np.flip(nc.variables[v][:,:,:,:], [1,2])
 for t in range(nt):
   c = np.array([var[t,:,idx[0][i], idx[1][i]] for i in range(mp)])
   for k in range(nlay1):
     zz[t,k,:,: ] = griddata(xyc, c[:,k], (x1, y1), method='linear')

時間(線性)內插

 for t in range(0,ntA,3):
   t3=int(t/3)
   var2[t+0,:,:,:]=zz[t3,:,:,:]
   var2[t+1,:,:,:]=zz[t3,:,:,:]*2/3+zz[t3+1,:,:,:]*1/3
   var2[t+2,:,:,:]=zz[t3,:,:,:]*1/3+zz[t3+1,:,:,:]*2/3

氣狀物單位轉換(重量混合比轉PPM)

   130	  if v in nms_gas:
   131	    nc1.variables[nms_gas[v]][:]=var2[:]*rate[v][0] * 28.E6/mws[dic[v]] #mixing ratio to ppm

粒狀物單位轉換(重量混合比轉 μ g/M³)

 else:
#    unit=1E9*dens[:] #28.E6/mvol #mixing ratio(kg/kg) to microgram/M3
   nms=nms_part[v]
   for nm in nms:
     nc1.variables[nm][:]+=var2[:] * rate[v][nms.index(nm)] * dens2[:]
 print(v)


給定全域屬性，存檔，離開。

nc1.SDATE,nc1.STIME=nc1.variables['TFLAG'][0,0,:]
nc1.NLAYS=nlay1
nc1.TSTEP=10000
nc1.close()


m3.nc檔案之後處理(combine)

複製一個完整批次的CCTM_ACONC檔案做為容器
擷取全月結果中該批次的日期部分，將濃度值倒入容器，即可進行CMAQ的後處理combine.exe，整併VOCs、PM₁₀、PM_2.5等項目。

結果檢視

2018/4/5～6 大陸沙塵暴之EAC4濃度水平分布mov
2018/4/5～6 大陸沙塵暴之EAC4濃度垂直分布（臺灣為中心）mov

程式下載

Download: grb2D1m3.py

Reference

ECMWF, EAC4 (ECMWF Atmospheric Composition Reanalysis 4), copernicus,record updated 2021-12-07 16:10:05 UTC
純淨天空, python numpy flip用法及代碼示例, vimsky -Python学习园, Scipy Tutorial-多维插值griddata, cpython