问题描述
我有两对纬度/经度(以十进制度表示)及其半径(以米为单位).我要实现的目标是确定这两个点之间是否存在相交(当然,显然这并不适用于此,但计划是在许多其他数据点中尝试该算法).为了检查这一点,我使用了Shapely的intersects()函数.但是我的问题是我应该如何处理不同的部门?我是否应该先进行某种转换\投影(纬度\经度和半径的单位相同)?
I have two pairs of lat/lon (expressed in decimal degrees) along with their radius (expressed in meters). What I am trying to achieve is to find if an intersect between these two points exits (of course, it is obvious that this doesn't hold here but the plan is to try this algorithm in many other data points). In order to check this I am using Shapely's intersects() function. My question however is how should I deal with the different units? Should I make some sort of transformation \ projection first (same units for both lat\lon and radius)?
48.180759,11.518950,19.0
47.180759,10.518950,10.0
我在这里找到了这个库( https://pypi.python.org/pypi/utm)似乎很有帮助.但是,我不能百分百确定我是否正确应用了它.有什么想法吗?
I found this library here (https://pypi.python.org/pypi/utm) which seems helpfull. However, I am not 100% sure if I apply it correctly. Any ideas?
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)
解决方案:
因此,我终于设法解决了我的问题.这是两个解决相同问题的不同实现:
So, I finally managed to solve my problem. Here are two different implementations that both solve the same problem:
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仅使用笛卡尔坐标系,因此为了了解公制距离,您需要:
Shapely only uses the Cartesian coordinate system, so in order to make sense of metric distances, you would need to either:
- 将坐标投影到使用米的距离单位(例如UTM区域)的本地投影系统中.
- 从(0,0)开始缓冲一点,并使用动态的方位角等距投影以纬度/经度点为中心,以投影到地理坐标.
- project the coordinates into a local projection system that uses distance units in metres, such as a UTM zone.
- buffer a point from (0,0), and use a dynamic azimuthal equidistant projection centered on the lat/lon point to project to geographic coords.
这是使用 shapely.ops.transform 和 pyproj
import pyproj
from shapely.geometry import Point
from shapely.ops import transform
from functools import partial
WGS84 = pyproj.Proj(init='epsg:4326')
def latlonbuffer(lat, lon, radius_m):
proj4str = '+proj=aeqd +lat_0=%s +lon_0=%s +x_0=0 +y_0=0' % (lat, lon)
AEQD = pyproj.Proj(proj4str)
project = partial(pyproj.transform, AEQD, WGS84)
return transform(project, Point(0, 0).buffer(radius_m))
A = latlonbuffer(48.180759, 11.518950, 19.0)
B = latlonbuffer(47.180759, 10.518950, 10.0)
print(A.intersects(B)) # False
您的两个缓冲点不相交.但是这些确实可以做到:
Your two buffered points don't intersect. But these do:
A = latlonbuffer(48.180759, 11.518950, 100000.0)
B = latlonbuffer(47.180759, 10.518950, 100000.0)
print(A.intersects(B)) # True
如绘制lon/lat坐标(扭曲圆)所示:
As shown by plotting the lon/lat coords (which distorts the circles):
这篇关于查找两个地理数据点之间的交点的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持!