znlgis 博客

GIS开发与技术分享

第13章：构造性几何操作

构造性几何操作是指从已有的几何对象中提取或构造新的几何特征的一系列方法。与集合运算不同，构造性操作通常作用于单个几何对象或几何集合，生成诸如边界、质心、凸包、三角剖分等派生几何。这些操作在空间分析、可视化和数据处理中有着广泛的应用。

13.1 构造性操作概述

13.1.1 什么是构造性几何操作

构造性几何操作（Constructive Geometric Operations）是指从一个或多个几何对象出发，通过数学变换或算法，生成新的几何对象的操作。这些操作不改变原始几何对象，而是返回新的几何结果。

13.1.2 操作分类

构造性操作可以分为以下几类：

类别 操作 说明
特征提取 boundary, exterior, interiors 提取几何的组成部分
特征点 centroid, representative_point 计算代表性点位
包络几何 envelope, convex_hull, concave_hull 构建包围几何
剖分与分区 delaunay_triangles, voronoi_polygons 空间分区
最值几何 minimum_bounding_circle, maximum_inscribed_circle, minimum_rotated_rectangle 极值包围几何
几何变换 line_merge, polygonize, build_area, normalize 几何类型转换
顶点操作 extract_unique_points 提取和操作顶点

13.1.3 在 GeoPandas 中的调用方式

构造性操作在 GeoPandas 中有两种调用方式：

# 1. 作为属性（property）访问
centroid = gdf.geometry.centroid
boundary = gdf.geometry.boundary

# 2. 作为方法（method）调用
hull = gdf.geometry.convex_hull
triangles = gdf.geometry.delaunay_triangles()

13.2 boundary - 提取边界

13.2.1 基本概念

boundary 属性返回几何对象的边界。不同类型几何对象的边界：

几何类型 边界
Point 空几何（GEOMETRYCOLLECTION EMPTY）
LineString 端点组成的 MultiPoint
Polygon 外环和内环组成的 MultiLineString 或 LineString
MultiPolygon 所有多边形边界的集合

13.2.2 代码示例

import geopandas as gpd
from shapely.geometry import Point, LineString, Polygon

# 不同几何类型的边界
gs = gpd.GeoSeries([
    Point(0, 0),
    LineString([(0, 0), (1, 1), (2, 0)]),
    Polygon([(0, 0), (2, 0), (2, 2), (0, 2)]),
    Polygon([(0, 0), (4, 0), (4, 4), (0, 4)],
            holes=[[(1, 1), (3, 1), (3, 3), (1, 3)]]),  # 带孔多边形
])

boundaries = gs.boundary
for i, b in enumerate(boundaries):
    print(f"几何 {i} ({gs.iloc[i].geom_type}): 边界类型={b.geom_type}")
# 几何 0 (Point): 边界类型=GeometryCollection  (空)
# 几何 1 (LineString): 边界类型=MultiPoint  (两个端点)
# 几何 2 (Polygon): 边界类型=LineString  (外环)
# 几何 3 (Polygon): 边界类型=MultiLineString  (外环+内环)

13.2.3 实际应用

# 提取行政区边界线
districts = gpd.read_file("districts.shp")
district_boundaries = districts.copy()
district_boundaries['geometry'] = districts.boundary

# 边界线可以用于：
# 1. 边界线绘图（避免面填充）
# 2. 边界长度计算
district_boundaries['perimeter'] = district_boundaries.geometry.length

13.3 centroid - 质心计算

13.3.1 基本概念

centroid 属性返回几何对象的几何质心（重心）。质心是几何对象的”重心”位置，对于均匀密度的平面图形，质心就是物理重心。

13.3.2 代码示例

import geopandas as gpd
from shapely.geometry import Polygon, LineString, MultiPoint

gs = gpd.GeoSeries([
    Polygon([(0, 0), (4, 0), (4, 2), (0, 2)]),   # 矩形
    Polygon([(0, 0), (4, 0), (2, 3)]),             # 三角形
    LineString([(0, 0), (4, 0), (4, 2)]),           # 折线
    MultiPoint([(0, 0), (4, 0), (2, 3)]),           # 多点
])

centroids = gs.centroid
for i, c in enumerate(centroids):
    print(f"几何 {i}: 质心 = ({c.x:.2f}, {c.y:.2f})")
