Dr. Kuang's Utilities

hsinda3G.py程式說明

 
  • Jul 15, 2022
  • views
>

背景

  • 這支程式讀取傳統的CAMx點源point_source格式檔案，從中抽取開發計畫的排放數據，填入CAMx6的nc檔案中。
  • CAMx6點源nc檔案的內容、以及與CAMx7的差異詳見比較表CAMx67。

程式之執行

檔案準備

CAMx6點源檔案模版

  • 檔案位置：/nas1/camxruns/2016_v7/ptse/XindaG3/template.nc
  • 檔案維度內容
kuang@master /nas1/camxruns/2016_v7/ptse/XindaG3
$ ncdump -h $nc|H
netcdf template {
dimensions:
        TSTEP = UNLIMITED ; // (1 currently)
        NSTK = 1 ;
        VAR = 41 ;
        DATE-TIME = 2 ;
variables:
        float ACET(TSTEP, NSTK) ;
                ACET:units = "mole/hr         " ;
                ACET:long_name = "ACET            " ;
...
  • 只有一支煙囪(開發計畫之集束煙囪)
  • 運用ncqdp來使時間成為可增加(unlimited)的維度

CAMx點源檔案

  • 檔名：fortBE.14.hsinda.1.h80.n5.09Mp
  • point_source格式，由fortran程式產生。為環評使用之檔案。

引數

  • 無
  • 程式內執行月份之迴圈。(不需同步多工進行)

執行程式

  • python hsinda3G.py

程式說明

IO模組

PseudoNetCDF模組

  • 使用PseudoNetCDF的pncopen指定格式’ppoint_source’來開啟二進位檔案
  • 並將煙囪參數陣列全部讀取出來
,4):
  exec('v'+str(j)+'=list(filter(lambda x:pt.variables[x].ndim=='+str(j)+', [i for i in pt.variables]))')
nhr,nvar,dt=pt.variables[v3[0]].shape
nt,nopts=pt.variables[v2[0]].shape
d={}
for v in 'XYHDTV':
  var=v+'STK'
  d.update({v:np.array(list(pt.variables[var][:]))})
d.update({'I':np.array([i for i in range(nopts)])})

讀取開發計畫之數據

  • 由於point_source([[2022-07-18-PTSRCE_format]]1)二進位檔案沒有煙囪編號，需藉由煙囪參數的特徵來辨識。
    • 位置及高度
idx=np.where(abs(d['X']-v1_xinda[0])<3000)
idy=np.where(abs(d['Y'][idx[0]]-v1_xinda[1])<3000)
idh=np.where(abs(d['H'][idx[0]][idy[0]]-80)<5)
I=d['I'][idx[0]][idy[0]][idh[0]]
parms={v:d[v][I] for v in 'XYHDTV'}
emiss={v:pt.variables[v][0,I] for v in v2 if pt.variables[v][0,I]>0}

對時間維度逐步開展

  • 時間展開必須在填入數據之前，以避免矩陣被遮罩(mask)
  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)]

將數據存到新檔案

  • 三集束
  d['D'][I]=d['D'][I]*np.sqrt(3.)
  for v in 'XYHDTV':
    var=v+'STK'
    nc.variables[var][0]=d[v][I]
  for v in emiss:
    nc.variables[v][:,0]=emiss[v][0]*3

程式下載

結果確認

  • 可以使用pt2em_d01.py以地面排放檔將點源nc檔網格化，再以VERDI檢視是否在開發計畫位置確有增加排放量。

  1. https://sinotec2.github.io/FAQ/2022/07/18/PTSRCE_format.html “ CAMx點源格式說明” ↩

Search