空品測站之Voronoi圖

Table of contents

背景
程式說明
分區之應用
- 公版模式範圍1公里解析度網格之分區
- 分區結果

背景

測站網絡內之個別站能代表甚麼範圍的地區，即測站代表性問題，是測站測值解釋的重要議題。此處要運用的技巧稱之為「Voronoi沃羅諾伊」圖。
- 詳細發展歷程見wiki
- 點狀資訊將會影響到著鄰近的資訊，所以「最靠近」、「距離最短」之類的問題，多半可以透過Voronoi Diagram 解決，如幾何、晶體學、建築學、地理學、氣象學、信息系統等許多領域有廣泛的應用。如空氣品質測站¹²³⁴。
Instance：
1. Sinica 勢力分佈圖 (Voronoi Diagram)@PM2.5 開放資料入口網站
2. Visualising air quality data with Voronoi diagrams@MODULO ERRORS
3. Mapping Air Quality Data with D3js - VORONOI by Ian Johnson(2020)

程式說明

程式運用到第3方模組geovoronoi⁵，綁定版本python 3.1 <= 3.6。過高版本將無法執行。

相依性

import numpy as np
import geopandas as gpd
import pandas as pd
import matplotlib.pyplot as plt
from shapely.ops import cascaded_union, unary_union
from shapely.geometry import Point
from geovoronoi.plotting import subplot_for_map, plot_voronoi_polys_with_points_in_area
from geovoronoi import voronoi_regions_from_coords, points_to_coords

安裝過程有些順序需要遵守

pip install --user pyproj
pip install --user fiona
pip install --user geovoronoi
pip install --user descartes

python3.6早期的netCDF安裝過程需指定HDF5及netCDF4的包括檔頭(header of include files)環境變數，否則裝不起來

export NETCDF4_DIR=/opt/netcdf/netcdf4_gcc
export HDF5_DIR=/opt/hdf/hdf5_gcc
pip install netCDF4

IO’s

讀取內政部之縣市界shp檔。從中讀取本島海岸線範圍(島內縣市範圍之聯集)。
- 採排他法：去除金門澎湖與連江範圍之物件
- 本島縣市也可能有離島(綠島、蘭嶼、龜山、小琉球等)，因此以最大面積之多邊形為代表。這個邏輯需要每個縣市一一檢討。

shp_fname="/home/kuang/bin/TWN_COUNTY.shp"
gdf = gpd.read_file(shp_fname)

CNTYNAM=set(gdf.COUNTYNAME)-{'金門縣','澎湖縣','連江縣'}
ifirst=1
for c in list(CNTYNAM)[:]:
    a=gdf.loc[gdf.COUNTYNAME==c].reset_index(drop=True)
    area=[i.area for i in a.geometry]
    imax=area.index(max(area))
    if len(a)==1:
        b=a.to_crs(epsg=4326)
    else:
        b=a.loc[a.index==imax].reset_index(drop=True).to_crs(epsg=4326)
    if ifirst==1:
        df0=b.to_crs(epsg=4326)
        ifirst=0
    else:
        df0=gpd.GeoDataFrame(pd.concat([df0,b],ignore_index=True))

讀取測站經緯度csv檔案
- 並將資料表寫成GeoDataFrame方便後續繪圖

stn=pd.read_csv('/nas1/cmaqruns/2016base/data/sites/sta_ll.csv')
stnpnt=[Point(i,j) for i,j in zip(stn.lon,stn.lat)]

for i in range(len(stn)):
    b=gpd.GeoDataFrame({'COUNTYSN':stn.loc[i,'ID'] ,'COUNTYNAME':stn.loc[i,'New'],'geometry':[stnpnt[i]]})
    df0=gpd.GeoDataFrame(pd.concat([df0,b],ignore_index=True))

df1=df0.loc[:21]
df1.to_file('mainisland.shp',mode='w')
df2=df0.loc[22:]
df1.to_file('stn.shp',mode='w')

取本島海岸線範圍內的測點

理論上plot_voronoi_polys_with_points_in_area似乎會直接繪製範圍界線內的點，但似乎程式不允許界線外存在有點。
使用within函數來判斷。(與邊界的距離長短似乎沒有作用)

boundary = gpd.read_file("mainisland.shp")
boundary = boundary.to_crs(epsg=4326)
boundary_shape = unary_union(boundary.geometry)