# 几何 0: 质心 = (2.00, 1.00)  矩形中心
# 几何 1: 质心 = (2.00, 1.00)  三角形重心
# 几何 2: 质心 = (2.67, 0.67)  线的质心
# 几何 3: 质心 = (2.00, 1.00)  点集的平均位置

13.3.3 注意事项

  • 质心可能不在几何对象内部（如凹多边形、多环多边形）
  • 如果需要保证点在几何内部，使用 representative_point()
  • 地理坐标系下的质心计算可能不准确，建议先投影
from shapely.geometry import Polygon

# L 形多边形，质心在外部
l_shape = Polygon([(0, 0), (2, 0), (2, 1), (1, 1), (1, 2), (0, 2)])
centroid = l_shape.centroid
print(f"L形质心: ({centroid.x:.2f}, {centroid.y:.2f})")
print(f"质心在多边形内: {l_shape.contains(centroid)}")
# 质心可能在多边形外部！

13.4 envelope - 最小外接矩形

13.4.1 基本概念

envelope 属性返回几何对象的最小外接矩形（Minimum Bounding Rectangle, MBR），也称为边界框（Bounding Box）。这个矩形的边与坐标轴平行。

13.4.2 代码示例

import geopandas as gpd
from shapely.geometry import Polygon, LineString, MultiPoint

gs = gpd.GeoSeries([
    Polygon([(0, 0), (3, 0), (2, 3), (1, 4)]),   # 不规则多边形
    LineString([(0, 0), (5, 3), (2, 5)]),           # 折线
    MultiPoint([(1, 1), (4, 2), (3, 5)]),           # 多点
])

envelopes = gs.envelope
for i, env in enumerate(envelopes):
    print(f"几何 {i}: 外接矩形 = {env}")

13.4.3 bounds 属性

除了 envelope，还可以使用 bounds 属性获取数值形式的边界框：

# bounds 返回 (minx, miny, maxx, maxy) 的 DataFrame
print(gs.bounds)
#    minx  miny  maxx  maxy
# 0   0.0   0.0   3.0   4.0
# 1   0.0   0.0   5.0   5.0
# 2   1.0   1.0   4.0   5.0

# total_bounds 返回整个 GeoSeries 的边界框
print(gs.total_bounds)
# [0. 0. 5. 5.]

13.4.4 实际应用

# 为空间查询创建边界框过滤
import geopandas as gpd

# 使用边界框进行快速的粗过滤
target = gpd.read_file("target.shp").geometry.iloc[0]
bbox = target.envelope

# 先用边界框过滤（快）
candidates = gpd.read_file("large_data.shp")
coarse_filter = candidates[candidates.intersects(bbox)]

# 再用精确几何过滤（慢但精确）
fine_filter = coarse_filter[coarse_filter.intersects(target)]

13.5 convex_hull - 凸包

13.5.1 基本概念

convex_hull 属性返回几何对象的凸包（Convex Hull）。凸包是包含几何对象的最小凸多边形，可以想象为用一根橡皮筋围绕几何对象时形成的形状。

13.5.2 代码示例

import geopandas as gpd
from shapely.geometry import MultiPoint, Polygon

# 从点集生成凸包
points = MultiPoint([(0, 0), (1, 3), (3, 3), (4, 0), (2, -1)])

gs = gpd.GeoSeries([points])
hull = gs.convex_hull
print(f"凸包类型: {hull.iloc[0].geom_type}")  # Polygon
print(f"凸包面积: {hull.iloc[0].area:.2f}")

# 凹多边形的凸包
concave = Polygon([(0, 0), (4, 0), (4, 4), (2, 2), (0, 4)])
gs2 = gpd.GeoSeries([concave])
hull2 = gs2.convex_hull
print(f"原始面积: {concave.area:.2f}")
print(f"凸包面积: {hull2.iloc[0].area:.2f}")  # 凸包面积 >= 原始面积

13.5.3 凸包的性质

  • 凸包是包含所有输入点的最小凸多边形
  • 凸包的面积 ≥ 原始几何的面积
  • 凸多边形的凸包就是其自身
  • 凸包的顶点都是原始几何的顶点

13.5.4 实际应用

# 计算一组 GPS 轨迹点的活动范围
track_points = gpd.read_file("gps_tracks.shp")

# 凸包表示活动范围的大致边界
activity_area = track_points.union_all().convex_hull
print(f"活动范围面积: {activity_area.area:.2f}")

