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

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

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

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

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

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

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

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

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) 面是否覆盖点,通常包含边界 业务上边界也算在区域内时更合适

最容易混淆的是 contains 和 intersects。如果点在多边形边界上,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 改写判断逻辑

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

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

如果你的业务含义是“点在面内或边界上都算属于该面”,优先考虑 covers 或 intersects。如果只是点和面之间的关系,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 适合大量数据,通常比手写循环更方便

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

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

  • 几何类型是否正确:点应该是 Point,面应该是 Polygon 或 MultiPolygon。
  • 坐标顺序是否正确: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。如果你只想判断严格内部,应该用 contains 或 within。

Shapely会自动处理坐标系吗?

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

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

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

MultiPolygon 怎么判断点在面内?

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

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 空间分析工作流。