由CWB數據計算軌跡

Table of contents

背景
內插程式說明
軌跡程式說明
- csv2bln
成果檢討
程式原始碼
Reference

背景

空間尺度之界定

逆(反)軌跡圖是空氣污染來源追蹤、溯源解析過程中經常使用的圖像工具。按照污染源與受體空氣品質關係的遠近，分析的尺度有：

近距離、~10公里範圍，高斯煙流擴散模式可以應用之範圍，時間約在小時解析度範圍。可以用代表性地面測站之風向直接進行(統計)解析研判。
城市~地區尺度，中距離約數10~百公里範圍，受海陸風影響的平坦~平緩地區，時間範圍約在1日~3日之間。可以用地面站網、氣象模式進行解析。
地區~長程傳輸現象，約數百~7千公里範圍，除前述現象外，也受到天氣現象的顯著影響，時間範圍約3日~週之間。須以氣象模式、HYSPLIT等模式進行三維風場與軌跡解析。此處要處理的是城市~地區尺度，因此需應用全臺自動站的風速風向數據。

軌跡正確性與風場模式

軌跡正確性的關鍵在於風場，因此有以：

高解析度觀測站數據內插、由於風場的正確性與測站所函蓋的範圍有關，密度較高的平地範圍，會有較高的正確性，而在海上或山區等測站密度較低範圍，可能有較低的正確性。
數值氣象預報模式產品、如WRF，有解析度與範圍的限制。一般WRF最高解析度為3公里
客觀分析等不同方式產生風場。(介於前2者之間)

策略方案

CODiS提供的數據至少有300個測站是有風速、風向數據
- 全台面積36,193平方公里，平均一站分攤119平方公里，約為11公里X11公里之解析度，平地還有更高的密度
- 用以內插建立風場，應有其充分性與正確性
CWB提供有WRF_3KM預報風場，亦足以代表海上、山區等測站較少範圍
建立內插方法
- 有測站範圍地區則以測站值為主
- 海上則以模式分析結果，做為3公里X3公里之虛擬測站，納入相同機制進行計算。
全區風場計算儲存的必要性
- 全區風場事先計算儲存看似方便，然佔用空間龐大，因電腦速度提高了，因此在實際計算時似乎也沒有提高太多效率。
- 主要由於軌跡線經過網格數有限，實在沒有必要進行全臺風場之計算或儲存。
發展路徑
- 直接內插觀測值
- 直接內插模擬值、FDDA模擬值
- 以WRFDA或其他模式合併觀測及模式

內插程式說明

見計算網格點與CODiS測站間距離反比加權值

軌跡程式說明

軌跡程式碼可以由githup下載，此處分段說明如下
引用模組

#!/opt/anaconda3/envs/py27/bin/python
import numpy as np
from pandas import *
import os, sys, subprocess, time, json
from scipy.io import FortranFile
from datetime import datetime, timedelta
import twd97
from pyproj import Proj


讀取引數的副程式
- -t 測站名稱，可以是環保署測站（單站）、經緯度組合、twd97座標組合，組合以逗號,區隔。
- -d 年月日時共10碼
- -b 是否是反軌跡（True or False）

def getarg():
 """ read time period and station name from argument(std input)
 traj2kml.py -t daliao -d 2017123101 """
 import argparse
 ap = argparse.ArgumentParser()
 ap.add_argument("-t", "--STNAM", required=True, type=str, help="station name(num),sep by ';' ,or Lat,Lon")
 ap.add_argument("-d", "--DATE", required=True, type=str, help="yyyymmddhh")
 ap.add_argument("-b", "--BACK", required=True, type=str, help="True or False")
 args = vars(ap.parse_args())
 return [args['STNAM'], args['DATE'],args['BACK']]


解析引數中的判別(布林值)，使用者可以輸入任何大小寫的yes, true, t, y, 1等等。

def str2bool(v):
   if isinstance(v, bool):
      return v
   if v.lower() in ('yes', 'true', 't', 'y', '1'):
       return True
   elif v.lower() in ('no', 'false', 'f', 'n', '0'):
       return False
   else:
       raise argparse.ArgumentTypeError('Boolean value expected.')


讀取環保署空品測站站名與編號對照表

def nstnam():
 import json
 fn = open(path+'sta_list.json')
 d_nstnam = json.load(fn)
 d_namnst = {v: k for k, v in d_nstnam.items()}
 return (d_nstnam, d_namnst)



由風速、風向計算U、V值

