目的
仕事(?)で使えればと昨年copilotと一緒に作成した。一人で勝手にやっているだけのため可読性・保守性は度外視でgit管理すらしていない。記事にするまで半年以上経過しておりまともに覚えていないため、最小限GPTにリファクタリングしてもらった上での雑な投稿である。
実現したいこと
GPTによると「線分データをトポロジー的に整理して閉じた面を生成し、それらを統合して最終的な外形を取得した後、指定した距離だけ外側へ拡張する」と表現するのが正しいらしい。
以下画像で示す、座標の点データで得た青線から赤線を求める感じ。シンプルそうだが、直接得られそうなものはなかった。
コード（雑）
import系は省略。また、当記事用に前後にあるIO関係処理含めて大きく省略している。plot系の関数は処理の途中で視覚的に確認できるようにするためのもののため、実際は不要。定数値は扱うデータのスケールや精度によって要調整。
SNAP_TOLERANCE = 0.01
SIMPLIFY_TOLERANCE = 0.01
THICKNESS = 5
MITRE_LIMIT = 99999
data = [
[-389, -388, -349, -388],
[-349, -437, -349, -388],
[-329, -291, -329, -281],
[-229, -291, -329, -291],
[-449, -487, -449, -437],
[-279, -487, -449, -487],
[-229, -437, -229, -291],
[-279, -487, -229, -437],
[-329, -281, -389, -281],
[-349, -437, -388, -437],
[-388, -437, -449, -437],
[-229, -291, -329, -281],
[-389, -437, -389, -388],
[-389, -388, -389, -281],
[-254, -353, -204, -353],
[-204, -403, -204, -353],
[-254, -403, -204, -403],
[-254, -403, -254, -353],
[-496, -292, -446, -292],
[-446, -342, -446, -292],
[-496, -342, -446, -342],
[-496, -342, -496, -292]
]
def plotLineString(lines):
fig, ax = plt.subplots()
for line in lines:
x, y = line.xy
ax.plot(x, y, color="blue")
ax.set_aspect("equal")
plt.show()
def plot(polygon):
x, y = polygon.exterior.xy
plt.figure()
plt.plot(x, y, label="Polygon")
plt.axis("equal")
plt.legend()
plt.show()
def plots(geometries):
plt.figure()
for geometry in geometries:
if isinstance(geometry, MultiPolygon):
polygons = geometry.geoms
elif isinstance(geometry, Polygon):
polygons = [geometry]
else:
polygons = geometry
for polygon in polygons:
x, y = polygon.exterior.xy
plt.plot(x, y, label="Polygon")
plt.axis("equal")
plt.legend()
plt.show()
def plot_vertices(geometry, title="Vertices"):
if geometry.is_empty:
print("The geometry is empty.")
return
if isinstance(geometry, Polygon):
for x, y in geometry.exterior.coords:
plt.plot(x, y, "o", label=f"Vertex ({x:.2f}, {y:.2f})")
elif hasattr(geometry, "geoms"):
for geom in geometry.geoms:
plot_vertices(geom, title)
plt.axis("equal")
plt.legend()
plt.title(title)
plt.show()
lines = [
LineString([
(coords[0], coords[1]),
(coords[2], coords[3]),
])
for coords in data
]
# プロット
# plotLineString(lines)
# 線同士の微小なズレを補正するため
line_geometry = MultiLineString(lines)
snapped_lines = snap(
line_geometry,
line_geometry,
tolerance=SNAP_TOLERANCE,
)
# 線同士の交差・接続関係を整理し、連続する線をまとめる
merged_lines = linemerge(
unary_union(snapped_lines)
)
# 閉じた線からPolygonを生成する
polygons = list(
polygonize(
unary_union(merged_lines)
)
)
# プロット
# plots(polygons)
# 複数の面を統合して外郭を作る
outer = unary_union(polygons)
# 統合後の形状を簡略化する
simplified = outer.simplify(
tolerance=SIMPLIFY_TOLERANCE,
preserve_topology=True,
).normalize()
# プロット
# plots([simplified])
# 外側へ指定距離だけ拡大する
appearance = simplified.buffer(
distance=THICKNESS,
join_style="mitre",
mitre_limit=MITRE_LIMIT,
)
# プロット
# plots([appearance])
ここで扱う座標データ
実在のデータではなく当記事用に適当に作成したもの。本来はテキストファイルを読み込んだものをdataに格納していたが、ここでは直接変数に入れている。配列の各要素(1行)が1本の線を表し、前2項目が線の始点の座標、後ろ2項目が線の終点の座標である。
※変にマイナス値を使用しているのは特に意味はない
data = [
[-389, -388, -349, -388],
...座標の点から線集合にする
座標データを格納している変数dataは、あくまで点の情報でしかない。これを後ほどShapelyで処理するために線データに変換する必要があるため、LineStringに変換。
lines = [
LineString([
(coords[0], coords[1]),
(coords[2], coords[3]),
])
for coords in data
]▼プロット結果
図のとおり、座標の始点と終点を結んだ線分の集合がプロットされる。この時点では点から得た線データのため、閉じた図形としては扱うことができない。
各線の端点同士をスナップ（必要な場合のみ）
取り込むデータの線の端点同士が完全に結ばれているなら、この処理は不要
取り込んだ線の端点同士が完全に結ばれておらず極僅かに離れているケースがあったため、この処理を加えた。後で面として扱う必要があるが、線によって完全に閉じられていなければ面にできないためだ。
snapメソッドの第3引数toleranceでは、どの程度までスナップするかを数値で指定するものであり、こちらも扱うデータのスケールによって調整のこと。
# 線同士の微小なズレを補正するため
line_geometry = MultiLineString(lines)
snapped_lines = snap(
line_geometry,
line_geometry,
tolerance=SNAP_TOLERANCE,
)なお、扱うデータのスケールに対してtoleranceの値が大きすぎると、以下のように想定していない線同士が結ばれて崩壊してしまう。（※画像は同データをtolerance=100で実行した場合）
線の交差・重複を整理してつなぐ
# 線同士の交差・接続関係を整理し、連続する線をまとめる
merged_lines = linemerge(
unary_union(snapped_lines)
)unary_union()はShapely公式ドキュメントによると“複数のジオメトリの和集合を返します”とのこと。
和集合は適当に画像検索して、ベン図で見るとわかりやすい。
今回の場合は重複部分をまとめ、交差する線を交点で分割してくれる。その結果をlinemerge()の引数で渡すと、端点でつながる線分を連続した線としてマージさせている。このとき、分岐点は、複数の線として残ることになる。
面（ポリゴン）データにする
# 閉じた線から Polygon を生成する
polygons = list(
polygonize(
unary_union(merged_lines)
)
)polygonize()は、閉じた線網から囲まれた領域をPolygonとして取得する。polygonize()の引数に再度unary_union()してる理由は忘れた。取り込むデータによってはないとダメだったかもしれないが、バージョン管理も何もしていないため、覚えていない。（不要かもしれないがそのままにしている）
これで得たpolygonsをプロットすると以下のようになる。単一の線だったものが、それぞれの閉じた領域で単一の図形としてみなされ、複数のポリゴン（図形）がプロットされている状態だ。一部の線は重複しているが、それぞれ閉じた領域として各ポリゴン（図形）の辺として機能している。
最外郭の多角形を得る
# 複数の面を統合して外郭を作る
outer = unary_union(polygons)
# 統合後の形状を簡略化する
simplified = outer.simplify(
tolerance=SIMPLIFY_TOLERANCE,
preserve_topology=True,
).normalize()複数のポリゴンををunary_union()で統合すると、その和集合を得ることでそれぞれが最外郭のポリゴンが取得できた。この処理(125行目)だけで済みそうなのだが、一部データで頂点のjoin_styleがbevel(※後述)になることがあったため128～131行目の処理をいれたようだ。（覚えていないが、当時のコード上のコメントブロックにFIXMEで記載していた）
得た最外郭（多角形）を外方向に拡大
前段で得た変数simplifiedにbuffer()を適用すると、境界から指定した数値分拡大できる。負の値なら狭めることもできるようだが、今回は広げることが目的のため正の値にしている。広げる大きさは引数distanceに渡し、拡大した際の頂点の処理方法はjoin_styleで指定する。
round：頂点を丸くつなぐ
デフォルト値
appearance = simplified.buffer(
distance=THICKNESS,
join_style="round",
)mitre：辺を延長して尖らせる
鋭角をとがらせる場合はこれ。引数が1つ増え、mitre_limitを指定している。
鋭角部では、元の外郭から頂点までの距離はdistanceよりも長くなるため、この限界値を指定するもの。これを超えた場合は後述するbevelのように面取りされた形状になる。
appearance = simplified.buffer(
distance=THICKNESS,
join_style="mitre",
mitre_limit=MITRE_LIMIT,
)bevel：頂点を斜めに落とす
bevelは辺の延長を頂点で止め、その間を直線で結ぶ
appearance = simplified.buffer(
distance=THICKNESS,
join_style="bevel",
)元コード（join_style=”mitre”)でプロットした結果、以下になる。少し分かりづらいが、前回のsimplifiedのプロット結果より拡大し、欲しい結果が得られた。