# 计算凸度（紧凑度）
# 凸度 = 原始面积 / 凸包面积，值越接近 1 形状越凸
original_area = track_points.union_all().area
convexity = original_area / activity_area.area if activity_area.area > 0 else 0
print(f"凸度: {convexity:.4f}")

13.6 concave_hull() - 凹包

13.6.1 基本概念

concave_hull() 计算几何对象的凹包（Concave Hull），也称为 alpha 形状。凹包比凸包更紧密地包围点集，可以呈现非凸的形状。

13.6.2 参数说明

GeoSeries.concave_hull(ratio=0.0, allow_holes=False)
参数 说明
ratio 0.0-1.0，控制凹陷程度。0.0 最凹，1.0 等同于凸包
allow_holes 是否允许结果中包含孔洞

13.6.3 代码示例

import geopandas as gpd
from shapely.geometry import MultiPoint
import numpy as np

# 生成 C 形分布的点
theta = np.linspace(0.5, 2 * np.pi - 0.5, 50)
x = np.cos(theta) * 10
y = np.sin(theta) * 10
points = MultiPoint(list(zip(x, y)))

gs = gpd.GeoSeries([points])

# 不同 ratio 的凹包
hull_convex = gs.convex_hull           # 凸包
hull_concave_1 = gs.concave_hull(ratio=0.3)   # 中等凹陷
hull_concave_2 = gs.concave_hull(ratio=0.0)   # 最大凹陷

print(f"凸包面积: {hull_convex.iloc[0].area:.2f}")
print(f"凹包(0.3)面积: {hull_concave_1.iloc[0].area:.2f}")
print(f"凹包(0.0)面积: {hull_concave_2.iloc[0].area:.2f}")

13.6.4 ratio 参数的影响

ratio 值 效果
0.0 最紧密的包围，可能非常不规则
0.3 较紧密的包围
0.5 中等凹陷
0.7 接近凸包
1.0 等同于凸包

13.6.5 实际应用

# 从建筑物点位推断建成区范围
buildings = gpd.read_file("buildings.shp")

# 合并所有建筑物点
all_buildings = buildings.union_all()

# 使用凹包获取建成区边界
# ratio 调整为适当值以获取合理的边界
urban_boundary = gpd.GeoSeries([all_buildings]).concave_hull(ratio=0.2)
print(f"建成区面积: {urban_boundary.iloc[0].area:.2f}")

13.7 exterior 与 interiors - 多边形外环与内环

13.7.1 exterior - 外环

exterior 属性返回多边形的外部环（外边界），类型为 LinearRing。

import geopandas as gpd
from shapely.geometry import Polygon

# 带孔多边形
poly = Polygon(
    [(0, 0), (10, 0), (10, 10), (0, 10)],  # 外环
    [[(2, 2), (4, 2), (4, 4), (2, 4)],      # 内环1（孔1）
     [(6, 6), (8, 6), (8, 8), (6, 8)]]      # 内环2（孔2）
)

gs = gpd.GeoSeries([poly])
exterior = gs.exterior
print(f"外环类型: {exterior.iloc[0].geom_type}")  # LinearRing
print(f"外环坐标数: {len(exterior.iloc[0].coords)}")
print(f"外环坐标: {list(exterior.iloc[0].coords)}")

13.7.2 interiors - 内环

interiors 属性返回多边形的内部环（孔洞边界）列表。

# 获取内环
interiors = gs.interiors
print(f"内环数量: {len(interiors.iloc[0])}")

for i, ring in enumerate(interiors.iloc[0]):
    print(f"内环 {i}: {ring}")

13.7.3 实际应用

# 分析多边形的孔洞
import geopandas as gpd

polygons = gpd.read_file("buildings_with_courtyards.shp")

# 统计每个多边形的孔洞数量
polygons['hole_count'] = polygons.geometry.apply(
    lambda g: len(list(g.interiors)) if g.geom_type == 'Polygon' else 0
)

# 提取所有孔洞作为独立的多边形
holes = []
for idx, row in polygons.iterrows():
    if row.geometry.geom_type == 'Polygon':
        for interior in row.geometry.interiors:
            holes.append(Polygon(interior))

holes_gdf = gpd.GeoDataFrame(geometry=holes)
print(f"共有 {len(holes_gdf)} 个庭院/孔洞")

13.8 representative_point() - 代表点

13.8.1 基本概念

