第8章：空间度量与计算

空间度量是 GIS 分析的基础。Shapely 提供了丰富的度量函数，涵盖面积、长度、距离、形状相似度等多种计算。这些函数基于平面几何（欧几里得空间），能够高效地计算几何对象的各种度量属性。本章将全面介绍 Shapely 中的所有度量函数，包括数学含义、应用场景和使用方法。

area：面积计算

基本用法

area(geometry) 返回几何对象的面积。只有面状几何（Polygon / MultiPolygon）才有非零面积。

from shapely import area
from shapely.geometry import Polygon, Point, LineString

# 矩形面积
rect = Polygon([(0, 0), (4, 0), (4, 3), (0, 3)])
print(f"矩形面积: {area(rect)}")  # 12.0

# 三角形面积
tri = Polygon([(0, 0), (6, 0), (3, 4)])
print(f"三角形面积: {area(tri)}")  # 12.0

# 线和点的面积为 0
line = LineString([(0, 0), (1, 1)])
point = Point(0, 0)
print(f"线面积: {area(line)}")    # 0.0
print(f"点面积: {area(point)}")   # 0.0

带洞多边形的面积

from shapely import area
from shapely.geometry import Polygon

# 外环面积 100，内环面积 4
poly_with_hole = Polygon(
    [(0, 0), (10, 0), (10, 10), (0, 10)],
    [[(3, 3), (5, 3), (5, 5), (3, 5)]]
)
print(f"带洞多边形面积: {area(poly_with_hole)}")  # 96.0

MultiPolygon 面积

from shapely import area
from shapely.geometry import MultiPolygon, Polygon

mp = MultiPolygon([
    Polygon([(0, 0), (1, 0), (1, 1), (0, 1)]),
    Polygon([(2, 2), (4, 2), (4, 4), (2, 4)]),
])
print(f"MultiPolygon 面积: {area(mp)}")  # 1 + 4 = 5.0

批量面积计算

from shapely import area, box
import numpy as np

# 利用向量化进行批量计算
boxes = np.array([box(0, 0, i, i) for i in range(1, 6)])
areas = area(boxes)
print(f"批量面积: {areas}")  # [1. 4. 9. 16. 25.]

length：长度 / 周长

基本用法

length(geometry) 返回几何对象的长度。对于线状几何返回线长度，对于面状几何返回周长。

from shapely import length
from shapely.geometry import LineString, Polygon, Point

# 线长度
line = LineString([(0, 0), (3, 4)])
print(f"线长度: {length(line)}")  # 5.0

# 折线长度
polyline = LineString([(0, 0), (3, 0), (3, 4)])
print(f"折线长度: {length(polyline)}")  # 7.0

# 多边形周长
poly = Polygon([(0, 0), (4, 0), (4, 3), (0, 3)])
print(f"矩形周长: {length(poly)}")  # 14.0

# 点没有长度
point = Point(0, 0)
print(f"点长度: {length(point)}")  # 0.0

圆的周长

from shapely import length
from shapely.geometry import Point
import math

circle = Point(0, 0).buffer(5)
print(f"近似圆周长: {length(circle):.4f}")
print(f"理论圆周长: {2 * math.pi * 5:.4f}")

bounds：边界框

基本用法

bounds(geometry) 返回几何对象的边界框，格式为 (minx, miny, maxx, maxy)。

from shapely import bounds
from shapely.geometry import Polygon, LineString, Point

poly = Polygon([(1, 2), (5, 1), (4, 6), (2, 5)])
print(f"多边形边界框: {bounds(poly)}")  # (1.0, 1.0, 5.0, 6.0)

line = LineString([(0, 0), (3, 4)])
print(f"线边界框: {bounds(line)}")  # (0.0, 0.0, 3.0, 4.0)

point = Point(3, 7)
print(f"点边界框: {bounds(point)}")  # (3.0, 7.0, 3.0, 7.0)

批量边界框

from shapely import bounds, box
import numpy as np

geoms = np.array([
    box(0, 0, 1, 1),
    box(2, 3, 5, 6),
    box(-1, -1, 4, 4),
])

