🌐
geom 0 is INVALID: Self-intersection 自己交差エラーとその回避
ポリゴンの扱いをしている中で以下のようなエラーが発生します。
GEOSException: TopologyException:
Input geom 0 is invalid:
Self-intersection at 182.71451535122969 3137.9381128435944
エラーの原因は Self-intersection の通りでポリゴン自体が 自己交差 してしまっていることが原因です。
ArcGIS の資料には 自己交差 の例が記載されています。
回避方法
自動処理の方が楽ですが対応できない複雑なポリゴンも存在するので
自動 を試してダメなら 手動 修正が順当です。
1. 自動処理
1. バッファ近所処理
バッファをある程度見て丸め込みます。一番簡単ですが、変形か異なるエラーが多いです。
polygon = polygon.buffer(0)
をするだけに修正できます。
バッファ関数
2. ポリゴン分割
交差点がある場合はそこを境に分割します。
intersections = polygon.buffer(0).intersection(polygon)
non_intersections = polygon.difference(intersections)
fixed_polygon = non_intersections.difference(intersections)
交差点の抽出関数
3. ポリゴンを合成する
交差領域を合わせて１つのポリゴンに結合します。
gdf = gpd.GeoSeries([polygon1, polygon2])
unioned_geometry = gdf.unary_union
gdf_unioned = gpd.GeoDataFrame(geometry=[unioned_geometry])
両方の領域の結合関数
2. 手動処理
1. QGIS の ジオメトリ修復機能
良しなに修正してくれます。しかし六芒星のような内外が定義できないものは修正できません。
2. 交差点を見つけて手動修正
諦めて QGIS を開いてポリゴンを手動で編集しましょう。
嬉しいことに QGIS には 有効性チェックと言う自動で交差点を点ポリゴンで検出してくれる機能があります。
すると、有効なジオメトリ と 不正なジオメトリ と エラー部分 をレイヤーで追加して確認することが可能です。
エラー部分の確認をしましょう。
エラー対処とその可視化
それぞれのエラー回避がどのようになっているか可視化します。
実際に無効なポリゴンを作成します。
from shapely.geometry import Polygon
import matplotlib.pyplot as plt
from matplotlib.patches import Polygon as MatplotlibPolygon
# ポリゴンを作成
polygon = Polygon([(0, 0), (1, 1), (1, 0), (0, 1), (0, 0)])
# Matplotlib の Figure と Axes を作成
fig, ax = plt.subplots()
# Shapely のポリゴンを Matplotlib の Polygon に変換
mpl_polygon = MatplotlibPolygon(list(polygon.exterior.coords), edgecolor='blue', facecolor='none')
# Axes に Polygon を追加
ax.add_patch(mpl_polygon)
# 描画範囲の設定
ax.set_xlim(-0.5, 1.5)
ax.set_ylim(-0.5, 1.5)
ax.set_aspect('equal')
# グリッドを表示
ax.grid(False)
# グラフを表示
plt.show()
このポリゴンに対して処理をかけていきます。
1. バッファ近所処理
実際に無効な
from shapely.geometry import Polygon
# 無効なポリゴン（自己交差が含まれる）を作成
invalid_polygon = Polygon([(0, 0), (1, 1), (1, 0), (0, 1), (0, 0)])
invalid_polygon = invalid_polygon.buffer(0)
fig, ax = plt.subplots()
# Axes に Polygon を追加
mpl_polygon = MatplotlibPolygon(
list(invalid_polygon.exterior.coords),
edgecolor='blue',
facecolor='none'
)
ax.add_patch(mpl_polygon)
# 描画範囲の設定
ax.set_xlim(-0.5, 1.5)
ax.set_ylim(-0.5, 1.5)
ax.set_aspect('equal')
# グリッドを表示
ax.grid(False)
# グラフを表示
plt.show();
因みに、 buffer(0) で処理しましたが、 buffer(0.1) でのミンコフスキー和の丸め込みでは以下のようになります。
2. ポリゴン分割
from shapely.geometry import Polygon
def split_self_intersecting_polygon(polygon):
# 交差している箇所と交差していない箇所に分割
intersections = polygon.buffer(0).intersection(polygon)
non_intersections = polygon.difference(intersections)
# 交差している箇所を元のポリゴンから差し引くことで、解消されたポリゴンを得る
fixed_polygon = non_intersections.difference(intersections)
return fixed_polygon
# 無効なポリゴン（自己交差が含まれる）を作成
invalid_polygon = Polygon([(0, 0), (1, 1), (1, 0), (0, 1), (0, 0)])
invalid_polygon = invalid_polygon.buffer(0)
# 自己交差を解消してポリゴンを分割
fixed_and_split_polygon = split_self_intersecting_polygon(invalid_polygon)
print(fixed_and_split_polygon.is_valid)
fixed_and_split_polygon
片方のみの状態から差し引かれて消えています。（例が良くないですね）
3. ポリゴンを合成する
交差点でぶつかる場合の新しいポリゴンへの合成です。
import geopandas as gpd
from shapely.geometry import Polygon
import matplotlib.pyplot as plt
# 2つのポリゴンを作成
polygon1 = Polygon([(0, 0), (1, 0), (1, 1), (0, 1)])
polygon2 = Polygon([(1, 1), (2, 1), (2, 2), (1, 2)])
for pol, col in zip([polygon1, polygon2], ['blue', 'green']):
# Matplotlib の Figure と Axes を作成
fig, ax = plt.subplots()
# Shapely のポリゴンを Matplotlib の Polygon に変換
mpl_polygon = MatplotlibPolygon(list(pol.exterior.coords), edgecolor=col, facecolor='none')
# Axes に Polygon を追加
ax.add_patch(mpl_polygon)
# 描画範囲の設定
ax.set_xlim(-0.5, 2.5)
ax.set_ylim(-0.5, 2.5)
ax.set_aspect('equal')
# グリッドを表示
ax.grid(False)
# グラフを表示
plt.show();
# GeoSeriesを作成
gdf = gpd.GeoSeries([polygon1, polygon2])
# GeoSeriesを結合
unioned_geometry = gdf.unary_union
# MultiPolygonを含むGeoDataFrameを作成
gdf_unioned = gpd.GeoDataFrame(geometry=[unioned_geometry])
# crsを指定 (例: EPSG:4326 - WGS 84)
gdf_unioned.crs = "EPSG:4326"
# 可視化
fig, ax = plt.subplots(figsize=(6, 6))
# 結合したジオメトリをプロット
gdf_unioned.plot(ax=ax, facecolor='none', edgecolor='blue');
# 原点のポリゴンもプロット
gdf.plot(ax=ax, facecolor='none', edgecolor='red')
# グラフの設定
ax.set_aspect('equal')
plt.show()
2. 手動処理
1. QGIS の ジオメトリ修復機能
QGIS での自動修復によって順序変換をすることで交差点回避や分割処理をしてくれます。
では、大井町埠頭をポリゴンで囲ってみます。
俯瞰してみると大丈夫そうですが、
あるあるなんですが、マウスで近くに点を打ったつもりが微妙にずれていて エラー部分 が出現します。
ジオメトリ修正をしてもらいましょう。
形的には変化はありません。
では何が起こったのでしょうか？
geojson で確認します。
- 自己交差あり
{
"type": "FeatureCollection",
"name": "shinagawa_oihuto",
"crs": {
"type": "name",
"properties": {
"name": "urn:ogc:def:crs:OGC:1.3:CRS84"
}
},
"features": [
{
"type": "Feature",
"properties": {
"id": 0
},
"geometry": {
"type": "MultiPolygon",
"coordinates": [
[
[
[
139.754688153336019,
35.633517459592603
],
[
139.759442611643806,
35.633640133591882
],
[
139.759669014420382,
35.618795212890909
],
[
139.752801463531341,
35.619040605883036
],
[
139.754159880190684,
35.632352047209459
],
[
139.754686957223441,
35.633516148877163
],
[
139.75468840138663,
35.63351757855736
],
[
139.754688153336019,
35.633517459592603
]
]
]
]
}
}
]
}
- ジオメトリ修復後
{
"type": "FeatureCollection",
"name": "fix_giometry",
"crs": {
"type": "name",
"properties": {
"name": "urn:ogc:def:crs:OGC:1.3:CRS84"
}
},
"features": [
{
"type": "Feature",
"properties": {
"id": 0
},
"geometry": {
"type": "MultiPolygon",
"coordinates": [
[
[
[
139.759669014420382,
35.618795212890909
],
[
139.752801463531341,
35.619040605883036
],
[
139.754159880190684,
35.632352047209459
],
[
139.754686957223441,
35.633516148877163
],
[
139.754688284638917,
35.633517462980464
],
[
139.759442611643806,
35.633640133591882
],
[
139.759669014420382,
35.618795212890909
]
]
],
[
[
[
139.75468840138663,
35.63351757855736
],
[
139.754688284638917,
35.633517462980464
],
[
139.754688153336019,
35.633517459592603
],
[
139.75468840138663,
35.63351757855736
]
]
]
]
}
}
]
}
２つに分割されていますね。
図で見てみましょう。
親のオレンジ ポリゴンと分離した 子の青色 ポリゴンのマルチポリゴンに再構成されています。
2. 交差点を見つけて手動修正
geojson などで直接編集してもいいですが、、、、
これが交差していて、
{
"type": "FeatureCollection",
"name": "self_intersection",
"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:OGC:1.3:CRS84" } },
"features": [
{ "type": "Feature", "properties": { "id": 0 }, "geometry": { "type": "MultiPolygon", "coordinates": [ [ [ [ 139.667870107444486, 35.319946493090285 ], [ 139.780789166767903, 35.303126493037304 ], [ 139.896396775122838, 35.409837537404144 ], [ 139.80498610805148, 35.547038422170466 ], [ 139.817532670198517, 35.615551006139647 ], [ 140.005731102404212, 35.649785312537226 ], [ 140.070256279160418, 35.619193650518945 ], [ 140.058605900023878, 35.541204857982713 ], [ 139.954648670805511, 35.458760450465988 ], [ 139.764657872578823, 35.505465019389057 ], [ 139.69027468270707, 35.439778968847648 ], [ 139.667870107444486, 35.319946493090285 ] ] ] ] } }
]
}
これが交差していないです。
{
"type": "FeatureCollection",
"name": "no_intersection",
"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:OGC:1.3:CRS84" } },
"features": [
{ "type": "Feature", "properties": { "id": 0 }, "geometry": { "type": "MultiPolygon", "coordinates": [ [ [ [ 139.685818312939801, 35.437104012515036 ], [ 139.564524675119202, 35.50098007005613 ], [ 139.68308236622201, 35.461620591682895 ], [ 139.799816092846356, 35.51137369516529 ], [ 139.828087542263177, 35.602631268155847 ], [ 139.960324966954772, 35.672302756233229 ], [ 140.07705869357909, 35.59744053631848 ], [ 140.029635617137956, 35.51582769412245 ], [ 139.904694050360405, 35.467562879948488 ], [ 139.892838281250107, 35.369459005325893 ], [ 139.68308236622201, 35.299524367255387 ], [ 139.685818312939801, 35.437104012515036 ] ] ] ] } }
]
}
はい、無理ですね。
QGIS で可視化してみましょう。
東京湾を囲おうとしてたんですね（雑すぎです）
手動修正的には、交差点を消しましょう。
それよりも根本的には、自己交差するようなポリゴンを 作成しないように意識 することが大切です。
小技
ポリゴン検証
is_valid という関数で以下のようなエラーが出るかどうか事前にチェックすることが可能です。
GEOSException: TopologyException: Input geom 1 is invalid: Self-intersection at 0.5 0.5
from shapely.geometry import MultiPolygon
MultiPolygon([Point(0, 0).buffer(2.0), Point(1, 1).buffer(2.0)]).is_valid
# False
交差点抽出
make_valid という関数では以下のように交差点を検出してくれます。
from shapely.validation import make_valid
coords = [(0, 0), (0, 2), (1, 1), (2, 2), (2, 0), (1, 1), (0, 0)]
p = Polygon(coords)
make_valid(p)
# <MULTIPOLYGON (((1 1, 0 0, 0 2, 1 1)), ((2 0, 1 1, 2 2, 2 0)))>
謝辞
発端は私が Instance segmentation の処理時に助けを求めたら以下の方々のように助けて頂いたことです。大変感謝しております。ありがとうございます！
最後に
調べてみたら意外と解決方法や様々あり、それに検索した時に日本語資料が少ないので、
根本的、網羅的に勉強する必要があると思ったのでまとめました。
皆様のご参考になればと幸いです。
Discussion