x,y直角坐标系转经纬度WGS-84坐标系

x,y直角坐标系转经纬度WGS-84坐标系

  坐标系的转换采用了白塞尔大地主题反算算法，需要指明x,y坐标系的中心点所对应的现实世界的经纬度，代码如下：

//获取方位角
double getAngle(double x, double y) {
	double angleTemp;
	if (x == 0 && y > 0)		//y轴正向
	{
		angleTemp = 0;
	}
	if (y >= 0 && x > 0)		//第一象限及x轴正向
	{
		angleTemp = 0.5*PI - atan(y/x);
	}
	else if (y < 0 && x > 0)	//第四象限
	{
		angleTemp = 0.5*PI + atan(abs(y)/x);
	}
	else if (x == 0 && y < 0)	//y轴负向
	{
		angleTemp = PI;
	}
	else if (y <= 0 && x < 0)	//第三象限及x轴负向
	{
		angleTemp = 1.5*PI - atan(y/x);
	}
	else						//第二象限
	{
		angleTemp = 1.5*PI + atan(y/abs(x));
	}
	return angleTemp;
}

//获取距离
double getDistance(double x, double y) {
	return sqrt(pow(x, 2) + pow(y, 2));
}

/*度分秒转换为弧度*/
double hudu(double x, double y, double z)
{
	double A0;
	A0 = (x + y / 60 + z / 3600)*PI / 180;
	return A0;
}

/*弧度转换为度*/
double du(double B0)
{
	double x0;
	x0 = (int)(B0 * 180 / PI);
	return x0;
}

/*弧度转换为分*/
double fen(double C0)
{
	double _y, y0;
	_y = (int)(C0 * 180 / PI);
	y0 = (fabs)((int)((C0 * 180 / PI - _y) * 60));
	return y0;
}
/*弧度转换为秒*/
double miao(double D0)
{
	double _z1, _z2, z0;
	_z1 = (int)(D0 * 180 / PI);
	_z2 = (int)((D0 * 180 / PI - _z1) * 60);
	z0 = (fabs)((double)(((D0 * 180 / PI - _z1) * 60 - _z2) * 60));
	return z0;
}

//日期:2021年1月2日
//码农:Wp
//描述:VTD(x,y)-> lat,lon double
void getConvertPosition(double x,double y, double &lat, double &lon){
    double ax, ay, az, bx, by, bz, cx, cy, cz, S, dz, ez, fz, B1, B2, L1, L2, A1, A2;
	int dx, dy, ex, ey, fx, fy;
	double W1, sinu1, cosu1, sinA0, coto1, sin2o1, cos2o1, sin2o, cos2o, A, B, C, r, t, o0, o, g, sinu2, q;
	double e2 = 0.006693421622966;
	A1 = getAngle(x, y);
	S = getDistance(x,y);
	/*调用函数*/
	B1 = hudu(30, 0, 0);
	L1 = hudu(120, 0, 0);
	//白塞尔大地主题解算
	W1 = sqrt(1 - e2*sin(B1)*sin(B1));
	sinu1 = sin(B1)*(sqrt(1 - e2)) / W1;
	cosu1 = cos(B1) / W1;
	sinA0 = cosu1*sin(A1);
	coto1 = cosu1*cos(A1) / sinu1;
	sin2o1 = 2 * coto1 / (coto1*coto1 + 1);
	cos2o1 = (coto1*coto1 - 1) / (coto1*coto1 + 1);
	A = 6356863.020 + (10718.949 - 13.474*(1 - sinA0*sinA0))*(1 - sinA0*sinA0);
	B = (5354.469 - 8.978*(1 - sinA0*sinA0))*(1 - sinA0*sinA0);
	C = (2.238*(1 - sinA0*sinA0))*(1 - sinA0*sinA0) + 0.006;
	r = 691.46768 - (0.58143 - 0.00144*(1 - sinA0*sinA0))*(1 - sinA0*sinA0);
	t = (0.2907 - 0.0010*(1 - sinA0*sinA0))*(1 - sinA0*sinA0);
	o0 = (S - (B + C*cos2o1)*sin2o1) / A;
	sin2o = sin2o1*cos(2 * o0) + cos2o1*sin(2 * o0);
	cos2o = cos2o1*cos(2 * o0) - sin2o1*sin(2 * o0);
	o = o0 + (B + 5 * C*cos2o)*sin2o / A;
	g = (r*o + t*(sin2o - sin2o1))*sinA0;
	//A1 = hudu(cx, cy, cz);
	/*求B2*/
	sinu2 = sinu1*cos(o) + cosu1*cos(A1)*sin(o);
	B2 = atan(sinu2 / (sqrt(1 - e2)*sqrt(1 - sinu2*sinu2)));

	/*求L2*/
	q = atan(sin(A1)*sin(o) / (cosu1*cos(o) - sinu1*sin(o)*cos(A1)));

	/*判断q*/
	if (sin(A1)>0 && tan(q)>0)
		q = fabs(q);
	else if (sin(A1)>0 && tan(q)<0)
		q = PI - fabs(q);
	else if (sin(A1)<0 && tan(q)<0)
		q = -fabs(q);
	else
		q = fabs(q) - PI;
	L2 = L1 + q - g / 3600 / 180 * PI;

	/*求A2*/
	A2 = atan(cosu1*sin(A1) / (cosu1*cos(o)*cos(A1) - sinu1*sin(o)));

	/*判断A2*/
	if (sin(A1)<0 && tan(A2)>0)
		A2 = fabs(A2);
	else if (sin(A1)<0 && tan(A2)<0)
		A2 = PI - fabs(A2);
	else if (sin(A1)>0 && tan(A2)>0)
		A2 = PI + fabs(A2);
	else
		A2 = 2 * PI - fabs(A2);

	/*调用函数*/
	dx = (int)(du(B2));
	dy = (int)(fen(B2));
	dz = miao(B2);
	ex = (int)(du(L2));
	ey = (int)(fen(L2));
	ez = miao(L2);
	fx = (int)(du(A2));
	fy = (int)(fen(A2));
	fz = miao(A2);

    lat = B2*180/PI;
    lon = L2*180/PI;
}

  如上，调用时直接调用getConvertPosition(x, y, &lat, &lon)即可，注意这里lat和lon参数传的是引用，这样在函数的修改的值会一直生效，引用是C++用法，如果是C的话，可以参照我之前的博客C 函数引用参数传递error
  除此之外需要注意X，Y坐标系的方位角与真实世界方位角的转换。（X,Y直角坐标系是以x轴正向为0°，逆时针旋转增加；现实世界方位角是以正北方向为0°，顺时针旋转增加）