地理空间坐标和公里数距离

14

这是对这个问题的后续。

我好像陷入了困境。基本上,我需要能够在标准度系统和通过沿着国际日期线测量南极点向北的距离以及从该点开始向东的距离之间来回转换坐标。为此(以及进行一些更通用的测距),我有一个方法来确定两个lat/lon点之间的距离,以及另一个方法,它接受一个lat/lon点,一个方向和一个距离,并返回该课程结束时的lat/lon点。

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

/* 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;
}
“结束点”应该是约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;

1
作为参考资料,您可能需要查看 http://www.movable-type.co.uk/scripts/latlong-vincenty.html 和 http://www.movable-type.co.uk/scripts/latlong-vincenty-direct.html(包括 JavaScript 实现)。 - Brendan Cashman
在你的第一个函数中,最好使用两个变量,“double angle = ...1st expression; double dist = ...2nd expression modified to reference angle;”。实际上,第一个值根本不是距离。这只是吹毛求疵的人。 - Jonathan Leffler
5个回答

64

你的代码中存在大量魔法数字。该表达式:

 (60 * 1.1515 * 1.609344 * 1000)

这段内容中出现了一些数字,下面为您进行解释:1.609344是一英里对应的公里数;60是一度对应的分钟数;1000是一公里对应的米数;1.1515是一个海里对应的英里数(感谢DanM)。一个海里的长度等于赤道上一纬度的长度。

我猜测你使用的是球体地球模型,而非椭圆地球模型?因为这个代数式并不足够复杂以应用于椭圆地球。

第一个公式——两组纬度和经度之间的转换——有点奇怪。您需要同时知道Δλ和Δφ才能得到答案。此外,这两组坐标之间的距离:

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

应该是相同的-但我非常确定你的代码产生了不同的答案。

因此,我认为您需要回到您的球面三角学参考资料中,看看您做错了什么。(我需要一段时间才能找到关于该主题的书籍,它需要从箱子里拆开.)

给定一个球面三角形,顶点处的角度为ABC,对应顶点的边长分别为abc(也就是说,边长a是从BC,等等),余弦公式为:

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

将此应用于问题,我们可以称给定的两个点为BC,并创建一个在A处呈直角的右球面三角形。

糟糕的 ASCII 艺术:

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

a的方程为:a = sqrt((b * 111.12)^2 + (c * 111.12 * cos(lat1))^2)

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

在我给出的方程中应用两次这个公式可能会得到Aviation Formulary中的公式。

(我的参考资料:《天文学:原理与实践,第四版》,作者为A E Roy和D Clarke(2003年),我手头的是1977年的第一版,Adam Hilger出版,ISBN 0-85274-346-7。)


NB 查看(谷歌)“define:‘nautical mile’”;现在一个海里的定义是1852米(1.852公里)。乘数1.1515对应于旧的海里定义,大约为6080英尺。使用带有10为小数精度的bc,我得到:

(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,以及AB,我们可以应用余弦公式来得到c。请注意,我们不能简单地重新应用正弦公式来获得c,因为我们没有C的值,并且因为我们正在处理球面三角学,所以没有方便的规则表明C = 90° - B(球面三角形的角度之和可能大于180°;考虑一个等边球面三角形,其所有角度都等于90°,这是完全可行的)。

1.1515是从英里转换为海里的系数。 - DanM

6

请查看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.cos1) * Math.cos2) *
        Math.sin(Δλ/2) * Math.sin(Δλ/2);
var c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a));

var d = R * c; 

尽情享受吧!


2

你的公里和弧度之间的转换是错误的。一个海里是一度的1/60,所以假设1.15是你从英里到海里的转换,而1.6是你从公里到英里的转换。

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

换句话说,我认为你的数据偏差了1000倍。

我真不知道当时 Distance 方法是怎么想的,它返回的是米而非千米。 - DanM

1
关于您更新的问题:不应该。
double lon1 = endOfCourse(lat1, 0.0, 90.0, distEast)[0];

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

0

我找出了这些公式的一个大问题,除了其他答案和更新中提到的实现错误。

大问题在于:距离方法(用于计算两点之间的距离)正在计算大圆距离。当然,这是有道理的 - 那是两点之间的最短路径。然而,同一纬线上的两点之间的大圆距离与沿着纬线直接行驶时这两点之间的距离不相同,除非你在赤道。

因此:这些函数是正确的;但是,我在原始问题中提出的替代坐标系统要求我们仅查看沿着IDL向北的距离,然后沿着得到的纬度上的平行线向东的距离。计算沿特定纬线的距离与计算沿大圆的距离非常不同!

总之,就是这样。


网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接