Shapely判断点在面内?几何关系咋算?

GIS基础理论
Dr.GIS
wowwwai GIS研习社 · 工具流程与项目排障

很多 GIS 初学者在用 Python 做空间分析时,都会遇到一个看似简单的问题:Shapely判断点在面内?几何关系咋算? 如果你已经有一个点坐标和一个行政区、多边形地块或缓冲区范围,想判断这个点是否落在面内,Shapely 是最常用、最轻量的 Python 几何计算库之一。

本文围绕一个具体任务展开:用 Shapely 判断点在面内,并顺带讲清楚常见的几何关系计算方法,包括 containswithinintersectstouches 等。重点不是背 API,而是知道什么时候该用哪个判断,以及边界点为什么经常让结果“看起来不对”。

引言:Shapely判断点在面内的典型场景

在 GIS 项目里,点面关系判断非常常见。例如:

  • 判断门店点是否位于某个商圈多边形内。
  • 判断采样点是否落在某个行政区范围内。
  • 判断 GPS 轨迹点是否进入电子围栏。
  • 判断事故点是否位于缓冲区分析结果内。
  • 批量给点数据匹配所属地块、街道或网格。

这些任务的核心都是一个问题:点和面之间是什么空间关系。Shapely 判断点在面内时,通常会用到 PointPolygon 以及 containswithin 方法。

Shapely判断点在面内与Shapely几何关系计算示意图
点在面内、点在边界上、点在面外时,Shapely 几何关系判断结果并不完全相同。

背景:为什么点在面内判断会出错

很多人第一次用 Shapely 判断点在面内,会写出类似代码:

from shapely.geometry import Point, Polygon

point = Point(1, 1)
polygon = Polygon([(0, 0), (3, 0), (3, 3), (0, 3)])

print(polygon.contains(point))

如果输出是 True,说明点在多边形内部。这个例子很顺利,但实际项目里经常会出现以下问题:

  • 点明明看起来在面内,结果返回 False
  • 点落在边界线上,contains 返回 False
  • 经纬度坐标和投影坐标混用,导致判断结果完全不可信。
  • 多边形无效,例如自相交,导致空间关系计算异常。
  • 坐标顺序写反,把经度纬度写成了纬度经度。

所以,Shapely 判断点在面内不是只记住一个函数名就够了,还要理解几何关系的定义、坐标系统一致性,以及边界点如何处理。

原理:Shapely几何关系咋算

Shapely 是基于 GEOS 几何引擎的 Python 库。GEOS 是很多开源 GIS 软件和数据库都在使用的几何计算核心,例如 PostGIS 也使用 GEOS 做大量空间关系判断。

在 Shapely 中,几何对象包括点、线、面等:

  • Point:点对象,例如一个 GPS 位置。
  • LineString:线对象,例如道路、河流或轨迹。
  • Polygon:面对象,例如行政区、地块、缓冲区。
  • MultiPolygon:多个面组成的集合,例如包含多个岛屿的行政区。

判断点在面内,本质上是在计算两个几何对象之间的拓扑关系。常用方法如下:

方法 含义 适合场景
polygon.contains(point) 面是否包含点,不包含边界点 严格判断点是否在面内部
point.within(polygon) 点是否位于面内部,不包含边界点 和 contains 方向相反,语义更符合点查询
polygon.intersects(point) 面是否与点有交集,包含内部和边界 点在面内或边界上都算命中
polygon.touches(point) 点是否接触面边界 专门识别边界点
polygon.covers(point) 面是否覆盖点,通常包含边界 业务上边界也算在区域内时更合适

最容易混淆的是 containsintersects。如果点在多边形边界上,contains 通常返回 False,但 intersects 会返回 True。这就是很多“点明明在线上,为什么不算在面内”的根源。

步骤:用Shapely判断点在面内

步骤1:安装 Shapely

如果你的 Python 环境还没有安装 Shapely,可以使用 pip 安装:

pip install shapely

如果你使用 Anaconda,也可以在对应环境中安装:

conda install shapely

建议在项目虚拟环境中安装,避免和 QGIS、ArcGIS Pro 自带 Python 环境混在一起。

步骤2:创建点和多边形

下面先用一个简单矩形演示 Shapely 判断点在面内:

from shapely.geometry import Point, Polygon

polygon = Polygon([
    (0, 0),
    (4, 0),
    (4, 3),
    (0, 3)
])

point_inside = Point(2, 1)
point_outside = Point(5, 1)
point_boundary = Point(4, 1)

print(polygon.contains(point_inside))
print(polygon.contains(point_outside))
print(polygon.contains(point_boundary))

输出结果通常是:

True
False
False

这里要注意,边界点 Point(4, 1) 不被 contains 认为是在面内部。

步骤3:用 within 改写判断逻辑

withincontains 是一组方向相反的关系。对于点在面内判断,很多人更喜欢写成:

print(point_inside.within(polygon))
print(point_outside.within(polygon))
print(point_boundary.within(polygon))

它的结果和 polygon.contains(point) 基本对应:

True
False
False

如果你写的是点查询,point.within(polygon) 可读性更强;如果你写的是区域筛选,polygon.contains(point) 也很自然。

步骤4:如果边界也算在区域内,用 covers 或 intersects

很多业务规则里,点落在行政区边界、地块边界、电子围栏边界上,也应该算作命中。这时不要只用 contains

print(polygon.intersects(point_boundary))
print(polygon.covers(point_boundary))

通常会得到:

True
True

如果你的业务含义是“点在面内或边界上都算属于该面”,优先考虑 coversintersects。如果只是点和面之间的关系,covers 的语义通常比 intersects 更精确。

步骤5:识别边界点

如果你需要单独判断点是否刚好落在面边界上,可以使用 touches

print(polygon.touches(point_inside))
print(polygon.touches(point_boundary))

输出通常是:

False
True

在数据质检中,这个判断很有用。例如你要检查采样点是否落在网格边界上,避免后续归属到多个网格。

步骤6:封装一个实用函数

实际项目中,建议把点在面内判断封装成函数,并明确是否包含边界:

from shapely.geometry import Point, Polygon

def point_in_polygon(x, y, polygon, include_boundary=True):
    point = Point(x, y)

    if include_boundary:
        return polygon.covers(point)

    return polygon.contains(point)


polygon = Polygon([
    (0, 0),
    (4, 0),
    (4, 3),
    (0, 3)
])

print(point_in_polygon(2, 1, polygon, include_boundary=True))
print(point_in_polygon(4, 1, polygon, include_boundary=True))
print(point_in_polygon(4, 1, polygon, include_boundary=False))

这个函数的好处是业务规则清楚:include_boundary=True 表示边界算在面内,include_boundary=False 表示只算严格内部。

步骤:从GeoJSON或WKT读取几何后判断

从 WKT 读取点和面

很多 GIS 数据库或接口会返回 WKT。WKT 是 Well-Known Text 的缩写,是一种用文本表达几何对象的格式。

from shapely import wkt

polygon_wkt = "POLYGON ((0 0, 4 0, 4 3, 0 3, 0 0))"
point_wkt = "POINT (2 1)"

polygon = wkt.loads(polygon_wkt)
point = wkt.loads(point_wkt)

print(polygon.contains(point))

如果你从 PostGIS、CSV 字段或接口响应中拿到 WKT,可以直接用这种方式转换成 Shapely 几何对象。

从 GeoJSON 读取几何

如果数据来自 WebGIS 或前端地图,常见格式是 GeoJSON。Shapely 可以通过 shape 把 GeoJSON 几何字典转成对象。

from shapely.geometry import shape, Point

geojson_polygon = {
    "type": "Polygon",
    "coordinates": [[
        [0, 0],
        [4, 0],
        [4, 3],
        [0, 3],
        [0, 0]
    ]]
}

polygon = shape(geojson_polygon)
point = Point(2, 1)

print(polygon.covers(point))

这里要特别注意 GeoJSON 坐标顺序是 [经度, 纬度],也就是 [x, y]。不要写成 [纬度, 经度]

常见坑:Shapely判断点在面内结果不对怎么办

坑1:边界点不被 contains 认为在面内

这是最常见的问题。contains 判断的是严格包含,边界不算内部。如果业务上边界也算命中,请改用:

polygon.covers(point)

或者在只需要判断是否有交集时使用:

polygon.intersects(point)

坑2:坐标系不一致

Shapely 本身不理解坐标系。它只把坐标当作普通的 x、y 数值计算。如果点是 WGS84 经纬度,而面是 CGCS2000 高斯投影或 Web Mercator 米制坐标,结果肯定不可靠。

判断前要确认:

  • 点和面是否来自同一个坐标参考系统。
  • 经纬度是否写成了正确的 x、y 顺序。
  • 是否需要先用 GeoPandas、pyproj 或 GIS 软件统一投影。

使用 GeoPandas 时,可以这样统一坐标系:

points = points.to_crs(polygons.crs)

如果只是 Shapely 对象,没有 CRS 信息,你需要自己保证坐标已经统一。

坑3:多边形无效

如果多边形自相交、环方向混乱、洞区不合法,Shapely 几何关系计算可能出现异常或不符合预期。

可以先检查几何是否有效:

print(polygon.is_valid)

如果返回 False,需要修复几何。常见方式是使用 Shapely 的 make_valid,或在 QGIS 中运行“修复几何”工具。

from shapely.validation import make_valid

fixed_polygon = make_valid(polygon)

坑4:把经纬度顺序写反

Shapely 的 Point(x, y) 中,x 通常对应经度或横坐标,y 对应纬度或纵坐标。很多地图 API 返回的是 lat, lon,如果直接传入 Point(lat, lon),点的位置会完全错误。

