勻稱和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數組做一個大的調用，這就是爲什麼它是有效的。 請注意，您可以傳遞實際上使多邊形膨脹或降低的半徑，也可以在執行檢查之前變換（投影...）。