計算時間變異係數對照表

Table of contents

背景
時間變異係數檔案前處理(prep_dfAdmw.py)分段說明
副程式(prep_df.py)分段說明
nc_fac.json檔案產生程式(prep_json.py)分段說明
檔案下載
Reference

背景

為消除資料庫在時間與空間的相依性，需要消除資料庫類別-縣市之維度，此處即針對面源的時間變異係數與類別-縣市之對照表(nc_fac.json)進行計算，以備未來套用。主要問題出現在環保署提供的時變係數檔，空間資料為中文名稱，與主資料庫是鄉鎮區代碼2碼編號不能對應，分3個程式說明：
- 時間變異係數(csv)檔案前處理
- csv檔案之產生
- 將csv檔案應用到面源排放量資料庫，並展開至全年逐時之序列，存成nc_fac.json。
排放量整體處理原則參見處理程序總綱、針對面源之處理及龐大.dbf檔案之讀取與重新計算網格座標，為此處之前處理。

時間變異係數檔案前處理(prep_dfAdmw.py)分段說明

此程式為可執行的主程式，主要處理縣市名稱與編號的對照，中文改成英文。

引用模組

會用到中文轉英文的pypinyin模組，以使程式可以在py27/py37都可以使用

引用副程式prep_df，說明如後。

kuang@node03 /nas1/TEDS/teds11/area
$ cat -n prep_dfAdmw.py
 1  #!coding=utf8
 2  import numpy as np
 3  from pandas import *
 4  from pypinyin import pinyin, lazy_pinyin
 5  import subprocess
 6  import json
 7  from datetime import datetime, timedelta
 8  from prep_df import PrepDf
 9  import sys, os
10
11

從中文檔案cnty.csv開始處理

補充連江縣

#input and prepare the cnty names
df_cnty=read_csv('cnty.csv',encoding='big5')
for i in range(len(df_cnty)):
  cha=df_cnty.loc[i,'cnty']
  ll=lazy_pinyin(cha)
  s=''
  for l in ll:
    s=s+l
  df_cnty.loc[i,'cnty']=s
df_cnty['no']=[str(s) for s in df_cnty.no]
df_cnty['no']=['0'*(2-len(s))+s for s in df_cnty.no]
df_cnty.loc[len(df_cnty),'no']=51
df_cnty.loc[len(df_cnty)-1,'cnty']='lianjiangxian'

補充五都名稱、機場名稱

  25
#dictionary of region name to county number
d_cnty={x:str(int(y)) for x,y in zip(df_cnty['cnty'],df_cnty['no'])}
d_cnty.update({'xinbeishi':'31','taoyuanshi':'32',\
'hualianjichang':'45','taidongjichang':'46','lvdaojichang':'46','lanyujichang':'46','magongjichang':'44','jinmenjichang':'50', \
'taizhongjichang':'36','tainanjichang':'21','jiayijichang':'40','qimeijichang':'44','wanganjichang':'44','beiganjichang':'51',\
'nanganjichang':'51','songshanjichang':'1','xiaogangjichang':'2','taoyuanjichang':'32'})
  32

新的直轄市與舊編號的對照

  33  #combination of 5 Municipalities
  34  NewCities={'tainanshi':{'21','41'},'taizhongshi':{'17','36'},'gaoxiongshi':{'2','42'}}
  35

補充港口名稱與縣市編號之對照

#including the harbors
gangko='gaoxionggang jilonggang anpinggang shenaozhuanyongyougang taibeigang shalunwaihaixieyoufutong suaogang taizhonggang mailiaogang budaigang xingdadianchangxiemeimatou yonganyehuatianranqijieshouzhan magonggang hualiangang hepinggang jinmengang mazugang'.split()
gangkoC=['2','11','21','31','31','36','34','36','39','40','42','42','44','45','45','50','51']
d_gangko={i:j for i,j in zip(gangko,gangkoC)}
d_cnty.update(d_gangko)
  41

環保署提供的時變係數檔案亦有以空品區為單位之情形，需另建對照表

# air quality basins
d_kpq={'beibukongpinqu':['taibeishi','taoyuanshi','xinbeishi','jilongshi'], \
'gaopingkongpinqu':['gaoxiongshi','gaoxiongxian','pingdongxian'],\
'huadongkongpinqu':['hualianxian','taidongxian'], \
'yilankongpinqu':['yilanxian'],'yunjianankongpinqu':['yunlinxian','jiayixian','jiayishi','tainanshi'], \
'zhongbukongpinqu':['taizhongshi','taizhongxian','zhanghuaxian','nantouxian'], \
'zhumiaokongpinqu':['miaolixian','xinzhushi','xinzhuxian']}
d_kpq.update({'quanguo':list(df_cnty['cnty'])})
d_kpq.update({'lidao':['jinmenxian','lianjiangxian','penghuxian']})
  51

