shapely 和 matplotlib 多边形点在地理位置上不准确
Posted
技术标签:
【中文标题】shapely 和 matplotlib 多边形点在地理位置上不准确【英文标题】:shapely and matplotlib point-in-polygon not accurate with geolocation 【发布时间】:2014-02-15 05:19:23 【问题描述】:我正在使用 matplotlib 和 shapely 测试多边形中的点函数。
这是一个 map 包含一个百慕大三角形多边形。
Google 地图 的多边形内点功能清楚地显示 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
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
这里到底发生了什么?谷歌的算法比这两个更好吗?
谢谢
【问题讨论】:
【参考方案1】:我只是为了测试这些点是否真的在三角形内:
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 地图,并且将多边形映射到球坐标上时,三角形会发生变形,请记住这一点。
无论如何,在 Google 地球中使用 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稍微在三角形之外;在欧几里得二维坐标中绘制时。 对于地理坐标中的几何计算,请查看 QGIS Python 控制台或 GDAL/OGR 工具;或者您可以使用 Google Maps API,就像链接在 this page 上的示例一样,其中涵盖了 2D 几何 VS 测地几何主题。
【讨论】:
感谢您的解释。 matplotlib可以使用球坐标进行匹配吗? 代码中的“x”和“y”是什么?您的示例不完整。【参考方案2】:记住:世界不是平的!如果 Google 地图的投影是您想要的答案,您需要将地理坐标投影到 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
似乎得到了正确的答案。
【讨论】:
非常感谢这个想法,它确实解决了问题。顺便说一句,如果运行它,我发现 shapely 比 matplotlib 慢得多,比如 1000 次。 很高兴知道,我得看看是什么让 matplotlib 这么快。 有什么办法可以用matplotlib
来代替吗?【参考方案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
【讨论】:
【参考方案4】:要检查一个多边形是否包含多个点,我会使用 matplotlib contains_points
,记录在这里:http://matplotlib.org/api/path_api.html#matplotlib.path.Path.contains_points
这使用 numpy 数组进行了一次大调用,这就是它高效的原因。 请注意,您可以传递一个实际上使多边形膨胀或缩小的半径,您还可以在进行检查之前进行变换(投影...)。
【讨论】:
以上是关于shapely 和 matplotlib 多边形点在地理位置上不准确的主要内容,如果未能解决你的问题,请参考以下文章
Python使用matplotlib可视化聚类图使用encircle函数绘制多边形标定属于同一聚类簇的数据点并自定义每个聚类簇的背景色(Cluster Plot)