representative_point() 返回一个保证在几何对象内部的点。与质心（centroid）不同，代表点总是位于几何对象内部或边界上。

13.8.2 代码示例

import geopandas as gpd
from shapely.geometry import Polygon

# L 形多边形（质心可能在外部）
l_shape = Polygon([(0, 0), (3, 0), (3, 1), (1, 1), (1, 3), (0, 3)])

gs = gpd.GeoSeries([l_shape])

centroid = gs.centroid.iloc[0]
rep_point = gs.representative_point().iloc[0]

print(f"质心: ({centroid.x:.2f}, {centroid.y:.2f})")
print(f"质心在多边形内: {l_shape.contains(centroid)}")

print(f"代表点: ({rep_point.x:.2f}, {rep_point.y:.2f})")
print(f"代表点在多边形内: {l_shape.contains(rep_point)}")
# 代表点总是在多边形内！

13.8.3 centroid vs representative_point

特性 centroid representative_point
定义 几何重心 保证在几何内部的点
可能在几何外部 是（凹形状） 否
位置稳定性 稳定 可能变化
用途 重心位置、标注位置 标注位置（安全）

13.8.4 实际应用

# 在地图上为每个区域添加标注
districts = gpd.read_file("districts.shp")

# 使用代表点作为标注位置（保证在区域内）
label_points = districts.representative_point()
districts['label_x'] = label_points.x
districts['label_y'] = label_points.y

# 绑定到地图标注
import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(10, 10))
districts.plot(ax=ax, edgecolor='black', facecolor='lightblue')

for idx, row in districts.iterrows():
    ax.annotate(row['name'],
                xy=(row['label_x'], row['label_y']),
                ha='center', fontsize=8)

13.9 delaunay_triangles() - Delaunay 三角剖分

13.9.1 基本概念

Delaunay 三角剖分是一种将点集连接成三角形网络的方法，具有以下性质：

  • 没有点位于任何三角形的外接圆内部
  • 最大化最小角度（避免细长三角形）
  • 是 Voronoi 图的对偶

13.9.2 代码示例

import geopandas as gpd
from shapely.geometry import MultiPoint

# 创建点集
points = MultiPoint([
    (0, 0), (1, 0), (2, 0),
    (0.5, 1), (1.5, 1),
    (1, 2)
])

gs = gpd.GeoSeries([points])

# 生成 Delaunay 三角剖分
triangles = gs.delaunay_triangles()
print(f"结果类型: {triangles.iloc[0].geom_type}")
# GeometryCollection，包含多个三角形

# 提取各个三角形
for i, tri in enumerate(triangles.iloc[0].geoms):
    print(f"三角形 {i}: {tri}")
    print(f"  面积: {tri.area:.2f}")

13.9.3 参数选项

# tolerance 参数控制捕捉容差
# only_edges=True 返回边而不是三角形
edges = gs.delaunay_triangles(only_edges=True)
print(f"Delaunay 边: {edges.iloc[0].geom_type}")
# MultiLineString，包含所有三角剖分的边

13.9.4 实际应用

# 从离散采样点构建连续的地形模型
elevation_points = gpd.read_file("elevation_samples.shp")

# 合并所有点
all_points = elevation_points.union_all()

# 生成 TIN（不规则三角网）
tin = gpd.GeoSeries([all_points]).delaunay_triangles()

# 将三角形转换为 GeoDataFrame
triangles_list = list(tin.iloc[0].geoms)
tin_gdf = gpd.GeoDataFrame(geometry=triangles_list, crs=elevation_points.crs)

print(f"TIN 包含 {len(tin_gdf)} 个三角形")

13.10 constrained_delaunay_triangles() - 约束 Delaunay 三角剖分

13.10.1 基本概念

约束 Delaunay 三角剖分在标准 Delaunay 三角剖分的基础上，可以强制保留输入几何对象的边。这在处理多边形三角剖分时特别有用。

13.10.2 代码示例

import geopandas as gpd
from shapely.geometry import Polygon

# 对多边形进行约束三角剖分
poly = Polygon([(0, 0), (4, 0), (4, 4), (2, 2), (0, 4)])

gs = gpd.GeoSeries([poly])
constrained = gs.constrained_delaunay_triangles()

print(f"约束三角剖分结果: {constrained.iloc[0].geom_type}")
# 三角形会遵循多边形的边界

13.10.3 与标准 Delaunay 的区别

