第5章：空间关系与谓词判断

空间关系判断是 GIS 分析的核心能力。本章将深入介绍 DE-9IM 空间关系模型，以及 Shapely 提供的所有空间谓词——contains、within、intersects、touches、crosses、overlaps、covers、equals 等，每个谓词都配有图形化示意和详细代码示例，并介绍如何使用 Prepared Geometry 加速重复判断。

5.1 DE-9IM 模型详解

5.1.1 Interior、Boundary、Exterior 概念

DE-9IM（Dimensionally Extended 9-Intersection Model）是描述两个几何体空间关系的数学模型。每个几何体被分为三个部分：

┌─────────────────────────────────────────────────────────────┐
│                  几何体的三个组成部分                          │
├─────────────────────────────────────────────────────────────┤
│                                                             │
│  ┌─── Interior（内部）───┐                                  │
│  │  几何体的内部区域      │                                  │
│  │  · 点没有内部          │                                  │
│  │  · 线的内部是端点之间  │                                  │
│  │  · 面的内部是边界围成  │                                  │
│  │    的区域              │                                  │
│  └────────────────────────┘                                 │
│                                                             │
│  ┌─── Boundary（边界）───┐                                  │
│  │  几何体的边界          │                                  │
│  │  · 点的边界为空        │                                  │
│  │  · 线的边界是两个端点  │                                  │
│  │  · 面的边界是外环和    │                                  │
│  │    内环               │                                  │
│  └────────────────────────┘                                 │
│                                                             │
│  ┌─── Exterior（外部）───┐                                  │
│  │  几何体以外的所有空间  │                                  │
│  │  · 无限延伸的平面      │                                  │
│  │    减去内部和边界      │                                  │
│  └────────────────────────┘                                 │
│                                                             │
└─────────────────────────────────────────────────────────────┘

对于不同维度的几何体：

几何类型	Interior（内部）	Boundary（边界）	Exterior（外部）
Point	点本身	空集	点以外所有空间
LineString	端点之间的部分	两个端点	线以外所有空间
Polygon	边界围成的区域	外环 + 内环	面以外所有空间

5.1.2 9 交叉矩阵

DE-9IM 通过计算两个几何体 A 和 B 的 Interior、Boundary、Exterior 两两交集的维度，构成一个 3×3 的矩阵：

         B.Interior  B.Boundary  B.Exterior
        ┌───────────┬───────────┬───────────┐
A.Int   │  dim(I∩I) │  dim(I∩B) │  dim(I∩E) │
        ├───────────┼───────────┼───────────┤
A.Bdy   │  dim(B∩I) │  dim(B∩B) │  dim(B∩E) │
        ├───────────┼───────────┼───────────┤
A.Ext   │  dim(E∩I) │  dim(E∩B) │  dim(E∩E) │
        └───────────┴───────────┴───────────┘

矩阵中的值：

T：非空交集（维度 ≥ 0）
F：空交集
0：0 维交集（点）
1：1 维交集（线）
2：2 维交集（面）
*：任意值（通配符）

5.1.3 使用 `relate()` 获取 DE-9IM 字符串

from shapely import Point, LineString, Polygon

# 两个重叠的多边形
a = Polygon([(0, 0), (4, 0), (4, 4), (0, 4)])
b = Polygon([(2, 2), (6, 2), (6, 6), (2, 6)])

# 获取 DE-9IM 字符串（9 个字符）
matrix = a.relate(b)
print(f"DE-9IM: {matrix}")  # 212101212

# 解读：
# 位置 0 (I∩I) = 2: 内部交集是面（重叠区域）
# 位置 1 (I∩B) = 1: A内部与B边界交集是线
# 位置 2 (I∩E) = 2: A内部与B外部交集是面
# 位置 3 (B∩I) = 1: A边界与B内部交集是线
# 位置 4 (B∩B) = 0: A边界与B边界交集是点
# 位置 5 (B∩E) = 1: A边界与B外部交集是线
# 位置 6 (E∩I) = 2: A外部与B内部交集是面
# 位置 7 (E∩B) = 1: A外部与B边界交集是线
# 位置 8 (E∩E) = 2: A外部与B外部交集是面

5.1.4 使用 `relate_pattern()` 模式匹配

from shapely import Point, Polygon

