如何将 lon/lat 坐标转换为地球表面上 N-E 米的距离？

Posted 2023-02-23

【中文标题】如何将 lon/lat 坐标转换为地球表面上 N-E 米的距离？

【英文标题】：how do I translate of lon/lat coordinate by some N-E meters distance on earth surface?

【发布时间】：2011-01-13 12:14:45

【问题描述】：

如何在地球表面上进行一些平移（以米为单位）后，从参考点（大地测量）获得新的大地测量坐标（纬度/经度），并且我还需要使用真正的地球椭球体进行计算WGS84等模型。

例如：

假设我的参考点是 10.32E，-4.31N 然后我进行 (3000,-2000) 米的平移（即将点向东移动 3000 米，向南移动 2000 米

在地球表面

。 那么我需要新点的坐标

大地测量

。

谢谢

【问题讨论】：

你说的东移和南移，真的是指沿经纬线移动吗？或者您的意思是沿着其他一些网格的网格线移动，例如英国使用的军械测量网格。

它沿着纬度/经度线。这就是为什么我称它为北/南 - 东/西，但我不需要以球坐标的度数单位移动，而是以地球表面的笛卡尔米为单位移动，我需要再次以球坐标表示结果（不是球，它应该是椭圆）

按什么顺序？先东后南，还是先南后东？有区别。

先东/西再北/南，注意：“米”距离不是直线上的距离，而是“弧长”，因为它是球面上的距离。

【参考方案1】：

查看开源库PROJ.4，您可以使用它来准确地将地理坐标（纬度/经度）转换为投影坐标（米），然后再转换回来。在您的情况下，您可以投影到 WGS 84 / World Mercator (EPSG:3395)，以米为单位执行转换，然后取消投影回地理。

【讨论】：

目前我正在这样做，即将大地（纬度/经度）的当前参考点转换为 UTM（通用横向墨卡托）并使用 UTM 坐标进行移动，然后从 UTM 转换回新位置进入大地测量学，但我讨厌 UTM 的一件事是他们有区域/球体枚举，我正在寻找没有任何分区的墨卡托投影，你知道吗？

WGS84 世界墨卡托与 UTM 不同。 WGS84 世界墨卡托的预计边界为：-20037508.3428、-15496570.7397、20037508.3428、18764656.2314。这涵盖了整个世界。

【参考方案2】：

找到答案：

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

来自： Vincenty 直接公式 - T Vincenty，“测地线的正解和逆解 Ellipsoid with application of nested equations”，Survey Review, vol XXII no 176, 1975

http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf

【参考方案3】：

此代码计算给定纬度/经度坐标的两点之间的距离（N 和 E）。您可以轻松地将其反转为您的目的。

看看函数

u8 GPS_CalculateDeviation()

在

http://svn.mikrokopter.de/filedetails.php?repname=NaviCtrl&path=/tags/V0.15c/GPS.c

【讨论】：

如果是关于获取两个纬度/经度坐标之间的距离，有更好的文章和代码：movable-type.co.uk/scripts/latlong-vincenty.html，但扭转它不是一件容易的事，这就是我在这里问的原因跨度>

@uray：即使（可能很昂贵）数值反演也不行？

@alex：好吧，我不是测地线数学专家，因为它涉及椭球面，如果它是完美的球面，也许会更简单

@uray：它不依赖于底层模型，你有一个函数可以获取一个位置并给你一个位移。你用例如反转它。具有有限差分导数的牛顿法。

使用这个很简单。只需获取您所在位置以北 1 度的弧长，将您的平移北除以得到纬度的变化。然后从您的新位置获得向东一度的弧长，然后将您的东平移除以该数字以获得您的经度变化。

【参考方案4】：

你要么找一些地理图书馆，要么自己做三角学。

在任何情况下，您都应该更准确地表述您的问题。你特别说：

然后我进行 (3000,-2000) 米的平移（即在地球表面将点向东移动 3000 米，向南移动 2000 米。

您应该注意，向东移动 3 公里然后向南移动 2 公里不同与向南移动 2 公里然后向东移动 3 公里。这些不是交换动作。因此，通过偏移 (3000, -2000) 调用它是不正确的。

【讨论】：

“或者自己做三角学”这不是一项简单的任务：地球表面不是一个球体。

不，因为 Raedwald 表示这不是数学。信息不足。该问题在指出需要 WGS-84 模型时已经承认这一点。

@Raedwald：我没有说这是一项微不足道的任务。然而，它仍然存在于解析几何领域。有了适当的数学背景，这样做没有问题

