匀称和matplotlib点多边形不准确与地理定位

Question

我正在用matplotlib和匀称测试点多边形函数。

这里是一个map包含一个百慕大三角形多边形。

谷歌地图的点在多边形的功能清楚地示出testingPoint和testingPoint2是这是一个正确的结果多边形的内部。

如果我测试matplotlib两个点，并且很体面，只有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 

匀称示出了如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的算法比那两个更好吗？

感谢

Answer 1

切记：世界不平坦！如果Google Maps的投影是您想要的答案，则需要将地理坐标投影到spherical Mercator上以获取不同的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 

似乎得到正确的答案。

Answer 2

我只是这样做是为了测试是否点实际上是在三角形内：

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() 

现在当你使用谷歌地图和多边形映射到球面坐标，三角形变形了，需要牢记。

无论如何，在Gookle 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是slighlty三角形，在欧几里得2D坐标绘制时的外部。 对于地理坐标中的几何计算，请检查QGIS Python控制台或GDAL/OGR工具。或者，您可以使用谷歌地图API，就像在示例中那样，链接在this page上，其中主题2D几何与测地几何被覆盖。

Answer 3

虽然你已经接受了答案，但除了@ 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 

Answer 4

要检查一个多边形包含多个点，我将使用matplotlib contains_points，这里记录：http://matplotlib.org/api/path_api.html#matplotlib.path.Path.contains_points

这使用numpy数组做一个大的调用，这就是为什么它是有效的。 请注意，您可以传递实际上使多边形膨胀或降低的半径，也可以在执行检查之前变换（投影...）。