特性 标准 Delaunay 约束 Delaunay
输入 点集 几何对象（含边）
边约束 无 保留输入边
最优性 完全 Delaunay 尽可能 Delaunay
应用 点集插值 多边形剖分

13.11 voronoi_polygons() - Voronoi 图

13.11.1 基本概念

Voronoi 图（也称为泰森多边形）是一种空间分区方法。对于每个输入点，Voronoi 多边形包含所有距离该点最近的空间位置。

13.11.2 代码示例

import geopandas as gpd
from shapely.geometry import MultiPoint

# 创建点集
points = MultiPoint([
    (1, 1), (3, 1), (5, 1),
    (2, 3), (4, 3),
    (3, 5)
])

gs = gpd.GeoSeries([points])

# 生成 Voronoi 多边形
voronoi = gs.voronoi_polygons()
print(f"Voronoi 结果: {voronoi.iloc[0].geom_type}")
# GeometryCollection，包含多个多边形

# 提取各个 Voronoi 多边形
for i, poly in enumerate(voronoi.iloc[0].geoms):
    print(f"Voronoi 多边形 {i}: 面积={poly.area:.2f}")

13.11.3 参数选项

# extend_to 参数指定 Voronoi 图的范围
from shapely.geometry import box

# 限制 Voronoi 图的范围
extent = box(0, 0, 6, 6)
voronoi_clipped = gs.voronoi_polygons(extend_to=extent)

# only_edges=True 返回 Voronoi 边而不是多边形
voronoi_edges = gs.voronoi_polygons(only_edges=True)

13.11.4 Delaunay 与 Voronoi 的对偶关系

Delaunay 三角剖分和 Voronoi 图互为对偶：

  • Delaunay 三角形的外接圆圆心是 Voronoi 顶点
  • Delaunay 的边与 Voronoi 的边垂直相交
  • Delaunay 的三角形对应 Voronoi 的顶点

13.11.5 实际应用

# 使用 Voronoi 图划分服务区域
import geopandas as gpd
from shapely.geometry import MultiPoint

# 气象站位置
stations = gpd.read_file("weather_stations.shp")

# 生成 Voronoi 多边形（每个站的最近服务区）
all_stations = stations.union_all()
voronoi = gpd.GeoSeries([all_stations]).voronoi_polygons()

# 将 Voronoi 多边形转换为 GeoDataFrame
voronoi_polys = list(voronoi.iloc[0].geoms)
service_areas = gpd.GeoDataFrame(
    geometry=voronoi_polys,
    crs=stations.crs
)

# 裁剪到研究区域
study_area = gpd.read_file("study_area.shp").geometry.iloc[0]
service_areas['geometry'] = service_areas.intersection(study_area)
service_areas = service_areas[~service_areas.is_empty]

print(f"划分了 {len(service_areas)} 个服务区域")

13.12 minimum_bounding_circle() - 最小外接圆

13.12.1 基本概念

minimum_bounding_circle() 计算包含几何对象的最小面积的圆（最小外接圆）。

13.12.2 代码示例

import geopandas as gpd
from shapely.geometry import MultiPoint, Polygon

# 不规则多边形
poly = Polygon([(0, 0), (3, 0), (4, 2), (2, 4), (0, 3)])

gs = gpd.GeoSeries([poly])
circle = gs.minimum_bounding_circle()

print(f"最小外接圆: {circle.iloc[0].geom_type}")  # Polygon（近似圆）
print(f"圆面积: {circle.iloc[0].area:.2f}")

# 圆的中心和半径
circle_centroid = circle.iloc[0].centroid
print(f"圆心: ({circle_centroid.x:.2f}, {circle_centroid.y:.2f})")

13.12.3 实际应用

# 评估地块的紧凑度
parcels = gpd.read_file("parcels.shp")

# 紧凑度 = 面积 / 最小外接圆面积
circles = parcels.minimum_bounding_circle()
parcels['compactness'] = parcels.area / circles.area

# 紧凑度接近 1 表示形状接近圆形
print(parcels[['name', 'compactness']].sort_values('compactness', ascending=False))

13.13 maximum_inscribed_circle() - 最大内切圆

13.13.1 基本概念

maximum_inscribed_circle() 计算完全在几何对象内部的最大圆。返回的是一条从圆心到边界最近点的线段，线段长度即为半径。

13.13.2 代码示例

