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°,顺时针旋转增加)