def ws_uv(ws, wd):
 PAI = np.pi
 RAD = (270. - wd) * PAI / 180.0
 u = ws * np.cos(RAD)
 v = ws * np.sin(RAD)
 return u, v



判斷軌跡點是否超出範圍界線

def beyond(xpp, ypp):
 xp_km, yp_km = int(xpp // 1000), int(ypp // 1000)
 boo = not ((xp_km - x_mesh[0]) * (xp_km - x_mesh[-1]) < 0 and (yp_km - y_mesh[0]) * (yp_km - y_mesh[-1]) < 0)
 return [boo, (xp_km, yp_km)]



開啟CODiS全台自動氣象站數據檔。事先下載儲存在年代目錄下。

def opendf(pdate):
 ymd = pdate.strftime('%Y%m%d')
 fname= path+'../' + ymd[:4] + '/cwb' + ymd + '.csv'
 if not os.path.isfile(fname):
   dfT=DataFrame({})
   print ('no file for '+fname)
 else:
   try:
     dfT = read_csv(fname)
     dfT['stno'] = [i[:6] for i in dfT.stno_name]
     dfT = dfT.loc[dfT.stno.map(lambda x: x in stno)].reset_index(drop=True)
#	  (dfT.WS>0)& 
     dfT = dfT.fillna(0)
     ws, wd = np.array(dfT.WS), np.array(dfT.WD)
     uv = np.array([ws_uv(i, j) for i, j in zip(ws, wd)])
     dfT['u'], dfT['v'] = (uv[:, i] for i in [0, 1])
     dfT.ObsTime = [int(i)-1 for i in dfT.ObsTime]
   except:
     dfT=DataFrame({})
 return dfT, ymd


計算u,v之加權平均值

def uvb(r,u,v):
 dfuv=DataFrame({'R':r,'u':u,'v':v})
 dfuv=dfuv.sort_values('R',ascending=False).reset_index(drop=True)
 rr,uu,vv=np.array(dfuv.R)[:ns3],np.array(dfuv.u)[:ns3],np.array(dfuv.v)[:ns3]
 rr=rr/sum(rr)
 ub,vb=sum(rr*uu),sum(rr*vv)
 return ub,vb



主程式
- 設定Proj座標轉換工具pnyc。因本次使用CODiS數據，還是以臺灣本島為主，因此主要還是使用twd97系統。
- 設定網格數及測站數nx, ny, ns
- 讀取權重計算結果
- 讀取環保署空品測站站名與編號對照表
- 讀取並解析引數中的布林(BACK)與日期(DATE)，並計算起訖時間(bdate,next_date)。
- 解析(nam)站名

path='/Users/Data/cwb/e-service/surf_trj/'

Latitude_Pole, Longitude_Pole = 23.61000, 120.9900
Xcent, Ycent = twd97.fromwgs84(Latitude_Pole, Longitude_Pole)
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)


# restore the matrix
nx, ny, ns = 252, 414, 431
fnameO = path+'R%d_%d_%d.bin' % (ny, nx, ns)
with FortranFile(fnameO, 'r') as f:
 R2 = f.read_record(dtype=np.float64)
R2 = R2.reshape(ny, nx, ns)
with FortranFile(path+'x_mesh.bin', 'r') as f:
 x_mesh = list(f.read_record(dtype=np.int64))
with FortranFile(path+'y_mesh.bin', 'r') as f:
 y_mesh = list(f.read_record(dtype=np.int64))

(d_nstnam, d_namnst) = nstnam()
stnam, DATE, BACK = getarg()
BACK=str2bool(BACK)
BF=-1
if not BACK:BF=1
bdate = datetime(int(DATE[:4]), int(DATE[4:6]), int(DATE[6:8]), int(DATE[8:]))
next_date= bdate + timedelta(hours=96*BF)
nam = [i for i in stnam.split(';')]

解析軌跡起始點(twd97系統)
- 測站名稱可以是環保署測站、經緯度組合、twd97座標組合
- 環保署測站：由sta_ll.csv中讀取位置
- 經緯度組合：順序為(維度,經度)，
- twd97座標組合：直接使用