bbox = bounds(geoms)
print(f"批量边界框:\n{bbox}")

total_bounds：批量总边界框

total_bounds(geometries) 返回一组几何对象的总边界框。

from shapely import total_bounds, box
import numpy as np

geoms = np.array([
    box(0, 0, 2, 2),
    box(3, 1, 5, 4),
    box(1, 3, 4, 6),
])

tb = total_bounds(geoms)
print(f"总边界框: {tb}")  # (0.0, 0.0, 5.0, 6.0)

distance：欧几里得距离

基本用法

distance(a, b) 返回两个几何对象之间的最短欧几里得距离。

from shapely import distance
from shapely.geometry import Point, Polygon, LineString

# 两点之间的距离
p1 = Point(0, 0)
p2 = Point(3, 4)
print(f"两点距离: {distance(p1, p2)}")  # 5.0

# 点到线的距离
line = LineString([(0, 1), (4, 1)])
print(f"点到线距离: {distance(Point(2, 3), line)}")  # 2.0

# 点到多边形的距离
poly = Polygon([(2, 2), (6, 2), (6, 6), (2, 6)])
print(f"内部点距离: {distance(Point(3, 3), poly)}")  # 0.0（在内部）
print(f"外部点距离: {distance(Point(0, 0), poly)}")

相交几何体的距离

from shapely import distance
from shapely.geometry import Polygon

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

# 相交时距离为 0
print(f"相交几何体距离: {distance(a, b)}")  # 0.0

批量距离计算

from shapely import distance
from shapely.geometry import Point
import numpy as np

origin = Point(0, 0)
points = np.array([Point(i, i) for i in range(1, 6)])

# 向量化计算
distances = distance(origin, points)
print(f"批量距离: {distances}")

hausdorff_distance：Hausdorff 距离

概述

Hausdorff 距离是衡量两个几何对象形状相似度的指标。它定义为：对于 A 上的每个点，找到距离 B 最近的点的距离，取所有这些距离中的最大值；对 B 做同样的操作；取两者中的较大值。

直觉上，如果 Hausdorff 距离为 d，那么任何一个几何体上的每个点，到另一个几何体的距离都不超过 d。

基本用法

from shapely import hausdorff_distance
from shapely.geometry import LineString

line1 = LineString([(0, 0), (1, 0), (2, 0)])
line2 = LineString([(0, 1), (1, 1), (2, 1)])

hd = hausdorff_distance(line1, line2)
print(f"Hausdorff 距离: {hd}")  # 1.0

形状相似度比较

from shapely import hausdorff_distance
from shapely.geometry import LineString

# 原始曲线
original = LineString([(0, 0), (1, 1), (2, 0), (3, 1), (4, 0)])

# 近似曲线
approx1 = LineString([(0, 0), (2, 0), (4, 0)])
approx2 = LineString([(0, 0.1), (1, 1.1), (2, 0.1), (3, 1.1), (4, 0.1)])

hd1 = hausdorff_distance(original, approx1)
hd2 = hausdorff_distance(original, approx2)

print(f"粗略近似 Hausdorff: {hd1:.3f}")
print(f"精细近似 Hausdorff: {hd2:.3f}")

densify 参数

densify 参数用于在计算前对几何体进行加密，以获得更精确的结果：

from shapely import hausdorff_distance
from shapely.geometry import LineString

line1 = LineString([(0, 0), (10, 0)])
line2 = LineString([(0, 1), (5, 1), (10, 0)])

# 不加密
hd_raw = hausdorff_distance(line1, line2)
print(f"不加密: {hd_raw:.3f}")

# 加密（densify 表示每段长度与总长度的比例）
hd_dense = hausdorff_distance(line1, line2, densify=0.1)
print(f"加密后: {hd_dense:.3f}")

frechet_distance：Fréchet 距离

概述

Fréchet 距离（弗雷歇距离）是衡量两条曲线相似度的指标，直觉上类似于”遛狗距离”：一个人沿着一条曲线走，他的狗沿着另一条曲线走，两者都只能前进不能后退， Fréchet 距离就是在所有可能的行走方案中，最短狗绳的最小长度。

