我有两对纬度/经度(以十进制表示)及其半径(以米表示)。我想要实现的是找到这两点之间是否存在相交(当然,很明显这在这里不成立,但计划是在许多其他数据点中尝试该算法)。为了检查这一点,我使用 Shapely 的 intersects() 函数。然而我的问题是我应该如何处理不同的单位?我应该首先进行某种变换\投影(纬度\经度和半径的单位相同)吗?
48.180759,11.518950,19.0
47.180759,10.518950,10.0
编辑:
我在这里找到了这个库(https://pypi.python.org/pypi/utm),这似乎很有帮助。但是,我不能 100% 确定我是否正确应用它。有什么想法吗?
X = utm.from_latlon(38.636782, 21.414384)
A = geometry.Point(X[0], X[1]).buffer(30.777)
Y = utm.from_latlon(38.636800, 21.414488)
B = geometry.Point(Y[0], Y[1]).buffer(23.417)
A.intersects(B)
解决方案:
所以,我终于解决了我的问题。以下是解决同一问题的两种不同实现:
X = from_latlon(48.180759, 11.518950)
Y = from_latlon(47.180759, 10.518950)
print(latlonbuffer(48.180759, 11.518950, 19.0).intersects(latlonbuffer(47.180759, 10.518950, 19.0)))
print(latlonbuffer(48.180759, 11.518950, 100000.0).intersects(latlonbuffer(47.180759, 10.518950, 100000.0)))
X = from_latlon(48.180759, 11.518950)
Y = from_latlon(47.180759, 10.518950)
print(geometry.Point(X[0], X[1]).buffer(19.0).intersects(geometry.Point(Y[0], Y[1]).buffer(19.0)))
print(geometry.Point(X[0], X[1]).buffer(100000.0).intersects(geometry.Point(Y[0], Y[1]).buffer(100000.0)))
Shapely 仅使用 笛卡尔坐标系,因此为了理解公制距离,您需要:
这是如何执行#2,使用 shapely.ops.transform 和 pyproj
import pyproj
import shapely
from shapely.geometry import Point
WGS84 = pyproj.Proj(init='epsg:4326')
def latlonbuffer(lat, lon, radius_m):
proj4str = f'+proj=aeqd +lat_0={lat} +lon_0={lon} +x_0=0 +y_0=0 +units=m'
AEQD = pyproj.CRS.from_proj4(proj4str)
proj = pyproj.Transformer.from_crs(crs_from=AEQD, crs_to=WGS84).transform
p0 = Point(0, 0).buffer(radius_m)
return shapely.ops.transform(proj, p0)
A = latlonbuffer(48.180759, 11.518950, 19.0)
B = latlonbuffer(47.180759, 10.518950, 10.0)
print(A.intersects(B)) # False
您的两个缓冲点不相交。但这些确实:
A = latlonbuffer(48.180759, 11.518950, 100000.0)
B = latlonbuffer(47.180759, 10.518950, 100000.0)
print(A.intersects(B)) # True
如绘制经/纬度坐标所示(这会扭曲圆圈):