煙囪及建築物座標系統之旋轉計算

Table of contents

背景
輸入檔
- 輸入檔說明
CGI_Python程式設計
RESULTS
- CaaS
- 結果檢視
Reference

背景

地圖數位板作業的結果是經緯度座標、而作為煙流模式前處理BPIP的輸入檔，必須是沒有歪斜的XY直角座標系統，需要進一步轉換。
尤有進者，實際個案的工廠北或多或少都會與背景地圖的真北有所偏差(夾角D)，因此都必須進行平移及旋轉，才能適用該前處理程式，並且不會發生疊圖誤差。
煙流模式對於風向是非常敏感的，風向的0度是真北，因此夾角D對煙流地面濃度的走向、建築物的影響等等，都有絕對的影響。經由數位板與此處的計算，可以降低量測與計算造成誤差的可能性。
除此之外，程式還提供由背景資料中找到個案所在地的煙囪基地高程E(步驟5)、煙囪頂的離地高度H(步驟6)等功能，以便進行計算。

輸入檔

此處設定為KML格式檔案，範例如example.kml
使用KML格式的理由詳見討論

輸入檔說明

必須資訊
- 煙囪位置(點)：將做為模擬個案之中心座標
- 囪煙名稱：此名稱將做為下述應用的個案名稱，
  - 污染源名稱、
  - 模式輸出檔案名稱
  - 氣象檔案名稱（適用於同一3Km網格範圍的其他個案）
輔助資訊
- 建築物頂點位置(多邊形)：方向不拘
- 建築物名稱(含高度、名稱高度間可以空格、逗號、分號、斜槓等區格)
- 煙囪如有高度跟在名稱之後，將會優先使用，否則會從TEDS資料庫中選取位置最接近的煙囪。
可以使用數位板點取座標值（詳見地圖數位板、另存成KML檔案）。
檢核：
- rotate_kml→旋轉角度，輸出成bpip輸入檔格式
- iscParser→解讀bpip輸入檔，輸出成kml格式(反解)
KML內之順序：(不限)
高度附掛單位（m or M)不影響計算，一律設定是m

CGI_Python程式設計

旋轉平移

對某一點旋轉座標：參考網友之程式碼
- 同一副程式應用在isc_parser.py系列、rd_kmlFull.py系列。

  def rotate_about_a_point(target_point,center_point,angle_rs):
    cp=np.subtract(target_point,center_point)
    px=cp[0]*math.cos(math.radians(angle_rs))+cp[1]*-math.sin(math.radians(angle_rs))
    py=cp[0]*math.sin(math.radians(angle_rs))+cp[1]*math.cos(math.radians(angle_rs))
    return(np.add([px,py],center_point))

先審視是否太過偏斜：
- 絕大工廠的廠房都是矩形，此處設定偏斜的標準為各建築物4面之角度的標準偏差>10度 (第75行)

dir=np.zeros(shape=(nplgs,4)) 
for i in range(nplgs): 
  diri=[90-math.atan2((y[i,j+1]-y[i,j]),(x[i,j+1]-x[i,j]))*180/math.pi for j in range(4)] 
  diri.sort() 
#  from North and clockwise 
  dir[i,:]=np.array(diri) 
if max( [np.std(dir[:,j]) for j in range(4)])>10: 
  print ('wrong direction or skewed!</br>') 
  for i in range(nplgs): 
    print (('dir for building# {:d} is: {:f} {:f} {:f} {:f}</br>').format(i,*dir[i,:])) 
  print ('</body></html>') 
  sys.exit('wrong direction or skewed!') 

找到最小旋轉角度
針對原點(KML輸入檔案的第1點，第85行)、做旋轉(第86行)、再做平移(89~90行)

P=[(i,j) for i,j in zip(x[:,:4].flatten(),y[:,:4].flatten())]
angl= min([np.mean(dir[:,j]) for j in range(4)])
if angl<0:angl+=360.
orig=P[0]
Pn=[rotate_about_a_point(pnt,orig,angl) for pnt in P]
Pn=np.array(Pn).flatten().reshape(nplgs,4,2)
#mnx, mny=(np.min(Pn[:,:,i]) for i in range(2))
Pn[:,:,0]+=-orig[0] #-mnx
Pn[:,:,1]+=-orig[1] #-mny

重整點位順序

校正數位化造成的不準度
讓bpip.inp內建築物的頂點為西南開始，逆時針順序輸入
找到最接近0點的頂點，做為原點，並補足每一點的不準度偏差(dx,dy)
煙囪點也同樣旋轉、平移(但沒有校正不準度)

for i in range(nplgs):
  xm,ym=np.mean(Pn[i,:,0]),np.mean(Pn[i,:,1])
  x1=sum([Pn[i,j,0] for j in range(4) if Pn[i,j,0] < xm])/2
  x2=sum([Pn[i,j,0] for j in range(4) if Pn[i,j,0] > xm])/2
  y1=sum([Pn[i,j,1] for j in range(4) if Pn[i,j,1] < ym])/2
  y2=sum([Pn[i,j,1] for j in range(4) if Pn[i,j,1] > ym])/2
  Pn[i,0,0],Pn[i,3,0]=x1,x1
  Pn[i,1,0],Pn[i,2,0]=x2,x2
  Pn[i,0,1],Pn[i,1,1]=y1,y1
  Pn[i,2,1],Pn[i,3,1]=y2,y2
  if i==0:
    dist=(Pn[i,:,0])**2+(Pn[i,:,1])**2
    idx=np.where(dist==np.min(dist))
    dx,dy=-Pn[i,idx[0],0],-Pn[i,idx[0],1]