基本用法

from shapely import frechet_distance
from shapely.geometry import LineString

line1 = LineString([(0, 0), (1, 0), (2, 0), (3, 0)])
line2 = LineString([(0, 1), (1, 1), (2, 1), (3, 1)])

fd = frechet_distance(line1, line2)
print(f"Fréchet 距离: {fd}")  # 1.0

Fréchet vs Hausdorff

from shapely import frechet_distance, hausdorff_distance
from shapely.geometry import LineString

# 两条方向相反的线
line1 = LineString([(0, 0), (1, 0), (2, 0)])
line2 = LineString([(2, 1), (1, 1), (0, 1)])

hd = hausdorff_distance(line1, line2)
fd = frechet_distance(line1, line2)

print(f"Hausdorff 距离: {hd:.3f}")  # 只看最近距离的最大值
print(f"Fréchet 距离: {fd:.3f}")    # 考虑行进方向，通常更大

densify 参数

from shapely import frechet_distance
from shapely.geometry import LineString

line1 = LineString([(0, 0), (10, 0)])
line2 = LineString([(0, 0), (5, 5), (10, 0)])

fd = frechet_distance(line1, line2)
fd_dense = frechet_distance(line1, line2, densify=0.1)
print(f"Fréchet 距离: {fd:.3f}")
print(f"加密后 Fréchet: {fd_dense:.3f}")

minimum_bounding_radius：最小外接圆半径

from shapely import minimum_bounding_radius
from shapely.geometry import Polygon, MultiPoint

poly = Polygon([(0, 0), (4, 0), (4, 3), (0, 3)])
radius = minimum_bounding_radius(poly)
print(f"最小外接圆半径: {radius:.3f}")  # 对角线一半

points = MultiPoint([(0, 0), (3, 0), (0, 4)])
radius_pts = minimum_bounding_radius(points)
print(f"三角形外接圆半径: {radius_pts:.3f}")

minimum_clearance：最小间隙

概述

minimum_clearance(geometry) 返回几何对象中，移动任何顶点使几何体变为无效（如自相交）所需的最小距离。这是衡量几何体”稳定性”的指标。

基本用法

from shapely import minimum_clearance
from shapely.geometry import Polygon

# 接近退化的矩形
poly = Polygon([(0, 0), (10, 0), (10, 0.001), (0, 0.001)])
mc = minimum_clearance(poly)
print(f"最小间隙: {mc}")  # 接近 0.001

# 正方形
square = Polygon([(0, 0), (1, 0), (1, 1), (0, 1)])
mc_sq = minimum_clearance(square)
print(f"正方形最小间隙: {mc_sq}")  # 1.0

minimum_clearance_line：最小间隙线

返回表示最小间隙的线段。

from shapely import minimum_clearance, minimum_clearance_line
from shapely.geometry import Polygon

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

mc = minimum_clearance(poly)
mc_line = minimum_clearance_line(poly)

print(f"最小间隙: {mc:.3f}")
print(f"最小间隙线: {mc_line}")
print(f"线长度: {mc_line.length:.3f}")

minimum_width：最小宽度

返回几何对象的最小宽度。

from shapely import minimum_width
from shapely.geometry import Polygon

# 窄矩形
narrow = Polygon([(0, 0), (10, 0), (10, 1), (0, 1)])
mw = minimum_width(narrow)
print(f"最小宽度: {mw.length}")  # 1.0

# 三角形
tri = Polygon([(0, 0), (6, 0), (3, 2)])
mw_tri = minimum_width(tri)
print(f"三角形最小宽度: {mw_tri.length:.3f}")

maximum_inscribed_circle：最大内切圆

返回几何对象内部的最大内切圆。

from shapely import maximum_inscribed_circle
from shapely.geometry import Polygon

# 正方形
square = Polygon([(0, 0), (10, 0), (10, 10), (0, 10)])
mic = maximum_inscribed_circle(square, tolerance=0.01)
print(f"最大内切圆类型: {mic.geom_type}")  # LineString（圆心到边的线）
print(f"内切圆半径: {mic.length:.2f}")  # 5.0