poly = Polygon([(0, 0), (4, 0), (4, 4), (0, 4)])
point = Point(2, 2)  # 在多边形内部

# contains 对应的 DE-9IM 模式：T*****FF*
print(f"relate: {poly.relate(point)}")  # 0F20F1FF2

# 用模式匹配验证 contains 关系
# 模式 T*****FF* 意味着：
# - I∩I 非空 (T)
# - E∩I 为空 (F)
# - E∩B 为空 (F)
print(f"匹配 contains: {poly.relate_pattern(point, 'T*****FF*')}")  # True

# within 对应的模式：T*F**F***
print(f"匹配 within: {point.relate_pattern(poly, 'T*F**F***')}")  # True

5.2 contains / within

5.2.1 `contains(other)` — 包含

判断几何体 A 是否完全包含几何体 B。要求 B 的所有点都在 A 的内部或边界上，且至少有一个点在 A 的内部。

contains 示意：

  A ┌────────────────┐
    │                │
    │    ┌──┐        │
    │    │B │        │   A.contains(B) = True
    │    └──┘        │   B 完全在 A 内部
    │                │
    └────────────────┘

  A ┌────────────────┐
    │                │
    │            ┌───┼───┐
    │            │   │ B │   A.contains(B) = False
    │            └───┼───┘   B 部分在 A 外部
    │                │
    └────────────────┘

from shapely import Point, Polygon

poly = Polygon([(0, 0), (10, 0), (10, 10), (0, 10)])

# 内部点
p_inside = Point(5, 5)
print(f"包含内部点: {poly.contains(p_inside)}")    # True

# 边界上的点
p_boundary = Point(0, 5)
print(f"包含边界点: {poly.contains(p_boundary)}")  # True

# 外部点
p_outside = Point(15, 5)
print(f"包含外部点: {poly.contains(p_outside)}")   # False

# 更小的多边形
small_poly = Polygon([(2, 2), (4, 2), (4, 4), (2, 4)])
print(f"包含小多边形: {poly.contains(small_poly)}")  # True

5.2.2 `within(other)` — 被包含

within 是 contains 的逆操作：A.within(B) 等价于 B.contains(A)。

from shapely import Point, Polygon

poly = Polygon([(0, 0), (10, 0), (10, 10), (0, 10)])
p = Point(5, 5)

# within 和 contains 互为逆操作
print(f"p.within(poly): {p.within(poly)}")       # True
print(f"poly.contains(p): {poly.contains(p)}")   # True

# 完全等价
assert p.within(poly) == poly.contains(p)

5.3 contains_properly

5.3.1 严格包含

contains_properly 要求 B 的所有点都在 A 的内部（不包括边界）。这比 contains 更严格。

contains_properly 与 contains 的区别：

  A ┌────────────────┐
    │                │
    │    · P         │   contains(P) = True
    │                │   contains_properly(P) = True
    └────────────────┘   P 在 A 内部

  A ┌────────────────┐
    │                │
    · P              │   contains(P) = True
    │                │   contains_properly(P) = False
    └────────────────┘   P 在 A 边界上

from shapely import Point, Polygon

poly = Polygon([(0, 0), (10, 0), (10, 10), (0, 10)])

# 内部点
p_inside = Point(5, 5)
print(f"contains: {poly.contains(p_inside)}")            # True
print(f"contains_properly: {poly.contains_properly(p_inside)}")  # True

# 边界上的点
p_boundary = Point(0, 5)
print(f"contains: {poly.contains(p_boundary)}")            # True
print(f"contains_properly: {poly.contains_properly(p_boundary)}")  # False

# 顶点
p_vertex = Point(0, 0)
print(f"contains: {poly.contains(p_vertex)}")              # True
print(f"contains_properly: {poly.contains_properly(p_vertex)}")    # False

5.4 covers / covered_by

5.4.1 `covers(other)` — 覆盖

covers 比 contains 更宽松。只要求 B 的所有点都在 A 的内部或边界上，没有”至少一个点在内部”的限制。

from shapely import Point, LineString, Polygon

poly = Polygon([(0, 0), (10, 0), (10, 10), (0, 10)])

# 边界上的线段
boundary_line = LineString([(0, 0), (10, 0)])

