使用ceres库将经纬度坐标GCJ02到WGS84精确转换

在之前的一篇博客中，介绍了经纬度坐标系的关系，WGS84是世界大地坐标系，属于原始坐标系，在商用中，需要通过火星加密算法将经纬度做转换，转换之后的坐标，称为国测局坐标，也叫火星坐标系，简称GCJ02。这个过程称为坐标偏移算法。反过来，我们将GCJ02坐标系转为WGS84就是纠偏算法，严格来说，他们不是一对正向与反向的可逆操作，但是有理论公式可以推导。

ceres库是一个c++开源库，用来解决非线性优化问题，可以解决边界约束鲁棒非线性最小二乘法优化问题（Ceres Solver 1 is an open source C++ library for modeling and solving large, complicated optimization problems. It can be used to solve Non-linear Least Squares problems with bounds constraints and general unconstrained optimization problems），通过这个库，我们几乎可以将GCJ02精确转换为WGS84坐标。下面来看示例。

代码是从github上拿过来的。

c++头文件：

/***************************************************************************** This class WGStoGCJ implements geodetic coordinate transform between geodetic* coordinate in WGS84 coordinate system and geodetic coordinate in GCJ-02* coordinate system.** details: https://blog.csdn.net/gudufuyun/article/details/106738942** Copyright (c) 2020.  All rights reserved.**  kikkimo*  School of Remote Sensing  and Information Engineering,*  WuHan University,*  Wuhan, Hubei, P.R. China. 430079*  fywhu@outlook.com*****************************************************************************                                                                         **   This program is free software; you can redistribute it and/or modify  **   it under the terms of the GNU General Public License as published by  **   the Free Software Foundation; either version 2 of the License, or     **   (at your option) any later version.                                   **                                                                         ****************************************************************************/
#pragma once
#include <utility>#ifdef _WGSTOGCJ_API
#if (defined _WIN32 || defined WINCE || defined __CYGWIN__)
#define WGSTOGCJ_API_EXPORTS __declspec(dllexport)
#elif defined __GNUC__ && __GNUC__ >= 4
#define WGSTOGCJ_API_EXPORTS __attribute__((visibility("default")))
#endif
#endif  #ifndef WGSTOGCJ_API_EXPORTS
#define WGSTOGCJ_API_EXPORTS
#endif/***  \brief Covert geodetic coordinate in WGS84 coordinate system to geodetic coordinate*         in GCJ-02 coordinate system**  \param [in] wgs84lon: longitude in WGS84 coordinate system [unit:degree]*  \param [in] wgs84lat: latitude in WGS84 coordinate system [unit:degree]*  \return Returns geodetic coordinate in GCJ-02 coordinate system*  \time 15:47:38 2020/06/12*/
WGSTOGCJ_API_EXPORTS std::pair<double, double>
Wgs2Gcj(const double& wgs84lon, const double& wgs84lat);/***  \brief Covert geodetic coordinate in GCJ-02 coordinate system to geodetic coordinate*         in WGS84 coordinate system**  \param [in] gcj02lon: longitude in GCJ-02 coordinate system [unit:degree]*  \param [in] gcj02lat: latitude in GCJ-02 coordinate system [unit:degree]*  \return Returns geodetic coordinate in WGS84 coordinate system*  \remark simple linear iteration*  \detail*  \time 15:51:13 2020/06/12*/
std::pair<double, double> __declspec(dllexport)
Gcj2Wgs_SimpleIteration(const double& gcj02lon, const double& gcj02lat);/***  \brief Covert geodetic coordinate in GCJ-02 coordinate system to geodetic coordinate*         in WGS84 coordinate system**  \param [in] gcj02lon: longitude in GCJ-02 coordinate system [unit:degree]*  \param [in] gcj02lat: latitude in GCJ-02 coordinate system [unit:degree]*  \return Returns geodetic coordinate in WGS84 coordinate system*  \remark the encryption formula is known and use an auto-differentiable cost function*  \time 15:51:13 2020/06/12*/
WGSTOGCJ_API_EXPORTS std::pair<double, double>
Gcj2Wgs_AutoDiff(const double& gcj02lon, const double& gcj02lat);/***  \brief Covert geodetic coordinate in GCJ-02 coordinate system to geodetic coordinate*         in WGS84 coordinate system**  \param [in] gcj02lon: longitude in GCJ-02 coordinate system [unit:degree]*  \param [in] gcj02lat: latitude in GCJ-02 coordinate system [unit:degree]*  \return Returns geodetic coordinate in WGS84 coordinate system*  \remark the encryption formula is unknown but we can covert point in WGS84 to point*          in GCJ-02 with an API,then use the numerical derivation method to solve the*          problem* \detail Assuming the encryption formula is**            gcj02lon = Wgs2Gcj(wgs84lon, wgs84lat)*            gcj02lat = Wgs2Gcj(wgs84lon, wgs84lat)**    In the rectification process, (wgs84lon, wgs84lat) are unknown items. Obviously,*   this is a system of nonlinear equations.**   The linear formed error functions of forward intersection come from*   consideration of a Taylor series expansion.*           V = AX - b*    here:*    V: The residuals of the observed values*    A: The jacobian matrix*    X: The modification of the unknown items*    b: The constant terms of the error functions**    Then the error functions written in vector form are:*    | V_lon | = | dlongcj_dlonwgs  dlongcj_dlatwgs |  | d_lonwgs | - | l_lon |*    | V_lat | = |         0        dlatgcj_dlatwgs |  | d_latwgs | - | l_lat |*    here:*    l_lon = longcj - longcj'                 // the modification of longcj*    l_lat = latgcj - latgcj'                 // the modification of latgcj*    longcj : the observed longitude in GCJ-02*    latgcj : the observed latitude in GCJ-02*    longcj' = Wgs2Gcj(wgs84lon',wgs84lat')    // estimated longitude in GCJ-02*    latgcj' = Wgs2Gcj(wgs84lon',wgs84lat')    // estimated latitude in GCJ-02*    wgs84lon' : estimated longitude in WGS84*    wgs84lat' : estimated latitude in WGS84*    d_lonwgs : unknown items*    d_latwgs : unknown items*    wgs84lon = wgs84lon' + d_lonwgs                           // update*    wgs84lat = wgs84lat' + d_latwgs**     let V = [V_lon V_lat]T = 0, then*     d_latwgs = (l_lon * dlatgcj_dlonwgs - l_lat * dlongcj_dlonwgs) /*            (dlongcj_dlatwgs * dlatgcj_dlonwgs - dlatgcj_dlatwgs * dlongcj_dlonwgs)*      d_lonwgs = (l_lon - dlongcj_dlatwgs * d_latwgs) / dlongcj_dlonwgs**    This iterative procedure is repeated until X= [d_lonwgs d_latwgs]T are*    sufficiently small.*  \time 17:42:01 2020/06/12*/
WGSTOGCJ_API_EXPORTS std::pair<double, double>
Gcj2Wgs_NumbericDiff(const double& gcj02lon, const double& gcj02lat);/***  \brief Covert geodetic coordinate in GCJ-02 coordinate system to geodetic coordinate*         in WGS84 coordinate system**  \param [in] gcj02lon: longitude in GCJ-02 coordinate system [unit:degree]*  \param [in] gcj02lat: latitude in GCJ-02 coordinate system [unit:degree]*  \return Returns geodetic coordinate in WGS84 coordinate system*  \remark The encryption formula is known,and use the analytical derivation method to*         solve the problem with high precision.* \detail Assuming the encryption formula is**            gcj02lon = Wgs2Gcj(wgs84lon, wgs84lat)*            gcj02lat = Wgs2Gcj(wgs84lon, wgs84lat)**    In the rectification process, (wgs84lon, wgs84lat) are unknown items. Obviously,*   this is a system of nonlinear equations.**   The linear formed error functions of forward intersection come from*   consideration of a Taylor series expansion.*           V = AX - b*    here:*    V: The residuals of the observed values*    A: The jacobian matrix*    X: The modification of the unknown items*    b: The constant terms of the error functions**    Then the error functions written in vector form are:*    | V_lon | = | dlongcj_dlonwgs  dlongcj_dlatwgs |  | d_lonwgs | - | l_lon |*    | V_lat | = | dlatgcj_dlonwgs  dlatgcj_dlatwgs |  | d_latwgs | - | l_lat |*    here:*    l_lon = longcj - longcj'                 // the modification of longcj*    l_lat = latgcj - latgcj'                 // the modification of latgcj*    longcj : the observed longitude in GCJ-02*    latgcj : the observed latitude in GCJ-02*    longcj' = Wgs2Gcj(wgs84lon',wgs84lat')    // estimated longitude in GCJ-02*    latgcj' = Wgs2Gcj(wgs84lon',wgs84lat')    // estimated latitude in GCJ-02*    wgs84lon' : estimated longitude in WGS84*    wgs84lat' : estimated latitude in WGS84*    d_lonwgs : unknown items*    d_latwgs : unknown items*    wgs84lon = wgs84lon' + d_lonwgs                           // update*    wgs84lat = wgs84lat' + d_latwgs**     let V = [V_lon V_lat]T = 0, then*     d_latwgs = (l_lon * dlatgcj_dlonwgs - l_lat * dlongcj_dlonwgs) /*            (dlongcj_dlatwgs * dlatgcj_dlonwgs - dlatgcj_dlatwgs * dlongcj_dlonwgs)*      d_lonwgs = (l_lon - dlongcj_dlatwgs * d_latwgs) / dlongcj_dlonwgs**    This iterative procedure is repeated until X= [d_lonwgs d_latwgs]T are*    sufficiently small.*  \time 01:54:46 2020/06/13*/
WGSTOGCJ_API_EXPORTS std::pair<double, double>
Gcj2Wgs_AnalyticDiff(const double& gcj02lon, const double& gcj02lat);

