Link Search Menu Expand Document
Focus on Air Quality

pt_const程式說明

Table of contents
  • 目的：產生CMAQ所需的點源排放量中常數、不隨時間變化的部分

    執行

  • 引數：CAMx分月高空點源，MM=01~12
for MM in {01..12};
do python pt_constLL.py fortBE.413_teds11.ptse$MM.nc
done
  • 輸入檔(ptse檔案)：ptsEnHRBE.py結果
  • 結果：teds11.YYMM.const.nc

版本

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

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

項目EAsia_81KTWN_3X3說明
資料庫REAS+TEDSTEDS前者範圍較大需要REAS資料庫
引數及主要檔逐月CAMx高空點源nc檔案逐月CAMx高空點源nc檔案(same)
其他輸入檔point_reas16’+mm+’.csv’, stack_groups_ptnonipm_12US1_2016ff_16j.ncstack_groups_ptnonipm_12US1_2016ff_16j.ncnc檔模版來自於USEPA之範例檔。csv詳見https://sinotec2.github.io/FAQ/2022/07/06/REASPtRd.html
結果檔名teds11.YYMM.const.ncteds11.YYMM.const.nc(same)

分段說明

  • 模組及時間標籤轉換副程式(dt2jul及jul2dt)
     1  import numpy as np
     2  import netCDF4
     3  import PseudoNetCDF as pnc
     4  import os,sys,twd97
     5  from pyproj import Proj
     6  import subprocess
     7  from pandas import *
     8  import datetime
     9
    10  def dt2jul(dt):
    11    yr=dt.year
    12    deltaT=dt-datetime.datetime(yr,1,1)
    13    deltaH=int((deltaT.total_seconds()-deltaT.days*24*3600)/3600.)
    14    return (yr*1000+deltaT.days+1,deltaH*10000)
    15
    16  def jul2dt(jultm):
    17    jul,tm=jultm[:]
    18    yr=int(jul/1000)
    19    ih=int(tm/10000.)
    20    return datetime.datetime(yr,1,1)+datetime.timedelta(days=int(jul-yr*1000-1))+datetime.timedelta(hours=ih)
    21
    22
  • 讀入檔案名稱的月份
  • 準備座標轉換之套件
    23  mon=int(sys.argv[1][-5:-3])
    24  #join the pollutants for this month
    25  pth='/nas1/cmaqruns/2016base/data/ptse/'
    26  #prepare interceptions
    27  Latitude_Pole, Longitude_Pole = 23.61000, 120.9900
    28  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點源檔案

  • 讀取各維度變數名稱
  • 讀取全域變數
    30  fname1=sys.argv[1]
    31  pt=netCDF4.Dataset(fname1, 'r')
    32  v3=list(filter(lambda x:pt.variables[x].ndim==3, [i for i in pt.variables]))
    33  v2=list(filter(lambda x:pt.variables[x].ndim==2, [i for i in pt.variables]))
    34  v1=list(filter(lambda x:pt.variables[x].ndim==1, [i for i in pt.variables]))
    35  nhr,nvar,dt=pt.variables[v3[0]].shape
    36  nt,nopts=pt.variables[v2[0]].shape
    37  XCELL=pt.XCELL
    38  XORIG=pt.XORIG
    39  YORIG=pt.XORIG
    40
  • 讀取起始時間(0800 LST)
  • 準備ncks程式
    41  tb=pt.STIME-80000 #UTC
    42  bdate=jul2dt([pt.SDATE,pt.STIME])
    43  fname1=fname1.replace('fortBE.413_','').replace('ptse','19').replace('.nc','')
    44  fname=fname1+'.const.nc'
    45  fname0=pth+'stack_groups_ptnonipm_12US1_2016ff_16j.nc' #as template
    46
    47  #ncks path
    48  path={'114-32-164-198.HINET-IP.hinet.net':'/opt/anaconda3/bin/', 'node03':'/usr/bin/', \
    49          'master':'/cluster/netcdf/bin/','DEVP':'/usr/bin/'}
    50  hname=subprocess.check_output('echo $HOSTNAME',shell=True).decode('utf8').strip('\n')
    51  if hname not in path:
    52    sys.exit('wrong HOSTNAME')
    53  ncks=path[hname]+'ncks'

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

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

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

    59  nc = netCDF4.Dataset(fname,'r+')
    60  nc.NROWS=nopts
    61  nc.GDNAM='TWN_3X3'
    62  nc.P_ALP = np.array(10.)
    63  nc.P_BET = np.array(40.)
    64  nc.P_GAM = np.array(120.98999786377)
    65  nc.XCENT = np.array(120.98999786377)
    66  nc.YCENT = np.array(23.6100196838379)
    67  nc.XCELL=XCELL #27000.000
    68  nc.YCELL=XCELL #27000.000
    69  nc.XORIG=XORIG #-877500.0
    70  nc.YORIG=XORIG #-877500.0
    71  nc.SDATE,nc.STIME=dt2jul(bdate+datetime.timedelta(hours=-8))
    72
  • 點源變數名稱：CMAQ對照到CAMx
  • 單位轉換：每小時->每秒
  • 給定經緯度、網格位置、時間標籤等
    73  mp={'STKDM':'stkdiam','STKHT':'stkheight','STKTK':'stktemp','STKVE':'stkspeed','XLOCA':'xcoord', 'YLOCA':'ycoord',}
    74
    75  for v in mp:
    76    nc.variables[v][0,0,:,0]=pt.variables[mp[v]][:]
    77  #velocity is m/hr in CAMx pt, but in m/s in CMAQ_pt
    78  nc.variables['STKVE'][0,0,:,0]=pt.variables[mp[v]][:]/3600.
    79  nc.variables['IFIP'][0,0,:,0]=[1000+i for i in range(nopts)]
    80  nc.variables['ISTACK'][0,0,:,0]=[1+i for i in range(nopts)]
    81
    82  nc.variables['COL'][0,0,:,0]=[int((i-nc.XORIG)/nc.XCELL) for i in x]
    83  nc.variables['ROW'][0,0,:,0]=[int((i-nc.YORIG)/nc.YCELL) for i in y]
    84  nc.variables['TFLAG'][0,:,0]=[nc.SDATE for i in range(nc.NVARS)]
    85  nc.variables['TFLAG'][0,:,1]=[nc.STIME for i in range(nc.NVARS)]
    86  nc.variables['LATITUDE'][0,0,:,0]=lat
    87  nc.variables['LONGITUDE'][0,0,:,0]=lon
    88  nc.close()

程式下載

Reference