import geopandas as gpd
from shapely.geometry import Polygon

# 创建多边形
poly = Polygon([(0, 0), (4, 0), (4, 2), (0, 2)])

gs = gpd.GeoSeries([poly])
mic = gs.maximum_inscribed_circle(tolerance=0.01)

print(f"结果类型: {mic.iloc[0].geom_type}")  # LineString
print(f"半径: {mic.iloc[0].length:.2f}")       # 内切圆半径

# 圆心是线段的起点
center = mic.iloc[0].coords[0]
print(f"圆心: ({center[0]:.2f}, {center[1]:.2f})")

13.13.3 tolerance 参数

tolerance 参数控制计算精度。较小的容差值会得到更精确的结果，但计算时间更长。

# 不同容差的比较
for tol in [1.0, 0.1, 0.01, 0.001]:
    mic = gs.maximum_inscribed_circle(tolerance=tol)
    print(f"容差 {tol}: 半径 = {mic.iloc[0].length:.6f}")

13.13.4 实际应用

# 评估建筑地块可用空间
parcels = gpd.read_file("parcels.shp")

# 最大内切圆半径反映了地块的有效宽度
mic = parcels.maximum_inscribed_circle(tolerance=0.1)
parcels['max_radius'] = mic.length

# 筛选足够宽的地块（如半径 > 20 米）
suitable = parcels[parcels['max_radius'] > 20]
print(f"满足宽度要求的地块: {len(suitable)} 个")

13.14 minimum_rotated_rectangle() - 最小旋转外接矩形

13.14.1 基本概念

minimum_rotated_rectangle() 计算包含几何对象的最小面积矩形，这个矩形可以旋转（不必与坐标轴平行）。也称为最小面积外接矩形（Minimum Area Bounding Rectangle）。

13.14.2 与 envelope 的区别

import geopandas as gpd
from shapely.geometry import Polygon

# 斜长条形多边形
poly = Polygon([(0, 0), (4, 3), (3, 4), (-1, 1)])

gs = gpd.GeoSeries([poly])

# 与坐标轴平行的外接矩形
env = gs.envelope
print(f"envelope 面积: {env.iloc[0].area:.2f}")

# 最小旋转外接矩形
mrr = gs.minimum_rotated_rectangle()
print(f"minimum_rotated_rectangle 面积: {mrr.iloc[0].area:.2f}")

# 最小旋转矩形的面积 <= envelope 的面积

13.14.3 实际应用

# 分析建筑物的朝向
buildings = gpd.read_file("buildings.shp")

# 获取每个建筑的最小旋转外接矩形
mrr = buildings.minimum_rotated_rectangle()

# 计算矩形的朝向（长轴方向）
import numpy as np

def get_orientation(rect):
    """计算矩形的主轴方向角度"""
    coords = list(rect.exterior.coords)
    edge1_len = Point(coords[0]).distance(Point(coords[1]))
    edge2_len = Point(coords[1]).distance(Point(coords[2]))

    if edge1_len > edge2_len:
        dx = coords[1][0] - coords[0][0]
        dy = coords[1][1] - coords[0][1]
    else:
        dx = coords[2][0] - coords[1][0]
        dy = coords[2][1] - coords[1][1]

    angle = np.degrees(np.arctan2(dy, dx))
    return angle % 180  # 0-180度

buildings['orientation'] = mrr.apply(get_orientation)

13.15 extract_unique_points() - 提取唯一顶点

13.15.1 基本概念

extract_unique_points() 从几何对象中提取所有唯一的顶点坐标，返回 MultiPoint。

13.15.2 代码示例

import geopandas as gpd
from shapely.geometry import Polygon, LineString

gs = gpd.GeoSeries([
    Polygon([(0, 0), (1, 0), (1, 1), (0, 1)]),   # 4 个顶点
    LineString([(0, 0), (1, 1), (2, 0)]),           # 3 个顶点
])

unique_points = gs.extract_unique_points()
for i, pts in enumerate(unique_points):
    print(f"几何 {i}: {len(pts.geoms)} 个唯一点")
    for p in pts.geoms:
        print(f"  ({p.x}, {p.y})")

13.15.3 注意事项

# 多边形的第一个和最后一个坐标相同（闭合环），
# 但 extract_unique_points 只返回唯一点
poly = Polygon([(0, 0), (1, 0), (1, 1), (0, 1)])
print(f"坐标数（含闭合点）: {len(poly.exterior.coords)}")    # 5
print(f"唯一点数: {len(gpd.GeoSeries([poly]).extract_unique_points().iloc[0].geoms)}")  # 4