分別執行月變化、週變化與日變化檔案之整理與輸出

#execution of dataframe formings and savings
csvs={'m':'mon.csv','w':'week.csv','d':'day.csv'}
for t in 'mwd':
  df='df_A'+t
  try:
    exec(df+'=read_csv("'+df+'")')
  except:
    exec(df+'=PrepDf("'+csvs[t]+'")')
    exec(df+'.set_index("nsc2").to_csv("'+df+'")')
  61

由於分段執行，需將對照表寫出備用。如合併執行則無此需要。

# output the dictionary
for kc in ['cnty','kpq']:
  with open('d_'+kc+'.json', 'w', newline='') as jsonfile:
    exec('json.dump(d_'+kc+', jsonfile)')

副程式(prep_df.py)分段說明

主程式處理對照，副程式則處理時變係數檔內部的不一致性

引用模組

import numpy as np
from pandas import *
from pypinyin import pinyin, lazy_pinyi

先將中文改成英文，以便對照

#prepare the dataframe
def PrepDf(fname):
  path=''#'/home/sespub/teds10/08-時間分配權重/月週日時間權重/'
  df=read_csv(path+fname,encoding='big5')
  df.loc[df.NSC_SUB.map(lambda x: isna(x)),'NSC_SUB']='b'
  df['nsc2']=[str(x)+y for x,y in zip(df['NSC'],df['NSC_SUB'])]
  for i in range(len(df.REGION)):
    cha=df.loc[i,'REGION']
    if type(cha) in [int,float]:continue
    ll=lazy_pinyin(cha)
    if len(ll)==0:continue
    s=''
    for l in ll:
      s=s+l
    df.loc[i,'REGION']=s

時變係數檔中有以全形頓號「、」來表示縣市的聯集，需逐一拆解改成資料庫對應關係

  a=df.loc[df.DICT.map(lambda x:type(x)==str and '、' in x)]
  idx=a.index
  df=df.drop(idx).reset_index(drop=True)
  a=a.reset_index(drop=True)
  for i in range(len(a)):
    dct=a.loc[i,'DICT'].split('、')
    b=DataFrame({})
    for j in range(len(dct)):
      b=b.append(a.loc[i],ignore_index=True)
      b.loc[j,'DICT']=dct[j]
    df=df.append(b,ignore_index=True)

時變係數檔中有空白者，意即全國，需填以明確的文字

  idx=df.loc[df.REGION.map(lambda x:type(x)==float and np.isnan(x))].index
  df.loc[idx,'REGION']=['quanguo' for i in idx]
  for c in df.columns:
    if c in ['REGION','NSC','NSC_SUB','nsc2','DICT']:continue
    df[c]=[float(i) for i in list(df[c])]

時變係數檔中以51代表所有類別的船舶，與主資料檔的作法又有出入

#change 51b to 51A~51D, no need changing it manually
  snsc2=set(df.nsc2)
  if '51b' in snsc2:
    df51=df.loc[df.nsc2=='51b'].reset_index(drop=True)
    for s in 'ABCD':
      if '51A' not in snsc2:
        tmp=df51
        tmp.nsc2='51'+s
        df=df.append(tmp,ignore_index='True',sort=False)

確認是否漏了臺中機場的時間變化，如是，則以臺南機場變化形態代之。

# taizhongjichang is missing, fill it according tainanjichang
  if 'taizhongjichang' not in set(df.REGION):
    a=df.loc[df.REGION=='tainanjichang'].reset_index(drop=True)
    a.DICT=3605
    a.REGION='taizhongjichang'
    df=df.append(a,ignore_index=True)

刪去重覆記錄、重置順序、輸出DataFrame檔案

  df.drop_duplicates(inplace=True)
  df=df.reset_index(drop=True)
  return df

`nc_fac.json`檔案產生程式(prep_json.py)分段說明

引用模組

kuang@node03 /nas1/TEDS/teds11/area
$ cat -n prep_json.py
#!coding=utf8
import numpy as np
from pandas import *
import subprocess
import json
from datetime import datetime, timedelta
import sys, os
   8

