town_aqst程式說明

Table of contents

背景
town_aqst.py程式說明
- 模組及IO
中央氣象局測站版本

背景

town_aqst.py這支程式的基本構想，是將環保署測站與鄉鎮區進行關聯，不論是鄉鎮區內的平均值、或者是鄰近鄉鎮區之外插，都非常需要。程式計算的結果產出為town_aqstEnew.csv。
因台灣地區鄉鎮區的分界線，與水系、山系的自然分野有關，此一邏輯概念可以擴展到其他觀測項目的空間內外插過程。
這支程式的演進有好幾個版本，最先是用在面源的測站調整，讓測站附近的排放量可以反映該測站的濃度變化，過程會需要建立測站與周圍鄉鎮區的關聯。其後因台灣地區鄉鎮區代碼更新，因而舊版不敷使用而更新。
- 2018/05 master:/home/camxruns/2013/emis/8760wt/area/town_aqst.py
- 2019/10 /nas1/TEDS/teds10_camx/HourlyWeighted/area/town_aqst.py
town_aqstEnew.csv用在下列程式中

town_aqst.py程式說明

模組及IO

系列作業本來源自於TEDS面源的內容，但因行政區整併、局部名稱更換等等因素，因此有用到後來產生的結果檔案進行修正，或局部修正。
輸入檔
- 'adj_dict.json'這個鄉鎮區鄰近其他鄉鎮區的名稱對照，是adj_dict.py所產生的。
- town2.csv'：為鄉鎮區名稱、含縣市名稱
- "sta_dict.csv"：測站編號、所在網格座標位置(1公里解析度)與所在縣市鄉鎮區名稱(漢語拼音)

#kuang@master /nas1/TEDS/teds10_camx/HourlyWeighted/area
#$ cat town_aqst.py
from pandas import *
import json
df_town=read_csv('town2.csv')
df_st=read_csv("sta_dict.csv")
fname='adj_dict.json'
fn=open(fname,'r')
adj_dict=json.load(fn)
ddict=set(df_st['dict'])
df_cnty=df_town
df_cnty['aq_st']=['' for i in xrange(len(df_cnty))]

逐一鄉鎮區檢討
- 如果行政區範圍內有測站，則在資料庫中紀錄

for i in xrange(len(df_cnty)):
    dd=df_cnty.loc[i,'Name']
    if dd in ddict:
        s=str(list(df_st.loc[df_st['dict']==dd,'idx'])[0])+';'
        df_cnty.loc[i,'aq_st']=s

如果沒有。則檢討周圍行政區是否有測站
- 如果有，則記錄測站。紀錄時先檢查之前是否已經包含了其他測站。

    else:
        for jj in adj_dict[dd]:
            if jj in ddict:
                s=str(list(df_st.loc[df_st['dict']==jj,'idx'])[0])+';'
                if len(df_cnty.loc[i,'aq_st'])==0:
                    df_cnty.loc[i,'aq_st']=s
                else:
                    df_cnty.loc[i,'aq_st']=df_cnty.loc[i,'aq_st']+s

如果都沒有測站：則記錄’0;’

    if len(df_cnty.loc[i,'aq_st'])==0:
         df_cnty.loc[i,'aq_st']='0'+';'
df_cnty.set_index('code').to_csv('town_aqst.csv')

中央氣象局測站版本

逐點比對法

這個版本的town_cwb.py與前述town_aqst.py非常相似，只是應用在中央氣象局的自記式測站。
比較這2個版本的差異如下
檔案位置與名稱

kuang@master /home/backup/data/cwb/e-service/read_web
$ diff town_cwb.py /nas1/TEDS/teds10_camx/HourlyWeighted/area/town_aqst.py
1d0
<
4,7c3,5
< path='/nas1/TEDS/teds10_camx/HourlyWeighted/area'
< df_town=read_csv(path+'/town3.csv',encoding='big5')
< df_st=read_csv("stats_dict.csv")
< fname=path+'/adj_dict.json'
---
> df_town=read_csv('town2.csv')
> df_st=read_csv("sta_dict.csv")
> fname='adj_dict.json'

py27與py37的差異

12,13c10,11
< df_cnty['aq_st']=''
< for i in range(len(df_cnty)):
---
> df_cnty['aq_st']=['' for i in xrange(len(df_cnty))]
> for i in xrange(len(df_cnty)):

