给定一个点的经纬度、距离和方位角，如何获取新的经纬度？

Question

我在网上找到了一段代码。它可以根据给定的经纬度点和距离计算最小边界矩形。

private static void GetlatLon(double LAT, double LON, double distance, double angle, out double newLon, out double newLat)
        {     
            double dx = distance * 1000 * Math.Sin(angle * Math.PI / 180.0);
            double dy = distance * 1000 * Math.Cos(angle * Math.PI / 180.0);
            double ec = 6356725 + 21412 * (90.0 - LAT) / 90.0;
            double ed = ec * Math.Cos(LAT * Math.PI / 180);
            newLon = (dx / ed + LON * Math.PI / 180.0) * 180.0 / Math.PI;
            newLat = (dy / ec + LAT * Math.PI / 180.0) * 180.0 / Math.PI;
            
            }
            
    public static void GetRectRange(double centorlatitude, double centorLogitude, double distance,
                                  out double maxLatitude, out double minLatitude, out double maxLongitude, out double minLongitude)
{       
            GetlatLon(centorlatitude, centorLogitude, distance, 0, out temp, out maxLatitude);
            GetlatLon(centorlatitude, centorLogitude, distance, 180, out temp, out minLatitude);
            GetlatLon(centorlatitude, centorLogitude, distance, 90, out minLongitude, out temp);
            GetlatLon(centorlatitude, centorLogitude, distance, 270, out maxLongitude, out temp);
}



double ec = 6356725 + 21412 * (90.0 - LAT) / 90.0; //why?
double ed = ec * Math.Cos(LAT * Math.PI / 180);    // why?
dx / ed                                            //why?
dy / ec                                            //why?

6378137是赤道半径，6356725是极半径，21412=6378137-6356725。

从这个链接中，我知道了一些含义。但是这四行，我不知道为什么。请帮助提供更多信息。您能帮忙解释一下公式的推导过程吗？

从这个链接中，在“起始点给定距离和方位角的目标点”一节中，它给出了另一个公式来得出结果。这个公式的推导过程是怎样的呢？

从这个链接中，我了解了Haversine公式的推导过程，非常有用。我认为“起始点给定距离和方位角的目标点”一节中的公式并不只是Haversine公式的简单反转。

非常感谢！

- gfan

---

你有原始代码的链接吗？ - tmthydvnprt

double ec = 6356725 + 21412 * (90.0 - LAT) / 90.0 看起来像是一种加权平均值，用于估算赤道和极地之间的地球半径。因此，如果您在极点上，只需使用半径=6356725，如果您在赤道上，则半径=6378137（请注意，在赤道上更大）。因此，当您处于两者之间时，需要估算半径。 - Scott

double ed = ec * Math.Cos(LAT * Math.PI / 180) 将估计的地球半径乘以给定纬度（转换为弧度）的余弦值。因此，它只是一个校正项。由于 dx 和 dy 是在平面上计算的距离，因此 ec 旨在校正地球的曲率。 - Scott

Answer 1

这是一个典型例子，说明为什么在代码中添加注释可以使其更易读和可维护。从数学上讲，您要看的是：

double ec = 6356725 + 21412 * (90.0 - LAT) / 90.0; //为什么？

这是一种衡量偏心率以某种方式解释赤道膨胀的方法。21412 是赤道和极地之间的地球半径差异，6356725 是极半径。(90.0 - LAT) / 90.0 在赤道时为 1，在极地时为 0。该公式估计了在任何给定纬度出现的膨胀量。

double ed = ec * Math.Cos(LAT * Math.PI / 180); // 为什么？

(LAT * Math.PI / 180) 是纬度从度转换为弧度的值。cos(0) = 1，cos(1) = 0，因此在赤道上，您正在应用完整的偏心率，在极地上则没有应用任何内容。与前一行类似。

dx / ed //为什么？

dy / ec //为什么？

以上似乎是与任何给定纬度/经度中的凸出有关的距离分数增量，这些增量用于在newLon newLat计算中使用，以到达新位置。

我没有对您发现的代码片段进行任何研究，但从数学上讲，此处正在发生的情况。希望这会引导您朝正确的方向前进。

---

