Link Search Menu Expand Document
Focus on Air Quality

高空點源之時變係數

Table of contents

背景

  • 高空點源的時變係數骨幹是CEMS數據，然而同一工廠無數據、鄰近工業區其他廠無數據者，亦會參考CEMS設定其時變係數。
  • 排放量整體處理原則參見處理程序總綱、針對點源之處理及龐大.dbf檔案之讀取，為此處之前處理。程式也會呼叫到ptse_sub中的副程式

程式說明

程式執行

因排放物質類別與污染源製造程序的特徵有關，必須分開個別處理，此處則以個別污染項目執行ptseE_ONS.py，執行方式如下：

for spe in CO NMHC NOX PM SOX;do python ptseE_ONS.py $spe;done
  • 由於程式消耗記憶體非常大量，如要同時進行，需注意記憶體的使用情形。
  • 污染源個數與排放高度限值的設定、地面PM排放條件之給定、以及數據年代等等都有關係，需配套紀錄。

排放與CEMS資料檔之讀取及準備

  • 引用模組
    • 程式用到ptse_sub中的副程式CORRECT, add_PMS, check_nan, check_landsea, FillNan, WGS_TWD, Elev_YPM
kuang@node03 /nas1/TEDS/teds11/ptse
$ cat -n ptseE_ONS.py
     1
     2  #! crding = utf8
     3  from pandas import *
     4  import numpy as np
     5  import os, sys, subprocess
     6  import netCDF4
     7  import twd97
     8  import datetime
     9  from calendar import monthrange
    10  from scipy.io import FortranFile
    11
    12  from ptse_sub import CORRECT, add_PMS, check_nan, check_landsea, FillNan, WGS_TWD, Elev_YPM
    13
  • 從工作目錄讀取teds版本與年代
      14  #Main
  15  P=subprocess.check_output('pwd',shell=True).decode('utf8').strip('\n')+'/'
  16  teds=int(P.split('/')[3][-2:])
  17  yr=2016+(teds-10)*3
  18  ndays=365
  19  if yr%4==0:ndays=366
  20  s365=set([i*24 for i in range(ndays)])
  21  nhrs=ndays*24
  22
  • 檔案讀取與品質確認
    • 此個案將所有點源資料庫的數據都以「高空」方式處理(cutting height of stacksHs=10)
    23  Hs=10 #cutting height of stacks
    24  #Input the TEDS csv file
    25  try:
    26    df = read_csv('point.csv', encoding='big5')
    27  except:
    28    df = read_csv('point.csv')
    29  # check_NOPandSCC(0)
    30  df = check_nan(df)
    31  # check and correct the X coordinates for isolated islands
    32  df = check_landsea(df)
    33  df = WGS_TWD(df)
    34  df = Elev_YPM(df)
  • 使用Hs進行篩選「高空」點源
    35  df=df.loc[(df.HEI>=Hs) & (df.NO_S.map(lambda x:x[0]=='P'))].reset_index(drop=True)
  • 排除沒有SNCPV排放者
    36  df['SUM']=[i+j+k+l+m for i,j,k,l,m in zip(df.SOX_EMI,df.NOX_EMI,df.CO_EMI,df.PM_EMI,df.NMHC_EMI)]
    37  df=df.loc[df.SUM>0].reset_index(drop=True)
    38  df['CP_NO'] = [i + j for i, j in zip(list(df['C_NO']), list(df['NO_S']))]
    39  df['DY1']=[i*j for i,j in zip(df.DW1,df.WY1)]
    40  df['HY1']=[i*j for i,j in zip(df.HD1,df.DY1)]
    41