# contains vs covers
print(f"contains(boundary_line): {poly.contains(boundary_line)}")  # False
print(f"covers(boundary_line): {poly.covers(boundary_line)}")      # True
# contains 返回 False 因为 boundary_line 没有点在 poly 内部
# covers 返回 True 因为所有点都在 poly 的内部或边界上

# 内部点——两者都返回 True
p = Point(5, 5)
print(f"contains(p): {poly.contains(p)}")  # True
print(f"covers(p): {poly.covers(p)}")      # True

5.4.2 `covered_by(other)` — 被覆盖

covered_by 是 covers 的逆操作：

from shapely import Point, Polygon

poly = Polygon([(0, 0), (10, 0), (10, 10), (0, 10)])
p = Point(0, 0)  # 顶点

print(f"p.covered_by(poly): {p.covered_by(poly)}")  # True
print(f"poly.covers(p): {poly.covers(p)}")           # True

# 但 within/contains 对于边界点有区别
print(f"p.within(poly): {p.within(poly)}")        # True
print(f"poly.contains(p): {poly.contains(p)}")    # True

5.5 intersects / disjoint

5.5.1 `intersects(other)` — 相交

判断两个几何体是否有任何公共点。intersects 是最宽泛的空间关系——只要有任何接触就返回 True。

intersects 示意：

  ┌──────┐  ┌──────┐
  │  A   │  │  B   │   intersects = False (分离)
  └──────┘  └──────┘

  ┌──────┐
  │  A   ├──────┐
  │      │  B   │       intersects = True (重叠)
  └──────┤      │
         └──────┘

  ┌──────┬──────┐
  │  A   │  B   │       intersects = True (相邻)
  └──────┴──────┘

from shapely import Polygon, Point, LineString

poly1 = Polygon([(0, 0), (4, 0), (4, 4), (0, 4)])
poly2 = Polygon([(2, 2), (6, 2), (6, 6), (2, 6)])
poly3 = Polygon([(5, 5), (8, 5), (8, 8), (5, 8)])

# 重叠的多边形
print(f"poly1 ∩ poly2: {poly1.intersects(poly2)}")  # True

# 分离的多边形
print(f"poly1 ∩ poly3: {poly1.intersects(poly3)}")  # False

# 点与多边形
p = Point(2, 2)
print(f"point ∩ poly1: {p.intersects(poly1)}")  # True

# 线与多边形
line = LineString([(0, 5), (5, 0)])
print(f"line ∩ poly1: {line.intersects(poly1)}")  # True

5.5.2 `disjoint(other)` — 不相交

disjoint 是 intersects 的逻辑否定：

from shapely import Polygon

a = Polygon([(0, 0), (2, 0), (2, 2), (0, 2)])
b = Polygon([(3, 3), (5, 3), (5, 5), (3, 5)])
c = Polygon([(1, 1), (3, 1), (3, 3), (1, 3)])

# disjoint 与 intersects 完全互补
print(f"a.disjoint(b): {a.disjoint(b)}")       # True  (不相交)
print(f"a.intersects(b): {a.intersects(b)}")   # False

print(f"a.disjoint(c): {a.disjoint(c)}")       # False (相交)
print(f"a.intersects(c): {a.intersects(c)}")   # True

# 验证互补关系
assert a.disjoint(b) == (not a.intersects(b))
assert a.disjoint(c) == (not a.intersects(c))

5.6 touches

5.6.1 仅边界相接

touches 判断两个几何体是否仅在边界上接触，内部不相交。

touches 示意：

  ┌──────┐┌──────┐
  │  A   ││  B   │   touches = True  (共享边界)
  └──────┘└──────┘

  ┌──────┐
  │  A   ├───┐
  │      │ B │        touches = False (内部相交)
  └──────┤   │
         └───┘

  ┌──────┐
  │  A   │
  └──────┘
           · P         touches = False (不接触)

from shapely import Polygon, Point, LineString

# 共享边的两个矩形
a = Polygon([(0, 0), (4, 0), (4, 4), (0, 4)])
b = Polygon([(4, 0), (8, 0), (8, 4), (4, 4)])

print(f"共享边 touches: {a.touches(b)}")  # True
print(f"共享边 intersects: {a.intersects(b)}")  # True

