检查地理点是否在多边形内部或外部

61

我正在使用Python,已经定义了地图上多边形的纬度和经度(以度为单位)。我的目标是检查一个通用点P的坐标x,y是否在多边形内部。因此,我需要一个函数来检查这个条件,并返回TrueFalse,以指示点是否在多边形内部或外部。

enter image description here

在这个例子中,该点位于多边形外部,所以结果应该是False

问题:有没有库/包可以实现我的目标?如果有,你会推荐哪一个?你能否给出一个小例子来说明如何使用它?

以下是我到目前为止编写的代码:

import numpy as np

# Define vertices of polygon (lat/lon)
v0 = [7.5, -2.5] 
v1 = [2, 3.5]
v2 = [-2, 4]
v3 = [-5.5, -4]
v4 = [0, -10]
lats_vect = np.array([v0[0],v1[0],v2[0],v3[0],v4[0]])
lons_vect = np.array([v0[1],v1[1],v2[1],v3[1],v4[1]])

# Point of interest P
x, y = -6, 5 # x = Lat, y = Lon

## START MODIFYING FROM HERE; DO NOT MODIFY POLYGON VERTICES AND DATA TYPE
# Check if point of interest falls within polygon boundaries
# If yes, return True
# If no, return False
为了绘制多边形和感兴趣的点,我使用了Cartopy并编写了以下代码行:
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
ax = plt.axes(projection=ccrs.PlateCarree())
ax.stock_img() 

# Append first vertex to end of vector to close polygon when plotting
lats_vect = np.append(lats_vect, lats_vect[0])
lons_vect = np.append(lons_vect, lons_vect[0])
plt.plot([lons_vect[0:-1], lons_vect[1:]], [lats_vect[0:-1], lats_vect[1:]],
         color='black', linewidth=1, 
         transform=ccrs.Geodetic(),
         )   

plt.plot(y, x, 
        '*',          # marker shape
        color='blue',  # marker colour
        markersize=8  # marker size
        )  

plt.show()  

注意

  • 点之间是通过大圆连接的!
  • 我在互联网上进行了研究,并最终找到了一些类似的问题 像这个,但我没有成功,因为它们都使用了我没有的.shp文件。

尝试将此算法转换为Python https://wrf.ecse.rpi.edu//Research/Short_Notes/pnpoly.html#The C Code - arboreal84
Python没有可以做任何事情的包。它只有一些预先构建的模块。包通常由社区提供。 - Uriel
1
多边形是否总是凸的? - user2261062
一般来说不是,它也可能是凹的。 - Federico Gentile
看一下这个问题并实现这个算法。它不难。https://dev59.com/p3VC5IYBdhLWcg3wrDJd - user2261062
1
以防万一:您始终可以从您的点向任何多边形边的中点投射光线。如果您的光线穿过多边形的边偶数次,则该点在外部。适用于凸多边形和凹多边形;使用测地线进行光线投射,可在球面(以及可能的任何1连通表面)上工作。当光线恰好通过一个顶点时存在一个特殊情况:您需要检查与该顶点相邻的边是否在光线的同侧。 - 9000
4个回答

54

这是解决我的问题的可能方法:

  1. 正确地存储地理坐标。例如,np.array([[Lon_A, Lat_A], [Lon_B, Lat_B], [Lon_C, Lat_C]])
  2. 创建多边形
  3. 创建待测试点
  4. 使用polygon.contains(point)来测试点是否在多边形内(True)还是外部(False)。

这是缺失的代码部分:

from shapely.geometry import Point
from shapely.geometry.polygon import Polygon

lons_lats_vect = np.column_stack((lons_vect, lats_vect)) # Reshape coordinates
polygon = Polygon(lons_lats_vect) # create polygon
point = Point(y,x) # create point
print(polygon.contains(point)) # check if polygon contains point
print(point.within(polygon)) # check if a point is in the polygon 

注意: 多边形不考虑大圆,因此需要将边缘分成许多线段,从而增加顶点的数量。


特殊情况: 如果点位于多边形边界上

例如:print(Polygon([(0,0),(1,0),(1,1)]).contains(Point(0,0))) 会失败

所以可以使用

print(polygon.touches(point)) # check if point lies on border of polygon 

2
最好先写纬度,然后是经度。这里的逻辑没有问题,但小心谨慎总是没错的。 - Zahran
1
@FedericoGentile - Shapely 不使用大圆距离,它使用欧几里得距离。 - gansub
这对于由光滑曲线界定的区域也适用吗?不仅仅是多边形。 - ap21
1
谢谢!答案简洁明了,回答了我所有需要的问题!太棒了! - ptee
Shapely 在多边形与日期变更线相交时也会失败。 - Alexander Bauer

28

还有一个新兴的 Python 库 turfpy,用于地理空间分析。

PyPI

Github

示例:

from turfpy.measurement import boolean_point_in_polygon
from geojson import Point, Polygon, Feature

point = Feature(geometry=Point((-46.6318, -23.5523)))
polygon = Polygon(
    [
        [
            (-46.653, -23.543),
            (-46.634, -23.5346),
            (-46.613, -23.543),
            (-46.614, -23.559),
            (-46.631, -23.567),
            (-46.653, -23.560),
            (-46.653, -23.543),
        ]
    ]
)
boolean_point_in_polygon(point, polygon)

1
我一直在安装shapely时遇到很多麻烦。这是一个更干净的解决方案 - 谢谢! - mustang
虽然 boolean_point_in_polygon(point, polygon) 完美地解决了 OP 的使用情况,但是对于聚合使用情况,值得检查 points_within_polygon(points, polygon) - Skippy le Grand Gourou

7

4

你可以使用pygeodesy包,它不依赖于系统级别的地理包,并且使用了Kenneth Gade的n-vector方法。

https://github.com/mrJean1/PyGeodesy

只需 pip install pygeodesy

示例代码

from pygeodesy.sphericalNvector import LatLon

p = LatLon(45.1, 1.1)
b = LatLon(45, 1), LatLon(45, 2), LatLon(46, 2), LatLon(46, 1)
print (p.isenclosedBy(b))

这应该会输出True


网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接