C 中的 Haversine 示例

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

double m2ft (double l) {            /* convert meters to feet       */
    return l/(1200.0/3937.0);
}

double ft2smi (double l) {          /* convert feet to statute miles*/
    return l/5280.0;
}

double km2smi (double l) {          /* convert km to statute mi.    */
    return ft2smi(m2ft( l * 1000.0 ));
}

static const double deg2rad = 0.017453292519943295769236907684886;
static const double earth_rad_m = 6372797.560856;

typedef struct pointd {
    double lat;
    double lon;
} pointd;

/* Computes the arc, in radian, between two WGS-84 positions.
   The result is equal to Distance(from,to)/earth_rad_m
      = 2*asin(sqrt(h(d/earth_rad_m )))
   where:
      d is the distance in meters between 'from' and 'to' positions.
      h is the haversine function: h(x)=sin²(x/2)
   The haversine formula gives:
      h(d/R) = h(from.lat-to.lat)+h(from.lon-to.lon)+cos(from.lat)*cos(to.lat)
   http://en.wikipedia.org/wiki/Law_of_haversines
 */
double arcradians (const pointd *from, const pointd *to)
{
    double latitudeArc  = (from-> lat - to-> lat) * deg2rad;
    double longitudeArc = (from-> lon - to-> lon) * deg2rad;

    double latitudeH = sin (latitudeArc * 0.5);
    latitudeH *= latitudeH;

    double lontitudeH = sin (longitudeArc * 0.5);
    lontitudeH *= lontitudeH;

    double tmp = cos (from-> lat * deg2rad) * cos (to-> lat * deg2rad);

    return 2.0 * asin (sqrt (latitudeH + tmp*lontitudeH));
}

/* Computes the distance, in meters, between two WGS-84 positions.
   The result is equal to earth_rad_m*ArcInRadians(from,to)
 */
double dist_m (const pointd *from, const pointd *to) {
    return earth_rad_m * arcradians (from, to);
}

int main (int argc, char **argv) {

    if (argc < 5 ) {
        fprintf (stderr, "Error: insufficient input, usage: %s (lat,lon) (lat,lon)\n", argv[0]);
        return 1;
    }

    pointd points[2];

    points[0].lat = strtod (argv[1], NULL);
    points[0].lon = strtod (argv[2], NULL);

    points[1].lat = strtod (argv[3], NULL);
    points[1].lon = strtod (argv[4], NULL);

    printf ("\nThe distance in meters from 1 to 2 (smi): %lf\n\n", km2smi (dist_m (&points[0], &points[1])/1000.0) );

    return 0;
}

/* Results/Example.
    ./bin/gce 31.77 -94.61 31.44 -94.698
    The distance in miles from Nacogdoches to Lufkin, Texas (smi): 23.387997 miles
 */

- David C. Rankin

---

谢谢David。在链接(http://www.movable-type.co.uk/scripts/latlong.html)的“给定起点距离和方位角的终点”部分，它提供了另一种公式来得出结果。这个公式的推导是什么？是否有像链接(http://mathforum.org/library/drmath/view/51879.html)那样清晰地演示Haversine公式的网页？ - gfan

Haversine公式还不错。我认为我已经将其移植到了C语言中。基本上，在所有的纬度/经度问题中，你都在处理大圆距离（这只是一个指地球表面上球形坐标系的花哨术语）。这一切源于一个事实，即在球体上，点a和点b之间的连线根本不是一条直线，而是一条弧线。因此，你要找到的距离是从a到b的弧长距离。archlength = radius * angle。这就是为什么你会看到三角函数，而r由于凸起而发生变化，所以你会看到ec、ed微调的原因。 - David C. Rankin

同样的公式也适用于“从起点出发，给定距离和方位角度的终点”问题。不同的是，你的初始条件不再是point1和point2，而是point1和北向夹角及弧长。在这里，你已经知道了距离，所以不需要再去计算距离，而是利用已知的距离来求解point2。（我已经很多年没有推导球面坐标方程了，而且以后也不会再做了，祝你好运:p） - David C. Rankin

Answer 2

我猜6356725与地球半径有关。看看这个答案，还要了解一下Haversine公式。

- tmthydvnprt