Shapely和Matplotlib在地理位置上的点多边形判断不准确

Question

我正在使用matplotlib和shapely测试点与多边形的关系。

这里有一个包含百慕大三角形多边形的地图。

Google maps中的点与多边形函数清晰地显示testingPoint和testingPoint2在多边形内，这是正确的结果。

如果我在matplotlib和shapely中测试这两个点，只有point2通过测试。

In [1]: from matplotlib.path import Path

In [2]: p = Path([[25.774252, -80.190262], [18.466465, -66.118292], [32.321384, -64.75737]]) 

In [3]: p1=[27.254629577800088, -76.728515625]

In [4]: p2=[27.254629577800088, -74.928515625]

In [5]: p.contains_point(p1)
Out[5]: 0

In [6]: p.contains_point(p2)
Out[6]: 1

Shapely 显示的结果与 matplotlib 相同。

In [1]: from shapely.geometry import Polygon, Point

In [2]: poly = Polygon(([25.774252, -80.190262], [18.466465, -66.118292], [32.321384, -64.75737]))

In [3]: p1=Point(27.254629577800088, -76.728515625)

In [4]: p2=Point(27.254629577800088, -74.928515625)

In [5]: poly.contains(p1)
Out[5]: False

In [6]: poly.contains(p2)
Out[6]: True

这里到底发生了什么？Google的算法比这两个更好吗？

谢谢。

- Chung

Answer 1

请记住：世界并非平面！ 如果您想使用Google Maps的投影作为您所需的答案，您需要将地理坐标转换成球形墨卡托以获得不同的X和Y坐标系。 Pyproj可以帮助您完成这项任务，只需在进行投影之前确保翻转坐标轴即可（即：X、Y或经度、纬度）。

import pyproj
from shapely.geometry import Polygon, Point
from shapely.ops import transform
from functools import partial

project = partial(
    pyproj.transform,
    pyproj.Proj(init='epsg:4326'),
    pyproj.Proj('+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs'))

poly = Polygon(([-80.190262, 25.774252], [-66.118292, 18.466465], [-64.75737, 32.321384]))
p1 = Point(-76.728515625, 27.254629577800088)

# Old answer, using long/lat coordinates
poly.contains(p1)  # False
poly.distance(p1)  # 0.01085626429747994 degrees

# Translate to spherical Mercator or Google projection
poly_g = transform(project, poly)
p1_g = transform(project, p1)

poly_g.contains(p1_g)  # True
poly_g.distance(p1_g)  # 0.0 meters

似乎得到了正确答案。

- Mike T

---

非常感谢您的想法，它确实解决了问题。顺便说一下，我发现如果运行1000次，shapely比matplotlib慢得多。 - Chung

好的，我得看看是什么让matplotlib运行更快。 - Mike T

有没有办法我可以用 matplotlib 来进行转换呢？ - Mikael S.

Answer 2

虽然您已经接受了一个答案，但除了 @MikeT 的答案之外，我还想为可能希望在 mpl_toolkit 中使用 matplotlib 和 basemap 的未来访问者添加以下内容：

from mpl_toolkits.basemap import Basemap
from matplotlib.path import Path


# Mercator Projection
# http://matplotlib.org/basemap/users/merc.html
m = Basemap(projection='merc', llcrnrlat=-80, urcrnrlat=80,
            llcrnrlon=-180, urcrnrlon=180, lat_ts=20, resolution='c')

# Poly vertices
p = [[25.774252, -80.190262], [18.466465, -66.118292], [32.321384, -64.75737]]

# Projected vertices
p_projected = [m(x[1], x[0]) for x in p]

# Create the Path
p_path = Path(p_projected)

# Test points
p1 = [27.254629577800088, -76.728515625]
p2 = [27.254629577800088, -74.928515625]

# Test point projection
p1_projected = m(p1[1], p1[0])
p2_projected = m(p2[1], p2[0])

if __name__ == '__main__':
    print(p_path.contains_point(p1_projected))  # Prints 1
    print(p_path.contains_point(p2_projected))  # Prints 1

- Sam R.

Answer 3

我刚刚这样做是为了测试点是否真的在三角形内部：

from matplotlib import pylab as plt
poly = [[25.774252, -80.190262],
        [18.466465, -66.118292],
        [32.321384, -64.75737],
        [25.774252, -80.190262]]
x = [point[0] for point in poly]
y = [point[1] for point in poly]
p1 = [27.254629577800088, -76.728515625]
p2 = [27.254629577800088, -74.928515625]
plt.plot(x,y,p1[0],p1[1],'*r',p2[0],p2[1],'*b')
plt.show()

现在当你使用Google Maps时，多边形被映射到球面坐标系上，这会导致三角形变形，请注意。

无论如何，使用Google Earth中的kml绘制数据时，也能显示点在三角形外面？！

<kml>
<Document>
<Placemark><name>Point 1</name><Point>
<coordinates> -76.728515625, 27.254629577800088,0</coordinates></Point></Placemark>
<Placemark><name>Point 2</name><Point>
<coordinates>-74.928515625, 27.254629577800088,     0</coordinates></Point></Placemark>
<Placemark><name>Poly</name><Polygon>
<outerBoundaryIs><LinearRing>
<coordinates> -80.190262,25.774252 -66.118292,18.466465 -64.75737,32.321384 -80.190262,25.774252</coordinates>
</LinearRing></outerBoundaryIs>
</Polygon></Placemark>
</Document>
</kml>

与 matplotlib 图像中的外观相同，点1 稍微在三角形外面；当在欧几里得 2D 坐标中绘制时。 对于地理坐标系中的几何计算，请查看 QGIS Python Console 或 GDAL/OGR 工具；或者您可以使用 Google Maps API，就像在此页面链接的示例中一样，其中涵盖了 2D 几何与大地测量几何的主题。

- Schuh

---

谢谢您的解释。Matplotlib能否使用球面坐标进行匹配？ - Chung

你的代码中的 'x' 和 'y' 是什么？你的示例不完整。 - Spacedman

Answer 4

为了检查多个点是否在多边形内部，我会使用matplotlib中的contains_points方法。你可以在这里找到它的文档：http://matplotlib.org/api/path_api.html#matplotlib.path.Path.contains_points 这个方法只需一次对numpy数组的调用，所以效率很高。 注意，你可以传递一个半径来膨胀或收缩多边形，你还可以在进行检查之前进行变换（投影等）。

- Christophe Roussy