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

大綱:

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

gap and sliver polygon

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

檢查方法

為了測試，我自行畫了一個範例：

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

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

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

用Geopandas讀自己畫的Polygon

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

將polygon丟進去multypolygon

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將聯集的部分合併

union = test.geometry.unary_union
union

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

找gap:

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

找sliver

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

把成果疊一起示意：

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')

找出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/