如何将WGS84坐标转换成经纬度坐标

Posted

tags:

篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了如何将WGS84坐标转换成经纬度坐标相关的知识,希望对你有一定的参考价值。

WGS84坐标系属于大地坐标,就是我们常说的经纬度坐标.

搜集应用区域内GPS“B”级网三个以上网点WGS84坐标系B、L、H值及我国坐标系(BJ54或西安80)B、L、h、x值。

B、L、H分别为大地坐标系中的大地纬度、大地经度及大地高,h、x分别为大地坐标系中的高程及高程异常。各参数可以通过各省级测绘局或测绘院具有“A”级、“B”级网的单位获得。  

可以计算不同坐标系三维直角坐标值。计算公式如下: X=(N+H)cosBcosL Y=(N+H)cosBsinL Z=[N(1-e2)+H]sinB

扩展资料

参心大地坐标系:指经过定位与定向后,地球椭球的中心不与地球质心重合而是接近地球质心。区域性大地坐标系。是我国基本测图和常规大地测量的基础。如Beijing54、Xian80。

地心大地坐标系:指经过定位与定向后,地球椭球的中心与地球质心重合。如CGCS2000、WGS84。

地理坐标,就是用经线(子午线)、纬线、经度、纬度表示地面点位的球面坐标。

一般地理坐标可分为三种,天文经纬度,大地经纬度,地心经纬度。通常地图上使用的经纬度都为大地经纬度,

参考资料来源:百度百科-坐标转换

参考技术A

1、如果是三维空间坐标系的转换,涉及到 WGS84 空间直角坐标(X,Y,Z)和 大地坐标(B、L、H)的相互换算。公式如下

2、如果是二维平面坐标系的转换,涉及到WGS84 平面直角坐标(X,Y) 和 大地坐标 (B、L)的相互换算,即高斯投影正算和反算公式。

以上空间和平面坐标的换算较为复杂,需要一定的坐标系基础理论知识,最好使用相应的专业软件进行转换,且在转换前需要对WGS54 坐标系所采用的参考椭球参数、高斯投影的中央经线进行设置,才能进行坐标的转换。建议使用坐标转换软件COORD 完成坐标转换。

本回答被提问者采纳
参考技术B WGS84坐标本身就是经纬度坐标啊,你还要转换什么啊?

使用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 meters
	double 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 radian
	double 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-02
	double 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.1mm
		if (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 meters
		T 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 radian
		T 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-02
		T lon2 = T(kRAD2DEG) * lon1 / (N * cosB);
		T lat2 = T(kRAD2DEG) * (lat1 * W * W / (N * (1 - kKRASOVSKY_ECCSQ)));

		//residuals
		residuals[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 meters
		double 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);
		else
			if (lon0 < 0)
				g_lon0 = -0.05 / sqrt(-lon0);
			else
				g_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 radian
		double 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 function

	std::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 function

	std::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 ceres
	std::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++项目结合,可以移步这里

以上是关于如何将WGS84坐标转换成经纬度坐标的主要内容,如果未能解决你的问题,请参考以下文章

如何将WGS84坐标转换成经纬度坐标(度分秒的形式)?

wgs84和cgcs2000坐标转换

从UTM32N WGS84坐标系怎么转化成经纬度?拜谢

如何wgs84与当地坐标之间的转换参数

wgs84转换为百度坐标要如何操作呢?

WGS84坐标,如:26.1937457,104.101245,如何换算成X,Y坐标?是不是有啥软件