【参考方案5】：

以下是从original version from ETH Zurich 略微修改的 C++ 代码。该文件仅依赖于Eigen library（如果需要通过自己编写矩阵乘法函数，可以通过一些琐碎的工作来消除）。您可以使用 geodetic2ned() 函数将纬度、经度、高度转换为 NED 帧。

//GeodeticConverter.hpp
#ifndef air_GeodeticConverter_hpp
#define air_GeodeticConverter_hpp

#include <math>
#include <eigen3/Eigen/Dense>

class GeodeticConverter

 public:
  GeodeticConverter(double home_latitude = 0, double home_longitude = 0, double home_altitude = 0)
    : home_latitude_(home_latitude), home_longitude_(home_longitude)
  
    // Save NED origin
    home_latitude_rad_ = deg2Rad(latitude);
    home_longitude_rad_ = deg2Rad(longitude);
    home_altitude_ = altitude;

    // Compute ECEF of NED origin
    geodetic2Ecef(latitude, longitude, altitude, &home_ecef_x_, &home_ecef_y_, &home_ecef_z_);

    // Compute ECEF to NED and NED to ECEF matrices
    double phiP = atan2(home_ecef_z_, sqrt(pow(home_ecef_x_, 2) + pow(home_ecef_y_, 2)));

    ecef_to_ned_matrix_ = nRe(phiP, home_longitude_rad_);
    ned_to_ecef_matrix_ = nRe(home_latitude_rad_, home_longitude_rad_).transpose();    
  

  void getHome(double* latitude, double* longitude, double* altitude)
  
    *latitude = home_latitude_;
    *longitude = home_longitude_;
    *altitude = home_altitude_;
  

  void geodetic2Ecef(const double latitude, const double longitude, const double altitude, double* x,
                     double* y, double* z)
  
    // Convert geodetic coordinates to ECEF.
    // http://code.google.com/p/pysatel/source/browse/trunk/coord.py?r=22
    double lat_rad = deg2Rad(latitude);
    double lon_rad = deg2Rad(longitude);
    double xi = sqrt(1 - kFirstEccentricitySquared * sin(lat_rad) * sin(lat_rad));
    *x = (kSemimajorAxis / xi + altitude) * cos(lat_rad) * cos(lon_rad);
    *y = (kSemimajorAxis / xi + altitude) * cos(lat_rad) * sin(lon_rad);
    *z = (kSemimajorAxis / xi * (1 - kFirstEccentricitySquared) + altitude) * sin(lat_rad);
  

  void ecef2Geodetic(const double x, const double y, const double z, double* latitude,
                     double* longitude, double* altitude)
  
    // Convert ECEF coordinates to geodetic coordinates.
    // J. Zhu, "Conversion of Earth-centered Earth-fixed coordinates
    // to geodetic coordinates," IEEE Transactions on Aerospace and
    // Electronic Systems, vol. 30, pp. 957-961, 1994.

    double r = sqrt(x * x + y * y);
    double Esq = kSemimajorAxis * kSemimajorAxis - kSemiminorAxis * kSemiminorAxis;
    double F = 54 * kSemiminorAxis * kSemiminorAxis * z * z;
    double G = r * r + (1 - kFirstEccentricitySquared) * z * z - kFirstEccentricitySquared * Esq;
    double C = (kFirstEccentricitySquared * kFirstEccentricitySquared * F * r * r) / pow(G, 3);
    double S = cbrt(1 + C + sqrt(C * C + 2 * C));
    double P = F / (3 * pow((S + 1 / S + 1), 2) * G * G);
    double Q = sqrt(1 + 2 * kFirstEccentricitySquared * kFirstEccentricitySquared * P);
    double r_0 = -(P * kFirstEccentricitySquared * r) / (1 + Q)
        + sqrt(
            0.5 * kSemimajorAxis * kSemimajorAxis * (1 + 1.0 / Q)
                - P * (1 - kFirstEccentricitySquared) * z * z / (Q * (1 + Q)) - 0.5 * P * r * r);
    double U = sqrt(pow((r - kFirstEccentricitySquared * r_0), 2) + z * z);
    double V = sqrt(
        pow((r - kFirstEccentricitySquared * r_0), 2) + (1 - kFirstEccentricitySquared) * z * z);
    double Z_0 = kSemiminorAxis * kSemiminorAxis * z / (kSemimajorAxis * V);
    *altitude = U * (1 - kSemiminorAxis * kSemiminorAxis / (kSemimajorAxis * V));
    *latitude = rad2Deg(atan((z + kSecondEccentricitySquared * Z_0) / r));
    *longitude = rad2Deg(atan2(y, x));
  

  void ecef2Ned(const double x, const double y, const double z, double* north, double* east,
                double* down)
  
    // Converts ECEF coordinate position into local-tangent-plane NED.
    // Coordinates relative to given ECEF coordinate frame.

    Vector3d vect, ret;
    vect(0) = x - home_ecef_x_;
    vect(1) = y - home_ecef_y_;
    vect(2) = z - home_ecef_z_;
    ret = ecef_to_ned_matrix_ * vect;
    *north = ret(0);
    *east = ret(1);
    *down = -ret(2);
  

  void ned2Ecef(const double north, const double east, const double down, double* x, double* y,
                double* z)
  
    // NED (north/east/down) to ECEF coordinates
    Vector3d ned, ret;
    ned(0) = north;
    ned(1) = east;
    ned(2) = -down;
    ret = ned_to_ecef_matrix_ * ned;
    *x = ret(0) + home_ecef_x_;
    *y = ret(1) + home_ecef_y_;
    *z = ret(2) + home_ecef_z_;
  

  void geodetic2Ned(const double latitude, const double longitude, const double altitude,
                    double* north, double* east, double* down)
  
    // Geodetic position to local NED frame
    double x, y, z;
    geodetic2Ecef(latitude, longitude, altitude, &x, &y, &z);
    ecef2Ned(x, y, z, north, east, down);
  

  void ned2Geodetic(const double north, const double east, const double down, double* latitude,
                    double* longitude, double* altitude)
  
    // Local NED position to geodetic coordinates
    double x, y, z;
    ned2Ecef(north, east, down, &x, &y, &z);
    ecef2Geodetic(x, y, z, latitude, longitude, altitude);
  

  void geodetic2Enu(const double latitude, const double longitude, const double altitude,
                    double* east, double* north, double* up)
  
    // Geodetic position to local ENU frame
    double x, y, z;
    geodetic2Ecef(latitude, longitude, altitude, &x, &y, &z);

    double aux_north, aux_east, aux_down;
    ecef2Ned(x, y, z, &aux_north, &aux_east, &aux_down);

    *east = aux_east;
    *north = aux_north;
    *up = -aux_down;
  

  void enu2Geodetic(const double east, const double north, const double up, double* latitude,
                    double* longitude, double* altitude)
  
    // Local ENU position to geodetic coordinates

    const double aux_north = north;
    const double aux_east = east;
    const double aux_down = -up;
    double x, y, z;
    ned2Ecef(aux_north, aux_east, aux_down, &x, &y, &z);
    ecef2Geodetic(x, y, z, latitude, longitude, altitude);
  