# L形
L_shape = Polygon([
    (0, 0), (6, 0), (6, 2), (2, 2), (2, 6), (0, 6)
])
mic_L = maximum_inscribed_circle(L_shape, tolerance=0.01)
print(f"L形最大内切圆半径: {mic_L.length:.2f}")

shortest_line：最短连线

返回连接两个几何对象的最短线段。

from shapely import shortest_line
from shapely.geometry import Point, Polygon, LineString

p = Point(0, 0)
poly = Polygon([(3, 3), (6, 3), (6, 6), (3, 6)])

sl = shortest_line(p, poly)
print(f"最短连线: {sl}")
print(f"长度: {sl.length:.3f}")

from shapely import shortest_line
from shapely.geometry import LineString

line1 = LineString([(0, 0), (4, 0)])
line2 = LineString([(1, 3), (3, 3)])

sl = shortest_line(line1, line2)
print(f"最短连线: {sl}")
print(f"起点: ({sl.coords[0][0]:.1f}, {sl.coords[0][1]:.1f})")
print(f"终点: ({sl.coords[1][0]:.1f}, {sl.coords[1][1]:.1f})")

批量最短连线

from shapely import shortest_line
from shapely.geometry import Point, LineString
import numpy as np

road = LineString([(0, 0), (10, 0)])
houses = np.array([Point(i, j) for i, j in [(2, 3), (5, 1), (8, 4)]])

lines = shortest_line(houses, road)
for i, sl in enumerate(lines):
    print(f"房屋{i} 到道路最短连线长度: {sl.length:.2f}")

批量计算与向量化操作

NumPy 数组输入

Shapely 2.x 的度量函数支持 NumPy 数组输入，实现高效的批量计算：

from shapely import area, length, distance, bounds
from shapely.geometry import Point
import numpy as np

# 创建大量随机多边形
np.random.seed(42)
n = 1000
points = np.array([Point(x, y).buffer(r)
                    for x, y, r in zip(
                        np.random.uniform(0, 100, n),
                        np.random.uniform(0, 100, n),
                        np.random.uniform(0.5, 3, n),
                    )])

# 批量计算
areas = area(points)
perimeters = length(points)

print(f"面积统计: min={areas.min():.2f}, max={areas.max():.2f}, mean={areas.mean():.2f}")
print(f"周长统计: min={perimeters.min():.2f}, max={perimeters.max():.2f}")

距离矩阵

from shapely import distance
from shapely.geometry import Point
import numpy as np

# 计算点之间的距离矩阵
pts = np.array([Point(0, 0), Point(1, 0), Point(0, 1), Point(1, 1)])
n = len(pts)

dist_matrix = np.zeros((n, n))
for i in range(n):
    # 向量化计算一行
    dist_matrix[i] = distance(pts[i], pts)

print("距离矩阵:")
print(np.round(dist_matrix, 3))

度量函数应用场景

场景 1：最近设施分析

from shapely import distance
from shapely.geometry import Point
import numpy as np

# 设施位置
facilities = np.array([
    Point(10, 20), Point(30, 10), Point(50, 40),
    Point(20, 50), Point(40, 25),
])

# 查询点
query = Point(25, 30)

# 找到最近的设施
distances = distance(query, facilities)
nearest_idx = np.argmin(distances)
print(f"最近的设施索引: {nearest_idx}")
print(f"最近距离: {distances[nearest_idx]:.2f}")

场景 2：形状紧凑度分析

from shapely import area, length
from shapely.geometry import Polygon, Point
import math

shapes = {
    "正方形": Polygon([(0, 0), (1, 0), (1, 1), (0, 1)]),
    "长矩形": Polygon([(0, 0), (10, 0), (10, 1), (0, 1)]),
    "圆形": Point(0, 0).buffer(1),
    "L形": Polygon([(0, 0), (2, 0), (2, 1), (1, 1), (1, 2), (0, 2)]),
}

print("形状紧凑度分析（4π × 面积 / 周长²）:")
for name, shape in shapes.items():
    a = area(shape)
    p = length(shape)
    compactness = 4 * math.pi * a / (p * p)
    print(f"  {name}: 面积={a:.2f}, 周长={p:.2f}, 紧凑度={compactness:.4f}")

