geolocation - 地理空间坐标和距离(以公里为单位)

Question

这是 this question 的后续.

我似乎被困在这一点上。基本上，我需要能够来回转换以引用标准度数系统中的坐标，或者通过测量沿国际日期变更线从南极向北的距离，然后从日期的那个点开始向东测量距离线。要做到这一点(以及一些更一般的距离测量的东西)，我有一种方法来确定两个纬度/经度点之间的距离，另一种方法是采用纬度/经度点、航向和距离，然后返回该类(class)结束时的纬度/经度点。

这是我定义的两个静态方法:

/* Takes two lon/lat pairs and returns the distance between them in kilometers.
*/
public static double distance (double lat1, double lon1, double lat2, double lon2) {
    double theta = toRadians(lon1-lon2);
    lat1 = toRadians(lat1);
    lon1 = toRadians(lon1);
    lat2 = toRadians(lat2);
    lon2 = toRadians(lon2);

    double dist = sin(lat1)*sin(lat2) + cos(lat1)*cos(lat2)*cos(theta);
    dist = toDegrees(acos(dist)) * 60 * 1.1515 * 1.609344 * 1000;

    return dist;
}

/* endOfCourse takes a lat/lon pair, a heading (in degrees clockwise from north), and a distance (in kilometers), and returns
 * the lat/lon pair that would be reached by traveling that distance in that direction from the given point.
 */
public static double[] endOfCourse (double lat1, double lon1, double tc, double dist) {
    double pi = Math.PI;
    lat1 = toRadians(lat1);
    lon1 = toRadians(lon1);
    tc = toRadians(tc);
    double dist_radians = toRadians(dist / (60 * 1.1515 * 1.609344 * 1000));
    double lat = asin(sin(lat1) * cos(dist_radians) + cos(lat1) * sin(dist_radians) * cos(tc));
    double dlon = atan2(sin(tc) * sin(dist_radians) * cos(lat1), cos(dist_radians) - sin(lat1) * sin(lat));
    double lon = ((lon1-dlon + pi) % (2*pi)) - pi;
    double[] endPoint = new double[2];
    endPoint[0] = lat; endPoint[1] = lon;
    return endPoint;
}

这是我用来测试它的函数:

public static void main(String args[]) throws java.io.IOException, java.io.FileNotFoundException {
    double distNorth = distance(0.0, 0.0, 72.0, 0.0);
    double distEast = distance(72.0, 0.0, 72.0, 31.5);
    double lat1 = endOfCourse(0.0, 0.0, 0.0, distNorth)[0];
    double lon1 = endOfCourse(lat1, 0.0, 90.0, distEast)[1];
    System.out.println("end at: " + lat1 + " / " + lon1);
    return;
}

“结束于”值应为 appx。 72.0/31.5。但相反，我得到了大约 1.25/0.021。

我想我一定是遗漏了一些愚蠢的东西，忘记在某处转换单位，或者什么......任何帮助将不胜感激!

更新 1:

我已经(正确地)编写了距离函数来返回米，但错误地在评论中写了公里......这当然让我今天回到它时感到困惑。无论如何，现在已经解决了，我已经修复了 endOfCourse 方法中的因式分解错误，而且我还意识到我也忘记了在该方法中将弧度转换回度数。无论如何:虽然看起来我现在得到了正确的纬度数(71.99 ...)，但经度数相差甚远(我得到的是 3.54 而不是 11.5)。

更新 2:
我在测试中有一个错字，如下所述。它现在已在代码中修复。然而，经度数仍然是错误的:我现在得到的是 -11.34 而不是 11.5。我认为这些行肯定有问题:

double dlon = atan2(sin(tc) * sin(dist_radians) * cos(lat1), cos(dist_radians) - sin(lat1) * sin(lat));
double lon = ((lon1-dlon + pi) % (2*pi)) - pi;

Answer 1

你有一个严重的代码中的魔数(Magic Number)案例。表达方式:

 (60 * 1.1515 * 1.609344 * 1000)

出现了两次，但没有太多解释。在一些帮助下:1.609344 是一英里的公里数； 60 是一个度的分钟数； 1000 是一公里的米数； 1.1515 是一海里的法定英里数(谢谢，DanM)。一海里是赤道上一分钟纬度的长度。