13.15.4 实际应用

# 提取所有道路的交叉节点
roads = gpd.read_file("roads.shp")

# 提取所有道路的顶点
all_points = roads.extract_unique_points()

# 统计每个位置出现的次数
from collections import Counter

point_counter = Counter()
for mp in all_points:
    for p in mp.geoms:
        point_counter[(round(p.x, 6), round(p.y, 6))] += 1

# 出现多次的点可能是交叉节点
intersections = [Point(coord) for coord, count in point_counter.items() if count > 1]
print(f"发现 {len(intersections)} 个可能的道路交叉节点")

13.16 line_merge() - 线段合并

13.16.1 基本概念

line_merge() 将连续的线段（共享端点）合并为更长的线段。这对于处理分段存储的道路或河流数据非常有用。

13.16.2 代码示例

import geopandas as gpd
from shapely.geometry import MultiLineString, LineString

# 三段连续的线
multi_line = MultiLineString([
    [(0, 0), (1, 0)],
    [(1, 0), (2, 1)],
    [(2, 1), (3, 1)],
])

gs = gpd.GeoSeries([multi_line])
merged = gs.line_merge()

print(f"合并前类型: {multi_line.geom_type}")        # MultiLineString
print(f"合并前段数: {len(multi_line.geoms)}")         # 3
print(f"合并后类型: {merged.iloc[0].geom_type}")      # LineString
print(f"合并后坐标数: {len(merged.iloc[0].coords)}")  # 4

# 不连续的线不会被合并
disconnected = MultiLineString([
    [(0, 0), (1, 0)],
    [(2, 0), (3, 0)],   # 不与第一段相连
])
gs2 = gpd.GeoSeries([disconnected])
merged2 = gs2.line_merge()
print(f"不连续线合并: {merged2.iloc[0].geom_type}")  # MultiLineString

13.16.3 directed 参数

# directed=True 时只合并方向一致的线段
multi_line = MultiLineString([
    [(0, 0), (1, 0)],     # 向右
    [(2, 0), (1, 0)],     # 向左（方向相反但端点共享）
])

gs = gpd.GeoSeries([multi_line])

# 无方向合并
merged = gs.line_merge()
print(f"无方向合并: {merged.iloc[0].geom_type}")  # LineString

# 有方向合并
merged_directed = gs.line_merge(directed=True)
print(f"有方向合并: {merged_directed.iloc[0].geom_type}")  # 可能是 MultiLineString

13.17 polygonize() - 线转多边形

13.17.1 基本概念

polygonize() 将线段构成的网络转换为多边形。线段必须形成闭合环才能生成多边形。

13.17.2 代码示例

import geopandas as gpd
from shapely.geometry import MultiLineString

# 组成一个矩形的四条边
lines = MultiLineString([
    [(0, 0), (2, 0)],   # 下边
    [(2, 0), (2, 2)],   # 右边
    [(2, 2), (0, 2)],   # 上边
    [(0, 2), (0, 0)],   # 左边
])

gs = gpd.GeoSeries([lines])
polygons = gs.polygonize()

print(f"结果类型: {polygons.iloc[0].geom_type}")
# GeometryCollection 包含一个 Polygon
print(f"多边形面积: {list(polygons.iloc[0].geoms)[0].area}")  # 4.0

13.17.3 复杂线网

# 网格状线网会生成多个多边形
grid_lines = MultiLineString([
    # 水平线
    [(0, 0), (2, 0)], [(0, 1), (2, 1)], [(0, 2), (2, 2)],
    # 垂直线
    [(0, 0), (0, 2)], [(1, 0), (1, 2)], [(2, 0), (2, 2)],
])

gs = gpd.GeoSeries([grid_lines])
polygons = gs.polygonize()

polys = list(polygons.iloc[0].geoms)
print(f"生成了 {len(polys)} 个多边形")  # 4 个格子

for i, p in enumerate(polys):
    print(f"多边形 {i}: 面积={p.area}")

13.17.4 实际应用

# 从道路网络生成街区
roads = gpd.read_file("roads.shp")

# 合并所有道路线
all_roads = roads.union_all()

# 从道路网络生成街区多边形
gs = gpd.GeoSeries([all_roads])
blocks = gs.polygonize()

