WGS84、BD09、GCJ02坐标转换

名词解释

WGS84

此坐标系解释参考笔者另一篇博客GIS坐标系、投影与转换

GCJ02

GCJ-02是由中国国家测绘局（G表示Guojia国家，C表示Cehui测绘，J表示Ju局）制订的地理信息系统的坐标系统。 是中国大陆地区的地图数据使用的坐标系。基于 WGS84 进行加密后形成。

BD09

BD09 是中国的百度地图使用的一种坐标系，全称是百度坐标系（BD-09）。它是在 GCJ-02（国测局坐标系）基础上进行的进一步加密处理。百度地图使用 BD-09 坐标系，这种坐标系主要用于在百度地图上的定位和导航。

坐标转换

WGS84转GCJ02

需要经过一个偏移算法

double offset_gcj02_latitude(double x, double y) {
    double ret = -100.0 + 2.0 * x + 3.0 * y + 0.2 * y * y + 0.1 * x * y + 0.2 * std::sqrt(abs(x));
    ret += (20.0 * std::sin(6.0 * x * pi) + 20.0 * std::sin(2.0 * x * pi)) * 2.0 / 3.0;
    ret += (20.0 * std::sin(y * pi) + 40.0 * std::sin(y / 3.0 * pi)) * 2.0 / 3.0;
    ret += (160.0 * std::sin(y / 12.0 * pi) + 320 * std::sin(y * pi / 30.0)) * 2.0 / 3.0;
    return ret;
}

double offset_gcj02_longitude(double x, double y) {
    double ret = 300.0 + x + 2.0 * y + 0.1 * x * x + 0.1 * x * y + 0.1 * std::sqrt(abs(x));
    ret += (20.0 * std::sin(6.0 * x * pi) + 20.0 * std::sin(2.0 * x * pi)) * 2.0 / 3.0;
    ret += (20.0 * std::sin(x * pi) + 40.0 * std::sin(x / 3.0 * pi)) * 2.0 / 3.0;
    ret += (150.0 * std::sin(x / 12.0 * pi) + 300.0 * std::sin(x / 30.0 * pi)) * 2.0 / 3.0;
    return ret;
}

void gcj02_delta(double lat, double lon, double& dLat, double& dLon) {
    const double a = 6378245.0;                 // 长半轴
    const double ee = 0.00669342162296594323;   // 扁率

    double radLat = lat / 180.0 * pi;
    double magic = std::sin(radLat);
    magic = 1 - ee * magic * magic;
    double sqrtMagic = std::sqrt(magic);
    dLat = (offset_gcj02_latitude(lon - 105.0, lat - 35.0) * 180.0) / ((a * (1 - ee)) / (magic * sqrtMagic) * pi);
    dLon = (offset_gcj02_longitude(lon - 105.0, lat - 35.0) * 180.0) / (a / sqrtMagic * std::cos(radLat) * pi);
}

void wgs84_to_gcj02(double wgsLat, double wgsLon, double& gcjLat, double& gcjLon) {
    if (outOfChina(wgsLat, wgsLon)) {
        gcjLat = wgsLat;
        gcjLon = wgsLon;
        return;
    }

    double lat_offset=0, lon_offset=0;
    gcj02_delta(wgsLat, wgsLon, lat_offset, lon_offset);
    gcjLat = wgsLat + lat_offset;
    gcjLon = wgsLon + lon_offset;
}

如上述代码所示， 其转换仅需在WGS84坐标上加上对应偏移即可。 因为是中国国测局制定的， 所以在转换开始判断坐标是否位于中国大陆境内。

逆转的话， 减去对应的偏移即可

WGS84转BD09

WGS84转为BD09需首先将其转换为GCJ02, 然后再在GCJ02基础上作加密偏移。以下代码仅仅示例了如何从GCJ02转换到BD09

void gcj02_to_bd09(double gcjLat, double gcjLon, double& bdLat, double& bdLon) {
    double z = std::sqrt(gcjLon * gcjLon + gcjLat * gcjLat) + 0.00002 * std::sin(gcjLat * pi);
    double theta = std::atan2(gcjLat, gcjLon) + 0.000003 * std::cos(gcjLon * pi);
    bdLon = z * std::cos(theta) + 0.0065;
    bdLat = z * std::sin(theta) + 0.006;
}

其逆转也是一个反向过程

void bd09_to_gcj02(double bdLat, double bdLon, double& gcjLat, double& gcjLon) {
    double x = bdLon - 0.0065;
    double y = bdLat - 0.006;
    double z = std::sqrt(x * x + y * y) - 0.00002 * std::sin(y * pi);
    double theta = std::atan2(y, x) - 0.000003 * std::cos(x * pi);
    gcjLon = z * std::cos(theta);
    gcjLat = z * std::sin(theta);
}

转换为笛卡尔坐标系

由于GCJ02与BD09都是基于WGS84的偏移加密， 所以其大地坐标系的基准都是与WGS84一致。所以转换为笛卡尔坐标系需要通过WGS84进行转换。