if len(nam) > 1:
 try:
   lat = float(nam[0])
   lon = float(nam[1])
 except:
   sys.exit('more than two station, suggest executing iteratively')
 else:
   # in case of lat,lon
   if lat < 90.:
     xy0 = twd97.fromwgs84(lat,lon)
     x0, y0 =([xy0[i]] for i in [0,1])
     nam[0] = str(round(lat,2))+'_'+str(round(lon,2))+'_'
   #   in case of twd97_x,y
   else:
     # test the coordinate unit
     if lat>1000.:
       x0, y0 = [lat],[lon]
       nam[0] = str(int(lat/1000))+'+'+str(int(lon/1000))+'_'
     else:
       x0, y0 = [lat*1000],[lon*1000]
       nam[0] = str(int(lat))+'_'+str(int(lon))+'_'

# len(nam)==1, read the location from csv files
else:
 for stnam in nam:
   try:
     nst=[int(stnam)]
   except:
     astnam=stnam  
     if astnam not in d_namnst: sys.exit("station name not right: " + stnam)
     nst = [int(d_namnst[i]) for i in nam]
   else:	
     if stnam not in d_nstnam: sys.exit("station number not right: " + stnam)
     nam[0]=d_nstnam[stnam]	  
 # locations of air quality stations
 # read from the EPA web.sprx
 fname = path+'sta_ll.csv'
 sta_list = read_csv(fname)
 x0, y0 = [], []
 for s in nst:
   sta1 = sta_list.loc[sta_list.ID == s].reset_index(drop=True)
   xx,yy=pnyc(list(sta1['lon'])[0],list(sta1['lat'])[0], inverse=False)
   print xx,yy
   x0.append(xx+Xcent) #list(sta1['twd_x'])[0]-Xcent)
   y0.append(yy+Ycent) #list(sta1['twd_y'])[0]-Ycent)


軌跡計算初始化
- 讀取CODiS測站位置stat_wnd.csv
- 軌跡點間距時間delt：15秒
- 讀取當日CODiS測站數據opendf(pdate)

xp, yp = x0, y0
dfS = read_csv(path+'stat_wnd.csv')
if len(dfS) != ns: sys.exit('ns not right')
stno = list(dfS.stno)
pdate = bdate
df, ymd0 = opendf(pdate)
if len(df)==0:sys.exit('no cwb data for date of:'+ymd0)
delt = 15
s = 0
o_ymdh,o_time,o_xp,o_yp,l_xp,l_yp=[],[],[],[],[],[]
itime=0
ymdh=int(DATE)
o_ymdh.append(nam[0]+'@'+pdate.strftime('%Y/%m/%d/%H:00'))
o_time.append('hour='+str(itime))
o_xp.append(xp[s])
o_yp.append(yp[s])
l_xp.append(xp[s])
l_yp.append(yp[s])
ns3=int(ns)

開始計算軌跡點，知道超出範圍界線
- 由df中篩出個別小時數據，另存成df1，如有缺漏，則補以0。
- 由df1中提取u, v

while not beyond(xp[s], yp[s])[0] and len(df)!=0:
 boo=pdate>next_date	
 if not BACK:boo=pdate<next_date
 if not boo:break
 df1 = df.loc[(df.ObsTime == ymdh) & (df.stno.map(lambda x:x in stno))].reset_index(drop=True)
 df1 = df1.drop_duplicates()
 ldf1=len(df1)
 if ldf1 < ns:
   if boo:
     ns2 = set(df1.stno)
     miss = set(stno) - set(ns2)
     if len(miss)!=0:
       for m in miss:
         df2= DataFrame({'stno_name':[m],'ObsTime':[ymdh]})
         df1=df1.append(df2,ignore_index=True, sort=False)
   else:
     print 'df1 not right' + str(ymdh)
     break
 df1=df1.sort_values('stno_name').reset_index(drop=True)
 df1 = df1.fillna(0)
 u, v = np.array(df1.u), np.array(list(df1.v))

執行該小時軌跡點的計算
- 由位置找到加權R2[iy, ix, :]
- 進行內插ub, vb = uvb(R2[iy, ix, :],u,v)
- 計算下一點xp[s], yp[s] = xp[s]+BF*delt * ub, yp[s]+BF*delt * vb
- 儲存結果

 for sec in range(0, 3601, delt):
   boo, (xp_km, yp_km) = beyond(xp[s], yp[s])
   if boo: break
   ix, iy = x_mesh.index(xp_km), y_mesh.index(yp_km)
   if sec == 0:
     ix0, iy0 = ix, iy
     ub, vb = uvb(R2[iy, ix, :],u,v)
   else:
     if ix0 != ix or iy0 != iy:
       # ub, vb = sum(R2[iy, ix, :] * u), sum(R2[iy, ix, :] * v)
       ub, vb = uvb(R2[iy, ix, :],u,v)
       ix0, iy0 = ix, iy
   xp[s], yp[s] = xp[s]+BF*delt * ub, yp[s]+BF*delt * vb
   l_xp.append(xp[s])	
   l_yp.append(yp[s])	
 pdate = pdate + timedelta(hours=BF)
 ymdh = int(pdate.strftime('%Y%m%d%H'))
 itime+=1
 o_ymdh.append(pdate.strftime('%Y/%m/%d/%H:00'))
 o_time.append('hour='+str(itime))
 o_xp.append(xp[s])
 o_yp.append(yp[s])
 if pdate.strftime('%Y%m%d') != ymd0:
   df, ymd0 = opendf(pdate)
   if len(df)==0:break