block_list = list(blocks.iloc[0].geoms)
blocks_gdf = gpd.GeoDataFrame(geometry=block_list, crs=roads.crs)
print(f"生成了 {len(blocks_gdf)} 个街区")

13.18 build_area() - 从线构建面

13.18.1 基本概念

build_area() 与 polygonize() 类似，但处理方式略有不同。它从线性几何构建面几何，可以处理嵌套的环（孔洞）。

13.18.2 代码示例

import geopandas as gpd
from shapely.geometry import MultiLineString

# 外环和内环
lines = MultiLineString([
    # 外环
    [(0, 0), (4, 0), (4, 4), (0, 4), (0, 0)],
    # 内环（孔）
    [(1, 1), (3, 1), (3, 3), (1, 3), (1, 1)],
])

gs = gpd.GeoSeries([lines])
area = gs.build_area()

print(f"结果类型: {area.iloc[0].geom_type}")  # Polygon（带孔）
print(f"面积: {area.iloc[0].area}")  # 16 - 4 = 12

# 检查是否有孔洞
print(f"孔洞数: {len(list(area.iloc[0].interiors))}")  # 1

13.18.3 build_area vs polygonize

特性 build_area polygonize
处理孔洞 自动识别内外环 内外环单独成为多边形
返回类型 单个几何（可含孔洞） GeometryCollection
适用场景 有嵌套环的场景 一般线网转面

13.19 normalize() - 几何标准化

13.19.1 基本概念

normalize() 将几何对象转换为标准化形式。标准化包括：

  • 坐标环的方向统一（外环逆时针，内环顺时针）
  • 起始点的选择统一
  • 多几何对象中子几何的排序

13.19.2 代码示例

import geopandas as gpd
from shapely.geometry import Polygon

# 两个相同的多边形但坐标顺序不同
poly_a = Polygon([(0, 0), (1, 0), (1, 1), (0, 1)])
poly_b = Polygon([(1, 1), (0, 1), (0, 0), (1, 0)])

gs = gpd.GeoSeries([poly_a, poly_b])
normalized = gs.normalize()

print(f"标准化前:")
print(f"  A: {list(poly_a.exterior.coords)}")
print(f"  B: {list(poly_b.exterior.coords)}")

print(f"标准化后:")
print(f"  A: {list(normalized.iloc[0].exterior.coords)}")
print(f"  B: {list(normalized.iloc[1].exterior.coords)}")
# 标准化后两者的坐标表示应该相同

13.19.3 实际应用

# 用于比较两个几何对象是否相同
# 标准化后可以使用字符串比较
normalized = gs.normalize()
wkt_a = normalized.iloc[0].wkt
wkt_b = normalized.iloc[1].wkt
print(f"WKT 相等: {wkt_a == wkt_b}")

13.20 本章小结

本章全面介绍了 GeoPandas 中的构造性几何操作。主要内容回顾：

特征提取与代表点

操作 说明
boundary 提取几何边界
exterior 提取多边形外环
interiors 提取多边形内环
centroid 计算质心
representative_point() 保证在几何内部的代表点

包络与包围几何

操作 说明
envelope 最小外接矩形（轴对齐）
convex_hull 凸包
concave_hull() 凹包
minimum_bounding_circle() 最小外接圆
maximum_inscribed_circle() 最大内切圆
minimum_rotated_rectangle() 最小旋转外接矩形

空间剖分

操作 说明
delaunay_triangles() Delaunay 三角剖分
constrained_delaunay_triangles() 约束 Delaunay 三角剖分
voronoi_polygons() Voronoi 图/泰森多边形

几何变换

操作 说明
extract_unique_points() 提取唯一顶点
line_merge() 合并连续线段
polygonize() 线网转多边形
build_area() 从线构建面
normalize() 几何标准化

使用建议

  1. 质心 vs 代表点：需要标注位置时优先使用 representative_point()
  2. 凸包 vs 凹包：凸包用于范围估计，凹包用于更精确的边界
  3. 三角剖分选择：一般用 delaunay_triangles()，有边约束时用约束版本
  4. 线段合并：处理分段数据时先 line_merge() 再分析
  5. 性能考虑：minimum_bounding_circle() 和 maximum_inscribed_circle() 可能较慢

下一章我们将学习缓冲区与简化操作，探索空间分析中最常用的几何变换方法。

↑