多邊形求面積與地理中心/形心/幾何中心近似值

最近想要處理 osm 上面撈出來的便利商店資料。 其中有些以點標示， 有些以多邊形標示。 為了簡化問題， 想用多邊形的 地理中心 (點) 來取代多邊形 (面)。 以數學術語來說， 就是需要計算 幾何中心， 也就是 centroid， 亦稱為形心。 在 平面上的多邊形形心計算公式 裡面， 還需要順便計算面積。 所以我寫了 gj-area-centroid.py 來計算 geojson 簡單多邊形的面積與形心的近似值。

1. 假設你已有一個 geojson 檔， 例如 從 osm 裡面撈出來 的 university.geojson。
2. 首先用 json 轉檔瑞士刀 挑出當中所有的 Polygon 及 MultiPolygon： jq '.features | map(select(.geometry.type | contains("Polygon"))) | {type:"FeatureCollection",features:.}' university.geojson > univ.geojson 。
3. 然後用 gj-area-centroid.py 找出 centroid： ./gj-area-centroid.py univ.geojson > univ-area.geojson
4. 再度用 jq 挖出 centroid 的座標， 並且把原來的資料除了 geometry 的部分之外原封不動照抄過來： jq '.features | map(.geometry={type:"Point",coordinates:.geometry.centroid}) | {type:"FeatureCollection",features:.}' univ-area.geojson > univ-centroid.geojson 這就是原多邊形的形心， 以 Point 的形式呈現的 geojson 檔。
5. 最後把 univ.geojson 跟 univ-centroid.geojson 一起放進 potluckmap， 得到 臺灣部分大學校園範圍及地理中心圖。 注意： 形心不一定落在多邊形的內部。 凹多邊形的形心就有可能落在外部， 例如中原大學。 MultiPolygon 的形心也有可能落在外部， 例如世新大學。

我採用的算式是把球面當成平面來近似的。 如果你想要知道真實的面積， 那就要先用 pip install area 安裝 這個套件， 然後用 ./gj-area-centroid.py -e 100000 univ.geojson > univ-area.geojson 這樣 univ-area.geojson 裡面會多出 area_true 及 area_error 這兩個欄位， 前者是真實面積， 後者是誤差比例乘以 100000。 然後 jq '.features | map([.properties.name, .geometry.area, .geometry.area_error] | @csv)' univ-area.geojson | perl -pe 's/\\"//g; s/"//g' > error.txt 可以產生誤差比例統計檔 error.txt。 再用 sort -n -t , -k 3 error.txt 可以看到絕大多數的值都在正負 1 之間。 也就是說， 以臺灣的校園大小的面積 (幾平方公里) 來說， 採用平面公式的誤差小於十萬分之一。 更低的緯度或更小的面積的情況下， 誤差只會更小。

[11/8] 有了這個工具， 就可以找出全臺最大的一些超商。 假設你已從超商資料 (下 osmfilter 指令時用 --keep="shop=convenience" 選項) 撈出有面積的店面 (約290個)， 並產生內含面積和形心的 conv-area.geojson。 再來可用 jq '.features | map(select(.geometry.area>800)) | {type:"FeatureCollection",features:.}' conv-area.geojson > big-conv.geojson 撈出面積大於 800 平方公尺的超商。 再用 jq '.features | map(select(.geometry.area>800) | .properties.area=.geometry.area | .geometry={type:"Point",coordinates:.geometry.centroid}) | {type:"FeatureCollection",features:.}' conv-area.geojson > big-conv-centroid.geojson 把 area 欄位從 .geometry 搬到 .properties 去。 最後放進 potluckmap 得到： 臺灣的特大號超商旅遊地圖。

(這篇根本就是活用 jq 示範教學~~)