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

之前在網上找過一個火星坐標的轉換算法實現 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