儲存逐時軌跡點檔案(twd97座標值)

print('beyond:',beyond(xp[s], yp[s])[0],'len(df)=',len(df))
df=DataFrame({'ymdh':o_ymdh,'xp':o_xp,'yp':o_yp,'Hour':o_time})
col=['xp','yp','Hour','ymdh']
dr='b'
if BF!=-1:dr='f'
name='trj_results/'+dr+'trj'+nam[0]+DATE+'.csv'
df[col].set_index('xp').to_csv(name)

計算並寫出經緯度值、寫出測站名稱

#geodetic LL
x,y=np.array(o_xp)-Xcent,np.array(o_yp)-Ycent
lon, lat = pnyc(x, y, inverse=True)
dfg=DataFrame({'lon':lon,'lat':lat})
dfg.set_index('lon').to_csv(name.replace('.csv','_mark.csv'),header=None)
with open('trj_results/filename.txt','w') as f:
 f.write(name.split('/')[1])


儲存逐點結果

# output the line segments for each delta_t
dfL=DataFrame({'TWD97_x':l_xp,'TWD97_y':l_yp})
dfL.set_index('TWD97_x').to_csv(name.replace('.csv','L.csv'))
#geodetic LL
x,y=np.array(l_xp)-Xcent,np.array(l_yp)-Ycent
lon, lat = pnyc(x, y, inverse=True)
dfg=DataFrame({'lon':lon,'lat':lat})
dfg.set_index('lon').to_csv(name.replace('.csv','_line.csv'),header=None)


呼叫外部程式轉換成KML檔案及bln檔(for SURFER)
- csv2kml.py
- KML是Google地圖、Open Street Map等地圖系統相容的檔案格式

#make kml file
dir='NC'
if not BACK:dir='RC'
os.system('/opt/local/bin/csv2kml.py -f '+name+' -n '+dir+' -g TWD97')
os.system('/opt/local/bin/csv2bln.cs '+name)

csv2bln

csv2bln.cs為下列腳本。

echo $(( $(wc -l $1|/opt/local/bin/awkk 1) - 1 )) > $1.bln
sed 1d $1 >> $1.bln

alias awkk=awk '{print $'$1'}'

成果檢討

2018/10/27 12時林園vs當天2時北高雄某廠燃燒塔之正軌跡(紅色)及林園測站反軌跡(白色)
燃燒塔正軌跡分析證實了以北高雄燃燒塔對反應生成測站O₃高值有較大的可能性，此處分析其他事件的情況，並將測站O₃尖峰時間的反軌跡繪出以進行比較驗證。
- 圖中白色點線為林園測站O₃尖峰時間之逆軌跡，
- 北方工業區燃燒塔以當天凌晨CEMS記錄最大流量發生時間，做為之正軌跡線的起始時間，圖中為紅色點線表示。
由圖中可以發現，紅、白2條軌跡線在高雄市西南側沿海與近海地區大多呈現平行運動，間距約為2~4公里，約為一般網格模式解析度範圍，應為污染源的精確位置、或風場模式內插所造成的誤差。
然而就污染物受日夜海陸風的性質而言，圖中正、反軌跡線可以確認造成林園測站O₃高值的污染源，非常可能就是仁大工業區之燃燒塔排放。

程式原始碼

可以在github找到:

python程式

Download: traj2kml.py

Reference

MM5/WRF之little_r格式
NOAA, HYSPLIT
Jimy Dudhia， WRF Four Dimensional Data Assimilation (FDDA), documen.site, May 12, 2018
Tom.Chen, Python converter between TWD97 and WGS84, pypi.org, Oct 22, 2014
Wiki, Keyhole Markup Language, wikipedia.org,last edited on 24 October 2021