# 仅共享顶点
c = Polygon([(4, 4), (8, 4), (8, 8), (4, 8)])
print(f"共享顶点 touches: {a.touches(c)}")  # True

# 重叠的多边形——不是 touches
d = Polygon([(2, 2), (6, 2), (6, 6), (2, 6)])
print(f"重叠 touches: {a.touches(d)}")  # False

# 点在多边形边界上
p_on_boundary = Point(0, 2)
print(f"边界点 touches: {a.touches(p_on_boundary)}")  # True

# 点在多边形内部
p_inside = Point(2, 2)
print(f"内部点 touches: {a.touches(p_inside)}")  # False

5.7 crosses

5.7.1 交叉关系

crosses 判断两个几何体是否交叉。交叉意味着几何体的内部有交集，但交集的维度低于两个几何体中维度较大者。

crosses 适用场景：

  线与线交叉：
  ─────┐
       │ ──────        crosses = True (交叉点是 0 维)
       │

  线穿过面：
  ────────────
     ┌────┤
     │    │            crosses = True (线穿过面)
     └────┘

from shapely import LineString, Polygon

# 两条交叉的线
line1 = LineString([(0, 0), (4, 4)])
line2 = LineString([(0, 4), (4, 0)])
print(f"线-线交叉: {line1.crosses(line2)}")  # True

# 平行线不交叉
line3 = LineString([(0, 1), (4, 5)])
print(f"平行线交叉: {line1.crosses(line3)}")  # False

# 线穿过多边形
poly = Polygon([(1, 1), (3, 1), (3, 3), (1, 3)])
line_through = LineString([(0, 2), (4, 2)])
print(f"线穿过面: {line_through.crosses(poly)}")  # True

# 线完全在多边形内——不是 crosses
line_inside = LineString([(1.5, 2), (2.5, 2)])
print(f"线在面内: {line_inside.crosses(poly)}")  # False

5.8 overlaps

5.8.1 重叠关系

overlaps 判断两个同维度几何体是否有部分重叠但不完全包含。

overlaps 示意：

  ┌──────┐
  │  A   ├──────┐
  │      │ 重叠 │     overlaps = True (部分重叠)
  └──────┤  B   │
         └──────┘

  ┌──────────────┐
  │  A           │
  │   ┌──────┐   │
  │   │  B   │   │   overlaps = False (B 被 A 包含)
  │   └──────┘   │
  └──────────────┘

from shapely import Polygon, LineString, Point

# 部分重叠的多边形
a = Polygon([(0, 0), (4, 0), (4, 4), (0, 4)])
b = Polygon([(2, 2), (6, 2), (6, 6), (2, 6)])

print(f"部分重叠: {a.overlaps(b)}")  # True

# 完全包含——不是 overlaps
small = Polygon([(1, 1), (3, 1), (3, 3), (1, 3)])
print(f"完全包含: {a.overlaps(small)}")  # False

# 不同维度——不适用 overlaps
line = LineString([(0, 0), (4, 4)])
print(f"面与线: {a.overlaps(line)}")  # False

# 部分重叠的线段
line1 = LineString([(0, 0), (3, 0)])
line2 = LineString([(1, 0), (5, 0)])
print(f"线段重叠: {line1.overlaps(line2)}")  # True

5.9 equals

5.9.1 拓扑相等

equals 判断两个几何体是否在拓扑上等价（点集完全相同）：

from shapely import Polygon

# 坐标起始点不同，但拓扑相同
poly1 = Polygon([(0, 0), (4, 0), (4, 4), (0, 4)])
poly2 = Polygon([(4, 0), (4, 4), (0, 4), (0, 0)])

print(f"拓扑相等: {poly1.equals(poly2)}")  # True
print(f"== 运算符: {poly1 == poly2}")      # True

# 不同形状
poly3 = Polygon([(0, 0), (5, 0), (5, 5), (0, 5)])
print(f"不同形状: {poly1.equals(poly3)}")  # False

5.9.2 `equals_exact(other, tolerance)`

在指定容差范围内判断相等：

from shapely import Point

p1 = Point(0, 0)
p2 = Point(0.0001, 0.0001)

print(f"精确相等: {p1.equals(p2)}")                    # False
print(f"容差 0.001: {p1.equals_exact(p2, 0.001)}")     # True
print(f"容差 0.0001: {p1.equals_exact(p2, 0.0001)}")   # False
print(f"容差 0.0002: {p1.equals_exact(p2, 0.0002)}")   # True

