计算给定方位和距离的坐标
Posted
技术标签:
【中文标题】计算给定方位和距离的坐标【英文标题】:Calculating coordinates given a bearing and a distance 【发布时间】:2010-10-27 00:36:20 【问题描述】:我在实现here 此处描述的功能时遇到问题。
这是我的 Java 实现:
private static double[] pointRadialDistance(double lat1, double lon1,
double radianBearing, double radialDistance)
double lat = Math.asin(Math.sin(lat1)*Math.cos(radialDistance)+Math.cos(lat1)
*Math.sin(radialDistance)*Math.cos(radianBearing));
double lon;
if(Math.cos(lat) == 0) // Endpoint a pole
lon=lon1;
else
lon = ((lon1-Math.asin(Math.sin(radianBearing)*Math.sin(radialDistance)/Math.cos(lat))
+Math.PI) % (2*Math.PI)) - Math.PI;
return (new double[]lat, lon);
我在调用函数之前将方位角转换为弧度并将距离(km)转换为弧度距离 - 所以这不是问题。
但是,当我输入坐标时,例如: 纬度 = 49.25705; lon = -123.140259; 方位角225(西南),距离1km
我得到这个返回: 纬度:-1.0085434360125864 lon:-3.7595299668539504
显然不正确,谁能看出我做错了什么?
谢谢
【问题讨论】:
尝试一些非常简单的输入。例如,输入 lat == lon == 0(实际上是在非洲附近),方位角为零,距离为零。你找回你的出发点了吗?如果是这样,请扩展它以尝试其他经纬度。然后尝试添加一个范围:你得到了什么有意义的东西吗? 【参考方案1】:您的代码中似乎存在以下问题:
-
在调用函数之前,您需要将
lat1
和 lon1
转换为弧度。
您可能错误地缩放radialDistance
。
测试浮点数是否相等是危险的。在精确算术之后相等的两个数字在浮点算术之后可能不完全相等。因此,abs(x-y) < threshold
比 x == y
更安全地测试两个浮点数 x
和 y
是否相等。
我认为您想将 lat
和 lon
从弧度转换为度数。
这是我在 Python 中的代码实现:
#!/usr/bin/env python
from math import asin,cos,pi,sin
rEarth = 6371.01 # Earth's average radius in km
epsilon = 0.000001 # threshold for floating-point equality
def deg2rad(angle):
return angle*pi/180
def rad2deg(angle):
return angle*180/pi
def pointRadialDistance(lat1, lon1, bearing, distance):
"""
Return final coordinates (lat2,lon2) [in degrees] given initial coordinates
(lat1,lon1) [in degrees] and a bearing [in degrees] and distance [in km]
"""
rlat1 = deg2rad(lat1)
rlon1 = deg2rad(lon1)
rbearing = deg2rad(bearing)
rdistance = distance / rEarth # normalize linear distance to radian angle
rlat = asin( sin(rlat1) * cos(rdistance) + cos(rlat1) * sin(rdistance) * cos(rbearing) )
if cos(rlat) == 0 or abs(cos(rlat)) < epsilon: # Endpoint a pole
rlon=rlon1
else:
rlon = ( (rlon1 - asin( sin(rbearing)* sin(rdistance) / cos(rlat) ) + pi ) % (2*pi) ) - pi
lat = rad2deg(rlat)
lon = rad2deg(rlon)
return (lat, lon)
def main():
print "lat1 \t lon1 \t\t bear \t dist \t\t lat2 \t\t lon2"
testcases = []
testcases.append((0,0,0,1))
testcases.append((0,0,90,1))
testcases.append((0,0,0,100))
testcases.append((0,0,90,100))
testcases.append((49.25705,-123.140259,225,1))
testcases.append((49.25705,-123.140259,225,100))
testcases.append((49.25705,-123.140259,225,1000))
for lat1, lon1, bear, dist in testcases:
(lat,lon) = pointRadialDistance(lat1,lon1,bear,dist)
print "%6.2f \t %6.2f \t %4.1f \t %6.1f \t %6.2f \t %6.2f" % (lat1,lon1,bear,dist,lat,lon)
if __name__ == "__main__":
main()
这是输出:
lat1 lon1 bear dist lat2 lon2
0.00 0.00 0.0 1.0 0.01 0.00
0.00 0.00 90.0 1.0 0.00 -0.01
0.00 0.00 0.0 100.0 0.90 0.00
0.00 0.00 90.0 100.0 0.00 -0.90
49.26 -123.14 225.0 1.0 49.25 -123.13
49.26 -123.14 225.0 100.0 48.62 -122.18
49.26 -123.14 225.0 1000.0 42.55 -114.51
【讨论】:
谢谢,我忘记将 rad 转换回 deg!正如你在这里所做的那样: lat = rad2deg(rlat) lon = rad2deg(rlon) 谢谢你抽出时间来帮助我 您可以使用内置的 math.radians 和 math.degrees 来代替自定义函数。 另外我认为:rdistance = distance / rEarth # 将线性距离标准化为弧度角是个坏主意,如果距离为整数,可能会导致 rdistance 为零。 @las3rjock 这太棒了!试图自己实现算法,这有很大帮助;我认为它只需要一个小小的改变。rlon = ( (rlon1 - asin( sin(rbearing)* sin(rdistance) / cos(rlat) ) + pi ) % (2*pi) ) - pi
应该是 rlon = ( (rlon1 + asin( sin(rbearing)* sin(rdistance) / cos(rlat) ) + pi ) % (2*pi) ) - pi
。如果假设负 lons 为西,则需要添加 arcsin,否则提供的方位会给您错误的方向(例如,90 应该是东,但会根据您的方程给您西)。
警告:此代码产生不正确的结果。 pointRadialDistance(52.20472, 0.14056, 90, 15)
应该输出 (52.20451603459371, 0.3599721245292571)
但我们得到 )(52.20451523812389, -0.07955815309722394)
。 pointRadialDistance(-32.06, 115.74, 225, 20000)
应该输出 (32.11195529143165, -63.95925278363718)
但我们得到 (31.96383452118179, 115.85329213631711)
。可以在这里找到好的实现:***.com/questions/7222382/…【参考方案2】:
我认为消息 5 中提供的算法存在问题。
它有效,但仅适用于纬度,对于经度,由于标志存在问题。
数据不言自明:
49.26 -123.14 225.0 1.0 49.25 -123.13
如果您从 -123.14° 开始向西走,您应该在西边有 FAR。我们回到东 (-123.13) !
公式应该包含某处:
degreeBearing = ((360-degreeBearing)%360)
在弧度转换之前。
【讨论】:
【参考方案3】:从根本上说,您的问题似乎是您将纬度、经度和方位作为度数而不是弧度传递。试着确保你总是将弧度传递给你的函数,看看你得到了什么。
PS:请参阅here 和here 讨论的类似问题。
【讨论】:
不,正如我所说,我将方位角转换为弧度,然后再将其传递给函数。但是,我没有转换纬度和经度,它们是否必须转换为弧度,可以吗? 是的,您应该以弧度输入所有角度输入。【参考方案4】:当我实现这个时,我得到的纬度是正确的,但经度是错误的。 例如起点:36.9460678N 9.434807E,方位角45.03334,距离15.0083313km 结果是 37.0412865N 9.315302E 那比我的起点更西,而不是更东。事实上,就好像方位角是 315.03334 度。
更多网络搜索导致我:http://www.movable-type.co.uk/scripts/latlong.html 经度代码如下所示(在 C# 中,所有内容均以弧度表示)
if ((Math.Cos(rLat2) == 0) || (Math.Abs(Math.Cos(rLat2)) < EPSILON))
rLon2 = rLon1;
else
rLon2 = rLon1 + Math.Atan2(Math.Sin(rBearing) * Math.Sin(rDistance) * Math.Cos(rLat1), Math.Cos(rDistance) - Math.Sin(rLat1) * Math.Sin(rLat2));
这对我来说似乎很好用。希望对您有所帮助。
【讨论】:
【参考方案5】:感谢您的 Python 代码 我尝试在我的用例中设置它 我试图在其他两个之间找到一个点的纬度 距第一点一定距离,因此与您的代码非常相似 appart 我的方位是动态计算的
起点(lat1) lon1/lat1 = 55.625541,-21.142463
终点 (lat2) lon2/lat2 = 55.625792,-22.142248
我的结果应该是 lon3/lat3 这两者之间的一个点 不幸的是,我得到 lon3/lat3 = 0.0267695450609,0.0223553243666
我认为这可能是纬度的差异,但没有 当我添加或替换它时,它并不好
任何建议都会很棒 谢谢
这是我的实现
距离 = 0.001 ε = 0.000001
动态计算方位
y = math.sin(distance) * math.cos(lat2);
x = math.cos(lat1)*math.sin(lat2) - math.sin(lat1)*math.cos(lat2)*math.cos(distance);
bearing = math.atan2(y, x)
动态计算 lat3 lon3
rlat1 = (lat1 * 180) / math.pi
rlon1 = (lon1 * 180) / math.pi
rbearing = (bearing * 180) / math.pi
rdistance = distance / R # normalize linear distance to radian angle
rlat = math.asin( math.sin(rlat1) * math.cos(rdistance) + math.cos(rlat1) * math.sin(rdistance) * math.cos(rbearing) )
if math.cos(rlat) == 0 or abs(math.cos(rlat)) < epsilon: # Endpoint a pole
rlon=rlon1
else:
rlon = ( (rlon1 + math.asin( math.sin(rbearing)* math.sin(rdistance) / math.cos(rlat) ) + math.pi ) % (2*math.pi) ) - math.pi
lat3 = (rlat * math.pi)/ 180
lon3 = (rlon * math.pi)/ 180
【讨论】:
【参考方案6】:一切都按预期工作,但问题是您的数学假设地球是一个球体,而实际上它近似于一个椭球体。
快速搜索一下您最喜欢的“Vincenty Formula”搜索引擎有望证明是有用的。
【讨论】:
关于Vincenty formulae的更多信息请查找here。以上是关于计算给定方位和距离的坐标的主要内容,如果未能解决你的问题,请参考以下文章