正确写法通常是:

lon = 116.397
lat = 39.908

point = Point(lon, lat)

坑5:批量判断时速度慢

如果你有几十万个点和大量多边形,逐个循环判断会很慢。Shapely 可以配合空间索引或 GeoPandas 空间连接来提高效率。

GeoPandas 中更推荐使用 sjoin 做批量点面匹配:

import geopandas as gpd

points = gpd.read_file("points.geojson")
polygons = gpd.read_file("polygons.geojson")

points = points.to_crs(polygons.crs)

result = gpd.sjoin(
    points,
    polygons,
    predicate="within",
    how="left"
)

result.to_file("point_polygon_result.geojson", driver="GeoJSON")

如果边界点也要算,可以根据 GeoPandas 和 Shapely 版本支持情况选择合适的 predicate,或者对边界点做单独处理。

方法比较:contains、within、covers、intersects怎么选

Shapely 判断点在面内时,不建议所有情况都用 contains。不同方法代表不同业务含义。

业务问题 推荐方法 说明
点必须严格在多边形内部 polygon.contains(point)point.within(polygon) 边界点返回 False
点在内部或边界都算属于区域 polygon.covers(point) 更符合“覆盖范围”的业务语义
只要点和面有任何接触就算命中 polygon.intersects(point) 内部和边界都会返回 True
专门识别点是否落在边界上 polygon.touches(point) 适合数据质检和边界争议处理
批量点匹配多边形属性 geopandas.sjoin 适合大量数据,通常比手写循环更方便

简单记忆:严格内部用 containswithin;边界也算用 covers;只看是否相交用 intersects;专查边界用 touches

检查清单:写代码前先确认这些问题

  • 几何类型是否正确:点应该是 Point,面应该是 PolygonMultiPolygon
  • 坐标顺序是否正确:Shapely 使用 x, y,经纬度场景通常是 lon, lat
  • 坐标系是否一致:Shapely 不存储 CRS,必须由你自己保证点和面坐标系统一致。
  • 边界是否算命中:如果算,请用 covers 或合适的边界处理逻辑。
  • 多边形是否有效:polygon.is_valid 检查,必要时先修复几何。
  • 是否需要批量计算:大数据量建议用 GeoPandas 空间连接或空间索引。
  • 结果是否抽样验证:把几条结果加载到 QGIS 中查看,确认空间位置和判断逻辑一致。

FAQ:Shapely判断点在面内常见问题

Shapely判断点在面内应该用 contains 还是 within?

两者都可以。polygon.contains(point) 表示面包含点,point.within(polygon) 表示点位于面内。对于点面关系,它们方向不同但含义对应。需要注意的是,二者通常都不把边界点算作内部。

点在多边形边界上为什么 contains 返回 False?

因为 contains 是严格包含关系,边界不属于内部。如果业务规则认为边界点也属于该面,建议使用 polygon.covers(point)

intersects 能不能用来判断点在面内?

可以,但语义要明确。对于点和面来说,intersects 会在点位于面内部或边界上时返回 True。如果你只想判断严格内部,应该用 containswithin

Shapely会自动处理坐标系吗?

不会。Shapely 只做几何数值计算,不知道 EPSG 编码,也不会自动投影转换。判断前必须确保点和面的坐标系一致。

GeoJSON坐标传给 Shapely 时要注意什么?

GeoJSON 坐标顺序是 [x, y],在经纬度数据中通常就是 [经度, 纬度]。如果你从地图接口拿到的是 lat, lon,需要先调整顺序再创建 Point

MultiPolygon 怎么判断点在面内?

Shapely 的 MultiPolygon 也支持 containscoversintersects 等方法。只要点落在任意一个子多边形内部或边界上,对应方法就会返回相应结果。

from shapely.geometry import Point, MultiPolygon, Polygon

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

multi = MultiPolygon([poly1, poly2])
point = Point(5, 1)

print(multi.contains(point))

大量点和大量面如何提高判断速度?

不要直接双重循环逐个判断。可以使用 GeoPandas 的 sjoin,或者使用 Shapely 的空间索引能力先筛选候选几何,再做精确判断。对于 GIS 工程项目,GeoPandas 空间连接通常更容易维护。

结论:先明确边界规则,再选择Shapely几何关系方法

Shapely判断点在面内的核心并不复杂,但要把结果算对,必须先明确业务规则:边界点算不算在内、坐标系是否一致、多边形是否有效、是否需要批量处理。

如果你只需要严格判断点是否在面内部,用 polygon.contains(point)point.within(polygon)。如果边界也算在区域内,用 polygon.covers(point) 更合适。如果只是判断是否有任何接触,可以用 intersects。掌握这些几何关系后,Shapely 不仅能判断点在面内,也能支撑更复杂的 Python GIS 空间分析工作流。