讀入前面prep_dfAdmw.py的執行成果：df_Am、df_Aw、df_Ad三個檔案

# read the time variation factors
csvs={'m':'mon.csv','w':'week.csv','d':'day.csv'}
for t in 'mwd':
  df='df_A'+t
  try:
    exec(df+'=read_csv("'+df+'")')
  except:
    sys.exit('df_A? not found, please re-run prep_dfAdmw.py ')
  17

s_nsc2：所有面源分類項目形成的集合，有時間變化者

  18  #union of nsc2 in time variation files
  19  s_nsc2=set(df_Ad.nsc2) | set(df_Am.nsc2) | set(df_Aw.nsc2)
  20

開啟檔案、刪除重複記錄、產生nsc2、CNTY，此2項目為對照表的主要索引。

#open the TEDS area csv file
pwd=subprocess.check_output('pwd',shell=True).decode('utf8').strip('\n')
teds=pwd.split('/')[3][4:6]
fname='areagrid'+teds+'LL.csv'
df = read_csv(fname)
df.drop_duplicates(inplace=True)
df=df.reset_index(drop=True)
if 'nsc2' not in df.columns:
  df.loc[df['NSC_SUB'].map(lambda x: (type(x)==float and np.isnan(x)==True) or ( x==' ')),'NSC_SUB']='b'
  df['nsc2']=[str(x)+y for x,y in zip(df['NSC'],df['NSC_SUB'])]
df.loc[df['DICT'].map(lambda x:np.isnan(x)==True),'DICT']=5300.
if 'CNTY' not in df.columns:
  df['CNTY']=[str(int(s/100)) for s in list(df['DICT'])]
  34

定義資料庫的座標系統，去尾到1公里解析度

#definition the coordinates of database
minx,miny=min(df.UTME),min(df.UTMN)
df.UTME=round(df.UTME-minx,-3)
df.UTMN=round(df.UTMN-miny,-3)
df['YX']=np.array(df.UTMN+df.UTME/1000,dtype=int)
  40

同1網格、同一縣市別、同一類別，可能有好幾筆資料(不同鄉鎮區，由時變檔案角度無法區分)，排放量予以加總。縮減資料表。

#sum-up the grids which length maybe smaller than 1KM
cole=['EM_SOX','EM_NOX','EM_CO','EM_PM25','EM_PM','EM_NMHC']
coli=['CNTY', 'nsc2','YX']
df=pivot_table(df,index=coli,values=cole,aggfunc=np.sum).reset_index()

形成縣市代碼、面源類別、位置的序列、序列長度、標籤對照表

#nXXX:len of list XXX; dXXX:index dictionary of listXXX, iXXX:index of XXX
for c in ['CNTY','nsc2','YX']:
  exec(c+'=list(set(df.'+c+'))')
  exec(c+'.sort()')
  exec('n'+c+'=len('+c+')')
  exec('d'+c+'={'+c+'[i]:i for i in range(n'+c+')}')
  exec('df["i'+c+'"]=[d'+c+'[i] for i in df.'+c+']')
  52

時變係數檔之面源類別可能只有大類，而資料表中含細類，必須全都適用大類的係數

#same NSC in df_A? and df, but without NSC_SUB, all add to s_nsc2
nsc2b=set([i for i in s_nsc2 if i[-1]=='b'])
for ii in nsc2b-set(df.nsc2):
  i=int(ii[:-1])
  s_nsc2=s_nsc2|set(df.loc[df.nsc2.map(lambda x:x[:-1]==i),'nsc2'])

區域代碼組合(rgn)的個數n_rgn,、集合s_rgn、對照表d_rgn

#Tuple_of_Length ={0, 101, 10000, 10101, 131313, 171717, 202020, 212121, 250101}
n_rgn,s_rgn,d_rgn={},[],{}
for n in s_nsc2:
  tl_rgn=(len(df_Am[df_Am.nsc2==n])*100+len(df_Aw[df_Aw.nsc2==n]))*100+len(df_Ad[df_Ad.nsc2==n])
  s_rgn.append(tl_rgn)
  n_rgn[n]=tl_rgn
  d_rgn[tl_rgn]=n #recording last nsc2 for calling
s_rgn=set(s_rgn) #tuple of region numbers
  66

由teds版本計算年代(yr)

dts為全年逐時之時間標籤（datetime）

yr=2016+(int(teds)-10)*3
bdate=datetime(yr,1,1)-timedelta(days=1)
nd365=365+2
if yr%4==0:nd365=366+2
nty=nd365*24
dts=[bdate+timedelta(days=i/24.) for i in range(nty)]
73