5.9.3 `equals_identical()`（Shapely 2.0+）

要求 WKB 表示完全相同：

import shapely
from shapely import Polygon

poly1 = Polygon([(0, 0), (1, 0), (1, 1), (0, 1)])
poly2 = Polygon([(1, 0), (1, 1), (0, 1), (0, 0)])  # 不同起始点
poly3 = Polygon([(0, 0), (1, 0), (1, 1), (0, 1)])  # 相同起始点

# equals: 拓扑相等
print(f"poly1.equals(poly2): {poly1.equals(poly2)}")  # True

# equals_identical: WKB 完全相同
print(f"identical(poly1, poly2): {shapely.equals_identical(poly1, poly2)}")  # False
print(f"identical(poly1, poly3): {shapely.equals_identical(poly1, poly3)}")  # True

5.10 快速坐标测试

5.10.1 `contains_xy()` 和 `intersects_xy()`

Shapely 2.0 提供了快速的坐标级别空间测试函数，避免创建临时 Point 对象：

import shapely
import numpy as np
from shapely import Polygon

poly = Polygon([(0, 0), (10, 0), (10, 10), (0, 10)])

# 单点测试
print(f"contains_xy(5, 5): {shapely.contains_xy(poly, 5, 5)}")    # True
print(f"contains_xy(15, 5): {shapely.contains_xy(poly, 15, 5)}")  # False

# 批量测试——向量化操作
x = np.array([5, 15, 5, 0, 10])
y = np.array([5, 5, 15, 0, 10])
results = shapely.contains_xy(poly, x, y)
print(f"批量测试: {results}")
# [True, False, False, True, True]

# intersects_xy 类似，但包括边界
results_intersects = shapely.intersects_xy(poly, x, y)
print(f"intersects_xy: {results_intersects}")

5.10.2 性能优势

import shapely
import numpy as np
from shapely import Point, Polygon
import time

poly = Polygon([(0, 0), (10, 0), (10, 10), (0, 10)])
n = 100000
x = np.random.uniform(-5, 15, n)
y = np.random.uniform(-5, 15, n)

# 方式 1：逐点创建 Point 对象（慢）
start = time.time()
results1 = [poly.contains(Point(xi, yi)) for xi, yi in zip(x, y)]
t1 = time.time() - start

# 方式 2：contains_xy 向量化（快）
start = time.time()
results2 = shapely.contains_xy(poly, x, y)
t2 = time.time() - start

print(f"逐点判断: {t1:.3f} 秒")
print(f"向量化判断: {t2:.3f} 秒")
print(f"加速比: {t1/t2:.1f}x")
# 通常加速 50-200 倍

5.11 `dwithin()` — 距离阈值测试

5.11.1 基本用法

dwithin 判断两个几何体的距离是否在指定阈值内，比先计算 distance() 再比较更高效：

import shapely
from shapely import Point, Polygon

p1 = Point(0, 0)
p2 = Point(3, 4)

# 两点距离为 5.0
print(f"距离: {p1.distance(p2)}")  # 5.0

# dwithin 距离阈值测试
print(f"距离 <= 5: {shapely.dwithin(p1, p2, 5.0)}")    # True
print(f"距离 <= 4.9: {shapely.dwithin(p1, p2, 4.9)}")  # False
print(f"距离 <= 6: {shapely.dwithin(p1, p2, 6.0)}")    # True

# 点到多边形的距离测试
poly = Polygon([(10, 10), (20, 10), (20, 20), (10, 20)])
print(f"点到多边形距离: {p1.distance(poly):.4f}")
print(f"距离 <= 15: {shapely.dwithin(p1, poly, 15.0)}")  # True
print(f"距离 <= 10: {shapely.dwithin(p1, poly, 10.0)}")  # False

5.11.2 批量距离测试

import shapely
import numpy as np
from shapely import Polygon

# 一个多边形和多个点
poly = Polygon([(0, 0), (10, 0), (10, 10), (0, 10)])
points = shapely.points(
    np.random.uniform(-20, 30, 1000),
    np.random.uniform(-20, 30, 1000)
)

