shapely

Shapely - 平面幾何

ShapelyはGeoPandasやその他多くのGISツールを支えるエンジンです。ジオメトリそのものに特化しており、交差、和集合、距離の計算や、空間的な関係（「この点はこのポリゴンの内側にあるか？」など）の判定を行います。

使用する場面

• 2Dの幾何学的な形状を正確に操作する。
• 集合論的演算（交差、和集合、差分）を行う。
• 空間述語（包含、内包、交差、接触）を判定する。
• 「汚れた」ジオメトリを修正し、検証する（自己交差の修正など）。
• 幾何学的なプロパティ（面積、長さ、重心、境界ボックス）を計算する。
• バッファーを生成したり、複雑な線を単純化したりする。
• 線形参照（線上の位置を見つける）。

参考ドキュメント

公式ドキュメント：https://shapely.readthedocs.io/
GEOS（エンジン）：https://libgeos.org/
検索パターン：shapely.geometry、shapely.ops.unary_union、shapely.validation.make_valid

基本原則

幾何オブジェクト

オブジェクトはイミュータブルです。一度作成したオブジェクトは変更せず、演算を行って新しいオブジェクトを返します。

• ポイント：0次元です。
• ラインストリング：1次元の曲線です。
• ポリゴン：任意の穴を含む2次元の面です。

デカルト幾何

Shapelyはデカルト平面上で動作します。地球の曲率、緯度、経度については認識しません。距離はsqrt(dx² + dy²)です。

ベクトル化（Shapely 2.0+）

最新のShapelyでは、ジオメトリオブジェクトのNumPy配列に対するベクトル化演算がサポートされており、旧バージョンより大幅に高速です。

クイックリファレンス

インストール

pip install shapely numpy

標準的なインポート

import numpy as np
from shapely import Point, LineString, Polygon, MultiPoint, MultiPolygon
from shapely import ops, wkt, wkb
import shapely

基本パターン - 作成と分析

from shapely.geometry import Point, Polygon

# 1. Create objects
p = Point(0, 0)
poly = Polygon([(0, 0), (2, 0), (2, 2), (0, 2)])

# 2. Check relationships
is_inside = p.within(poly) # True
is_on_border = p.touches(poly) # False (interior counts as within)

# 3. Calculate
print(f"Area: {poly.area}")
print(f"Distance: {p.distance(Point(10, 10))}")

重要なルール

✅ 推奨事項

• 妥当性を確認する - 複雑な演算を行う前に.is_validを使用してください。無効なジオメトリ（自己交差するポリゴンなど）はエラーの原因になります。
• unary_unionを使用する - 多数のポリゴンを結合する場合、ops.unary_union([list])はp1.union(p2)をループで実行するよりも桁違いに高速です。
• ベクトル化関数を優先する - パフォーマンス向上のため、ループの代わりにshapely.intersects(array_a, array_b)を使用してください。
• prepare()を使用する - 同じポリゴンに対して多数の点を判定する場合、後続のクエリを高速化するためにshapely.prepare(poly)を使用してください。
• 分析用に単純化する - 高い精度が不要であれば、複雑な境界に.simplify(tolerance)を使用してパフォーマンスを向上させてください。

❌ 非推奨事項

• デカルト座標と球面座標を混在させる - 緯度・経度の点で距離を計算しないでください。結果は意味のない「度」になります。
• ポリゴンの向きを決めつける - Shapelyが処理しますが、外周リングは反時計回り、穴は時計回りが理想的であることを覚えておいてください。
• 3D計算に使用する - ShapelyはZ座標を格納できますが、面積や交差などの多くの演算ではZを無視し、XY平面に投影します。
• 大規模なコレクションをループ処理する - Shapely 2.0+で提供されるNumPy形式のベクトル化を使用してください。

アンチパターン（絶対に避ける）

from shapely.geometry import Point, Polygon
from shapely.ops import unary_union

# ❌ BAD: Merging geometries in a loop (O(n²) complexity)
result = geometries[0]
for g in geometries[1:]:
    result = result.union(g)

# ✅ GOOD: Use unary_union (O(n log n) complexity)
result = unary_union(geometries)

# ❌ BAD: Checking many points without preparing
for p in many_points:
    if complex_poly.contains(p): # Slow for complex shapes
        pass

# ✅ GOOD: Prepare the geometry (Builds a spatial index)
from shapely import prepare
prepare(complex_poly)
for p in many_points:
    if complex_poly.contains(p): # Much faster
        pass

ジオメトリの種類と作成

標準的な基本図形

# Point (x, y, z)
pt = Point(1.0, 2.0)

# LineString (Ordered sequence of points)
line = LineString([(0, 0), (1, 1), (2, 0)])

# Polygon (Shell, [Holes])
shell = [(0, 0), (10, 0), (10, 10), (0, 10)]
hole = [(2, 2), (2, 4), (4, 4), (4, 2)]
poly = Polygon(shell, [hole])

# Multi-Geometries (Collections)
points = MultiPoint([(0,0), (1,1)])

空間述語（関係）

オブジェクト同士の関係を確認する

a = Point(1, 1).buffer(1.5) # A circle
b = Polygon([(0,0), (2,0), (2,2), (0,2)]) # A square

print(a.intersects(b))  # Shared space?
print(a.contains(b))    # B entirely inside A?
print(a.disjoint(b))    # No shared space?
print(a.overlaps(b))    # Same dimension, shared space, but not within?
print(a.touches(b))     # Only boundaries share space?
print(a.crosses(b))     # Line crossing a polygon?

集合論的演算