讀取或產生nsc2及可能應用的縣市代碼CNTY對照表

#n_cnty: the dictionary of nsc2 vs cnties applied
if os.path.exists('n_cnty.json'):
  with open('n_cnty.json', 'r', newline='') as jsonfile:
    n_cnty=json.load(jsonfile)
else:
  for kc in ['cnty','kpq']:
    jsname='d_'+kc+'.json'
    if os.path.exists(jsname):
      with open(jsname, 'r', newline='') as jsonfile:
        exec('d_'+kc+'=json.load(jsonfile)')
    else:
      sys.exit('d_'+kc+'.json not found, please re-run prep_dfAdmw.py ')

對每種區域代碼組合進行迴圈

挑出符合對類別、建立二者的對照關係

存檔備用

  n_cntys={}
  for tl in s_rgn: #tuple of lengs
    n=d_rgn[tl]
    cntys=[]
    for t in 'mwd':
      exec('regs=list(df_A'+t+'.loc[df_A'+t+'.nsc2==n,"REGION"])')
      if len(regs)==0:continue
      for r in regs:
        if r in d_kpq:
          cntys+=[d_cnty[i] for i in d_kpq[r]]
        else:
          cntys+=[d_cnty[r]]
    cntys=list(set(cntys))
    for n in [i for i in n_rgn if n_rgn[i]==tl]:
      n_cntys[n]=cntys
  with open('n_cnty.json', 'w', newline='') as jsonfile:
    json.dump(n_cntys, jsonfile)
 103
 104

如果沒有nc_fac.json，建立逐時時間標籤的月份、星期、小時序列備用

if not os.path.exists('nc_fac.json'):
  mns=np.array([dts[i].month-1 for i in range(nty)])
  wks=np.array([dts[i].weekday() for i in range(nty)])
  hrs=np.array([dts[i].hour for i in range(nty)])
 109

將df_A?的2欄內容(df_A.nsc2,df_A.REGION)，轉變成單一組合(f_A[(n,c)])對照表

先將df_A拉長

  nts={'m':12,'w':7,'d':24}
  for t in 'mwd':
    exec('df_A=df_A'+t)
    df=DataFrame({})
    df_A['CNTY']=0
    for i in range(len(df_A)):
      cntys=[]
      a=df_A.loc[i]
      r=df_A.loc[i,'REGION']
      if r in d_kpq:
        cntys+=[d_cnty[i] for i in d_kpq[r]]
      else:
        cntys+=[d_cnty[r]]
      for c in cntys:
        a.CNTY=c
        df=df.append(a,ignore_index=True)

整併標籤與其後的時變係數欄位內容

    f_A={}
    for j in range(len(df)):
      n=df.loc[j,'nsc2']
      c=df.loc[j,'CNTY']
      f_A[(n,c)]=[df.loc[j,str(i)] for i in range(1,nts[t]+1)]
    exec('f_A'+t+'=f_A')

有月變化而無週變化者，週變化係數設為1.0、有月變化而無日變化者，日變化係數設為1.0、

  for i in set(f_Am)-set(f_Aw):
    f_Aw[i]=np.ones(shape=7)
  for i in set(f_Am)-set(f_Ad):
    f_Ad[i]=np.ones(shape=24)
 136

對每一個(nsc2,cnty)組合，計算全年的時變係數、存檔

  nc_fac={}
  for n in s_nsc2:
    for c in n_cntys[n]:
      tup=(n,c)
      tp=n+'_'+c
      if tup not in f_Am:continue
      lfac=np.array([f_Am[tup][m]*f_Aw[tup][w]*f_Ad[tup][d] for m,w,d in zip(mns,wks,hrs)])
      sfac=sum(lfac)
      nc_fac[tp]=list(lfac/sfac)
 146
  with open('nc_fac.json', 'w', newline='') as jsonfile:
    json.dump(nc_fac, jsonfile)
 149

檔案下載

Download: 環保署時變係數檔案:day.csv、mon.csv、week.csv、python程式：prep_dfAdmw.py、prep_df.py、prep_json.py

計算時間變異係數對照表

背景

時間變異係數檔案前處理(prep_dfAdmw.py)分段說明

副程式(prep_df.py)分段說明

nc_fac.json檔案產生程式(prep_json.py)分段說明

檔案下載

Reference

`nc_fac.json`檔案產生程式(prep_json.py)分段說明