# 找到距离多边形 5 个单位以内的所有点
within_5 = shapely.dwithin(points, poly, 5.0)
print(f"距离 ≤ 5 的点数: {within_5.sum()}")

5.12 完整谓词对比表

from shapely import Point, Polygon, LineString

# 准备测试几何体
poly = Polygon([(0, 0), (10, 0), (10, 10), (0, 10)])
p_inside = Point(5, 5)      # 内部
p_boundary = Point(0, 5)    # 边界
p_outside = Point(15, 5)    # 外部
p_vertex = Point(0, 0)      # 顶点

print("=" * 65)
print(f"{'谓词':25s} {'内部点':8s} {'边界点':8s} {'外部点':8s} {'顶点':8s}")
print("=" * 65)

predicates = [
    ('contains', lambda g: poly.contains(g)),
    ('contains_properly', lambda g: poly.contains_properly(g)),
    ('covers', lambda g: poly.covers(g)),
    ('intersects', lambda g: poly.intersects(g)),
    ('within (reversed)', lambda g: g.within(poly)),
    ('touches', lambda g: poly.touches(g)),
    ('disjoint', lambda g: poly.disjoint(g)),
]

for name, func in predicates:
    results = [func(p) for p in [p_inside, p_boundary, p_outside, p_vertex]]
    print(f"{name:25s} {str(results[0]):8s} {str(results[1]):8s} "
          f"{str(results[2]):8s} {str(results[3]):8s}")

# 输出：
# 谓词                       内部点    边界点    外部点    顶点
# ================================================================
# contains                  True     True     False    True
# contains_properly         True     False    False    False
# covers                    True     True     False    True
# intersects                True     True     False    True
# within (reversed)         True     True     False    True
# touches                   False    True     False    True
# disjoint                  False    False    True     False

5.13 谓词关系总结

5.13.1 谓词分类

谓词	说明	DE-9IM 模式
`equals`	拓扑相等	TF**FFF
`contains`	A 包含 B	T****FF
`contains_properly`	A 严格包含 B（不含边界）	—
`covers`	A 覆盖 B	T*FF 或 TFF 或 TFF*
`within`	A 在 B 内	TFF
`covered_by`	A 被 B 覆盖	TFF 或 TFF 或 FTF
`intersects`	A 与 B 有公共点	≠ FFFF***
`disjoint`	A 与 B 无公共点	FFFF***
`touches`	A 与 B 仅边界接触	FT*** 或 FT* 或 FT***
`crosses`	A 与 B 交叉	依维度组合而异
`overlaps`	A 与 B 部分重叠	依维度组合而异

5.13.2 互为逆操作的谓词

from shapely import Point, Polygon

poly = Polygon([(0, 0), (10, 0), (10, 10), (0, 10)])
p = Point(5, 5)

# 逆操作对
assert poly.contains(p) == p.within(poly)
assert poly.covers(p) == p.covered_by(poly)
assert poly.intersects(p) == p.intersects(poly)  # 对称
assert poly.touches(p) == p.touches(poly)         # 对称
assert poly.disjoint(p) == p.disjoint(poly)       # 对称

# 互为否定
assert poly.intersects(p) == (not poly.disjoint(p))

print("所有逆操作和否定关系验证通过 ✓")

5.14 使用 Prepared Geometry 加速

5.14.1 Prepared Geometry 概述

当需要用同一个几何体重复判断与多个其他几何体的空间关系时，使用 Prepared Geometry 可以显著提升性能。Prepared Geometry 会预先构建空间索引，加速后续的谓词判断。

5.14.2 使用方法

from shapely import prepare, Point, Polygon
import time
import numpy as np
import shapely

# 创建一个复杂多边形
poly = Polygon([(0, 0), (100, 0), (100, 100), (50, 80), (0, 100)])

# 准备几何体（预构建索引）
prepare(poly)

# 生成大量测试点
n = 100000
x = np.random.uniform(-10, 110, n)
y = np.random.uniform(-10, 110, n)
points = [Point(xi, yi) for xi, yi in zip(x, y)]

# Prepared 几何体的谓词判断自动加速
start = time.time()
results = [poly.contains(p) for p in points]
t = time.time() - start
print(f"Prepared 判断 {n} 个点: {t:.3f} 秒")
print(f"包含的点数: {sum(results)}")

