shapely

Shapely - 平面几何

Shapely 是 GeoPandas 和许多其他 GIS 工具的核心引擎。它专注于几何本身：计算交集、并集、距离，以及检查空间关系（例如“这个点是否在这个多边形内？”）。

适用场景

• 精确操作二维几何形状。
• 执行集合论运算（交集、并集、差集）。
• 检查空间谓词（包含、位于内部、相交、接触）。
• 清理和验证“脏”几何（修复自相交）。
• 计算几何属性（面积、长度、质心、边界）。
• 生成缓冲区或简化复杂线条。
• 线性参考（查找沿线的点）。

参考文档

官方文档：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 可以处理，但请记住，外环理想情况下应为逆时针方向，孔洞应为顺时针方向。
• 用于三维数学计算 - 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. 妥善处理 GeometryCollections——需要时按几何类型筛选
10. 对于复杂多边形，使用 representative_point() 获取保证位于内部的点

Shapely 是一种专用而锋利的工具。它不关心你的坐标系或文件格式——只关心形状之间纯粹的数学关系。掌握它是构建高级空间算法的关键。