private:
  // Geodetic system parameters
  static double kSemimajorAxis = 6378137;
  static double kSemiminorAxis = 6356752.3142;
  static double kFirstEccentricitySquared = 6.69437999014 * 0.001;
  static double kSecondEccentricitySquared = 6.73949674228 * 0.001;
  static double kFlattening = 1 / 298.257223563;

  typedef Eigen::Vector3d Vector3d;
  typedef Eigen::Matrix<double, 3, 3> Matrix3x3d;

  inline Matrix3x3d nRe(const double lat_radians, const double lon_radians)
  
    const double sLat = sin(lat_radians);
    const double sLon = sin(lon_radians);
    const double cLat = cos(lat_radians);
    const double cLon = cos(lon_radians);

    Matrix3x3d ret;
    ret(0, 0) = -sLat * cLon;
    ret(0, 1) = -sLat * sLon;
    ret(0, 2) = cLat;
    ret(1, 0) = -sLon;
    ret(1, 1) = cLon;
    ret(1, 2) = 0.0;
    ret(2, 0) = cLat * cLon;
    ret(2, 1) = cLat * sLon;
    ret(2, 2) = sLat;

    return ret;
  

  inline double rad2Deg(const double radians)
  
    return (radians / M_PI) * 180.0;
  

  inline double deg2Rad(const double degrees)
  
    return (degrees / 180.0) * M_PI;
  

  double home_latitude_rad_, home_latitude_;
  double home_longitude_rad_, home_longitude_;
  double home_altitude_;

  double home_ecef_x_;
  double home_ecef_y_;
  double home_ecef_z_;

  Matrix3x3d ecef_to_ned_matrix_;
  Matrix3x3d ned_to_ecef_matrix_;

; // class GeodeticConverter

#endif