既存のジオメトリから新しいジオメトリを作成する

poly1 = Point(0, 0).buffer(1)
poly2 = Point(1, 0).buffer(1)

# Intersection (Shared area)
inter = poly1.intersection(poly2)

# Union (Combined area)
union = poly1.union(poly2)

# Difference (Area in poly1 NOT in poly2)
diff = poly1.difference(poly2)

# Symmetric Difference (Area in either but NOT both)
sdiff = poly1.symmetric_difference(poly2)

構成的メソッド

バッファー処理、分割、単純化

# Buffer: Expand/shrink geometry
# cap_style: 1=Round, 2=Flat, 3=Square
line_thick = line.buffer(0.5, cap_style=2)

# Centroid: Geometric center
center = poly.centroid

# Representative Point: Guaranteed to be INSIDE the geometry
# Useful for label placement in U-shaped polygons
label_pt = poly.representative_point()

# Simplify: Reduce number of vertices
simple_line = complex_line.simplify(tolerance=0.1, preserve_topology=True)

# Convex Hull: Smallest convex box containing all points
hull = MultiPoint(points).convex_hull

線形参照

LineString上の位置を扱う

line = LineString([(0, 0), (0, 10), (10, 10)])

# Find distance along line to the point nearest to (5, 5)
dist = line.project(Point(5, 5)) # returns 5.0 (it's at (0, 5))

# Find the actual point at a specific distance along the line
pt = line.interpolate(15.0) # returns Point(5, 10)

I/O：WKT、WKB、NumPy

シリアライズとデータ交換

# WKT (Well-Known Text) - Human readable
text = "POINT (10 20)"
p = wkt.loads(text)
print(p.wkt)

# WKB (Well-Known Binary) - Fast and compact
binary = wkb.dumps(p)
p_new = wkb.loads(binary)

# NumPy Integration (Shapely 2.0)
points_array = np.array([Point(0,0), Point(1,1), Point(2,2)])
areas = shapely.area(points_array) # Returns array of zeros
dist_matrix = shapely.distance(points_array[:, np.newaxis], points_array)

実践的なワークフロー

1. 無効なジオメトリを修正する

from shapely.validation import make_valid

def safe_area(geom):
    """Calculates area even for invalid/self-intersecting polygons."""
    if not geom.is_valid:
        geom = make_valid(geom)
    
    # After make_valid, a Polygon might become a MultiPolygon or GeometryCollection
    return geom.area

2. ポイント・イン・ポリゴン検索（最適化）

from shapely import prepare

def find_points_in_poly(points, poly):
    """Efficiently filters points inside a complex polygon."""
    prepare(poly) # Builds internal STRtree or spatial index
    
    # Using vectorized intersection (much faster)
    mask = shapely.contains(poly, points)
    return points[mask]

3. 線でポリゴンを分割する

from shapely.ops import split

def divide_land(polygon, line):
    """Splits a polygon into multiple parts using a LineString."""
    result = split(polygon, line)
    # Returns a GeometryCollection of the resulting parts
    return list(result.geoms)

パフォーマンスの最適化

最近傍検索のためのSTRtree

数千のジオメトリがあり、ある点の近くにあるものを見つける必要がある場合は、STRtreeを使用してください。

from shapely import STRtree

tree = STRtree(geometries)

# Find indices of geometries whose bounding boxes intersect the point's buffer
indices = tree.query(Point(0,0).buffer(10))

# Find the single nearest geometry index
nearest_idx = tree.nearest(Point(0,0))

よくある問題と解決策

「スリバーポリゴン」の問題

交差などの計算では、浮動小数点誤差により、非常に小さくほとんど見えないポリゴンが生成されることがあります。

# ✅ Solution: Filter by area
intersection = p1.intersection(p2)
if intersection.area < 1e-9:
    intersection = None

緯度と経度の混同

ポイントは（x, y）です。GIS では、通常（経度、緯度）を意味します。

# ❌ Error: Point(Latitude, Longitude)
# This will plot your maps sideways!
# ✅ Solution: Always use (Lon, Lat) to match (X, Y)
nyc = Point(-74.006, 40.7128)

GeometryCollections

分割や交差などの操作では、GeometryCollection が返されることがあります。これは、異なる型を格納するコンテナーです。

# ❌ Problem: Calling .area on a collection with Lines and Polygons
# ✅ Solution: Filter for the type you want
polys = [g for g in collection.geoms if g.geom_type == 'Polygon']

ベストプラクティス

1. 複雑な操作の前に、必ず .is_valid でジオメトリを検証してください
2. 和集合をループで処理する代わりに、unary_union を使用してください
3. 同じ形状に対して多数のポイントをチェックする場合は、prepare() でジオメトリを準備してください
4. パフォーマンス向上のため、Shapely 2.0+ ではベクトル化操作を使用してください
5. Shapely はデカルト座標系を使用します。距離の計算に緯度・経度を直接使用しないでください
6. 幾何操作の後、面積のしきい値を使ってスリバーポリゴンを除外してください
7. 多数のジオメトリを扱う場合は、空間インデックスに STRtree を使用してください
8. 高精度が必要でない場合は、複雑なジオメトリを単純化してください
9. GeometryCollection は適切に処理してください。必要に応じて、ジオメトリの型でフィルタリングしてください
10. 複雑なポリゴン内に確実に収まるポイントには、representative_point() を使用してください

Shapely は、特化した鋭利なツールです。座標系やファイル形式には関知せず、形状間の純粋な数学的関係だけを扱います。これを習得することが、高度な空間アルゴリズムを構築する鍵となります。