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

Posted 2023-02-14

【中文标题】检查地理点是在多边形内部还是外部

【英文标题】：Check if geo-point is inside or outside of polygon

【发布时间】：2017-10-09 02:25:49

【问题描述】：

我正在使用 python，并且我已经在地图上定义了多边形的纬度和经度（以度为单位）。 我的目标是检查坐标x,y 的通用点P 是否落在这样的多边形内。因此，我希望有一个函数允许我检查这种情况，如果该点在多边形内部或外部，则返回 True 或 False。

在这个例子中，点在外面，所以结果是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()  

注意：

点通过大圆圈相互连接！ 我在互联网上进行了研究，最终发现了一些类似的问题like this one，但我没有成功，因为它们都使用了我没有的.shp 文件。

【问题讨论】：

尝试将此算法转换为 Python wrf.ecse.rpi.edu//Research/Short_Notes/pnpoly.html#The C 代码

python 没有可以做任何事情的包。它有少量的预构建模块。软件包通常由社区提供。

多边形总是凸的吗？

一般不会，也可以是凹的

以防万一：您始终可以将光线从您的点投射到多边形任何边的中点。如果您的光线穿过多边形的边偶数次，则该点位于外侧。适用于凸多边形和凹多边形；使用测地线作为射线在球体表面（可能是任何 1 连通表面）上工作。当光线完全通过一个顶点时有边缘情况：您需要检查入射到顶点的边缘是否在光线的同一侧。

【参考方案1】：

这是我的问题的可能解决方案。

必须正确存储地理坐标。示例np.array([[Lon_A, Lat_A], [Lon_B, Lat_B], [Lon_C, Lat_C]]) 创建多边形 创建要测试的点 使用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 

【讨论】：

最好先写纬度，后写经度。这里的逻辑没有错，但为了安全起见。

@FedericoGentile - Shapely 不使用大圆距离。它使用欧几里得

这是否也适用于以平滑曲线为界的区域？不仅仅是一个多边形。

谢谢！答案很紧凑，但回答了我需要的所有问题！酷！

当多边形与日期线相交时，Shapely 也会失败

【参考方案2】：

还有一个新兴的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)

【讨论】：

我在安装匀称时遇到了很多麻烦。这是一个更清洁的解决方案 - 谢谢！

【参考方案3】：

另一种方法是使用此链接中解释的奇偶算法https://wrf.ecse.rpi.edu//Research/Short_Notes/pnpoly.html ***中给出了python代码https://en.wikipedia.org/wiki/Even–odd_rule

伙计们，请记住使多边形很重要的点顺序！我的意思是，不同的顺序会产生不同的多边形。

【参考方案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