鐵人賽2019 Day23 資料品質:幾何檢查(續篇)

  1. 1. gap and sliver polygon
  2. 2. 檢查方法
  3. 3. 修正
  4. 4. 後記

今天繼續GIS的幾何檢查的部分,主要處理的是polygon的接邊問題,這是一種常見的topology的問題。

大綱:

  • gap and sliver polygon
  • 檢查方法
  • 修正

gap and sliver polygon

兩個相接的polygon可能因為繪製的錯誤,或是浮點數的位數問題,造成邊界些微的不吻合(shp檔沒有儲存topology),可能有兩種情況,第一種是中間有縫隙(gap),第二種則是中間有一些重疊(sliver)。

https://ithelp.ithome.com.tw/upload/images/20181107/20107816SzAZQoYlpv.png
https://ithelp.ithome.com.tw/upload/images/20181107/20107816Fc64O9SC1P.png

檢查方法

為了測試,我自行畫了一個範例:
https://ithelp.ithome.com.tw/upload/images/20181107/20107816euXaTBcYPC.png

總共有三個Polygon,紅色部分highlight起來凸顯邊界問題,另外一個Polygon則是分開的

要檢查是否有因邊界不吻合所造成gap或sliver,
最簡單的方法,可以將polygon組為multyPolygon,利用shapely的is_valid檢查:

A valid MultiPolygon may not collect any overlapping polygons. — Shapely

用Geopandas讀自己畫的Polygon

1
2
3
import geopandas as gpd
from shapely.geometry import MultiPolygon, Polygon
test=gpd.read_file('data/test.shp',encoding='utf-8')

https://ithelp.ithome.com.tw/upload/images/20181107/20107816a5xCXg1QjZ.png

將polygon丟進去multypolygon

1
2
3
4
5
6
7
from shapely.geometry import MultiPolygon, Polygon
data=[]
for i in range(len(test.geometry)):
data.append(test['geometry'][i])
mp = MultiPolygon(data)
mp_obj = gpd.GeoSeries(mp)
mp_obj.is_valid

回傳:False
由於很明顯有不一致的現象,所以回傳是False

修正

確定這分資料有邊界的問題後,接下來是修正
昨天提到,修正最好還是在GIS軟體中編修介面中操作,
但我們可以試著將要修正的區塊找出來

首先可以使用unary_union將聯集的部分合併

1
2
union = test.geometry.unary_union
union

https://ithelp.ithome.com.tw/upload/images/20181107/20107816nbMLpVpqm6.png

結果由三個polygon被合併為兩個
1號polygon是被合併的部分,我們試著找出邊界不一致的區塊
可以例用differenceintersection的運算找出gap或sliver

找gap:

1
2
3
# gap
gap = list(union)[1] #
gap=gap.difference(test['geometry'][2]).difference(test['geometry'][1])

https://ithelp.ithome.com.tw/upload/images/20181107/20107816orUuH7zKCb.png

找sliver

1
2
# sliver
sliver=test['geometry'][2].intersection(test['geometry'][1])

https://ithelp.ithome.com.tw/upload/images/20181107/20107816cWR4bVRJxM.png

把成果疊一起示意:

1
2
3
4
5
ax=mp_obj.plot()
gap_obj = gpd.GeoSeries(gap)
sliver_obj = gpd.GeoSeries(sliver)
ax=gap_obj.plot(ax=ax,color='yellow')
sliver_obj.plot(ax=ax,color='red')

https://ithelp.ithome.com.tw/upload/images/20181107/20107816BIYcx2Ldzr.png

找出gap及sliver的區塊以後,建議回到QGIS加以編修,利用正確的數化工具將錯誤修正。
個人認為找出錯誤的部分比自動修正還要重要,
因為自動修正時編號1號的區塊很難程式化的決定要採取左邊或右邊的邊界,
有時候會遇到硬條件或是其他因素,自動化會比較瑣碎
(工人智慧?XD)

後記

幾何資料的正確性對資料品質影響很大,錯誤的資料在幾何運算也會有問題
今天的案例是一個測試,跟幾何有關的情況還有一些,有機會再來研究。

參考資料
https://gis.stackexchange.com/questions/277334/shapely-polygon-union-results-in-strange-artifacts-of-tiny-non-overlapping-area
https://www.gislounge.com/digitizing-errors-in-gis/