地理坐标系：GCJ02和WGS84互转

 地理坐标系介绍

主流地理坐标系、投影坐标系和投影方法的区别和联系_ 一只博客-CSDN博客https://blog.csdn.net/qq_42276781/article/details/122597363

关于坐标系的转换，可以参考如下两个开源项目gcoord和coordtransform。

GitHub - hujiulong/gcoord: 地理坐标系转换工具https://github.com/hujiulong/gcoordGitHub - wandergis/coordtransform: 提供了百度坐标（BD09）、国测局坐标（火星坐标，GCJ02）、和WGS84坐标系之间的转换https://github.com/wandergis/coordtransform我将gcoord中的GCJ02和WGS84互转的代码改写成了C语言的形式，便于大家理解。

#include<stdio.h>
#include<math.h>
#include<stdlib.h>
#define PI 3.1415926535897932384626

double a = 6378245;
double ee = 0.006693421622965823;

// roughly check whether coordinates are in China.
int isInChinaBbox(double lon, double lat){
  	return lon >= 72.004 && lon <= 137.8347 && lat >= 0.8293 && lat <= 55.8271;
}

double transformLat(double x, double y){
	double ret = -100.0 + 2.0 * x + 3.0 * y + 0.2 * y * y + 0.1 * x * y + 0.2 * sqrt(fabs(x));
	ret += ((20.0 * sin(6.0 * x * PI) + 20.0 * sin(2.0 * x * PI)) * 2.0) / 3.0;
	ret += ((20.0 * sin(y * PI) + 40.0 * sin((y / 3.0) * PI)) * 2.0) / 3.0;
	ret += ((160.0 * sin((y / 12.0) * PI) + 320.0 * sin((y * PI) / 30.0)) * 2.0) / 3.0;
	return ret;
}

double transformLon(double x, double y){
	double ret = 300.0 + x + 2.0 * y + 0.1 * x * x + 0.1 * x * y + 0.1 * sqrt(fabs(x));
	ret += ((20.0 * sin(6.0 * x * PI) + 20.0 * sin(2.0 * x * PI)) * 2.0) / 3.0;
	ret += ((20.0 * sin(x * PI) + 40.0 * sin((x / 3.0) * PI)) * 2.0) / 3.0;
	ret += ((150.0 * sin((x / 12.0) * PI) + 300.0 * sin((x / 30.0) * PI)) * 2.0) / 3.0;
	return ret;
}

double* delta(double lon, double lat){
	double *d = (double *)malloc(2*sizeof(double));
	double dLon = transformLon(lon - 105.0, lat - 35.0);
    double dLat = transformLat(lon - 105.0, lat - 35.0);
	
	double radLat = (lat / 180.0) * PI;
	double magic = sin(radLat);
	
	magic = 1.0 - ee * magic * magic;
	
	double sqrtMagic = sqrt(magic);
	d[0] = (dLon * 180.0) / (a / sqrtMagic * cos(radLat) * PI);
	d[1] = (dLat * 180.0) / ((a * (1.0 - ee)) / (magic * sqrtMagic) * PI);
	return d;
}

double* WGS84ToGCJ02(double lon_WGS, double lat_WGS){
	double *GCJ = (double *)malloc(2*sizeof(double));
	if (!isInChinaBbox(lon_WGS, lat_WGS)){
		GCJ[0] = lon_WGS;
		GCJ[1] = lat_WGS;
		return GCJ;
	} 
	double *d = delta(lon_WGS, lat_WGS);
	GCJ[0] = lon_WGS + d[0];
	GCJ[1] = lat_WGS + d[1];
	return GCJ;
}

double* GCJ02ToWGS84(double lon_GCJ, double lat_GCJ) {
	double *WGS = (double *)malloc(2*sizeof(double));
	if (!isInChinaBbox(lon_GCJ, lat_GCJ)){
		WGS[0] = lon_GCJ;
		WGS[1] = lat_GCJ;
		return WGS;
	} 
	WGS[0] = lon_GCJ; WGS[1] = lat_GCJ;
	double *temp = WGS84ToGCJ02(WGS[0], WGS[1]);
	double dx = temp[0] - lon_GCJ;
	double dy = temp[1] - lat_GCJ;
	while (fabs(dx) > 1e-6 || fabs(dy) > 1e-6) {
		WGS[0] -= dx;
		WGS[1] -= dy;
		temp = WGS84ToGCJ02(WGS[0], WGS[1]);
		dx = temp[0] - lon_GCJ;
		dy = temp[1] - lat_GCJ;
	}
	return WGS;
}

int main(){
	double lon, lat;
	double x = 114.5048503, y = 37.0705911; 
	double *q = GCJ02ToWGS84(x, y);
	printf("%.10lf %.10lf\n", q[0], q[1]);
	getchar(); 
}