Pn[:,:,0]=Pn[:,:,0]+dx
Pn[:,:,1]=Pn[:,:,1]+dy
...
x,y=pnyc(lonp,latp, inverse=False)
x+=Xcent
y+=Ycent
P=[(i,j) for i,j in zip(x,y)]
Pp=[rotate_about_a_point(pnt,orig,angl) for pnt in P]
Pp=np.array(Pp).flatten().reshape(npnts,2)
Pp[:,0]+=-orig[0] #-mnx
Pp[:,1]+=-orig[1] #-mny

從既有資料庫中找到需要的數據

基地高程：使用UCAR 333M解析度檔案(WPS之geo_em格式)

geo_name='/Users/WRF4.1/WPS/geo_em.d04_333m.nc'
...
nc = netCDF4.Dataset(geo_name, 'r')
v='HGT_M'
c=np.array(nc.variables[v][0,:,:])
for v in ['CLAT','CLONG']:
    exec(v+'=nc.variables[v][0,:,:]')
xg,yg=pnyc(CLONG,CLAT, inverse=False)
xg+=Xcent
yg+=Ycent
base=[]
for ii in range(nplgs):
  i=ii*4
  d=(xg-P[i][0])*(xg-P[i][0])+(yg-P[i][1])*(yg-P[i][1])
  idx=np.where(d==np.min(d))
  base.append(c[idx[0][0],idx[1][0]])

從TEDS資料庫中找到煙囪高度

tedsp_name=WEB+'isc_results/point_QC.csv'
...
#the stack heights are read from TEDS database IF that hgts are not contained in the name strings
df=read_csv(tedsp_name)
df.UTM_E+=Xcent
df.UTM_N+=Ycent
a=[]
for ll in range(1,6):
  L=ll*1000
  a=df.loc[df.UTM_E.map(lambda s:abs(s-P[0][0])<L) & df.UTM_N.map(lambda s:abs(s-P[0][1])<L)]
  if len(a)>0:
    df=a
    break
if len(a)==0:
  print ('the point source seems not existing in database. </body></html>')
  sys.exit('fail')
cole=['CO_EMI', 'NMHC_EMI', 'NOX_EMI', 'PM25_EMI', 'PM_EMI', 'SOX_EMI']
c2m={'SOX':64,'NOX':46,'CO':28,'PM25':24.5,'PM':24.5,'NMHC':12*4+10}
unit={i:'ppb' for i in c2m if 'PM' not in i}
unit.update({i:'ug/m3' for i in c2m if 'PM' in i})
hdtv=[ 'HEI', 'DIA', 'TEMP', 'VEL']
tims=[ 'DY1', 'HD1', 'HY1']
for v in cole+hdtv+tims:
  exec(v+'=[]')
for k in range(npnts):
  df['dist']=[np.sqrt((i-P[k][0])**2+(j-P[k][1])**2) for i,j in zip(list(df.UTM_E),list(df.UTM_N))]
  idx=df.loc[df.dist==min(df.dist)].index
  if len(idx)>1:
    idx=df.loc[idx].sort_values('HEI',ascending=False).head(1).index
  for v in cole+hdtv+tims:
    exec(v+'.append(list(df.'+v+'[idx])[0])')
  if len(hgts)<nplgs+npnts:
    hgts.append(list(df.HEI[idx])[0])
for v in cole:
  exec(v+'=['+v+'[i]/HY1[i]/3.6*1000. for i in range(npnts)]')
TEMP=[i+273 for i in TEMP]

輸出到bpip.inp檔案

命名為fort.10，以符合BPIP程式的約定。以隨機產生的目錄名稱做為使用者個案間的辨識。

with open(pth+'fort.10','w') as f:
  f.write(("'BPIP input file with "+'{:2d}'+' bldg and '+'{:2d}'+" stacks,originated at [{:.1f},{:.1f}](TWD97m).'\n")
  .format(nplgs,npnts,orig[0],orig[1]))
  f.write(("'P'\n"+"'METERS' 1.00\n'UTMN',"+'{:5.0f}\n').format(angl))
  f.write(('{:2d}\n').format(nplgs))
  for i in range(nplgs):
    f.write(("'"+names[i]+"' 1"+'{:6.1f}\n').format(base[i]))
    f.write(('4 '+'{:5.0f}\n').format(hgts[i]))
    for j in range(4):
      f.write(('{:5.1f}  {:5.1f}\n').format(Pn[i,j,0],Pn[i,j,1]))
  f.write(('{:2d}\n').format(npnts))
  for i in range(npnts):
    ii=i+nplgs
    f.write(("'"+names[ii]+"' "+'{:6.1f} {:6.1f} {:6.1f} {:6.1f} \n').format(base[ii],hgts[ii],Pp[i,0],Pp[i,1]))

CGI_Python下載

rotate_kml.py

RESULTS

CaaS

Rotate a KML file

結果檢視

檢討由KML產生bpip.inp，再由iscParser檢查，結果如下(TSMC.kml)
- 並無系統性之錯誤
- 仍然存在描圖的不準度，然誤差已在可接受範圍。
- 建築物輪廓線尚可保持方向之一致性


數位化、且經旋轉、歪斜修正後的廠房位置圖

Reference

Mathematics Stack Exchange, calculate new positon of rectangle corners based on angle