火星坐标(GCJ02)高精度反算

Posted 2021-05-27 [e-mail]

之前在网上找过一个火星坐标的转换算法实现 https://github.com/googollee/eviltransform ，但是其在部分区域的精度较低，达不到我们使用的要求。因为其没有进行有效的迭代计算，所以误差比较大。

参考 从地球到火星　～ 论 GCJ-02 及其衍生 这篇文章，我简单的借鉴文章内的代码写了一个 C++ 版本的转换计算实现代码，GCJ02 转 WGS84 的实现中使用了迭代计算逼近的方法来保证精度。经过我的测试，大部分情况下在 2 次迭代计算的情况下就可以达到 1 米精度的要求。

代码及测试

转换计算代码实现如下：

#include <cmath>

// China Transformed
namespace CT
{
/// 克拉索夫斯基椭球参数
static const double GCJ_A = 6378245;
static const double GCJ_EE = 0.00669342162296594323; // f = 1/298.3; e^2 = 2*f - f**2
static const double PI = 3.14159265358979323846;

struct Coord
{
    double lon; // 经度
    double lat; // 纬度
};

bool OutOfChina(Coord coords)
{
    return coords.lat < 0.8293 || coords.lat > 55.8271 ||
           coords.lon < 72.004 || coords.lon > 137.8347;
}

Coord WGS84ToGCJ02(Coord wgs, bool checkChina = true)
{
    if (checkChina && OutOfChina(wgs)) { return wgs; }

    // 将经度减去 105°，纬度减去 35°，求偏移距离
    double x = wgs.lon - 105;
    double y = wgs.lat - 35;

    // 将偏移距离转化为在 SK-42 椭球下的经纬度大小
    double dLat_m = -100.0 + 2.0 * x + 3.0 * y + 0.2 * y * y + 0.1 * x * y + 0.2 * std::sqrt(std::abs(x)) +
                    (2.0 * std::sin(x * 6.0 * PI) + 2.0 * std::sin(x * 2.0 * PI) +
                     2.0 * std::sin(y * PI) + 4.0 * std::sin(y / 3.0 * PI) +
                     16.0 * std::sin(y / 12.0 * PI) + 32.0 * std::sin(y / 30.0 * PI)) *
                      20.0 / 3.0;
    double dLon_m = 300.0 + x + 2.0 * y + 0.1 * x * x + 0.1 * x * y + 0.1 * std::sqrt(std::abs(x)) +
                    (2.0 * std::sin(x * 6.0 * PI) + 2.0 * std::sin(x * 2.0 * PI) +
                     2.0 * std::sin(x * PI) + 4.0 * std::sin(x / 3.0 * PI) +
                     15.0 * std::sin(x / 12.0 * PI) + 30.0 * std::sin(x / 30.0 * PI)) *
                      20.0 / 3.0;

    double radLat = wgs.lat / 180.0 * PI;
    double magic = 1.0 - GCJ_EE * std::pow(std::sin(radLat), 2.0); // just a common expr

    double lat_deg_arclen = (PI / 180.0) * (GCJ_A * (1.0 - GCJ_EE)) / std::pow(magic, 1.5);
    double lon_deg_arclen = (PI / 180.0) * (GCJ_A * std::cos(radLat) / std::sqrt(magic));

    // 往原坐标加偏移量
    return Coord{ wgs.lon + (dLon_m / lon_deg_arclen), wgs.lat + (dLat_m / lat_deg_arclen)};
}

// 计算两个坐标值的偏差值
Coord DiffCoord(Coord a, Coord b)
{
    return Coord{a.lon - b.lon,a.lat - b.lat};
}

// 精度控制（最大误差）
static const double g2w_precision = 1.0/111391.0;

Coord GCJ02ToWGS84(Coord gcj, bool checkChina = true)
{
    // 计算输入 gcj 坐标与将其计算为84坐标的偏差
    // 用当前的 gcj 坐标减去这个偏差，其近似于 gcj 对应的 84 坐标
    // 使用这个近似坐标去计算火星坐标，与输入的 gcj 进行比较，看是否符合精度
    // 如果不符合精度，则将近似坐标加上上面得到的偏差，再进行计算一次

    Coord wgs = DiffCoord(gcj, DiffCoord(WGS84ToGCJ02(gcj, checkChina), gcj));
    Coord d = DiffCoord(gcj,WGS84ToGCJ02(wgs));
    
    int MaxIterations = 10; // 最大迭代次数
    
    while((MaxIterations-- > 0) &&
          (std::abs(d.lon) > g2w_precision || std::abs(d.lat) > g2w_precision)) {
        wgs.lon += d.lon;
        wgs.lat += d.lat;
        d = DiffCoord(gcj,WGS84ToGCJ02(wgs));
    };
    return wgs;
}
}

测试代码如下：

#include <iostream>

int main()
{
    std::cout.precision(12);

    CT::Coord a{ 120.34, 36.10 };
    std::cout << "wgs84 = " << a.lon << "    " << a.lat << std::endl;
    auto b = CT::WGS84ToGCJ02(a);
    std::cout << "gcj02 = " << b.lon << "    " << b.lat << std::endl;
    auto c = CT::DiffCoord(b,a);
    std::cout << "diff  = " << c.lon << "    " << c.lat << std::endl;
    a = CT::GCJ02ToWGS84(b);
    std::cout << "wgs84 = " << a.lon << "    " << a.lat << std::endl;
    return 0;
}

编译输出情况：

$ clang++ gcj02.cpp

o@x-pc MINGW64 /g/ctemp
$ ./a.exe
wgs84 = 120.34    36.1
gcj02 = 120.345088846    36.1002223937
diff  = 0.0050888458279    0.000222393684851
wgs84 = 120.340000014    36.1000000171