經緯座標轉行政區

最近接到了一個工作，需要將經緯座標轉換為行政區，以往的做法可能會去(偷)用 google 的地址轉換來取得，不過這次因為是關在封閉環境，加上是 ETL 工作的一部份，一筆一筆發 request 效能似乎也不是太好，所以試試新的做法。

由於目前在 opendata 中已經有台灣各行政區的地理資訊的 shapefile ，其中包成了由經緯座標組成的的多邊形(polygon)資訊，所以要取得某個點位所屬的行政區，只要去做 pip (point-in-polygon) 測試即可，而這個測試很多函式庫都有，以下是一個混合使用python 函式庫 pyshp/fiona 的範例，其中 VILLAGE_MOI_1070119.shp 是由 https://data.gov.tw/dataset/7438 取得，另外要稍微注意一下的可能只有有的行政區是由多個多邊形組成的，記得在看到 type 是 “MultiPolygon” 時逐一測試就好了。

geo.py

import shapefile
import uniout
from shapely.geometry import shape,mapping, Point, Polygon, MultiPolygon
import fiona

vills = fiona.open("VILLAGE_MOI_1070119.shp")
test_pt = Point(121.289761, 24.920151)

for v in vills:
        if v['geometry']['type'] == 'Polygon':
                poly = Polygon( v['geometry']['coordinates'][0] )
                if poly.contains ( test_pt ):
                        print v['properties']['TOWNNAME'], v['properties']['VILLNAME']
        elif v['geometry']['type'] == 'MultiPolygon':
                print v['properties']['TOWNNAME'], v['properties']['VILLNAME']
                for shape in v['geometry']['coordinates'][0]:
                        poly = Polygon( shape )
                        if poly.contains (test_pt):
                                print v['properties']['TOWNNAME'], v['properties']['VILLNAME']
                                break

輸出結果 shell> python geo.py

大溪區 仁善里

p.s 在 github (https://github.com/g0v/twgeojson)上有人另外以 geojson 的格式整理了同樣的資料