Shapely の便利な小技集
この記事は toggle holding の アドベントカレンダー 2024 (URL は分かり次第更新) の 1 日目の記事です。
概要
最近、Python の Shapely を使う機会が多いのですが、複雑な計算に少しテクニックが必要だったり、日本語での記事が少なかったりと、ちょっとしたハマりどころが多かったので、それらをまとめてみました。
そんなに複雑なコードは出てこないので、Shapely の知識がなくても Python が分かれば読める内容になっています。
Shapely の基本的な使い方
Shapely は Python のライブラリで、このライブラリを用いることで幾何学的なオブジェクトを簡単に扱えるようになります。
たとえば、ポリゴンや点の間の距離を計算したり、ポリゴン同士の共通部分や和集合を簡単に計算することができ、GIS データのポリゴンの処理などで活用できます。
具体的な計算例を見てみましょう。
たとえば、以下のコードは、
- 正方形のポリゴンを定義
- そのポリゴンを平行移動
- 元のポリゴンとの和集合を計算
- 描画
を行うコードです。
from shapely.geometry import Polygon
from shapely.affinity import translate
from shapely.plotting import plot_polygon
import matplotlib.pyplot as plt
# 正方形の polygon を定義
polygon = Polygon(((0, 0), (0, 10000), (10000, 10000), (10000, 0)))
# polygon を平行移動
translated = translate(polygon, xoff=5000, yoff=5000)
# 2つのポリゴンの和集合
union = polygon.union(translated)
# 描画
plot_polygon(polygon)
plot_polygon(translated, color="red")
plt.title("polygon and translated")
plt.show()
plot_polygon(union)
plt.title("union")
plt.show()
簡単なコードで集合の論理計算をやってくれて非常に便利です！
小技集
以下では、実際に業務で使っていて便利だった小技をいくつか紹介します。
また、以下のコードでは、描画部分は本質的ではないため省略しています。
Polygon だけを抜き出す
オブジェクトの intersection を計算したりすると、結果が GeometryCollection になってしまうことがあります。
GeometryCollection は Polygon や LineString などの複数のオブジェクトをまとめて持つオブジェクトです。
たとえば、以下のような2つのポリゴンの intersection (緑色で示した部分) を計算するとGeometryCollection が返されます。
polygon1 = Polygon(((0, 0), (0, 100), (100, 100), (100, 0)))
polygon2 = Polygon(((-100, 0), (-100, 100), (20, 100), (20, 80), (-20, 80), (-20, 20), (0, 20), (0, 0)))
intersection = polygon1.intersection(polygon2)
しかし、このような GeometryCollection から Polygon だけを抜き出して使いたい場合があります。
そのような場合は、以下のようにして Polygon だけを抜き出すことができます。
from shapely.geometry import Polygon, MultiPolygon, GeometryCollection
from shapely.geometry.base import BaseGeometry
from itertools import chain
def pickup_polygons(geom: BaseGeometry) -> tuple[Polygon, ...]:
"""
geom の中から polygon を抜き出して tuple にして返す関数
Args:
geom (BaseGeometry): 対象の shapely のインスタンス
Returns:
tuple[Polygon, ...]: 抜き出した polygon の tuple
"""
if isinstance(geom, Polygon) and not geom.is_empty:
return (geom, )
if isinstance(geom, MultiPolygon | GeometryCollection):
return tuple(chain.from_iterable(
map(pickup_polygons, geom.geoms)
))
return ()
polygon1 = Polygon(((0, 0), (0, 100), (100, 100), (100, 0)))
polygon2 = Polygon(((-100, 0), (-100, 100), (20, 100), (20, 80), (-20, 80), (-20, 20), (0, 20), (0, 0)))
intersection = polygon1.intersection(polygon2)
# Polygon だけを抜き出す
polygons = pickup_polygons(intersection)
自己交差の解消
複雑なポリゴンの計算をしたり、現実のデータを扱うと、ポリゴンが自己交差していることがあります。
たとえば、以下のようなポリゴンがあったとします。
polygon = Polygon(((0, 0), (0, 100), (100, 100), (100, 0), (0, 0)))
ここから2つのポリゴンを抜き出すには、以下のような関数を適用すれば良いです。
ポイントは、ポリゴンの外周だけを抜き出して、それを buffer で膨らませてから縮めることです。
さらに、その際に join_style を mitre にすることで、角が丸まるのを防ぎます。
def remove_self_intersection(
polygon: Polygon,
size: float = 0.1
) -> tuple[Polygon, ...]:
"""
ポリゴンの自己交差を解消する関数
Args:
polygon (Polygon): 自己交差を解消するポリゴン
size (float, optional): ポリゴンを拡大縮小するサイズ. Defaults to 0.1.
Returns:
tuple[Polygon, ...]: 自己交差を解消したポリゴンの tuple
"""
ring = polygon.exterior
buffered = ring.buffer(size, join_style=JOIN_STYLE.mitre)
outer_polygon = Polygon(buffered.exterior)
shirinked = outer_polygon.buffer(-size, join_style=JOIN_STYLE.mitre)
polygons = pickup_polygons(shirinked)
return polygons
polygon = Polygon(((0, 0), (100, 100), (0, 100), (100, 0)))
polygons = remove_self_intersection(polygon)
print(polygons) # (<POLYGON ((0 0, 50 50, 100 0, 0 0))>, <POLYGON ((50 50, 0 100, 100 100, 50 50))>)
オブジェクトの精度設定
小技ではないですが、日本語で書かれた記事がなかったので書いておきます。
Shapely のオブジェクトの演算では、たまに意図しない結果が返されることがあります。
(まれに発生するんですが、再現できなかったのでコードはないです。)
すでに closed ですがこのような issue が Shapely のリポジトリにありました。
このような問題は、grid_size を設定することで回避できます。
オブジェクトに対して grid_size を設定すると、各頂点が grid_size の精度で丸められます。
from shapely import set_precision
from shapely.geometry import Point
point = Point(0.05, 0.049)
precision_point = set_precision(point, 0.1)
print(precision_point) # POINT (0.1 0)
どのくらいの精度で丸めるかはコンテキストによりますが、可能ならば精度を設定しておくと想定外のエラーを防げそうです。
一方で、Shapely を使うたびに精度を設定するのは面倒なので、以下のような精度つき Shapely オブジェクトの生成クラスを作成しておくと便利です。
from shapely.geometry import (
Point as ShapelyPoint,
MultiPoint as ShapelyMultiPoint,
LineString as ShapelyLineString,
MultiLineString as ShapelyMultiLineString,
LinearRing as ShapelyLinearRing,
Polygon as ShapelyPolygon,
MultiPolygon as ShapelyMultiPolygon,
GeometryCollection as ShapelyGeometryCollection,
)
from shapely import set_precision
type PointData = tuple[float, float]
type LineData = tuple[PointData, PointData]
type PolygonData = tuple[PointData, ...]
class PShapely:
"""
Grid_size を自動設定した Shapely のオブジェクトを生成するクラス
"""
_GRID_SIZE = 0.1
@staticmethod
def Point(point: PointData | ShapelyPoint) -> ShapelyPoint:
return set_precision(
ShapelyPoint(point),
grid_size=PShapely._GRID_SIZE,
)
@staticmethod
def LineString(
points: Iterable[PointData | ShapelyPoint]
) -> ShapelyLineString:
return set_precision(
ShapelyLineString(points),
grid_size=PShapely._GRID_SIZE,
)
@staticmethod
def Polygon(
polygon: PolygonData | ShapelyPolygon
) -> ShapelyPolygon:
return set_precision(
ShapelyPolygon(polygon),
grid_size=PShapely._GRID_SIZE,
)
# 他にも MultiPoint, MultiLineString, MultiPolygon なども同様なので省略
polygon = PShapely.Polygon(((0.001, 0.001), (0, 100), (100, 100), (100, 0)))
print(polygon) # POLYGON ((0 100, 100 100, 100 0, 0 0, 0 100))
精度を設定する際には、ポリゴンに自己交差があるとエラーが発生するので注意が必要です。
polygon = PShapely.Polygon(((0.001, 0.001), (100, 100), (0, 100), (100, 0)))
# shapely.errors.GEOSException: TopologyException: side location conflict at 50 50. This can occur if the input geometry is invalid.
最後に
今後もいい感じの使い方見つけたら更新していきますー！
明日は、toggle の CTO 久森さんの記事(URLは分かり次第更新)です！
Discussion