WGS84 大地坐标系转换为 ECEF 空间直角坐标系的原理涉及几何学和椭球体的数学描述。简单来说，转换过程是通过将地球表面的点（经纬度和高度）转换为以地球中心为原点的三维直角坐标（X, Y, Z）。

转换的原理核心在于通过几何公式将椭球上的点投影到三维空间中的一个点。ECEF 坐标系统与地球同步旋转，因此适合用于全球定位系统（GPS）和其他涉及地球表面点的位置计算的系统。这种转换能够在全球范围内提供一个统一的三维坐标系，使得不同地点的测量数据能够准确地相互比较和计算。

以下是这个转换的详细原理：

1. 椭球模型描述

• WGS84 椭球模型：地球并非完美的球体，而是一个扁球体，赤道半径大于极半径。WGS84 使用一个椭球来近似地球的形状：
  • 长半轴 (a)：赤道的半径，大约为 6378137 米。
  • 短半轴 (b)：极点的半径，大约为 6356752 米。
  • 扁率 (f)：描述了椭球的扁平程度，计算公式为 $f=\frac{a-b}{a}$。

2. 大地坐标系

• 纬度 (lat)：从地心到椭球表面点的法线与赤道平面的夹角。
• 经度 (lon)：地球表面点所在的子午线与本初子午线（经度为 0° 的子午线）之间的夹角。
• 高度 (h)：地球表面点沿法线方向相对于椭球面的距离。

3. 空间直角坐标系 (ECEF)

• ECEF 系统：这个坐标系的原点位于地球中心，X 轴指向经度为 0° 的赤道上，Y 轴指向经度为 90° 的赤道上，Z 轴指向北极。

4. 转换公式推导

为了将大地坐标转换为 ECEF 坐标，需理解以下几个概念：

曲率半径 (N)

• 对于纬度为 $\text{lat}$ 的点，曲率半径 $N$ 表示的是椭球体表面上的一个点到椭球中心的距离，定义为：

$$N=\frac{a}{\sqrt{1-e^2{\sin }^2(\text{lat})}}$$

其中，$e^2=\frac{2f-f^2}{1}$ 是椭球的第一偏心率的平方，描述了椭球的形状。

ECEF 坐标公式

• 使用 $N$ 和地心坐标系的三维直角坐标，X、Y、Z 计算如下：
  $$X=(N+h)\cos (\text{lat})\cos (\text{lon})$$
  $$Y=(N+h)\cos (\text{lat})\sin (\text{lon})$$
  $$Z=[(1-e^2)N+h]\sin (\text{lat})$$

这些公式从几何上来讲是通过将地球表面的一个点投影到地心坐标系中来实现的。X 和 Y 是根据纬度和经度的余弦值来确定，Z 则根据纬度的正弦值计算。

5. 转换过程

1. 计算 N：根据纬度来计算曲率半径 $N$，这涉及到 WGS84 椭球的偏心率和纬度的三角函数。
2. 计算 ECEF 坐标：使用上面推导出的公式，通过纬度、经度和高度来计算该点在地心坐标系中的 X, Y, Z 坐标。

6. 代码示例

void wgs84_to_ecef(double lat, double lon, double alt, double& x, double& y, double& z) {
    const double a = 6378137.0;             // WGS-84 椭球体的长半轴（赤道半径）
    const double f = 1 / 298.257223563;     // 地球扁率
    const double e2 = 2 * f - f * f;        // 0.00669437999014;  WGS-84 椭球体的第一偏心率的平方

    double radLat = lat * M_PI / 180.0;     // 将纬度转换为弧度
    double radLon = lon * M_PI / 180.0;     // 将经度转换为弧度

    double N = a / std::sqrt(1 - e2 * std::sin(radLat) * std::sin(radLat));  // 椭球曲率半径

    x = (N + alt) * std::cos(radLat) * std::cos(radLon);
    y = (N + alt) * std::cos(radLat) * std::sin(radLon);
    z = (N * (1 - e2) + alt) * std::sin(radLat);
}

逆转:

void ecef_to_wgs84(double x, double y, double z, double& lat, double& lon, double& alt) {

    const double a = 6378137.0;                 // WGS-84 椭球体的长半轴（赤道半径）
    const double f = 1.0 / 298.257223563;       // 扁率
    const double e2 = 2 * f - f * f;            //  0.00669437999014;  WGS-84 椭球体的第一偏心率的平方
    const double b = a * (1 - f);               // 6356752.314245; WGS-84 椭球体的短半轴（极半径）


    double eps = 1e-12;   // 精度
    double p = std::sqrt(x * x + y * y);
    lon = std::atan2(y, x);    // 计算经度

    // 迭代计算纬度
    double theta = std::atan2(z * a, p * b);
    double sinTheta = std::sin(theta);
    double cosTheta = std::cos(theta);

    lat = std::atan2(z + e2 * b * sinTheta * sinTheta * sinTheta,
        p - e2 * a * cosTheta * cosTheta * cosTheta);

    double N = a / std::sqrt(1 - e2 * std::sin(lat) * std::sin(lat));
    alt = p / std::cos(lat) - N;

    lat = lat * 180.0 / M_PI;   // 转换为度数
    lon = lon * 180.0 / M_PI;   // 转换为度数
}