地理空间坐标和距离（以公里为单位）

Posted 2023-04-14

【中文标题】地理空间坐标和距离（以公里为单位）

【英文标题】：Geospatial coordinates and distance in kilometers

【发布时间】：2010-09-28 04:49:44

【问题描述】：

这是this question 的后续。

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

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

/* 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：

我已经（正确地）写了距离函数来返回米，但是错误地在 cmets 中写了公里......当我今天回到它时，这当然让我感到困惑。无论如何，现在已经解决了，并且我已经修复了 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;

【问题讨论】：

作为参考资料，您可能需要查看movable-type.co.uk/scripts/latlong-vincenty.html 和movable-type.co.uk/scripts/latlong-vincenty-direct.html（包括javascript 实现。

在你的第一个函数中，最好使用两个变量，'double angle = ...1st expression; double dist = ...第二个表达式修改为参考角度；'。事实上，第一个值根本不是距离。只是迂腐的大惊小怪。

【参考方案1】：

代码中的幻数很严重。表达式：

 (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 处有角度 A、B、C，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 然后通过乘以地球的半径转换为距离。或者，给定 a 度数（以及度数的小数），那么 60 海里对应一个度数，因此 60 * 1.1515 法定英里数和 60 * 1.1515 * 1.609344 公里对应一个度数。除非您想要以米为单位的距离，否则我认为不需要 1000 倍。

Paul Tomblin 指出 Aviation Formulary v1.44 是方程式的来源 - 事实上，它就在那里，还有一个数值更稳定的版本，适用于位置差异很小的情况。

转到基本的三角学，我们也知道：

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

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

（我的参考："Astronomy: Principles and Practice, Fourth Edition" A E Roy 和 D Clarke（2003 年）；我的副本是 1977 年的第一版，Adam Hilger，ISBN 0-85274-346-7。）

---

注意查看（谷歌）'define：“海里”'；根据定义，现在一海里似乎是 1852 米（1.852 公里）。乘数 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°的等边球面三角形，这是完全可行的）。

---

【讨论】：

1.1515 是法定英里到海里的转换。

【参考方案2】：

查看http://www.movable-type.co.uk/scripts/latlong.html

该网站有很多不同的公式和 Javascript 代码可以帮助您。我已经成功地将它翻译成 C# 和 SQL Server UDF，而且我到处都在使用它们。

例如 Javascript 计算距离：

var R = 6371; // km
var φ1 = lat1.toRadians();
var φ2 = lat2.toRadians();
var Δφ = (lat2-lat1).toRadians();
var Δλ = (lon2-lon1).toRadians();

var a = Math.sin(Δφ/2) * Math.sin(Δφ/2) +
        Math.cos(φ1) * Math.cos(φ2) *
        Math.sin(Δλ/2) * Math.sin(Δλ/2);
var c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a));

var d = R * c; 

享受吧！

【参考方案3】：

您在公里和弧度之间的转换是错误的。海里是 1/60 度，所以假设 1.15... 是您从英里到海里的转换，1.6... 是您从公里到法定英里的转换，

   nm = km /  (1.1515 * 1.609344);
   deg = nm / 60;
   rad = toRadians(deg);

换句话说，我认为你的差距是 1000 倍。

【讨论】：

呃。我不知道我在用距离方法在想什么；返回以米为单位的距离，而不是公里。

【参考方案4】：

关于您更新的问题：不应该

double lon1 = endOfCourse(lat1, 0.0, 90.0, distEast)[0];

成为

double lon1 = endOfCourse(lat1, 0.0, 90.0, distEast)[1];

【参考方案5】：

除了其他答案和更新中提到的实现错误之外，我发现了这些公式的大问题。

最大的问题是：距离方法（用于计算两点之间的距离）是计算大圆距离。当然，这是有道理的——这是两点之间的最短路径。 但是，位于同一平行线（纬线）上的两点之间的大圆距离与直接沿纬线行驶时这两点之间的距离不同，除非您'在赤道。

所以：功能正常工作；但是，我在原始问题中提出的替代坐标系要求我们仅查看沿 IDL 向北的距离，然后查看沿所得纬度平行线向东的距离。而且计算沿特定平行线的距离与计算沿大圆的距离完全不同！

不管怎样，你有它。