查找地球坐标（纬度、经度）、距离（米）和方位（角度）

Question

我需要以各种方式处理地球坐标。 C/C++ 中没有函数可以直接执行此操作。
参考以下问题：

从第一个和movable type scripts website，我发现下面是公式：

求2个坐标之间的方位角（角度）

x = cos(lat1Rad)*sin(lat2Rad) - sin(lat1Rad)*cos(lat2Rad)*cos(lon2Rad-lon1Rad);
y = sin(lon2Rad-lon1Rad) * cos(lat2Rad);
bearing = atan2(y, x);  // In radians; 
// Convert to degrees and for -ve add 360

求2个坐标之间的距离（米）

PI = 3.14159265358979323846, earthDiameterMeters = 2*6371*1000;
x = sin((lat2Rad-lat1Rad) / 2);
y = sin((lon2Rad-lon1Rad) / 2);
meters = earthDiameterMeters * asin(sqrt(x*x + y*y*cos(lat1Rad)*cos(lat2Rad)));

坐标+距离+角度求坐标

meters *= 2 / earthDiameterMeters;
lat2Rad = asin(sin(lat1Rad)*cos(meters) + cos(lat1Rad)*sin(meters)*cos(bearing));
lon2Rad = lon1Rad + atan2(sin(bearing)*sin(meters)*cos(lat1Rad), 
                          cos(meters) - sin(lat1Rad)*sin(lat2Rad));

下面的伪代码应该相互验证以上 3 个等式：

struct Coordinate { double lat, lon; } c1, c2;  
auto degree = FindBearing(c1, c2);
auto meters = FindDistance(c1, c2);
auto cX = FindCoordiante(c1, degree, meters);

现在实际上答案几乎接近但不正确。即 cX 不等于 c2!
经度值总是相差0.0005。例如

c1 = (12.968460,77.641308)  
c2 = (12.967862,77.653130)  
angle = 92.97         ^^^
distance = 1282.74  
cX = (12.967862,77.653613)
                      ^^^

我的数学知识不多'Havesine Forumla. But what I know is that from the website fcc.gov，答案总是正确的。

我做错了什么？

Code仅供参考

虽然语法是在 C++ 中，但所有数学函数都来自 C，并且在 C 中也很容易移植（因此标记为两者）

#include<iostream>
#include<iomanip>
#include<cmath>


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

static const double PI = 3.14159265358979323846, earthDiameterMeters = 6371.0 * 2 * 1000;

double degreeToRadian (const double degree) { return (degree * PI / 180); };
double radianToDegree (const double radian) { return (radian * 180 / PI); };

double CoordinatesToAngle (double latitude1,
                           const double longitude1,
                           double latitude2,
                           const double longitude2)
{
  const auto longitudeDifference = degreeToRadian(longitude2 - longitude1);
  latitude1 = degreeToRadian(latitude1);
  latitude2 = degreeToRadian(latitude2);

  using namespace std;
  const auto x = (cos(latitude1) * sin(latitude2)) -
                 (sin(latitude1) * cos(latitude2) * cos(longitudeDifference));
  const auto y = sin(longitudeDifference) * cos(latitude2);

  const auto degree = radianToDegree(atan2(y, x));
  return (degree >= 0)? degree : (degree + 360);
}

double CoordinatesToMeters (double latitude1,
                            double longitude1,
                            double latitude2,
                            double longitude2)
{
  latitude1 = degreeToRadian(latitude1);
  longitude1 = degreeToRadian(longitude1);
  latitude2 = degreeToRadian(latitude2);
  longitude2 = degreeToRadian(longitude2);

  using namespace std;
  auto x = sin((latitude2 - latitude1) / 2), y = sin((longitude2 - longitude1) / 2);
#if 1
  return earthDiameterMeters * asin(sqrt((x * x) + (cos(latitude1) * cos(latitude2) * y * y)));
#else
  auto value = (x * x) + (cos(latitude1) * cos(latitude2) * y * y);
  return earthDiameterMeters * atan2(sqrt(value), sqrt(1 - value));
#endif
}

std::pair<double,double> CoordinateToCoordinate (double latitude,
                                                 double longitude,
                                                 double angle,
                                                 double meters)
{
  latitude = degreeToRadian(latitude);
  longitude = degreeToRadian(longitude);
  angle = degreeToRadian(angle);
  meters *= 2 / earthDiameterMeters;

  using namespace std;
  pair<double,double> coordinate;

  coordinate.first = radianToDegree(asin((sin(latitude) * cos(meters))
                             + (cos(latitude) * sin(meters) * cos(angle))));
  coordinate.second = radianToDegree(longitude
                    + atan2((sin(angle) * sin(meters) * cos(latitude)),
                    cos(meters) - (sin(latitude) * sin(coordinate.first))));

  return coordinate;
}

int main ()
{
  using namespace std;
  const auto latitude1 = 12.968460, longitude1 = 77.641308,
             latitude2 = 12.967862, longitude2 = 77.653130;

  cout << std::setprecision(10);
  cout << "(" << latitude1 << "," << longitude1 << ") --- "
          "(" << latitude2 << "," << longitude2 << ")\n";

  auto angle = CoordinatesToAngle(latitude1, longitude1, latitude2, longitude2);
  cout << "Angle =  " << angle << endl;

  auto meters = CoordinatesToMeters(latitude1, longitude1, latitude2, longitude2);
  cout << "Meters = " << meters << endl;

  auto coordinate = CoordinateToCoordinate(latitude1, longitude1, angle, meters);
  cout << "Destination = (" << coordinate.first << "," << coordinate.second << ")\n";
}

Answer 1

部分回答

角度 92.97° 然后转换为弧度，调用 sin/cos/tan 将有效地将角度更改为 2.97°。由于周期减少发生在度到弧度转换后和在三角函数调用中，仅此步骤就失去了 6 位精度。

可以提高大角度三角函数的精度。使用幸运相位，一圈有正好 360.0 度。执行 "modulo 45°"，可能使用 remquo(angle, 45, &octant)，然后然后转换为弧度，然后再使用弧度参数调用三角函数。

示例

您对 77.641308 和 77.653130 的回答大约相差 6600 分之一（~13 位精度）。 This answer可能无法完全解释这一点，但应该有所帮助。（如果在某处出现了 float 的某些用法，则应设为 double。）

Answer 2

在 CoordinateToCoordinate 中，您使用的 sin(coordinate.first) 已经是度数了。使用 sin(degreeToRadian(coordinate.first))。

或者更清洁：

... CoordinateToCoordinate (...)
{
  ...
  coordinate.first = asin((sin(latitude) * cos(meters))
                        + (cos(latitude) * sin(meters) * cos(angle)));
  coordinate.second = longitude + atan2((sin(angle) * sin(meters) * cos(latitude)), 
         cos(meters) - (sin(latitude) * sin(coordinate.first)));

  coordinate.first = radianToDegree(coordinate.first);
  coordinate.second = radianToDegree(coordinate.second);

  return coordinate;
}

这解决了问题。 Live Demo.

查找地球坐标（纬度、经度）、距离（米）和方位（角度）

Finding Earth Coordinates (latitude,longitude), Distance (meters) and Bearing (angle)

c

c++

geometry

coordinates

earthdistance

Code仅供参考