5.14.3 向量化 + Prepared

import shapely
import numpy as np
from shapely import Polygon

poly = Polygon([(0, 0), (100, 0), (100, 100), (0, 100)])

# 准备几何体
shapely.prepare(poly)

# 向量化操作自动利用 prepared 几何体
n = 1000000
points = shapely.points(
    np.random.uniform(-10, 110, n),
    np.random.uniform(-10, 110, n)
)

import time
start = time.time()
results = shapely.contains(poly, points)
t = time.time() - start
print(f"向量化 + Prepared 判断 {n} 个点: {t:.3f} 秒")
print(f"包含的点数: {results.sum()}")

5.14.4 性能对比总结

方法	适用场景	相对性能
逐点 `contains()`	少量点，简单几何	1x（基准）
Prepared + 逐点	大量点，复杂几何	5-20x
`contains_xy()`	大量坐标，避免创建 Point	50-200x
向量化 + Prepared	大规模批量操作	100-500x

5.15 实战示例

5.15.1 地理围栏检测系统

from shapely import Point, Polygon, prepare
import shapely
import numpy as np

# 定义多个地理围栏区域
zones = {
    '办公区': Polygon([(0, 0), (100, 0), (100, 50), (0, 50)]),
    '生活区': Polygon([(0, 50), (100, 50), (100, 100), (0, 100)]),
    '禁止区': Polygon([(40, 40), (60, 40), (60, 60), (40, 60)]),
}

# 预先准备所有围栏几何体
for zone in zones.values():
    prepare(zone)

# 模拟设备位置
devices = {
    '设备A': Point(50, 25),
    '设备B': Point(50, 75),
    '设备C': Point(50, 50),
    '设备D': Point(150, 50),
}

# 检测每个设备所在的区域
for dev_name, dev_loc in devices.items():
    in_zones = []
    for zone_name, zone_poly in zones.items():
        if zone_poly.contains(dev_loc):
            in_zones.append(zone_name)
    
    if not in_zones:
        print(f"{dev_name}: 不在任何已知区域内 ⚠️")
    else:
        zone_str = ', '.join(in_zones)
        is_forbidden = '禁止区' in in_zones
        status = '🚫 警告！' if is_forbidden else '✅'
        print(f"{dev_name}: 位于 [{zone_str}] {status}")

5.15.2 路径与区域的关系分析

from shapely import LineString, Polygon

# 定义路径和区域
road = LineString([(0, 5), (20, 5)])
park = Polygon([(3, 2), (8, 2), (8, 8), (3, 8)])
lake = Polygon([(12, 3), (16, 3), (16, 7), (12, 7)])
building = Polygon([(18, 4), (20, 4), (20, 6), (18, 6)])

areas = {'公园': park, '湖泊': lake, '建筑': building}

print("道路与各区域的空间关系分析：")
print("-" * 50)

for name, area in areas.items():
    relations = []
    
    if road.intersects(area):
        relations.append("相交")
    if road.crosses(area):
        relations.append("穿越")
    if road.touches(area):
        relations.append("接触")
    if road.within(area):
        relations.append("完全在内")
    if road.disjoint(area):
        relations.append("不相交")
    
    # 计算穿过区域的路段长度
    if road.intersects(area):
        intersection = road.intersection(area)
        print(f"  {name}: {', '.join(relations)}")
        print(f"    穿过长度: {intersection.length:.2f}")
        print(f"    DE-9IM: {road.relate(area)}")
    else:
        print(f"  {name}: {', '.join(relations)}")
        print(f"    最近距离: {road.distance(area):.2f}")

5.16 本章小结

本章我们学习了 Shapely 中完整的空间关系与谓词判断体系：

DE-9IM 模型：理解 Interior、Boundary、Exterior 概念和 9 交叉矩阵
relate() 和 relate_pattern()：获取和匹配 DE-9IM 字符串
包含关系：contains、within、contains_properly、covers、covered_by
相交关系：intersects、disjoint
边界关系：touches
交叉/重叠：crosses、overlaps
相等判断：equals、equals_exact、equals_identical
快速测试：contains_xy、intersects_xy、dwithin
性能优化：Prepared Geometry 预构建索引加速重复判断

掌握这些空间谓词是进行空间分析的基础，下一章我们将学习几何集合运算。