第11章拓扑关系与空间谓词是什么？

第11章：拓扑关系与空间谓词 空间拓扑关系是地理信息系统（GIS）中最核心的概念之一。它描述了两个几何对象在空间中的相对位置关系，是空间查询、空间分析和空间推理的理论基础。GeoPandas 基于 Shapely 库提供了丰富的拓扑关系判断方法，本章将全面介绍这些空间谓词的原理与用法。 11.1 空间关系概述 - DE-9IM 模型基础 11.1.1 什么是拓扑关系 拓扑关系（Topological Relationship）描述的是两个几何对象之间在空间上的定性关系，例如"是否相交"、"是否包含"、"是否接触"等。与度量关系（如距离、面积）不同，拓扑关系不依赖于坐标系的选择，具有在连续变换下保持不变的特性。 在 GIS 中，拓扑关系被广泛用于： 空间查询：查找与某个区域相交的所有要素 空间连接：基于空间关系将两个数据集关联起来 空间分析：确定要素之间的空间关联 数据验证：检查几何数据的拓扑一致性 11.1.2 DE-9IM 模型 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） * 任意值（不关心） 11.1.3 DE-9IM 矩阵示例 以两个重叠的多边形为例： 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维） 以此类推... 11.1.4 GeoPandas 中的空间谓词 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 11.2 intersects - 边界或内部相交判断 11.2.1 基本概念 intersects() 是最常用的空间谓词之一，用于判断两个几何对象是否在空间上有任何交集。只要两个几何对象的内部或边界存在公共部分，intersects() 就返回 True。 intersects() 实际上是 disjoint() 的逆操作： A.intersects(B) == NOT A.disjoint(B) 11.2.2 GeoSeries 与单个几何对象 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 11.2.3 GeoSeries 与 GeoSeries 逐元素比较 # 创建两个 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 11.2.4 实际应用：查找与区域相交的要素 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)}") 11.2.5 注意事项 intersects() 是一个非常宽松的条件，只要有任何接触就返回 True 对于大数据集，建议配合空间索引使用以提高性能 如果需要更精确的关系判断，可以使用 contains()、within() 等更严格的谓词 11.3 disjoint - 完全分离判断 11.3.1 基本概念 disjoint() 用于判断两个几何对象是否完全没有公共部分，即两个几何对象在空间上完全分离。它是 intersects() 的逆操作。 A.disjoint(B) == NOT A.intersects(B) DE-9IM 模式为 FF*FF****，表示 A 和 B 的内部不相交、A 和 B 的边界不相交、A 的内部与 B 的边界不相交、A 的边界与 B 的内部不相交。 11.3.2 代码示例 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 - 接触边界，不分离 11.3.3 实际应用：排除指定区域内的要素 # 查找不在禁区内的设施 facilities = gpd.read_file("facilities.shp") restricted_area = Polygon([(...)]) # 禁区多边形 # 筛选完全在禁区外的设施 safe_facilities = facilities[facilities.disjoint(restricted_area)] 11.4 touches - 边界接触（内部不相交） 11.4.1 基本概念 touches() 判断两个几何对象是否仅在边界上接触，而内部没有公共部分。这是一种比 intersects() 更精确的关系。 满足 touches() 的条件是： 两个几何对象的内部没有交集（I(A)∩I(B) = ∅） 但边界之间、或一个的内部与另一个的边界存在交集 11.4.2 代码示例 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 - 点在内部 11.4.3 touches 的几何类型组合 几何类型组合 touches 的含义 多边形-多边形 共享边或共享顶点，但内部不重叠 多边形-线 线的端点在多边形边界上，线不穿过内部 多边形-点 点在多边形边界上 线-线 线段端点相触，但不交叉 线-点 点是线段的端点 11.4.4 实际应用：查找相邻的行政区 # 查找与某省相邻的所有省份 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()}") 11.5 contains 与 within - 包含与被包含关系 11.5.1 contains 基本概念 contains() 判断一个几何对象是否完全包含另一个几何对象。如果 A 包含 B，则 B 的所有点都在 A 的内部或边界上，且 B 至少有一个点在 A 的内部。 DE-9IM 模式为 T*****FF*。 11.5.2 within 基本概念 within() 是 contains() 的逆操作。如果 B 在 A 内，则等价于 A 包含 B： A.contains(B) == B.within(A) 11.5.3 代码示例 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 11.5.4 边界点的特殊情况 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()。 11.5.5 实际应用：点在多边形查询 # 判断哪些点落在哪个行政区内 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') 11.6 contains_properly - 不含边界的严格包含 11.6.1 基本概念 contains_properly() 是比 contains() 更严格的包含关系。它要求被包含的几何对象的所有点都在包含几何对象的内部（不包括边界）。 A.contains_properly(B) 意味着: - B 的所有点都在 A 的内部 - B 与 A 的边界没有交集 11.6.2 contains 与 contains_properly 的区别 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 11.6.3 总结对比 情况 contains() contains_properly() 点在内部 True True 点在边界上 True False 几何完全在内部 True True 几何接触边界 True False 几何部分在外部 False False 11.6.4 实际应用 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)] 11.7 crosses - 交叉关系（维度低于最大值） 11.7.1 基本概念 crosses() 判断两个几何对象是否交叉。交叉关系的特点是： 两个几何对象的内部存在交集 但交集的维度低于两个几何对象中维度较大者的维度 简单来说，crosses() 通常用于： 线与线的交叉（交于点，不重叠） 线与面的交叉（线穿过面） 11.7.2 线与线的交叉 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（重叠不算交叉） 11.7.3 线与面的交叉 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（完全包含不算交叉） 11.7.4 crosses 与其他谓词的关系 如果 A.crosses(B) 为 True，则 A.intersects(B) 一定为 True 如果 A.crosses(B) 为 True，则 A.contains(B) 和 A.within(B) 一定为 False crosses() 是对称的：A.crosses(B) == B.crosses(A)（对于相同维度的几何对象） 11.8 overlaps - 重叠关系 11.8.1 基本概念 overlaps() 判断两个几何对象是否部分重叠。重叠关系的特点是： 两个几何对象具有相同的维度 交集的维度等于两个几何对象的维度 交集既不等于 A 也不等于 B（即两者都有不在交集中的部分） 11.8.2 代码示例 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 11.8.3 不同维度的几何对象 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 11.8.4 overlaps vs intersects 关系 说明 overlaps 两个同维度几何对象部分重叠（有公共部分，但双方都有独立部分） intersects 两个几何对象有任何公共部分（包括接触、重叠、包含等） 11.9 covers 与 covered_by - 覆盖关系 11.9.1 covers 基本概念 covers() 判断几何对象 A 是否覆盖几何对象 B。覆盖关系与 contains() 类似，但处理边界的方式不同： A.covers(B)：B 的每个点都在 A 的内部或边界上 A.contains(B)：B 的内部与 A 的外部不相交，且 B 与 A 的内部相交 在大多数情况下，covers() 和 contains() 的结果相同，但在某些边界情况下会有差异。 11.9.2 covered_by 基本概念 covered_by() 是 covers() 的逆操作： A.covers(B) == B.covered_by(A) 11.9.3 代码示例 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 11.9.4 covers vs contains 在某些边缘情况下，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() 在处理低维几何对象（如点在线段边界上）时行为更加直观。 11.10 dwithin - 距离内判断 11.10.1 基本概念 dwithin() 判断两个几何对象之间的距离是否在指定阈值之内。这是一种基于距离的空间谓词。 A.dwithin(B, distance) == (A.distance(B) <= distance) 但 dwithin() 通常比先计算距离再比较更高效，因为它可以利用空间索引进行优化。 11.10.2 代码示例 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 11.10.3 使用场景 # 查找某个位置 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)} 所") 11.10.4 注意事项 dwithin() 的距离单位取决于坐标参考系（CRS） 对于地理坐标系（如 EPSG:4326），距离单位为度，不是米 建议在使用 dwithin() 前先投影到合适的投影坐标系 11.11 距离计算 - distance()、hausdorff_distance()、frechet_distance() 11.11.1 distance() - 最短距离 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 11.11.2 hausdorff_distance() - Hausdorff 距离 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 距离常用于： 形状匹配 轨迹相似度比较 几何简化的质量评估 11.11.3 frechet_distance() - Fréchet 距离 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 距离不考虑顺序。 11.11.4 距离计算的注意事项 注意事项 说明 坐标系 距离单位取决于 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} 米") 11.12 几何相等判断 - geom_equals()、geom_equals_exact()、geom_equals_identical() 11.12.1 geom_equals() - 几何相等 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 11.12.2 geom_equals_exact() - 近似相等 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 11.12.3 geom_equals_identical() - 完全一致 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 11.12.4 三种相等判断的比较 方法 坐标顺序 浮点容差 起始点 严格程度 geom_equals() 不要求相同 无容差 不要求相同 中 geom_equals_exact() 要求相同 可指定容差 要求相同 低-高（取决于容差） geom_equals_identical() 要求完全相同 无容差 要求完全相同 最高 11.13 relate 与 DE-9IM - relate()、relate_pattern() 详解 11.13.1 relate() - 获取 DE-9IM 矩阵 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 11.13.2 DE-9IM 矩阵的解读 以内部点的 0F2FF1FF2 为例： poly内部 poly边界 poly外部 point内部 [ 0 ] [ F ] [ F ] point边界 [ F ] [ F ] [ F ] point外部 [ 2 ] [ 1 ] [ 2 ] 解读： 点的内部与多边形内部交集为 0 维（点本身）→ 点在多边形内部 点的外部与多边形内部交集为 2 维（面）→ 多边形还有大部分不与点相交 点没有边界（0维几何），所以边界行全为 F 11.13.3 relate_pattern() - 模式匹配 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 11.13.4 常用 DE-9IM 模式 空间关系 DE-9IM 模式 相等（Equals） TF**FFF 分离（Disjoint） FFFF*** 接触（Touches） FT******* 或 FT*** 或 FT* 包含（Contains） T*****FF* 被包含（Within） TFF 交叉-线/面（Crosses） TT***** 重叠-面/面（Overlaps） TTT 11.13.5 自定义空间关系 通过 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 11.14 二元谓词的 align 参数 11.14.1 align 参数的作用 GeoPandas 中的二元空间谓词方法（如 intersects()、contains() 等）在两个 GeoSeries 之间操作时，支持 align 参数。该参数控制是否在比较前对齐两个 GeoSeries 的索引。 # align=True（默认）：根据索引对齐后逐元素比较 # align=False：按位置逐元素比较，忽略索引 11.14.2 align=True（默认行为） 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)) 11.14.3 align=False # 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)) 11.14.4 使用建议 当两个 GeoSeries 的索引有意义且需要匹配时，使用 align=True 当两个 GeoSeries 按位置对应时，使用 align=False 对于与单个几何对象的比较，align 参数不起作用 11.15 实际应用案例 11.15.1 案例一：城市服务设施覆盖分析 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']} 未被任何医院覆盖！") 11.15.2 案例二：河流与道路交叉分析 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("这些交叉点可能需要建设桥梁") 11.15.3 案例三：多重空间条件筛选 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 11.15.4 案例四：空间关系矩阵 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) 11.16 本章小结 本章全面介绍了 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() 自定义 下一章我们将学习集合运算与几何操作，探索如何对几何对象进行布尔运算。