我假设您使用的是球形地球模型，而不是球形地球？代数不够复杂，不足以成为球体。

第一个公式——两个纬度和经度对之间的转换——是奇怪的。您需要 delta-lat (Δλ) 和 delta-lon (Δφ) 来整理答案。此外，对之间的距离:

(60° N, 30° W), (60° N, 60° W)
(60° N, 60° W), (60° N, 90° W)

应该是一样的 - 但我很确定你的代码会产生不同的答案。

所以，我认为你需要回到你的球面三角学引用资料，看看你做错了什么。 (我需要一段时间才能找到关于这个主题的书——无论它在哪个盒子里，它都需要打开包装。)

[...时间过去了...开箱完成...]

给定一个球面三角形，其顶点为 A、B、C，与这些顶点相对的边 a、b、c(即边 a 是从 B 到 C 等)，余弦公式为:

cos a = cos b . cos c + sin b . sin c . cos A

将此应用到问题中，我们可以调用给定 B 和 C 的两个点，并创建一个在 A 处为直角的直角球面三角形。

最糟糕的 ASCII 艺术:

                  + C
                 /|
                / |
            a  /  | b
              /   |
             /    |
            /     |
         B +------+ A
              c

c边等于经度差； b边等于纬度差；角度 A 是 90°，所以 cos A = 0。因此，我相信 a 的方程是:

cos a = cos Δλ . cos Δφ + sin Δλ . sin Δφ . cos 90°

a = arccos (cos Δλ . cos Δφ)

然后通过乘以地球半径将以弧度为单位的角度 a 转换为距离。或者，给定以度数(和度数的分数)，则有 60 海里到一个度，因此有 60 * 1.1515 法定英里，60 * 1.1515 * 1.609344 公里到一个度。除非您想要以米为单位的距离，否则我认为不需要 1000 的系数。

保罗·汤布林指向 Aviation Formulary v1.44作为等式的来源 - 事实上，当位置差异很小时，它与数值上更稳定的版本一起存在。

谈到基本的三角学，我们也知道:

cos (A - B) = cos A . cos B + sin A . sin B

在我给出的等式中应用两次很可能最终会出现在航空公式中的公式中。

(我的引用文献:"Astronomy: Principles and Practice, Fourth Edition" 作者:AE Roy 和 D Clarke (2003)；我的副本是 1977 年的第一版，Adam Hilger，ISBN 0-85274-346-7。)

备注 查看(谷歌)'定义:“海里”'；根据定义，海里现在似乎是 1852 m (1.852 km)。乘数 1.1515 对应于海里的旧定义，大约为 6080 英尺。使用 bc以 10 为尺度，我得到:

(1852/(3*0.3048))/1760
1.1507794480

哪个因素适合你取决于你的基础是什么。

从第一性原理看第二个问题，我们有一个稍微不同的设置，我们需要“另一个”球面三角方程，正弦公式:

sin A   sin B   sin C
----- = ----- = -----
sin a   sin b   sin c

改编上一张图:

                  + C
                 /|
                / |
            a  /  | b
           |  /   |
           |X/    |
           |/     |
         B +------+ A
              c

给定起点 B，角度 X = 90º - B，长度(角度)a，角度 A = 90°。您所追求的是 b(纬度增量)和 c(经度增量)。

所以，我们有:

sin a   sin b
----- = ----
sin A   sin B

或者

        sin a . sin B
sin b = -------------
            sin A

或者，由于 A = 90°，sin A = 1，sin B = sin (90° - X) = cos X:

sin b = sin a . cos X

这意味着您将行进的距离转换为角度 a，取其正弦，乘以航向的余弦，然后取结果的反正弦。

给定a、b(刚刚计算)和A、B，我们可以应用余弦公式得到c。请注意，我们不能简单地重新应用正弦公式来得到 c，因为我们没有 C 的值，而且因为我们在玩球面三角学，没有方便的规则 C = 90° - B(总和球面三角形中的角数可以大于 180°；考虑一个所有角都等于 90° 的等边球面三角形，这是完全可行的)。