空间拓扑关系是地理信息系统（GIS）中最核心的概念之一。它描述了两个几何对象在空间中的相对位置关系，是空间查询、空间分析和空间推理的理论基础。GeoPandas 基于 Shapely 库提供了丰富的拓扑关系判断方法，本章将全面介绍这些空间谓词的原理与用法。
拓扑关系（Topological Relationship）描述的是两个几何对象之间在空间上的定性关系，例如"是否相交"、"是否包含"、"是否接触"等。与度量关系（如距离、面积）不同，拓扑关系不依赖于坐标系的选择，具有在连续变换下保持不变的特性。
在 GIS 中，拓扑关系被广泛用于：
DE-9IM（Dimensionally Extended Nine-Intersection Model，维度扩展九交模型）是描述两个几何对象拓扑关系的数学框架。该模型由 Clementini 和 Egenhofer 在 1991 年提出，是 OGC（开放地理空间联盟）标准的基础。
对于任意两个几何对象 A 和 B，每个几何对象都有三个组成部分：
| 组成部分 | 符号 | 说明 |
|---|---|---|
| 内部（Interior） | I(A), I(B) | 几何对象的内部区域 |
| 边界（Boundary） | B(A), B(B) | 几何对象的边界 |
| 外部（Exterior） | E(A), E(B) | 几何对象之外的空间 |
DE-9IM 矩阵通过计算这三个部分两两之间的交集维度，构建一个 3×3 的矩阵：
B的内部 B的边界 B的外部
A的内部 [ dim(I(A)∩I(B)) dim(I(A)∩B(B)) dim(I(A)∩E(B)) ]
A的边界 [ dim(B(A)∩I(B)) dim(B(A)∩B(B)) dim(B(A)∩E(B)) ]
A的外部 [ dim(E(A)∩I(B)) dim(E(A)∩B(B)) dim(E(A)∩E(B)) ]
矩阵中每个元素的值可以是：
| 值 | 含义 |
|---|---|
| -1 或 F | 交集为空集 |
| 0 | 交集的维度为 0（点） |
| 1 | 交集的维度为 1（线） |
| 2 | 交集的维度为 2（面） |
| T | 交集非空（维度 ≥ 0） |
| * | 任意值（不关心） |
以两个重叠的多边形为例：
import geopandas as gpd
from shapely.geometry import Polygon
# 创建两个重叠的多边形
poly_a = Polygon([(0, 0), (2, 0), (2, 2), (0, 2)])
poly_b = Polygon([(1, 1), (3, 1), (3, 3), (1, 3)])
# 计算 DE-9IM 矩阵
matrix = poly_a.relate(poly_b)
print(f"DE-9IM 矩阵字符串: {matrix}")
# 输出: 212101212
这个 9 位字符串 212101212 按行排列表示：
2 1 2
1 0 1
2 1 2
解读：
I(A)∩I(B) = 2：两个多边形的内部相交，交集是一个面（2维）I(A)∩B(B) = 1：A 的内部与 B 的边界相交，交集是一条线（1维）GeoPandas 提供了以下空间谓词方法，它们都是基于 DE-9IM 模型定义的：
| 方法 | 说明 | DE-9IM 模式 |
|---|---|---|
intersects() |
几何对象是否相交 | T******** 或 T****** 或 T** 或 T |
disjoint() |
几何对象是否分离 | FFFF*** |
touches() |
几何对象是否接触 | FT******* 或 FT*** 或 FT* |
contains() |
是否包含另一个几何对象 | T*****FF* |
within() |
是否在另一个几何对象内 | TFF |
crosses() |
几何对象是否交叉 | 取决于几何类型 |
overlaps() |
几何对象是否重叠 | 取决于几何类型 |
covers() |
是否覆盖另一个几何对象 | T*****FF* 或 T**FF 或 TFF 或 **TFF |
covered_by() |
是否被另一个几何对象覆盖 | TFF 或 TFF 或 FTF 或 FTF |
intersects() 是最常用的空间谓词之一，用于判断两个几何对象是否在空间上有任何交集。只要两个几何对象的内部或边界存在公共部分，intersects() 就返回 True。
intersects() 实际上是 disjoint() 的逆操作：
A.intersects(B) == NOT A.disjoint(B)
import geopandas as gpd
from shapely.geometry import Point, Polygon, LineString
import pandas as pd
# 创建 GeoDataFrame
gdf = gpd.GeoDataFrame(geometry=[
Polygon([(0, 0), (2, 0), (2, 2), (0, 2)]),
Polygon([(3, 3), (5, 3), (5, 5), (3, 5)]),
Polygon([(1, 1), (3, 1), (3, 3), (1, 3)]),
])
# 创建查询几何
query_point = Point(1, 1)
# 判断每个几何对象是否与查询点相交
result = gdf.intersects(query_point)
print(result)
# 0 True
# 1 False
# 2 True
# dtype: bool
# 创建两个 GeoSeries
gs1 = gpd.GeoSeries([
Polygon([(0, 0), (1, 0), (1, 1), (0, 1)]),
Polygon([(2, 2), (3, 2), (3, 3), (2, 3)]),
Point(5, 5),
])
gs2 = gpd.GeoSeries([
Point(0.5, 0.5), # 在第一个多边形内部
Point(4, 4), # 不在第二个多边形内
LineString([(4, 4), (6, 6)]), # 与第三个点相交
])
# 逐元素判断
result = gs1.intersects(gs2)
print(result)
# 0 True
# 1 False
# 2 True
# dtype: bool
import geopandas as gpd
from shapely.geometry import box
# 读取数据
cities = gpd.read_file("cities.shp")
# 定义查询区域
query_area = box(100, 20, 120, 40) # 经度100-120, 纬度20-40
# 查找与查询区域相交的城市
intersecting_cities = cities[cities.intersects(query_area)]
print(f"与查询区域相交的城市数量: {len(intersecting_cities)}")
intersects() 是一个非常宽松的条件，只要有任何接触就返回 Truecontains()、within() 等更严格的谓词disjoint() 用于判断两个几何对象是否完全没有公共部分，即两个几何对象在空间上完全分离。它是 intersects() 的逆操作。
A.disjoint(B) == NOT A.intersects(B)
DE-9IM 模式为 FF*FF****，表示 A 和 B 的内部不相交、A 和 B 的边界不相交、A 的内部与 B 的边界不相交、A 的边界与 B 的内部不相交。
from shapely.geometry import Point, Polygon
# 创建几何对象
poly = Polygon([(0, 0), (2, 0), (2, 2), (0, 2)])
point_inside = Point(1, 1) # 在多边形内部
point_outside = Point(5, 5) # 在多边形外部
point_on_boundary = Point(2, 1) # 在多边形边界上
gdf = gpd.GeoDataFrame(geometry=[poly, poly, poly])
# 分别判断
print(gpd.GeoSeries([poly]).disjoint(point_inside)[0]) # False - 点在内部，不分离
print(gpd.GeoSeries([poly]).disjoint(point_outside)[0]) # True - 完全分离
print(gpd.GeoSeries([poly]).disjoint(point_on_boundary)[0]) # False - 接触边界，不分离
# 查找不在禁区内的设施
facilities = gpd.read_file("facilities.shp")
restricted_area = Polygon([(...)]) # 禁区多边形
# 筛选完全在禁区外的设施
safe_facilities = facilities[facilities.disjoint(restricted_area)]
touches() 判断两个几何对象是否仅在边界上接触，而内部没有公共部分。这是一种比 intersects() 更精确的关系。
满足 touches() 的条件是：
I(A)∩I(B) = ∅）from shapely.geometry import Polygon, Point, LineString
# 两个相邻但不重叠的多边形
poly_a = Polygon([(0, 0), (1, 0), (1, 1), (0, 1)])
poly_b = Polygon([(1, 0), (2, 0), (2, 1), (1, 1)]) # 共享边 (1,0)-(1,1)
gs = gpd.GeoSeries([poly_a])
print(gs.touches(poly_b)[0]) # True - 边界接触
# 两个重叠的多边形
poly_c = Polygon([(0.5, 0), (1.5, 0), (1.5, 1), (0.5, 1)])
print(gs.touches(poly_c)[0]) # False - 内部有交集
# 点在多边形边界上
point_on_edge = Point(1, 0.5)
print(gs.touches(point_on_edge)[0]) # True - 点在边界上
# 点在多边形内部
point_inside = Point(0.5, 0.5)
print(gs.touches(point_inside)[0]) # False - 点在内部
| 几何类型组合 | touches 的含义 |
|---|---|
| 多边形-多边形 | 共享边或共享顶点，但内部不重叠 |
| 多边形-线 | 线的端点在多边形边界上，线不穿过内部 |
| 多边形-点 | 点在多边形边界上 |
| 线-线 | 线段端点相触，但不交叉 |
| 线-点 | 点是线段的端点 |
# 查找与某省相邻的所有省份
provinces = gpd.read_file("provinces.shp")
target_province = provinces[provinces['name'] == '湖北省'].geometry.iloc[0]
# 使用 touches 找到所有相邻省份
adjacent = provinces[provinces.touches(target_province)]
print(f"与湖北省相邻的省份: {adjacent['name'].tolist()}")
contains() 判断一个几何对象是否完全包含另一个几何对象。如果 A 包含 B，则 B 的所有点都在 A 的内部或边界上，且 B 至少有一个点在 A 的内部。
DE-9IM 模式为 T*****FF*。
within() 是 contains() 的逆操作。如果 B 在 A 内，则等价于 A 包含 B：
A.contains(B) == B.within(A)
from shapely.geometry import Polygon, Point
# 创建几何对象
outer_poly = Polygon([(0, 0), (4, 0), (4, 4), (0, 4)])
inner_poly = Polygon([(1, 1), (3, 1), (3, 3), (1, 3)])
point_inside = Point(2, 2)
point_on_boundary = Point(0, 0)
point_outside = Point(5, 5)
gs = gpd.GeoSeries([outer_poly])
# contains 测试
print(f"外部多边形包含内部多边形: {gs.contains(inner_poly)[0]}") # True
print(f"外部多边形包含内部点: {gs.contains(point_inside)[0]}") # True
print(f"外部多边形包含边界点: {gs.contains(point_on_boundary)[0]}") # True（注意！）
print(f"外部多边形包含外部点: {gs.contains(point_outside)[0]}") # False
# within 测试
gs_inner = gpd.GeoSeries([inner_poly])
print(f"内部多边形在外部多边形内: {gs_inner.within(outer_poly)[0]}") # True
contains() 对于边界点的处理需要特别注意：
from shapely.geometry import Polygon, Point
poly = Polygon([(0, 0), (1, 0), (1, 1), (0, 1)])
# 边界上的点
point_on_edge = Point(0.5, 0) # 在底边中点
gs = gpd.GeoSeries([poly])
print(f"contains 边界点: {gs.contains(point_on_edge)[0]}") # True（边界也算包含）
如果需要排除边界上的点，应使用 contains_properly()。
# 判断哪些点落在哪个行政区内
points = gpd.read_file("sample_points.shp")
districts = gpd.read_file("districts.shp")
# 方法一：逐个判断
for idx, district in districts.iterrows():
mask = points.within(district.geometry)
points_in_district = points[mask]
print(f"{district['name']}: {len(points_in_district)} 个点")
# 方法二：使用空间连接（推荐，更高效）
result = gpd.sjoin(points, districts, predicate='within')
contains_properly() 是比 contains() 更严格的包含关系。它要求被包含的几何对象的所有点都在包含几何对象的内部（不包括边界）。
A.contains_properly(B) 意味着:
- B 的所有点都在 A 的内部
- B 与 A 的边界没有交集
from shapely.geometry import Polygon, Point, LineString
poly = Polygon([(0, 0), (2, 0), (2, 2), (0, 2)])
# 内部点
point_inside = Point(1, 1)
# 边界上的点
point_on_boundary = Point(0, 1)
# 接触边界的线
line_touching = LineString([(0.5, 0.5), (2, 1)]) # 一端在内部，一端在边界上
gs = gpd.GeoSeries([poly])
print(f"contains(内部点): {gs.contains(point_inside)[0]}") # True
print(f"contains_properly(内部点): {gs.contains_properly(point_inside)[0]}") # True
print(f"contains(边界点): {gs.contains(point_on_boundary)[0]}") # True
print(f"contains_properly(边界点): {gs.contains_properly(point_on_boundary)[0]}") # False
print(f"contains(接触边界的线): {gs.contains(line_touching)[0]}") # True
print(f"contains_properly(接触边界的线): {gs.contains_properly(line_touching)[0]}") # False
| 情况 | contains() | contains_properly() |
|---|---|---|
| 点在内部 | True | True |
| 点在边界上 | True | False |
| 几何完全在内部 | True | True |
| 几何接触边界 | True | False |
| 几何部分在外部 | False | False |
contains_properly() 常用于以下场景：
# 只选择严格位于保护区内部的监测站
stations = gpd.read_file("monitoring_stations.shp")
reserve = gpd.read_file("nature_reserve.shp").geometry.iloc[0]
# 严格在内部的站点
strict_inside = stations[gpd.GeoSeries([reserve] * len(stations), index=stations.index)
.contains_properly(stations.geometry)]
crosses() 判断两个几何对象是否交叉。交叉关系的特点是：
简单来说，crosses() 通常用于：
from shapely.geometry import LineString
# 十字交叉的两条线
line_a = LineString([(0, 0), (2, 2)])
line_b = LineString([(0, 2), (2, 0)])
gs = gpd.GeoSeries([line_a])
print(f"两条线交叉: {gs.crosses(line_b)[0]}") # True
# 不交叉的两条平行线
line_c = LineString([(0, 1), (2, 3)])
print(f"平行线交叉: {gs.crosses(line_c)[0]}") # False
# 部分重叠的线
line_d = LineString([(1, 1), (3, 3)])
print(f"重叠线交叉: {gs.crosses(line_d)[0]}") # False（重叠不算交叉）
from shapely.geometry import LineString, Polygon
poly = Polygon([(0, 0), (2, 0), (2, 2), (0, 2)])
# 穿过多边形的线
line_through = LineString([(-1, 1), (3, 1)]) # 线穿过多边形
gs = gpd.GeoSeries([line_through])
print(f"线穿过多边形: {gs.crosses(poly)[0]}") # True
# 完全在多边形内的线
line_inside = LineString([(0.5, 0.5), (1.5, 1.5)])
gs2 = gpd.GeoSeries([line_inside])
print(f"线在多边形内: {gs2.crosses(poly)[0]}") # False（完全包含不算交叉）
A.crosses(B) 为 True，则 A.intersects(B) 一定为 TrueA.crosses(B) 为 True，则 A.contains(B) 和 A.within(B) 一定为 Falsecrosses() 是对称的：A.crosses(B) == B.crosses(A)（对于相同维度的几何对象）overlaps() 判断两个几何对象是否部分重叠。重叠关系的特点是：
from shapely.geometry import Polygon, Point, LineString
# 部分重叠的两个多边形
poly_a = Polygon([(0, 0), (2, 0), (2, 2), (0, 2)])
poly_b = Polygon([(1, 1), (3, 1), (3, 3), (1, 3)])
gs = gpd.GeoSeries([poly_a])
print(f"部分重叠: {gs.overlaps(poly_b)[0]}") # True
# 完全包含的多边形（不算重叠）
poly_c = Polygon([(0.5, 0.5), (1.5, 0.5), (1.5, 1.5), (0.5, 1.5)])
print(f"完全包含: {gs.overlaps(poly_c)[0]}") # False
# 相同的多边形（不算重叠）
print(f"完全相同: {gs.overlaps(poly_a)[0]}") # False
# 仅接触的多边形（不算重叠）
poly_d = Polygon([(2, 0), (4, 0), (4, 2), (2, 2)])
print(f"仅接触: {gs.overlaps(poly_d)[0]}") # False
overlaps() 要求两个几何对象具有相同的维度。因此：
# 点与多边形不能 overlap
poly = Polygon([(0, 0), (2, 0), (2, 2), (0, 2)])
point = Point(1, 1)
gs = gpd.GeoSeries([poly])
print(f"点与多边形 overlaps: {gs.overlaps(point)[0]}") # False
# 两条部分重叠的线可以 overlap
line_a = LineString([(0, 0), (2, 0)])
line_b = LineString([(1, 0), (3, 0)])
gs2 = gpd.GeoSeries([line_a])
print(f"两条线 overlaps: {gs2.overlaps(line_b)[0]}") # True
| 关系 | 说明 |
|---|---|
overlaps |
两个同维度几何对象部分重叠（有公共部分，但双方都有独立部分） |
intersects |
两个几何对象有任何公共部分（包括接触、重叠、包含等） |
covers() 判断几何对象 A 是否覆盖几何对象 B。覆盖关系与 contains() 类似，但处理边界的方式不同：
A.covers(B)：B 的每个点都在 A 的内部或边界上A.contains(B)：B 的内部与 A 的外部不相交，且 B 与 A 的内部相交在大多数情况下，covers() 和 contains() 的结果相同，但在某些边界情况下会有差异。
covered_by() 是 covers() 的逆操作：
A.covers(B) == B.covered_by(A)
from shapely.geometry import Polygon, Point, LineString
poly = Polygon([(0, 0), (2, 0), (2, 2), (0, 2)])
# 内部点
point_inside = Point(1, 1)
# 边界点
point_on_boundary = Point(0, 1)
# 外部点
point_outside = Point(3, 3)
gs = gpd.GeoSeries([poly])
# covers 测试
print(f"covers 内部点: {gs.covers(point_inside)[0]}") # True
print(f"covers 边界点: {gs.covers(point_on_boundary)[0]}") # True
print(f"covers 外部点: {gs.covers(point_outside)[0]}") # False
# covered_by 测试
gs_point = gpd.GeoSeries([point_inside])
print(f"内部点 covered_by 多边形: {gs_point.covered_by(poly)[0]}") # True
在某些边缘情况下，covers 和 contains 的结果会不同：
from shapely.geometry import LineString, Point
# 线段
line = LineString([(0, 0), (1, 0)])
# 线段的端点
endpoint = Point(0, 0)
gs_line = gpd.GeoSeries([line])
# 线段 covers 其端点
print(f"线 covers 端点: {gs_line.covers(endpoint)[0]}") # True
# 线段 contains 其端点
print(f"线 contains 端点: {gs_line.contains(endpoint)[0]}") # True（注意：Shapely 2.0 行为）
一般来说，covers() 在处理低维几何对象（如点在线段边界上）时行为更加直观。
dwithin() 判断两个几何对象之间的距离是否在指定阈值之内。这是一种基于距离的空间谓词。
A.dwithin(B, distance) == (A.distance(B) <= distance)
但 dwithin() 通常比先计算距离再比较更高效，因为它可以利用空间索引进行优化。
from shapely.geometry import Point, Polygon
import geopandas as gpd
# 创建 GeoSeries
gs = gpd.GeoSeries([
Point(0, 0),
Point(3, 0),
Point(5, 0),
Point(10, 0),
])
# 查询点
query = Point(0, 0)
# 判断哪些点在距离 4 以内
result = gs.dwithin(query, 4)
print(result)
# 0 True (距离 0)
# 1 True (距离 3)
# 2 False (距离 5)
# 3 False (距离 10)
# dtype: bool
# 查找某个位置 5 公里范围内的所有学校
schools = gpd.read_file("schools.shp")
# 注意：确保使用投影坐标系（单位为米）
schools_proj = schools.to_crs(epsg=32650) # UTM Zone 50N
target = Point(500000, 3500000) # 投影坐标
nearby = schools_proj[schools_proj.dwithin(target, 5000)] # 5000 米
print(f"5公里范围内的学校: {len(nearby)} 所")
dwithin() 的距离单位取决于坐标参考系（CRS）dwithin() 前先投影到合适的投影坐标系distance() 计算两个几何对象之间的最短距离（欧几里得距离）。
from shapely.geometry import Point, Polygon, LineString
import geopandas as gpd
# 点到点的距离
gs = gpd.GeoSeries([Point(0, 0), Point(3, 0)])
other = Point(0, 4)
print(gs.distance(other))
# 0 4.0
# 1 5.0
# dtype: float64
# 点到多边形的距离
poly = Polygon([(2, 2), (4, 2), (4, 4), (2, 4)])
gs2 = gpd.GeoSeries([Point(0, 0)])
print(gs2.distance(poly))
# 0 2.828427 (约 2√2)
# dtype: float64
# 点在多边形内部，距离为 0
gs3 = gpd.GeoSeries([Point(3, 3)])
print(gs3.distance(poly))
# 0 0.0
# dtype: float64
Hausdorff 距离衡量两个几何对象之间的最大最近距离。直观地说，它是两个几何对象之间"最差匹配"的度量。
from shapely.geometry import LineString
import geopandas as gpd
# 两条相似但不完全相同的线
line_a = LineString([(0, 0), (1, 0), (2, 0)])
line_b = LineString([(0, 0.5), (1, 0.3), (2, 0.5)])
gs = gpd.GeoSeries([line_a])
print(f"Hausdorff 距离: {gs.hausdorff_distance(line_b)[0]:.4f}")
# 输出: Hausdorff 距离: 0.5000
Hausdorff 距离常用于：
Fréchet 距离（也称为"遛狗距离"）考虑了沿曲线行进的顺序，更适合比较有方向的曲线。
from shapely.geometry import LineString
import geopandas as gpd
# 两条路径
path_a = LineString([(0, 0), (1, 1), (2, 0)])
path_b = LineString([(0, 0), (1, 2), (2, 0)])
gs = gpd.GeoSeries([path_a])
print(f"Fréchet 距离: {gs.frechet_distance(path_b)[0]:.4f}")
Fréchet 距离与 Hausdorff 距离的区别在于，Fréchet 距离考虑了曲线上点的顺序关系，而 Hausdorff 距离不考虑顺序。
| 注意事项 | 说明 |
|---|---|
| 坐标系 | 距离单位取决于 CRS，地理坐标系下为度 |
| 投影选择 | 建议使用等距投影进行距离计算 |
| 性能 | 大规模距离矩阵计算可使用空间索引优化 |
| 精度 | 对于大范围地理距离，考虑使用测地线距离 |
# 正确的距离计算流程
import geopandas as gpd
# 1. 读取数据（通常为 WGS84）
points = gpd.read_file("points.shp")
# 2. 投影到合适的坐标系
points_proj = points.to_crs(epsg=32650)
# 3. 计算距离（单位为米）
distances = points_proj.distance(Point(500000, 3500000))
print(f"平均距离: {distances.mean():.2f} 米")
geom_equals() 判断两个几何对象在拓扑意义上是否相等，即它们表示相同的点集。坐标顺序和起始点可以不同。
from shapely.geometry import Polygon
import geopandas as gpd
# 相同的多边形，但顶点起始位置不同
poly_a = Polygon([(0, 0), (1, 0), (1, 1), (0, 1)])
poly_b = Polygon([(1, 0), (1, 1), (0, 1), (0, 0)]) # 不同的起始点
gs = gpd.GeoSeries([poly_a])
print(f"拓扑相等: {gs.geom_equals(poly_b)[0]}") # True
geom_equals_exact() 判断两个几何对象在给定容差范围内是否相等。这对于处理浮点数精度问题非常有用。
from shapely.geometry import Point
import geopandas as gpd
point_a = Point(0, 0)
point_b = Point(0.0000001, 0.0000001) # 微小差异
gs = gpd.GeoSeries([point_a])
# 精确比较
print(f"精确相等: {gs.geom_equals(point_b)[0]}") # False
# 带容差比较
print(f"近似相等(容差1e-6): {gs.geom_equals_exact(point_b, tolerance=1e-6)[0]}") # True
print(f"近似相等(容差1e-8): {gs.geom_equals_exact(point_b, tolerance=1e-8)[0]}") # False
geom_equals_identical() 是最严格的相等判断，要求两个几何对象在所有方面完全一致，包括：
from shapely.geometry import Polygon
import geopandas as gpd
poly_a = Polygon([(0, 0), (1, 0), (1, 1), (0, 1)])
poly_b = Polygon([(0, 0), (1, 0), (1, 1), (0, 1)]) # 完全相同
poly_c = Polygon([(1, 0), (1, 1), (0, 1), (0, 0)]) # 起始点不同
gs = gpd.GeoSeries([poly_a])
print(f"identical(完全相同): {gs.geom_equals_identical(poly_b)[0]}") # True
print(f"identical(起始点不同): {gs.geom_equals_identical(poly_c)[0]}") # False
# 对比 geom_equals
print(f"geom_equals(起始点不同): {gs.geom_equals(poly_c)[0]}") # True
| 方法 | 坐标顺序 | 浮点容差 | 起始点 | 严格程度 |
|---|---|---|---|---|
geom_equals() |
不要求相同 | 无容差 | 不要求相同 | 中 |
geom_equals_exact() |
要求相同 | 可指定容差 | 要求相同 | 低-高（取决于容差） |
geom_equals_identical() |
要求完全相同 | 无容差 | 要求完全相同 | 最高 |
relate() 方法返回两个几何对象之间的 DE-9IM 矩阵字符串，这是最底层的拓扑关系描述。
from shapely.geometry import Point, Polygon, LineString
import geopandas as gpd
# 创建几何对象
poly = Polygon([(0, 0), (2, 0), (2, 2), (0, 2)])
# 内部点
point_inside = Point(1, 1)
gs = gpd.GeoSeries([poly])
matrix = gs.relate(point_inside)
print(f"内部点的 DE-9IM: {matrix[0]}") # 0F2FF1FF2
# 边界点
point_boundary = Point(0, 0)
matrix = gs.relate(point_boundary)
print(f"边界点的 DE-9IM: {matrix[0]}") # FF20F1FF2
# 外部点
point_outside = Point(5, 5)
matrix = gs.relate(point_outside)
print(f"外部点的 DE-9IM: {matrix[0]}") # FF2FF10F2
以内部点的 0F2FF1FF2 为例：
poly内部 poly边界 poly外部
point内部 [ 0 ] [ F ] [ F ]
point边界 [ F ] [ F ] [ F ]
point外部 [ 2 ] [ 1 ] [ 2 ]
解读：
relate_pattern() 允许使用 DE-9IM 模式字符串来自定义空间关系的判断。模式字符串中可以使用以下通配符：
T：交集非空（维度 ≥ 0）F：交集为空*：任意值0、1、2：指定维度from shapely.geometry import Polygon, Point
import geopandas as gpd
poly = Polygon([(0, 0), (2, 0), (2, 2), (0, 2)])
point = Point(1, 1)
gs = gpd.GeoSeries([poly])
# 使用 contains 的 DE-9IM 模式
contains_pattern = "T*****FF*"
print(f"matches contains: {gs.relate_pattern(point, contains_pattern)[0]}") # True
# 使用 within 的 DE-9IM 模式
within_pattern = "T*F**F***"
gs_point = gpd.GeoSeries([point])
print(f"matches within: {gs_point.relate_pattern(poly, within_pattern)[0]}") # True
# 自定义模式：内部与内部相交，维度为 2（面）
custom_pattern = "2********"
poly2 = Polygon([(1, 1), (3, 1), (3, 3), (1, 3)])
print(f"内部相交为面: {gs.relate_pattern(poly2, custom_pattern)[0]}") # True
| 空间关系 | DE-9IM 模式 |
|---|---|
| 相等（Equals） | TF**FFF |
| 分离（Disjoint） | FFFF*** |
| 接触（Touches） | FT******* 或 FT*** 或 FT* |
| 包含（Contains） | T*****FF* |
| 被包含（Within） | TFF |
| 交叉-线/面（Crosses） | TT***** |
| 重叠-面/面（Overlaps） | TTT |
通过 relate_pattern()，可以定义标准谓词无法表达的复杂空间关系：
# 例如：判断两个多边形是否共享一条边（而非仅一个点）
# 条件：边界的交集维度为 1（线）
shared_edge_pattern = "****1****"
poly_a = Polygon([(0, 0), (1, 0), (1, 1), (0, 1)])
poly_b = Polygon([(1, 0), (2, 0), (2, 1), (1, 1)]) # 共享边
gs = gpd.GeoSeries([poly_a])
print(f"共享边: {gs.relate_pattern(poly_b, shared_edge_pattern)[0]}") # True
# 仅共享一个点的情况
poly_c = Polygon([(1, 1), (2, 1), (2, 2), (1, 2)]) # 仅共享点(1,1)
print(f"仅共享点: {gs.relate_pattern(poly_c, shared_edge_pattern)[0]}") # False
GeoPandas 中的二元空间谓词方法（如 intersects()、contains() 等）在两个 GeoSeries 之间操作时，支持 align 参数。该参数控制是否在比较前对齐两个 GeoSeries 的索引。
# align=True（默认）：根据索引对齐后逐元素比较
# align=False：按位置逐元素比较，忽略索引
import geopandas as gpd
from shapely.geometry import Point, Polygon
# 创建两个索引不同的 GeoSeries
gs1 = gpd.GeoSeries(
[Polygon([(0,0),(1,0),(1,1),(0,1)]), Polygon([(2,2),(3,2),(3,3),(2,3)])],
index=[0, 1]
)
gs2 = gpd.GeoSeries(
[Point(2.5, 2.5), Point(0.5, 0.5)],
index=[1, 0]
)
# align=True：按索引对齐
result = gs1.intersects(gs2, align=True)
print(result)
# 0 True (gs1[0] vs gs2[0]，即第一个多边形 vs Point(0.5,0.5))
# 1 True (gs1[1] vs gs2[1]，即第二个多边形 vs Point(2.5,2.5))
# align=False：按位置对齐
result = gs1.intersects(gs2, align=False)
print(result)
# 0 False (gs1位置0 vs gs2位置0，即第一个多边形 vs Point(2.5,2.5))
# 1 False (gs1位置1 vs gs2位置1，即第二个多边形 vs Point(0.5,0.5))
align=Truealign=Falsealign 参数不起作用import geopandas as gpd
from shapely.geometry import Point
import pandas as pd
# 模拟城市数据
# 创建医院位置
hospitals = gpd.GeoDataFrame({
'name': ['第一医院', '第二医院', '第三医院', '第四医院'],
'geometry': [
Point(116.40, 39.90),
Point(116.45, 39.92),
Point(116.38, 39.88),
Point(116.42, 39.95),
]
}, crs='EPSG:4326')
# 创建居民区位置
residential = gpd.GeoDataFrame({
'name': ['小区A', '小区B', '小区C', '小区D', '小区E'],
'geometry': [
Point(116.41, 39.91),
Point(116.46, 39.93),
Point(116.39, 39.87),
Point(116.50, 39.98),
Point(116.35, 39.85),
]
}, crs='EPSG:4326')
# 转换到投影坐标系（米为单位）
hospitals_proj = hospitals.to_crs(epsg=32650)
residential_proj = residential.to_crs(epsg=32650)
# 分析每个医院 3000 米服务范围
service_radius = 3000 # 米
hospital_buffers = hospitals_proj.copy()
hospital_buffers['geometry'] = hospitals_proj.buffer(service_radius)
# 判断每个小区是否在某个医院的服务范围内
for idx, community in residential_proj.iterrows():
covered = hospital_buffers[hospital_buffers.contains(community.geometry)]
if len(covered) > 0:
print(f"{residential.iloc[idx]['name']} 被 {covered['name'].tolist()} 覆盖")
else:
print(f"{residential.iloc[idx]['name']} 未被任何医院覆盖！")
import geopandas as gpd
# 读取数据
rivers = gpd.read_file("rivers.shp")
roads = gpd.read_file("roads.shp")
# 确保坐标系一致
roads = roads.to_crs(rivers.crs)
# 找到所有与河流交叉的道路
crossing_pairs = []
for r_idx, river in rivers.iterrows():
crossing_roads = roads[roads.crosses(river.geometry)]
for rd_idx, road in crossing_roads.iterrows():
crossing_pairs.append({
'river_name': river['name'],
'road_name': road['name'],
'intersection': river.geometry.intersection(road.geometry)
})
crossing_df = pd.DataFrame(crossing_pairs)
print(f"共发现 {len(crossing_df)} 个河流-道路交叉点")
print("这些交叉点可能需要建设桥梁")
import geopandas as gpd
from shapely.geometry import Point
# 综合使用多种空间谓词进行复杂筛选
# 场景：选择新建学校的位置
# 条件1：必须在规划区域内（within）
# 条件2：不能在已有学校 1000 米范围内（disjoint + buffer）
# 条件3：必须在主干道 500 米范围内（dwithin 或 intersects + buffer）
# 条件4：不能在工厂 200 米范围内（disjoint + buffer）
def find_suitable_sites(candidate_sites, planning_area, existing_schools,
main_roads, factories):
"""综合筛选适合建校的位置"""
# 条件1：在规划区域内
mask1 = candidate_sites.within(planning_area)
# 条件2：不在已有学校 1000 米范围内
school_buffers = existing_schools.buffer(1000).union_all()
mask2 = candidate_sites.disjoint(school_buffers)
# 条件3：在主干道 500 米范围内
road_buffers = main_roads.buffer(500).union_all()
mask3 = candidate_sites.intersects(road_buffers)
# 条件4：不在工厂 200 米范围内
factory_buffers = factories.buffer(200).union_all()
mask4 = candidate_sites.disjoint(factory_buffers)
# 综合所有条件
suitable = candidate_sites[mask1 & mask2 & mask3 & mask4]
return suitable
import geopandas as gpd
import pandas as pd
# 计算所有区域两两之间的空间关系
regions = gpd.read_file("regions.shp")
n = len(regions)
relation_matrix = pd.DataFrame(index=regions['name'], columns=regions['name'])
for i in range(n):
for j in range(n):
if i == j:
relation_matrix.iloc[i, j] = 'self'
else:
geom_i = regions.geometry.iloc[i]
geom_j = regions.geometry.iloc[j]
if geom_i.contains(geom_j):
relation_matrix.iloc[i, j] = 'contains'
elif geom_i.within(geom_j):
relation_matrix.iloc[i, j] = 'within'
elif geom_i.overlaps(geom_j):
relation_matrix.iloc[i, j] = 'overlaps'
elif geom_i.touches(geom_j):
relation_matrix.iloc[i, j] = 'touches'
elif geom_i.disjoint(geom_j):
relation_matrix.iloc[i, j] = 'disjoint'
else:
relation_matrix.iloc[i, j] = geom_i.relate(geom_j)
print(relation_matrix)
本章全面介绍了 GeoPandas 中的拓扑关系与空间谓词。主要内容回顾：
| 概念 | 说明 |
|---|---|
| DE-9IM 模型 | 描述两个几何对象拓扑关系的数学框架 |
| 空间谓词 | 基于 DE-9IM 定义的布尔空间关系判断 |
| 距离度量 | 欧几里得距离、Hausdorff 距离、Fréchet 距离 |
| 谓词 | 含义 | 对称性 |
|---|---|---|
intersects |
有任何公共部分 | 对称 |
disjoint |
完全分离 | 对称 |
touches |
仅边界接触 | 对称 |
contains |
完全包含 | 非对称（逆为 within） |
within |
完全被包含 | 非对称（逆为 contains） |
contains_properly |
严格包含（不含边界） | 非对称 |
crosses |
交叉 | 对称 |
overlaps |
部分重叠 | 对称 |
covers |
覆盖 | 非对称（逆为 covered_by） |
covered_by |
被覆盖 | 非对称（逆为 covers） |
dwithin |
距离范围内 | 对称 |
sjoin 等方法可自动利用空间索引contains 与 contains_properly 在边界处理上的差异relate_pattern() 自定义下一章我们将学习集合运算与几何操作，探索如何对几何对象进行布尔运算。
此内容由惯性聚合(RSS阅读器)自动聚合整理，仅供阅读参考。 原文来自 — 版权归原作者所有。