根据方位角和距离计算坐标

Question

根据方位角和距离计算坐标

6

我在实现这里描述的函数时遇到了问题。

这是我的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});
}

我将角度转换为弧度，并将距离（公里）转换为弧度距离后再调用函数 - 所以这不是问题。

但是，当我输入以下坐标：纬度=49.25705；经度=-123.140259；方位角为225（西南），距离为1公里

我得到如下结果：纬度：-1.0085434360125864 经度：-3.7595299668539504

显然不正确，有人能看出我做错了什么吗？

谢谢

- user106996

尝试一些非常简单的输入。例如，输入lat == lon == 0（实际上是靠近非洲）的方向为零和距离为零。你是否得到了你的起点？如果是，请尝试扩展到其他纬度和经度。然后尝试添加一个范围：你是否得到了合理的结果？ - Bob Cross

7个回答

4

我认为消息5中提供的算法存在问题。

它可以处理纬度，但对于经度会有问题，因为符号不正确。

数据说明问题：

49.26 -123.14 225.0 1.0 49.25 -123.13

如果您从-123.14°开始向西行驶，应该找到一个很远的东部地区。但是我们却回到了东部（-123.13）！

该公式应该包含以下内容：

degreeBearing = ((360-degreeBearing)%360)

在弧度转换之前进行操作。

- Stéphane

3

基本上，你的问题在于你传递的纬度、经度和方位角是以度为单位而不是弧度。请尝试确保始终将弧度作为参数传递给函数，并查看返回结果。

PS：请参见类似的讨论：这里和这里。

- Bob Cross

不，正如我所说的那样，在将方位角传递给函数之前，我会将其转换为弧度。但是，我还没有转换纬度和经度，它们需要转换为弧度吗？这可行吗？ - user106996

1

是的，您应该将所有角度输入转换为弧度。 - Bob Cross

1

当我实现这个功能时，我的纬度是正确的，但经度是错误的。例如起点：36.9460678N 9.434807E，方位角45.03334，距离15.0083313公里结果是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));
        }

这对我来说似乎很好用。希望它能帮到您。

- user164104

0

# -*- coding: utf-8 -*-
from math import asin, cos, pi, sin
import pandas as pd`enter code here`

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):
    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():
    source_path = "Source_File.xlsx"
    destination_path = "Destination_File.xlsx"
    
    df = pd.read_excel(source_path)
    for index, row in df.iterrows():
        latitude = row['Latitude']
        longitude = row['Longitude']
        bearing = 60
        distance = 0.060
        (lat,lon) = pointRadialDistance(latitude, longitude, bearing, distance)
        df.loc[index, 'New Latitude'] = lat
        df.loc[index, 'New Longitude'] = lon
        
    df.to_excel(destination_path, index=False)


if __name__ == "__main__":
    main()

- Vishal Gaddam

1

请详细说明您对这个12年前的问题的回答。提供示例输入/输出也会很有帮助。 - AcK

0

感谢您的Python代码，我尝试在我的用例中设置它，其中我正在尝试找到两个点之间的一个点的纬度和经度，该点距离第一个点的一定距离，因此它与您的代码非常相似，除了我的方位角是动态计算的。

起始点（lat1）lon1 / lat1 = 55.625541，-21.142463

终点（lat2）lon2 / lat2 = 55.625792，-22.142248

我的结果应该是这两个点之间的一个点，其经度/纬度为lon3 / lat3，但不幸的是，我得到的lon3 / lat3 = 0.0267695450609,0.0223553243666。

我认为这可能是纬度和经度的差异，但实际上并不是这样，当我加或减时结果并不正确。

任何建议都将非常棒，谢谢。

这是我的实现：

距离= 0.001 epsilon = 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

- Tibor

0

一切都按照预期运作，但问题在于你的数学假设地球是一个球体，而实际上它近似于一个椭球体。

在你喜欢的搜索引擎上快速搜索“Vincenty公式”将会有所帮助。

- Ximines

请在此处查看更多关于Vincenty公式的信息。 - menjaraz

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- las3rjock · Accepted Answer

您的代码似乎存在以下问题：

在调用函数之前，您需要将lat1和lon1转换为弧度。
您可能错误地缩放了radialDistance。
对浮点数进行相等性测试是危险的。即使在精确算术后相等的两个数字，在浮点算术后也可能不完全相等。因此，对于测试两个浮点数x和y是否相等，abs(x-y) < threshold比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