场景 3：多边形面积与周长统计

from shapely import area, length
from shapely.geometry import Polygon
import numpy as np

# 随机多边形
np.random.seed(0)
polygons = []
for _ in range(20):
    cx, cy = np.random.uniform(0, 100, 2)
    r = np.random.uniform(1, 10)
    angles = np.sort(np.random.uniform(0, 2 * np.pi, 6))
    coords = [(cx + r * np.cos(a), cy + r * np.sin(a)) for a in angles]
    coords.append(coords[0])
    try:
        p = Polygon(coords)
        if p.is_valid and not p.is_empty:
            polygons.append(p)
    except Exception:
        pass

geoms = np.array(polygons)
areas = area(geoms)
perimeters = length(geoms)

print(f"有效多边形数: {len(geoms)}")
print(f"平均面积: {areas.mean():.2f}")
print(f"平均周长: {perimeters.mean():.2f}")
print(f"面积标准差: {areas.std():.2f}")

注意事项

平面计算 vs 球面计算

Shapely 只进行平面（笛卡尔坐标系）计算，不支持球面（地理坐标系）计算。

from shapely import distance
from shapely.geometry import Point

# 这是坐标差值，不是真实的地球表面距离！
beijing = Point(116.4, 39.9)
shanghai = Point(121.5, 31.2)

# 这个结果的单位是"度"，不是"米"或"千米"
d = distance(beijing, shanghai)
print(f"坐标距离: {d:.2f}（度，不是实际距离）")

# 如需计算真实地球距离，应该：
# 1. 将坐标投影到合适的平面坐标系（如 UTM）
# 2. 使用 pyproj 或 geopy 进行球面距离计算

浮点精度

from shapely import area
from shapely.geometry import Polygon

# 极小多边形可能有精度问题
tiny = Polygon([
    (0, 0),
    (1e-15, 0),
    (1e-15, 1e-15),
    (0, 1e-15),
])
print(f"极小面积: {area(tiny)}")

# 极大坐标也可能有精度问题
huge = Polygon([
    (1e15, 1e15),
    (1e15 + 1, 1e15),
    (1e15 + 1, 1e15 + 1),
    (1e15, 1e15 + 1),
])
print(f"极大坐标面积: {area(huge)}")

空几何体的度量

from shapely import area, length, distance, bounds
from shapely import wkt

empty = wkt.loads("GEOMETRYCOLLECTION EMPTY")

print(f"空几何面积: {area(empty)}")     # 0.0
print(f"空几何长度: {length(empty)}")   # 0.0
print(f"空几何边界框: {bounds(empty)}")

度量函数速查表

函数	说明	适用类型
`area(geom)`	面积	Polygon, MultiPolygon
`length(geom)`	长度/周长	Line*, Polygon
`bounds(geom)`	边界框	所有类型
`total_bounds(geoms)`	总边界框	几何数组
`distance(a, b)`	最短距离	所有类型
`hausdorff_distance(a, b)`	Hausdorff 距离	所有类型
`frechet_distance(a, b)`	Fréchet 距离	LineString
`minimum_bounding_radius(geom)`	最小外接圆半径	所有类型
`minimum_clearance(geom)`	最小间隙	所有类型
`minimum_clearance_line(geom)`	最小间隙线	所有类型
`minimum_width(geom)`	最小宽度	Polygon
`maximum_inscribed_circle(geom)`	最大内切圆	Polygon
`shortest_line(a, b)`	最短连线	所有类型

本章小结

本章全面介绍了 Shapely 的空间度量函数：

基础度量：area、length、bounds
距离度量：distance、hausdorff_distance、frechet_distance
形状度量：minimum_clearance、minimum_width、maximum_inscribed_circle
连接度量：shortest_line

关键要点：

Shapely 只做平面计算，地理坐标需先投影
度量函数支持 NumPy 数组输入，可高效批量计算
注意浮点精度问题和空几何体的特殊行为

下一章将介绍仿射变换，学习如何对几何对象进行平移、旋转、缩放等操作。