讀取並填滿CEMS資料檔

  • 填滿資料表
    • 程式運作需要每筆管煙(管編+煙編)、每個小時都要有數值。將DataFrame轉成矩陣，再轉回DataFrame即可。
    42  #71 factories with CEMS will emit (at ground) when stacks are operating
    43  fname=P+'point_cems.csv'
    44  cems=read_csv(fname)
    45  val='SOX PM NOX FLOW X_BLANK1 X_BLANK2'.split()
    46  nval=len(val)
  • 管煙欄位為處理過後檔案的特徵
    47  if 'CP_NO' not in cems.columns: #pre-process
    48    cems=cems.drop(cems.loc[cems.C_NO=='C_NO'].index).reset_index(drop=True)
  • 新增管編+煙編(管煙)之新標籤、新增PM(設為SN之平均值)
    49    cems['CP_NO'] = [i + j for i, j in zip(list(cems['C_NO']), list(cems['NO_S']))]
    50    cems['PM']=[(i+j)/2 for i,j in zip(cems.SOX,cems.NOX)]
  • 新增時間標籤MDH。有的年度提供的CEMS檔案小時標記為0000~2300，有的是0~23，因此需要判斷一下
    51    if max(cems.HOUR)>100:
    52      cems['MDH']=[int(i*10000+j*100+k/100) for i,j,k in zip(cems.MONTH,cems.DATE,cems.HOUR)]
    53    else:
    54      cems['MDH']=[int(i*10000+j*100+k) for i,j,k in zip(cems.MONTH,cems.DATE,cems.HOUR)]
  • 一個(管煙,時籤)組合只應對應一組CEMS數據，以pivot_table sum進行整併
    55    cems=pivot_table(cems,index=['CP_NO','MDH'],values=val,aggfunc=sum).reset_index()
  • 維度標籤之計算。如果使用標籤=序列.index(值)指令將會非常耗時，直接使用dict{值:標籤}，會快很多。
    56    #cems(df) convert to cemsM(matrix)
    57    for MC in ['CP_NO','MDH']:
    58      mc=MC.lower()
    59      exec(mc+'=list(set(cems.'+MC+'))');exec(mc+'.sort()')
    60      exec('n'+MC+'=len('+mc+')')
    61      exec('d'+MC+'={'+mc+'[i]:i for i in range(n'+MC+')}')
    62      exec('cems["i'+MC+'"]=[d'+MC+'[i] for i in cems.'+MC+']')
    63    if len(mdh)!=ndays*24:sys.exit('mdh coverage not enough!')
  • DataFrame轉成Array，因所有值預設為0，如此就補滿空缺值了。有數據的部分再填入Array。
    64    cemsM=np.zeros(shape=(nMDH,nCP_NO,nval))
    65    for i in range(nval):
    66      cemsM[cems.iMDH[:],cems.iCP_NO[:],i]=cems[val[i]]
  • Array再轉回DataFrame。因後面CEMS數據使用邏輯太多樣了，Array形式不敷應用。
    67    DD={}
    68    for i in range(nval):
    69      DD[val[i]]=cemsM[:,:,i].flatten()
    70    DD['MDH']  =[i for i in mdh for j in cp_no]
    71    DD['CP_NO']=[j for i in mdh for j in cp_no]
    72    cems=DataFrame(DD)
    73    cems['C_NO']=[i[:8] for i in cems.CP_NO]
    74    cems['MD']=[i//100 for i in cems.MDH]
    75    cems.set_index('CP_NO').to_csv(fname)

整理CEMS各廠運作時間模式

  • 讀取資料庫欄位(管編C_NO, 管煙CP_NO, 月日時MDH, 月日MD,)
    76  for MC in ['CP_NO','MDH','MD','C_NO']:
    77    mc=MC.lower()
    78    exec(mc+'=list(set(cems.'+MC+'))');exec(mc+'.sort()')
    79    exec('n'+MC+'=len('+mc+')')
    80
  • 沒有CEMS之同家工廠其它煙道、或鄰近工廠，可能具有上下游連動關係，因此適用該CEMS的操作特性，此處先行備妥。
    • 個別廠的日操作特性，分析全年SOX排放量之小時變化，按排放量排序，代表最有可能運作之小時。
    • 先執行pivot_table加總
    81  #Hour of Day pattern
    82  cems['HR']=[i%100 for i in cems.MDH]
    83  pv_cems1=pivot_table(cems,index=['C_NO','HR'],values='SOX',aggfunc=sum).reset_index()
    84
  • 以管編為索引之新的DataFrame，其內容為小時之序位
    85  cems_HROD=DataFrame({'C_NO':c_no})
    86  cems_HROD['SOX_HR_ODER']=0
    87  for ic in cems_HROD.index:
    88    pv1=pv_cems1.loc[pv_cems1.C_NO==c_no[ic]]
    89    pv3=pv1.sort_values('SOX',ascending=False).reset_index(drop=True)
    90    cems_HROD.loc[ic,'SOX_HR_ODER']=''.join(['{:d} '.format(i) for i in pv3.HR])
  • 同樣原理應用在全年的逐日變化，內容為序列mdh的標籤。排序依據全天的流量總計，值越高者表當天最可能運作。
    91  #orders for DY1
    92  pv_cems2=pivot_table(cems,index=['C_NO','MD'],values='FLOW',aggfunc=sum).reset_index()
  • MD欄位由月日轉變成mdh的標籤，轉變需要2步驟，否則會出錯
    93  #Indexing is an exhaustive process.
    94  iMD=[mdh.index(i*100) for i in pv_cems2.MD] #change the MMDD into index sequence among MMDD00's
    95  pv_cems2.MD=iMD
  • 以管編為索引之新的DataFrame，其內容為日期標籤之序位
    96  cems_DAOD=DataFrame({'C_NO':c_no})
    97  cems_DAOD['FLOW_DA_ODER']=0
    98  for ic in cems_DAOD.index:
    99    pv1=pv_cems2.loc[pv_cems2.C_NO==c_no[ic]]
   100    pv3=pv1.sort_values('FLOW',ascending=False).reset_index(drop=True)
   101    cems_DAOD.loc[ic,'FLOW_DA_ODER']=''.join(['{:d} '.format(i) for i in pv3.MD])
   102

套用CEMS或運作時間模式

  • 建立各管編座標值之資料表
   103  dfxy=pivot_table(df,index='C_NO',values=['UTM_E','UTM_N'],aggfunc=np.mean).reset_index()
   104
  • BLS為污染項目與布林值的對照表。以直接提取資料庫中符合布林值、要處理的筆數。
    • {'NMHC':'PM', 'NMHC':'PM'}因為NMHC及CO沒有CEMS數據，假設與PM、NOX一樣。
   105  #booleans for pollutant selection
   106  c2v={'NMHC':'PM','SOX':'SOX','NOX':'NOX','PM':'PM','NMHC':'PM'} #point.csv vs cems.csv
   107  BLS={c:df[c+'_EMI']>0 for c in c2v}
   108  colT=['HD1','DY1','HY1']
   109  col=['C_NO','CP_NO','HD1','DY1','HY1']+[i for i in df.columns if 'EMI' in i]
  • s為物質種類，須由引數讀取，且限定在BLS的索引範圍(c2v)
   110  for spe in [s for s in [sys.argv[1]] if s in BLS]:
   111    dfV=df[col].loc[BLS[spe]].reset_index(drop=True)
  • 對全廠加總並形成新的資料庫dfV
   112    dfV1=pivot_table(dfV,index='CP_NO',values=spe+'_EMI',aggfunc=sum).reset_index()
   113    dfV2=pivot_table(dfV,index='CP_NO',values=colT,aggfunc=np.mean).reset_index()
   114    dfV=merge(dfV1,dfV2,on='CP_NO')
   115    dfV['C_NO']=[i[:8] for i in dfV.CP_NO]
   116    for c in colT:
   117      dfV[c]=np.array(dfV[c],dtype=int)
  • 對照資料庫(a)與cems(b)的管編cp
    • ab為二者交集
    • b1為b-a，有CEMS數據卻不在資料庫中、沒有座標
    • c1為ab-b1(就等於ab)
    • 列出所有c1的座標序列
   118    a,b=list(set(dfV.C_NO)),list(set(cems.C_NO));a.sort();b.sort()
   119    ab=[i for i in a if i in b]
   120    cp=list(set(dfV.CP_NO))
   121    cp.sort()
   122    ons=np.zeros(shape=(len(cp),nMDH))#,dtype=int)
   123    #other fatories without CEMS, take the nearest one
   124    b1=set(b)-set(dfxy.C_NO) #cems factory but without UTM location
   125    c1=[c for c in b if c not in b1 and c in a] #cems plant with X,Y
   126    cemsX=np.array([list(dfxy.loc[dfxy.C_NO==c,'UTM_E'])[0] for c in c1])
   127    cemsY=np.array([list(dfxy.loc[dfxy.C_NO==c,'UTM_N'])[0] for c in c1])
  • 逐一管編進行迴圈
    • 如果工廠沒有CEMS(c not in ab)，則以最近的一家廠的日變化特徵為其時變係數，令c_cems為該廠管編
    • 如果該廠日期標籤之序位不足全年日數(ndays)，也將補足缺漏日之標籤list(s365-set(pv2MD))
   128    #loop for every factories
   129    for c in [i for i in a if i not in b1]:
   130      c_cems=c
   131      if c not in ab:
   132        x0,y0=list(dfxy.loc[dfxy.C_NO==c,'UTM_E'])[0],list(dfxy.loc[dfxy.C_NO==c,'UTM_N'])[0]
   133        dist=(cemsX-x0)**2+(cemsY-y0)**2
   134        idx=list(dist).index(min(dist))
   135        c_cems=c1[idx]
   136      pv2MD=np.array(list(cems_DAOD.loc[cems_DAOD.C_NO==c_cems,'FLOW_DA_ODER'])[0].split(),dtype=int)
   137      if len(pv2MD)<ndays: pv2MD=np.array(list(pv2MD)+list(s365-set(pv2MD)))
   138      df_cp=dfV.loc[dfV.C_NO==c].reset_index(drop=True)
  • 該廠所有煙編逐一進行
   139      #loop for every NO_S in this factory
   140      for p in set(df_cp.CP_NO):
   141        ip=cp.index(p)
  • 如果是CEMS煙道，其全年變化即為CEMS資料表中數值
  • point_cems.csv檔案內容須事先確認其單位，此處要求其全年總合為1.0
   142        if p in set(cems.CP_NO):
   143          ons[ip,:]=cems.loc[cems.CP_NO==p,c2v[spe]]*nhrs
  • 將工廠操作矩陣(維度為(日數、時數))先設定為0
  • 如果不是，則由資料庫中讀取工作日數及小時數
    • 如果是全天連續操作，則逐時都需標籤
       144        else:
 145          dy1=dfV.DY1[ip]
 146          hd1=dfV.HD1[ip]
 147          md3=pv2MD[:dy1]
 148          days=np.zeros(shape=(dy1,hd1),dtype=int)
 149          if hd1==24:
 150            hrs=np.array([i for i in range(24)],dtype=int)
  • 不是24小時，最有可能的小時，則由前面準備好的cems_HROD.SOX_HR_ODER第1小時開始，依序填入工作時數。
    • 工作日則由為前述最可能工作的日期標籤之序位依序填入，形成days(月、日、時標籤)
    • 將days壓平成為一維矩陣
    • 這些有運作的時間，其ons值為1(其餘內設為0)
   151          else:
   152            first=np.array(list(cems_HROD.loc[cems_HROD.C_NO==c_cems,'SOX_HR_ODER'])[0].split(),dtype=int)[0]
   153            hrs=np.array([(first+ih)%24 for ih in range(hd1)])
   154          for id in range(dy1):
   155            days[id,:]=md3[id]+hrs[:]
   156          idx=days.flatten()
   157          ons[ip,idx]=1.

輸出結果

   158  #other sources
   159    fnameO=spe+'_ECP'+str(len(cp))+'_MDH'+str(len(mdh))+'_ONS.bin'
   160    with FortranFile(fnameO, 'w') as f:
   161      f.write_record(cp)
   162      f.write_record(mdh)
   163      f.write_record(ons)

檔案下載

Reference