2個資料表的欄位名稱(測站編號、行政區名稱)有點差異

10c8
< ddict=set(df_st['name'])
---
> ddict=set(df_st['dict'])
16c14
<         s=str(list(df_st.loc[df_st['name']==dd,'stno'])[0])+';'
---
>         s=str(list(df_st.loc[df_st['dict']==dd,'idx'])[0])+';'
21c19
<                 s=str(list(df_st.loc[df_st['name']==jj,'stno'])[0])+';'
---
>                 s=str(list(df_st.loc[df_st['dict']==jj,'idx'])[0])+';'

去除沒有測站的鄉鎮區。存檔

28,29c26,30
< df_cnty=df_cnty.loc[df_cnty.aq_st.map(lambda x:x!='0;')].reset_index(drop=True)
< df_cnty.set_index('code').to_csv('town_cwb.csv')
---
> df_cnty.set_index('code').to_csv('town_aqst.csv')
>

新版鄉鎮區代碼

此版本直接使用內政部鄉鎮區邊界之多邊形、以及鄉鎮區代碼來計算
使用shapely Point.within()函數來判斷測站落在哪一個行政區

kuang@master /home/backup/data/cwb/e-service/read_web
$ cat stats_dict2.py
from shapely.geometry import Polygon,Point
from pandas import *
import numpy as np
import sys
import os
import json

fname='/home/QGIS/Data/TWN_town/polygons.csv'
df=read_csv(fname,encoding='big5')
df['lonlats']=[j.replace(',','').replace(')','').replace('(','').replace('[','').replace(']','').split() for j in df.lonlats]
df['lonlats']=[[float(i) for i in j] for j in df.lonlats]
df['lonlats']=[[(j[i],j[i+1]) for i in range(0,len(j),2)] for j in df.lonlats]
df['polygon']=[Polygon(i) for i in df.lonlats]

fname='stats_tab.csv'
st=read_csv(fname)
pnts=[Point(i,j) for i,j in zip(st.LON,st.LAT)]

TOWNENG,twnid=[],[]
for p in pnts:
  sea=True
  for j in range(len(df)):
    if p.within(df['polygon'][j]):
      TOWNENG.append(df.TOWNENG[j])
      twnid.append(df.twnid[j])
      sea=False
      break
  if sea:
    TOWNENG.append('sea')
    twnid.append(5300)
st['TOWNENG']=TOWNENG
st['twnid']=twnid
cols=['stno','stat_name','TOWNENG','twnid']
st[cols].set_index('stno').to_csv('stats_dict.csv')

town_cwb.py

這個最後版本使用intersects與buffer來擴展鄉鎮區的代表範圍(詳見附近行政區之定位#intersects與buffer(最後版本))，其結果檔為adj_dict.json。
需要前面計算出來的'stats_dict.csv'檔案、以及內政部的原始檔TOWN_MOI_1090727E.csv
將會產出每個鄉鎮區與其測站名稱之對照關係(town_cwb.csv)，以便進行內積計算。
迴圈的設計與環保署測站版本略有差異，後者在對照不到測站名稱時，才在外圈尋找。此處不論是否找得到都搜尋一遍，而將是否重複勘列，放在每一次加入之前來做判斷。

...
for i in range(len(df_cnty)):
  dd=str(df_cnty.loc[i,'TOWNCODE'])
  if dd in ddict:
    stnos=list(df_st.loc[df_st['twnid']==int(dd),'stno'])
    if len(stnos)==0:sys.exit('fail')
    for s in stnos:
      if s not in df_cnty.loc[i,'aq_st']:df_cnty.loc[i,'aq_st']+=s+';'
for i in range(len(df_cnty)):
  dd=str(df_cnty.loc[i,'TOWNCODE'])
  if dd in adj_dict:
    for jj in adj_dict[dd.decode('utf-8')].split(';'):
      j=str(jj)
      if len(j)==0:continue
      if j in ddict:
        stnos=list(df_st.loc[df_st['twnid']==int(j),'stno'])
        if len(stnos)==0:sys.exit('fail')
        for s in stnos:
          if s not in df_cnty.loc[i,'aq_st']:df_cnty.loc[i,'aq_st']+=s+';'
...

程式下載

Download: town_cwb2.py