c++源文件：

#include <iostream>
#include <cmath>
#include "WGS84GCJ02.h"
#include <ceres/ceres.h>constexpr double kKRASOVSKY_A = 6378245.0;               // equatorial radius [unit: meter]
constexpr double kKRASOVSKY_B = 6356863.0187730473;    // polar radius
constexpr double kKRASOVSKY_ECCSQ = 6.6934216229659332e-3; // first eccentricity squared
constexpr double kKRASOVSKY_ECC2SQ = 6.7385254146834811e-3; // second eccentricity squared
constexpr double PI = 3.14159265358979323846;   //πconstexpr double kDEG2RAD = PI / 180.0;
constexpr double kRAD2DEG = 180.0 / PI;constexpr inline double Deg2Rad(const double deg) {return deg * kDEG2RAD;
}constexpr inline double Rad2Deg(const double rad) {return rad * kRAD2DEG;
}std::pair<double, double> GetGeodeticOffset(const double& wgs84lon, const double& wgs84lat)
{//get geodetic offset relative to 'center china'double lon0 = wgs84lon - 105.0;double lat0 = wgs84lat - 35.0;//generate an pair offset roughly in metersdouble lon1 = 300.0 + lon0 + 2.0 * lat0 + 0.1 * lon0 * lon0 + 0.1 * lon0 * lat0 + 0.1 * sqrt(fabs(lon0));lon1 = lon1 + (20.0 * sin(6.0 * lon0 * PI) + 20.0 * sin(2.0 * lon0 * PI)) * 2.0 / 3.0;lon1 = lon1 + (20.0 * sin(lon0 * PI) + 40.0 * sin(lon0 / 3.0 * PI)) * 2.0 / 3.0;lon1 = lon1 + (150.0 * sin(lon0 / 12.0 * PI) + 300.0 * sin(lon0 * PI / 30.0)) * 2.0 / 3.0;double lat1 = -100.0 + 2.0 * lon0 + 3.0 * lat0 + 0.2 * lat0 * lat0 + 0.1 * lon0 * lat0 + 0.2 * sqrt(fabs(lon0));lat1 = lat1 + (20.0 * sin(6.0 * lon0 * PI) + 20.0 * sin(2.0 * lon0 * PI)) * 2.0 / 3.0;lat1 = lat1 + (20.0 * sin(lat0 * PI) + 40.0 * sin(lat0 / 3.0 * PI)) * 2.0 / 3.0;lat1 = lat1 + (160.0 * sin(lat0 / 12.0 * PI) + 320.0 * sin(lat0 * PI / 30.0)) * 2.0 / 3.0;//latitude in radiandouble B = Deg2Rad(wgs84lat);double sinB = sin(B), cosB = cos(B);double W = sqrt(1 - kKRASOVSKY_ECCSQ * sinB * sinB);double N = kKRASOVSKY_A / W;//geodetic offset used by GCJ-02double lon2 = Rad2Deg(lon1 / (N * cosB));double lat2 = Rad2Deg(lat1 * W * W / (N * (1 - kKRASOVSKY_ECCSQ)));return { lon2, lat2 };
}std::pair<double, double> Wgs2Gcj(const double& wgs84lon, const double& wgs84lat)
{std::pair<double, double> dlonlat = GetGeodeticOffset(wgs84lon, wgs84lat);double gcj02lon = wgs84lon + dlonlat.first;double gcj02lat = wgs84lat + dlonlat.second;return { gcj02lon, gcj02lat };
}std::pair<double, double> Gcj2Wgs_SimpleIteration(const double& gcj02lon,const double& gcj02lat)
{std::pair<double, double> lonlat = Wgs2Gcj(gcj02lon, gcj02lat);double lon0 = lonlat.first;double lat0 = lonlat.second;int iterCounts = 0;while (++iterCounts < 1000){std::pair<double, double> lonlat1 = Wgs2Gcj(lon0, lat0);double lon1 = lonlat1.first;double lat1 = lonlat1.second;double dlon = lon1 - gcj02lon;double dlat = lat1 - gcj02lat;lon1 = lon0 - dlon;lat1 = lat0 - dlat;//1.0e-9 degree corresponding to 0.1mmif (fabs(dlon) < 1.0e-9 && fabs(dlat) < 1.0e-9)break;lon0 = lon1;lat0 = lat1;}return { lon0 , lat0 };
}class AutoDiffCostFunc
{
public:AutoDiffCostFunc(const double lon, const double lat) :mlonGcj(lon), mlatGcj(lat) {}template <typename T>bool operator() (const T* const lonWgs, const T* const latWgs, T* residuals) const{//get geodetic offset relative to 'center china'T lon0 = lonWgs[0] - T(105.0);T lat0 = latWgs[0] - T(35.0);//generate an pair offset roughly in metersT lon1 = T(300.0) + lon0 + T(2.0) * lat0 + T(0.1) * lon0 * lon0 + T(0.1) * lon0 * lat0+ T(0.1) * ceres::sqrt(ceres::abs(lon0));lon1 = lon1 + (T(20.0) * ceres::sin(T(6.0) * lon0 * T(PI)) + T(20.0) * ceres::sin(T(2.0) * lon0 * T(PI))) * T(2.0) / T(3.0);lon1 = lon1 + (T(20.0) * ceres::sin(lon0 * T(PI)) + T(40.0) * ceres::sin(lon0 / T(3.0) * T(PI))) * T(2.0) / T(3.0);lon1 = lon1 + (T(150.0) * ceres::sin(lon0 / T(12.0) * T(PI)) + T(300.0) * ceres::sin(lon0 * T(PI) / T(30.0))) * T(2.0) / T(3.0);T lat1 = T(-100.0) + T(2.0) * lon0 + T(3.0) * lat0 + T(0.2) * lat0 * lat0 + T(0.1) * lon0 * lat0+ T(0.2) * ceres::sqrt(ceres::abs(lon0));lat1 = lat1 + (T(20.0) * ceres::sin(T(6.0) * lon0 * T(PI)) + T(20.0) * ceres::sin(T(2.0) * lon0 * T(PI))) * T(2.0) / T(3.0);lat1 = lat1 + (T(20.0) * ceres::sin(lat0 * T(PI)) + T(40.0) * ceres::sin(lat0 / T(3.0) * T(PI))) * T(2.0) / T(3.0);lat1 = lat1 + (T(160.0) * ceres::sin(lat0 / T(12.0) * T(PI)) + T(320.0) * ceres::sin(lat0 * T(PI) / T(30.0))) * T(2.0) / T(3.0);//latitude in radianT B = latWgs[0] * T(kDEG2RAD);T sinB = ceres::sin(B), cosB = ceres::cos(B);T W = ceres::sqrt(T(1) - T(kKRASOVSKY_ECCSQ) * sinB * sinB);T N = T(kKRASOVSKY_A) / W;//geodetic offset used by GCJ-02T lon2 = T(kRAD2DEG) * lon1 / (N * cosB);T lat2 = T(kRAD2DEG) * (lat1 * W * W / (N * (1 - kKRASOVSKY_ECCSQ)));//residualsresiduals[0] = lonWgs[0] + lon2 - mlonGcj;residuals[1] = latWgs[0] + lat2 - mlatGcj;return true;}private:double mlonGcj;double mlatGcj;
};std::pair<double, double> Gcj2Wgs_AutoDiff(const double& gcj02lon,const double& gcj02lat)
{ceres::Problem * poProblem = new ceres::Problem;AutoDiffCostFunc* pCostFunc = new AutoDiffCostFunc(gcj02lon, gcj02lat);double wgslon = gcj02lon, wgslat = gcj02lat;poProblem->AddResidualBlock(new ceres::AutoDiffCostFunction<AutoDiffCostFunc, 2, 1, 1>(pCostFunc),nullptr,&wgslon,&wgslat);ceres::Solver::Options options;options.max_num_iterations = 30;options.linear_solver_type = ceres::DENSE_QR;options.minimizer_progress_to_stdout = false;options.gradient_tolerance = 1e-16;options.function_tolerance = 1e-12;options.parameter_tolerance = 1e-14;ceres::Solver::Summary summary;ceres::Solve(options, poProblem, &summary);delete poProblem;        //auto free memory of cost function "pCostFunc"return { wgslon, wgslat };
}std::pair<double, double> GetPartialDerivative_Lon(const double& wgs84lon, const double& wgs84lat, const double& dlon)
{double lonBk = wgs84lon + dlon;double lonFw = wgs84lon - dlon;std::pair<double, double > gcjlatlonBk = Wgs2Gcj(lonBk, wgs84lat);double gcjlonBk = gcjlatlonBk.first;double gcjlatBk = gcjlatlonBk.second;std::pair<double, double > gcjlatlonFw = Wgs2Gcj(lonFw, wgs84lat);double gcjlonFw = gcjlatlonFw.first;double gcjlatFw = gcjlatlonFw.second;double dlongcj_dlonwgs = (gcjlonBk - gcjlonFw) / (dlon * 2.0);double dlatgcj_dlonwgs = (gcjlatBk - gcjlatFw) / (dlon * 2.0);return { dlongcj_dlonwgs , dlatgcj_dlonwgs };
}std::pair<double, double> GetPartialDerivative_Lat(const double& wgs84lon, const double& wgs84lat, const double& dlat)
{double latBk = wgs84lat + dlat;double latFw = wgs84lat - dlat;std::pair<double, double> gcjlonlatBk = Wgs2Gcj(wgs84lon, latBk);double gcjlonBk = gcjlonlatBk.first;double gcjlatBk = gcjlonlatBk.second;std::pair<double, double> gcjlonlatFw = Wgs2Gcj(wgs84lon, latFw);double gcjlonFw = gcjlonlatFw.first;double gcjlatFw = gcjlonlatFw.second;double dlongcj_dlatwgs = (gcjlonBk - gcjlonFw) / (dlat * 2.0);double dlatgcj_dlatwgs = (gcjlatBk - gcjlatFw) / (dlat * 2.0);return { dlongcj_dlatwgs , dlatgcj_dlatwgs };
}std::pair<double, double> Gcj2Wgs_NumbericDiff(const double& gcj02lon,const double& gcj02lat)
{double wgs84lon = gcj02lon, wgs84lat = gcj02lat;int nIterCount = 0;double tol = 1e-9;while (++nIterCount < 1000){std::pair<double, double > data1 = GetPartialDerivative_Lon(wgs84lon, wgs84lat, tol);double dlongcj_dlonwgs = data1.first;double dlatgcj_dlonwgs = data1.second;std::pair<double, double > data2 = GetPartialDerivative_Lat(wgs84lon, wgs84lat, tol);double dlongcj_dlatwgs = data2.first;double dlatgcj_dlatwgs = data2.second;std::pair<double, double > data3 = Wgs2Gcj(wgs84lon, wgs84lat);double gcj02lonEst = data3.first;double gcj02latEst = data3.second;double l_lon = gcj02lon - gcj02lonEst;double l_lat = gcj02lat - gcj02latEst;double d_latwgs = (l_lon * dlatgcj_dlonwgs - l_lat * dlongcj_dlonwgs) /(dlongcj_dlatwgs * dlatgcj_dlonwgs - dlatgcj_dlatwgs * dlongcj_dlonwgs);double d_lonwgs = (l_lon - dlongcj_dlatwgs * d_latwgs) / dlongcj_dlonwgs;if (fabs(d_latwgs) < tol && fabs(d_lonwgs) < tol)break;wgs84lon = wgs84lon + d_lonwgs;wgs84lat = wgs84lat + d_latwgs;}return { wgs84lon, wgs84lat };
}std::pair<double, double> Gcj2Wgs_AnalyticDiff(const double& gcj02lon,const double& gcj02lat)
{double wgs84lon = gcj02lon, wgs84lat = gcj02lat;int nIterCount = 0;while (++nIterCount < 1000){//get geodetic offset relative to 'center china'double lon0 = wgs84lon - 105.0;double lat0 = wgs84lat - 35.0;//generate an pair offset roughly in metersdouble lon1 = 300.0 + lon0 + 2.0 * lat0 + 0.1 * lon0 * lon0 + 0.1 * lon0 * lat0 + 0.1 * sqrt(fabs(lon0));lon1 = lon1 + (20.0 * sin(6.0 * lon0 * PI) + 20.0 * sin(2.0 * lon0 * PI)) * 2.0 / 3.0;lon1 = lon1 + (20.0 * sin(lon0 * PI) + 40.0 * sin(lon0 / 3.0 * PI)) * 2.0 / 3.0;lon1 = lon1 + (150.0 * sin(lon0 / 12.0 * PI) + 300.0 * sin(lon0 * PI / 30.0)) * 2.0 / 3.0;double lat1 = -100.0 + 2.0 * lon0 + 3.0 * lat0 + 0.2 * lat0 * lat0 + 0.1 * lon0 * lat0 + 0.2 * sqrt(fabs(lon0));lat1 = lat1 + (20.0 * sin(6.0 * lon0 * PI) + 20.0 * sin(2.0 * lon0 * PI)) * 2.0 / 3.0;lat1 = lat1 + (20.0 * sin(lat0 * PI) + 40.0 * sin(lat0 / 3.0 * PI)) * 2.0 / 3.0;lat1 = lat1 + (160.0 * sin(lat0 / 12.0 * PI) + 320.0 * sin(lat0 * PI / 30.0)) * 2.0 / 3.0;double g_lon0 = 0;if (lon0 > 0)g_lon0 = 0.05 / sqrt(lon0);elseif (lon0 < 0)g_lon0 = -0.05 / sqrt(-lon0);elseg_lon0 = 0;double PIlon0 = PI * lon0, PIlat0 = PI * lat0;double dlon1_dlonwgs = 1 + 0.2 * lon0 + 0.1 * lat0 + g_lon0+ ((120 * PI * cos(6 * PIlon0) + 40 * PI * cos(2 * PIlon0))+ (20 * PI * cos(PIlon0) + 40 * PI / 3.0 * cos(PIlon0 / 3.0))+ (12.5 * PI * cos(PIlon0 / 12.0) + 10 * PI * cos(PIlon0 / 30.0))) * 2.0 / 3.0;double dlon1_dlatwgs = 2 + 0.1 * lon0;double dlat1_dlonwgs = 2 + 0.1 * lat0 + 2 * g_lon0+ (120 * PI * cos(6 * PIlon0) + 40 * PI * cos(2 * PIlon0)) * 2.0 / 3.0;double dlat1_dlatwgs = 3 + 0.4 * lat0 + 0.1 * lon0+ ((20 * PI * cos(PIlat0) + 40.0 * PI / 3.0 * cos(PIlat0 / 3.0))+ (40 * PI / 3.0 * cos(PIlat0 / 12.0) + 32.0 * PI / 3.0 * cos(PIlat0 / 30.0))) * 2.0 / 3.0;//latitude in radiandouble B = Deg2Rad(wgs84lat);double sinB = sin(B), cosB = cos(B);double WSQ = 1 - kKRASOVSKY_ECCSQ * sinB * sinB;double W = sqrt(WSQ);double N = kKRASOVSKY_A / W;double dW_dlatwgs = -PI * kKRASOVSKY_ECCSQ * sinB * cosB / (180.0 * W);double dN_dlatwgs = -kKRASOVSKY_A * dW_dlatwgs / WSQ;double PIxNxCosB = PI * N * cosB;double dlongcj_dlonwgs = 1.0 + 180.0 * dlon1_dlonwgs / PIxNxCosB;double dlongcj_dlatwgs = 180 * dlon1_dlatwgs / PIxNxCosB -180 * lon1 * PI * (dN_dlatwgs * cosB - PI * N * sinB / 180.0) / (PIxNxCosB * PIxNxCosB);double PIxNxSubECCSQ = PI * N * (1 - kKRASOVSKY_ECCSQ);double dlatgcj_dlonwgs = 180 * WSQ * dlat1_dlonwgs / PIxNxSubECCSQ;double dlatgcj_dlatwgs = 1.0 + 180 * (N * (dlat1_dlatwgs * WSQ + 2.0 * lat1 * W * dW_dlatwgs) - lat1 * WSQ * dN_dlatwgs) /(N * PIxNxSubECCSQ);std::pair<double, double> gcj02lonlatEst = Wgs2Gcj(wgs84lon, wgs84lat);double gcj02lonEst = gcj02lonlatEst.first;double gcj02latEst = gcj02lonlatEst.second;double l_lon = gcj02lon - gcj02lonEst;double l_lat = gcj02lat - gcj02latEst;double d_latwgs = (l_lon * dlatgcj_dlonwgs - l_lat * dlongcj_dlonwgs) /(dlongcj_dlatwgs * dlatgcj_dlonwgs - dlatgcj_dlatwgs * dlongcj_dlonwgs);double d_lonwgs = (l_lon - dlongcj_dlatwgs * d_latwgs) / dlongcj_dlonwgs;if (fabs(d_latwgs) < 1.0e-9 && fabs(d_lonwgs) < 1.0e-9)break;wgs84lon = wgs84lon + d_lonwgs;wgs84lat = wgs84lat + d_latwgs;}return { wgs84lon, wgs84lat };
}void TestCase1()
{double lonWgs = 116.39123343289631;double latWgs = 39.9072885060602;std::cout.precision(16);std::cout << "WGS84 Point: (" << latWgs << ", " <<  lonWgs << ")\n";std::pair<double, double> lonlatGcj = Wgs2Gcj(lonWgs, latWgs);double lonGcj = lonlatGcj.first; double latGcj = lonlatGcj.second;std::cout << "GCJ-02 Point [Wgs2Gcj]: (" << latGcj << ", " <<  lonGcj << ")\n";//simple linear iteration std::pair<double, double> lonlatWgsNS = Gcj2Wgs_SimpleIteration(lonGcj, latGcj);double lonWgsNS = lonlatWgsNS.first;double latWgsNS = lonlatWgsNS.second;std::cout << "WGS84 Point [simple linear iteration]: (" << latWgsNS << ", " <<  lonWgsNS << ")\n";//numerically differentiated cost functionstd::pair<double, double> lonlatWgsND = Gcj2Wgs_NumbericDiff(lonGcj, latGcj);double lonWgsND = lonlatWgsND.first;double latWgsND = lonlatWgsND.second;std::cout << "WGS84 Point [numerically derivation]: (" << latWgsND << ", " <<  lonWgsND << ")\n";//analytical differentiated cost functionstd::pair<double, double> lonlatWgsNAn = Gcj2Wgs_AnalyticDiff(lonGcj, latGcj);double lonWgsNAn = lonlatWgsNAn.first;double latWgsNAn = lonlatWgsNAn.second;std::cout << "WGS84 Point [analytical derivation]: (" << latWgsNAn << ", " <<  lonWgsNAn << ")\n";//auto differentiable, use ceresstd::pair<double, double> lonlatWgsNA = Gcj2Wgs_AutoDiff(lonGcj, latGcj);double lonWgsNA = lonlatWgsNA.first;double latWgsNA = lonlatWgsNA.second;std::cout << "WGS84 Point [auto differentiable]: (" << latWgsNA << ", " <<  lonWgsNA << ")\n";}int main() {TestCase1();return 0;
}

运行结果：

从打印结果来看，普通的转换，虽然精度上可以接近原始gps坐标，但是都无法完全还原，只有ceres代码执行的结果才能达到完全还原的效果。

有关ceres库的编译以及与c++项目结合，可以移步这里。

使用ceres库将经纬度坐标GCJ02到WGS84精确转换相关推荐

【地图转换工具类】：GCJ02与WGS84标准转换
[地图转换工具类]:GCJ02与WGS84标准转换 public class GCJ02_WGS84 {public static double pi = 3.14159265358979323846 ...
2020FME博客大赛——解放大脑经纬度坐标自动重投影至常用投影坐标系
作者:崔欣单位:中国石油天然气管道工程有限公司摘要:非测绘专业以及学艺不精的测绘人员对经纬度.投影带.带号.假东.假北.比例因子.高斯克吕格3度分带投影.高斯克吕格6度分带投影.墨卡托投影.通用横 ...
经纬度坐标转化为XYZ坐标的理解
前言经纬度坐标和XYZ笛卡尔坐标的转换常常应用在有关全景图的研究方面,在看了几篇提案和论文后才理解了坐标的转换方法正文通过一张图来说明按照我们直观的理解,会认为φ角是蓝色实线和y轴的夹角,θ角 ...
php 经纬度坐标转换 WGS84、火星坐标 (GCJ-02)、百度坐标 (BD-09)
项目有gps上报的功能, 由于前端插件问题导致大量gps定位数据转换百度坐标(BD-09)时产生极大偏移, 故需要后端做经纬度坐标转换, 看到一篇java的相关技术帖, 拿来做了修改 Ps: 坐标转换 ...
WGS84 / BD09 / GCJ02 / MapBar 经纬度坐标互转
WGS84 / BD09 / GCJ02 / MapBar 经纬度坐标互转 Geolocataion converting between WGS84, BD09 and GCJ02. WGS84 / ...
java 经纬度坐标转换 WGS84、火星坐标 (GCJ-02)、百度坐标 (BD-09)
会有偏移,但是还能接受 WGS84 国际标准,从 GPS 设备中取出的数据的坐标系国际地图提供商使用的坐标系火星坐标 (GCJ-02) 中国标准,从国行移动设备中定位获取的坐标数据使用这个坐标系 ...
根据经纬度坐标获得省市区县行政区划城市名称，自建数据库 java python php c# .net 均适用
文章目录步骤一.下载省市区边界数据步骤二.解析CSV文件导入数据库步骤三.在程序中根据坐标解析获得城市在LBS应用中,根据坐标来解析获得对应是哪个城市是一个很常见的功能,比如App里面通过手机 ...
【GUI界面】基于Python的WSG84三点定位系统（经纬度坐标与平面坐标转换法求解）
[GUI界面]基于Python的WSG84三点定位系统(经纬度坐标与平面坐标转换法求解) 方法汇总: blog.csdn.net/weixin_53403301/article/details/128 ...
【Python】利用Python实现精准三点定位（经纬度坐标与平面坐标转换法求解）
[Python]利用Python实现精准三点定位(经纬度坐标与平面坐标转换法求解) 众所周知,如果已知三个点的坐标,到一个未知点的距离,则可以利用以距离为半径画圆的方式来求得未知点坐标. 如果只有两个 ...
【精准三点定位求解汇总】利用Python或JavaScript高德地图开放平台实现精准三点定位（经纬度坐标与平面坐标转换法求解、几何绘图法求解）
[精准三点定位求解汇总]利用Python或JavaScript高德地图开放平台实现精准三点定位(经纬度坐标与平面坐标转换法求解.几何绘图法求解) 众所周知,如果已知三个点的坐标,到一个未知点的距离,则 ...

使用ceres库将经纬度坐标GCJ02到WGS84精确转换

使用ceres库将经纬度坐标GCJ02到WGS84精确转换相关推荐

最新文章

热门文章