df2=df2.reset_index(drop=True)
for i in range(len(df2)):
    p=df2.loc[i,'geometry']
    if not p.within(boundary_shape): # or boundary_shape.exterior.distance(p) < 0.01:
        df2=df2.drop(i)
df2=df2.reset_index(drop=True)

執行第3方模組

points_to_coords這個函數是否有必要值得挑戰，檢討原始碼發現也沒有太大的作用，就是將geometry形式的物件改成array而已。
此處主要呼叫voronoi_regions_from_coords，會輸出2個字典(dict)，編號為其內定順序，並無意義，只是用來串聯點及多邊形。
- 如果需要測點原來的ID，還是需要從測站經緯度csv檔案來連結。

gdf_proj = df2.to_crs(boundary.crs)
coords = points_to_coords(gdf_proj.geometry)
region_polys, region_pts = voronoi_regions_from_coords(coords, boundary_shape)

繪製本島之邊框範圍

做為測試邊框及內部測點

boundary = gpd.read_file("mainisland.shp")
fig, ax = plt.subplots(figsize=(12, 10))
boundary.plot(ax=ax, color="gray")
df1.plot(ax=ax, markersize=3.5)#, color="brown")
df2.plot(ax=ax, markersize=3.5, color="brown")
ax.axis("off")
plt.axis("equal")
plt.show()

繪製最終結果

fig, ax = subplot_for_map(figsize=(12, 10))
plot_voronoi_polys_with_points_in_area(ax, boundary_shape, region_polys, coords, region_pts)
ax.set_title('Voronoi regions of Air Quality Station Networks')
plt.tight_layout()
plt.show()


環保署測站	環保署+微型感測

大致還能保持方、圓形為主的分布，表示測站在範圍內還算均勻分配。
細長條情況仍然存在(如屏東站)，然較Sinica少很多。

程式下載

Download: 從shp檔繪製空品測站的Voronoi分區圖：Voronoi.py

分區之應用

公版模式範圍1公里解析度網格之分區

網格座標詳mk_gridLL
同樣使用shapely.geometry.Point的內設函數within來判斷。

ll=pd.read_csv('/nas2/cmaqruns/2019TZPP/output/Annual/aTZPP/LGHAP.PM25.D001/gridLL.csv')
ll['AQID']=0
for i in range(len(ll)):
    p=ll.Point[i]
    p=Point([float(i) for i in p.replace('(','').strip(')').split()[1:]])
    if not p.within(boundary_shape):continue
    j=0
    for b in dfv.geometry:
        if p.within(b):
            ll.AQID[i]=dfv.COUNTYSN[j]
            break
        j+=1

import netCDF4
fname='tempTW.nc'
nc = netCDF4.Dataset(fname,'r+')
nc['PM25_TOT'][0,0,:,:]=np.array(ll.AQID).reshape(393,276)

分區結果

由於縣市界多半是以山稜分水嶺等自然界限為準，北宜、竹宜、中投、投花等交界似乎還能符合。


測站Voronoi分區圖	鄉鎮區範圍平均後之分布

Ditsuhi Iskandaryan(2023) Study and Prediction of Air Quality in Smart Cities through Machine Learning Techniques Considering Spatiotemporal Components, A dissertation presented for the degree of Doctor of Computer Science, Universitat Jaume I.(pdf) ↩
Deligiorgi, Despina, 及Kostas Philippopoulos. 「Spatial Interpolation Methodologies in Urban Air Pollution Modeling: Application for the Greater Area of Metropolitan Athens, Greece」, 2011. https://doi.org/10.5772/17734. ↩
Chen, Ling-Jyh, Yao Ho, Hu-Cheng Lee, Hsuan-Cho Wu, Hao Min Liu, Hsin-Hung Hsieh, Yu-Te Huang and Shih-Chun Candice Lung. 「An Open Framework for Participatory PM2.5 Monitoring in Smart Cities」. IEEE Access PP (2017年7月6日): 1–1. .。 ↩
Liu, Xiaohong, Ying Zhu, Weili Wang及Fengmin Liu. 「3D GIS modeling of air pollution effects」. 收入 2010 3rd International Congress on Image and Signal Processing, 6:2714–17, 2010. https://doi.org/10.1109/CISP.2010.5647463. ↩
Markus Konrad markus.konrad@wzb.eu / post@mkonrad.net, March 2022, geovoronoi – a package to create and plot Voronoi regions inside geographic areas. ↩