2017年11月7日 星期二

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

幾所大學的校園範圍及地理中心 最近想要處理 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. MultiPolygon 的形心不一定落在內部, 例如世新大學 最後把 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 示範教學~~)

1 則留言:

  1. 1. 更正 gj-area-centroid.py 的數值問題: 把所有座標減去大約中心座標,以避免數值誤差 (numerical instability)

    2. 補上 「臺灣特大號超商旅遊地圖」。

    回覆刪除