📄 coordtrans.cpp
字号:
N = BJ54a / sqrt (1 - BJ54e2 * sin(geodetic->dLatitude) * sin(geodetic->dLatitude));
geodetic->dHeight = p / cos (geodetic->dLatitude)-N;
do {
oldgeodetic.dLatitude = geodetic->dLatitude;
oldgeodetic.dHeight = geodetic->dHeight;
q = sqrt (1 - BJ54e2 * sin(geodetic->dLatitude) * sin(geodetic->dLatitude));
N = BJ54a / q;
geodetic->dLatitude = atan2 (cartesian2.dZ*(N + geodetic->dHeight),p*(N * (1 - BJ54e2)+geodetic->dHeight));
geodetic->dHeight = cartesian2.dZ / sin(geodetic->dLatitude) - N * (1 - BJ54e2);
}while(fabs (oldgeodetic.dLatitude - geodetic->dLatitude) > EPSILON &&
fabs (oldgeodetic.dHeight - geodetic->dHeight) > EPSILON);
*/
Matrix x1(3,1);
Matrix zhuan(3,3);
Matrix X2(3,1);
Matrix deltax(3,1);
double m = -6.052572e-6;
zhuan(0,0) = 1;
zhuan(0,1) = -2.147541 / 206265.0000;
zhuan(0,2) = -4.955782 / 206265.0000;
zhuan(1,0) = 2.147541/206265.000;
zhuan(1,1) = 1;
zhuan(1,2) = -1.959986 / 206265.000;
zhuan(2,0) = 4.955782 / 206265.00;
zhuan(2,1) = 1.959986 / 206265.000;
zhuan(2,2) = 1;
deltax(0,0) = 131.9082;
deltax(1,0) = 204.8094;
deltax(2,0) = 87.7259;
X2(0,0) = cartesian2.dX;
X2(1,0) = cartesian2.dY;
X2(2,0) = cartesian2.dZ;
x1 = (1+m)*zhuan*X2 + deltax;
cartesian1->dX = x1 (0,0);
cartesian1->dY = x1 (1,0);
cartesian1->dZ = x1 (2,0);
}
void CCoordTrans::PZ90toWGS(COORDINATECARTESIAN cartesian2, COORDINATECARTESIAN *cartesian1)
{
Matrix x1(3,1);
Matrix zhuan(3,3);
Matrix X2(3,1);
Matrix deltax(3,1);
// double m = -6.052572e-6;
zhuan(0,0) = 1;
zhuan(0,1) = -1.9e-6 / 206265.0000;
zhuan(0,2) = 0;
zhuan(1,0) = 1.9e-6/206265.000;
zhuan(1,1) = 1;
zhuan(1,2) = 0;
zhuan(2,0) = 0;
zhuan(2,1) = 0;
zhuan(2,2) = 1;
deltax(0,0) = 0;
deltax(1,0) = 2.5;
deltax(2,0) = 0;
X2(0,0) = cartesian2.dX;
X2(1,0) = cartesian2.dY;
X2(2,0) = cartesian2.dZ;
x1 = zhuan*X2 + deltax;
cartesian1->dX = x1 (0,0);
cartesian1->dY = x1 (1,0);
cartesian1->dZ = x1 (2,0);
}
void CCoordTrans::WGSTOPLANE(COORDINATECARTESIAN cartesian,COORDINATEGEODETIC *geodetic,COORDINATEPLANE* plane,int *MiddleOfLat,int projection)
{
double M0;
double M2;
double M4;
double M6;
double M8;
double a0;
double a2;
double a4;
double a6;
double deltaL;
double X;
double t;
double YITA2;
double w;
double we;
long delt;
double B1;
double W;
double N;
double WGSep2;
geodetic->dLongtitude = atan(cartesian.dY/cartesian.dX);
if(geodetic->dLongtitude < 0)
geodetic->dLongtitude = geodetic->dLongtitude + PI;
geodetic->dLatitude = atan(cartesian.dZ / sqrt(cartesian.dX * cartesian.dX+cartesian.dY*cartesian.dY));
W = sqrt (1 - WGSe2 * sin(geodetic->dLatitude) * sin(geodetic->dLatitude));
N = WGSa / W;
do
{
B1 = geodetic->dLatitude;
W = sqrt (1 - WGSe2 * sin(B1) * sin(B1));
N = WGSa / W;
geodetic->dLatitude = atan((cartesian.dZ+N*WGSe2*sin(B1))/sqrt(cartesian.dX*cartesian.dX+cartesian.dY*cartesian.dY));
}while(fabs(geodetic->dLatitude-B1) > EPSILON);
geodetic->dHeight = sqrt(cartesian.dX * cartesian.dX+cartesian.dY*cartesian.dY)/
cos(geodetic->dLatitude)-WGSa/sqrt (1 - WGSe2 * sin(geodetic->dLatitude) * sin(geodetic->dLatitude));
if (projection == 3) {
*MiddleOfLat = floor ((geodetic->dLongtitude * 180.0 / PI + 1.5) / 3.0) * 3.0;
delt = (long) (*MiddleOfLat / 3.0) * 1000000;
}
if (projection == 6) {
*MiddleOfLat = (floor (geodetic->dLongtitude * 180.0 / (6.0 * PI)) + 1 ) * 6.0 - 3.0;
delt = (long) (floor (geodetic->dLongtitude * 180.0 / (6.0 * PI)) + 1 ) * 1000000;
}
M0 = WGSa * (1 - WGSe2);
M2 = 3 * WGSe2 * M0 / 2;
M4 = 5 * WGSe2 * M2 / 4;
M6 = 7 * WGSe2 * M4 / 6;
M8 = 9 * WGSe2 * M6 / 8;
a0 = M0 + M2 / 2.0 + 3.0 * M4 / 8.0 + 5.0 * M6 / 16.0 + 35.0 * M8 / 128.0;
a2 = M2 / 2.0 + M4 / 2.0 + 15.0 * M6 / 32.0 + 7.0 * M8 / 16.0;
a4 = M4 / 8.0 + 3.0 * M6 / 16.0 + 7.0 * M8 / 32.0;
a6 = M6 / 32.0 + M8 / 16.0;
X = a0 * geodetic->dLatitude - a2 * sin (2 * geodetic->dLatitude) / 2.0
+ a4 * sin (4 * geodetic->dLatitude) / 4.0 - a6 * sin (6 * geodetic->dLatitude) / 6.0;
WGSep2 = WGSe2 / (1 - WGSe2);
deltaL = geodetic->dLongtitude - *MiddleOfLat * PI / 180.0;
t = tan (geodetic->dLatitude);
YITA2 = WGSep2 * cos (geodetic->dLatitude) * cos (geodetic->dLatitude);
W = sqrt (1 - WGSe2 * sin(geodetic->dLatitude) * sin(geodetic->dLatitude));
N = WGSa / W;
w = cos (geodetic->dLatitude) * cos (geodetic->dLatitude) * cos (geodetic->dLatitude);
we = w * cos (geodetic->dLatitude) * cos (geodetic->dLatitude);
plane->dX = X + N * sin (geodetic->dLatitude) * cos (geodetic->dLatitude)
* deltaL * deltaL / 2.0 + N * sin (geodetic->dLatitude) * w * (5 - t * t
+ 9 * YITA2 + 4 * YITA2 * YITA2) * deltaL * deltaL * deltaL * deltaL / 24.0
+ N * sin (geodetic->dLatitude) * we * (61 - 58 * t * t + t * t * t * t) * deltaL * deltaL * deltaL * deltaL * deltaL * deltaL / 720.0;
plane->dY = N * cos (geodetic->dLatitude) * deltaL + N * w * (1 - t * t + YITA2) * deltaL * deltaL * deltaL / 6.0
+ N * we * (5 - 18 * t * t + t * t * t * t + 14 * YITA2 - 58 * YITA2 * t * t) * deltaL * deltaL * deltaL * deltaL * deltaL / 120.0 + 500000;// +delt;
}
void CCoordTrans::WGSTOPZ(COORDINATEGEODETIC wgdt, COORDINATECARTESIAN *wcar, COORDINATECARTESIAN *pcar, COORDINATEGEODETIC *pgdt)
{
double Nw;
double Np;
double PZ90b;
double pB0,pH0;
Nw = WGSa / sqrt (1- WGSe2 * sin(wgdt.dLatitude) * sin(wgdt.dLatitude));
wcar->dX = (Nw + wgdt.dHeight) * cos(wgdt.dLatitude) * cos(wgdt.dLongtitude);
wcar->dY = (Nw + wgdt.dHeight) * cos(wgdt.dLatitude) * sin(wgdt.dLongtitude);
wcar->dZ = (Nw * (1 - WGSe2) + wgdt.dHeight) * sin(wgdt.dLatitude);
pcar->dX = (wcar->dX + 1.9e-6 * (wcar->dY - 2.5)) / (1 + 1.9e-6 * 1.9e-6);
pcar->dY = (-1.9e-6 * wcar->dX + wcar->dY - 2.5) / (1 + 1.9e-6 * 1.9e-6);
pcar->dZ = wcar->dZ;
PZ90b = PZ90a * sqrt(1 - PZ90e2);
pgdt->dLongtitude = atan(pcar->dY / pcar->dX);
Np = PZ90a;
pgdt->dHeight = sqrt(pcar->dX * pcar->dX + pcar->dY * pcar->dY + pcar->dZ * pcar->dZ) - sqrt(PZ90a * PZ90b);
pgdt->dLatitude = atan(pcar->dZ / ((1 - PZ90e2 * Np / (Np + pgdt->dHeight)) * sqrt(pcar->dX * pcar->dX + pcar->dY * pcar->dY)));
do
{
pB0 = pgdt->dLatitude;
Np = PZ90a / sqrt(1 - PZ90e2 * sin(pB0) * sin(pB0));
pH0 = pgdt->dHeight;
pgdt->dHeight = sqrt(pcar->dX * pcar->dX + pcar->dY * pcar->dY) / cos(pB0) - Np;
pgdt->dLatitude = atan(pcar->dZ / ((1 - PZ90e2 * Np / (Np + pgdt->dHeight)) * sqrt(pcar->dX * pcar->dX + pcar->dY * pcar->dY)));
}while(fabs(pgdt->dLatitude - pB0) > 0.00001 * PI / (3600.0 *180.0) || fabs(pgdt->dHeight - pH0) > 0.001 );
}
////////////////////////////////大地坐标九参数转换////////////////////////////////////////////
void CCoordTrans::ninetrans(COMMANCOORD coord, NINEPARAMETER *np)
{
int i;
double rou = 206265.0;
double M;
double N;
double W;
Matrix x1(12,1);
Matrix x2(12,1);
Matrix B(12,7);
Matrix dx(7,1);
Matrix l(12,1);
Matrix v(12,1);
double dl,db,dh;
for(i = 0; i < 4; i ++)
{
x1(i*3,0) = coord.B1[i];//WGS84坐标系大地坐标
x1(i*3+1,0) = coord.L1[i];
x1(i*3+2,0) = coord.H1[i];
x2(i*3,0) = coord.B2[i];//北京54坐标系大地坐标
x2(i*3+1,0) = coord.L2[i];
x2(i*3+2,0) = coord.H2[i];
}
W = sqrt(1 - WGSe2 * sin(coord.B1[i]) * sin(coord.B1[i]));
M = WGSa * (1 - WGSe2) / (W * W * W);
N = WGSa / W;
for(i = 0; i < 4; i ++)
{
l(i*3,0) = coord.L2[i] - coord.L1[i];
l(i*3+1,0) = coord.B2[i] - coord.B1[i] + (WGSa - BJ54a) * N * WGSe2 * sin(coord.B1[i]) * cos(coord.B1[i]) * rou / ((M + coord.H1[i]) * WGSa) + (WGSF - BJ54F) * M * (2 - WGSe2 * sin(coord.B1[i]) * sin(coord.B1[i])) * sin(coord.B1[i]) * cos(coord.B1[i]) * rou / ((M + coord.H1[i]) * (1 - WGSF));
l(i*3+2,0) = coord.H2[i] - coord.H1[i] + (WGSa - BJ54a) * (-N) * (1 - WGSe2 * sin(coord.B1[i]) * sin(coord.B1[i])) / WGSa + (WGSF - BJ54F) * M * (1 - WGSe2 * sin(coord.B1[i]) * sin(coord.B1[i])) * sin(coord.B1[i]) * sin(coord.B1[i]) / (1 - WGSF);
}
for(i = 0; i < 4; i ++)
{
B(i*3,0) = -sin(coord.L1[i]) * rou / ((N + coord.H1[i]) * cos(coord.B1[i])); //WGS84下的直角坐标
B(i*3,1) = cos(coord.L1[i]) * rou / ((N + coord.H1[i]) * cos(coord.B1[i]));
B(i*3,2) = 0;
B(i*3,3) = tan(coord.B1[i]) * cos(coord.L1[i]);
B(i*3,4) = tan(coord.B1[i]) * sin(coord.L1[i]);
B(i*3,5) = -1;
B(i*3,6) = 0;
B(i*3+1,0) = -sin(coord.B1[i]) * cos(coord.L1[i]) * rou / (M + coord.H1[i]);
B(i*3+1,1) = -sin(coord.B1[i]) * sin(coord.L1[i]) * rou / (M + coord.H1[i]);
B(i*3+1,2) = cos(coord.B1[i]) * rou / (M + coord.H1[i]);
B(i*3+1,3) = -sin(coord.L1[i]);
B(i*3+1,4) = cos(coord.L1[i]);
B(i*3+1,5) = 0;
B(i*3+1,6) = -N * WGSe2 * sin(coord.B1[i]) * cos(coord.B1[i]) * rou / (M + coord.H1[i]);
B(i*3+2,0) = cos(coord.B1[i]) * cos(coord.L1[i]);
B(i*3+2,1) = cos(coord.B1[i]) * sin(coord.L1[i]);
B(i*3+2,2) = sin(coord.B1[i]);
B(i*3+2,3) = -N * WGSe2 * sin(coord.B1[i]) * cos(coord.B1[i]) * sin(coord.L1[i]) / rou;
B(i*3+2,4) = N * WGSe2 * sin(coord.B1[i]) * cos(coord.B1[i]) * cos(coord.L1[i]) / rou;
B(i*3+2,5) = 0;
B(i*3+2,6) = N * (1 - WGSe2 * sin(coord.B1[i]) * sin(coord.B1[i])) + coord.H1[i];
}
dx = (!(~B*B))*(~B)*l;
np->dx = dx(0,0);
np->dy = dx(1,0);
np->dz = dx(2,0);
np->ex = dx(3,0);
np->ey = dx(4,0);
np->ez = dx(5,0);
np->m = dx(6,0);
np->deltaa = WGSa - BJ54a;
np->deltaalpha = WGSF - BJ54F;;
v = B * dx + l;
dl = v(0,0);
db = v(1,0);
dh = v(2,0);
coord.L1[0] = coord.L1[0] - dl;
coord.B1[0] = coord.B1[0] - db;
coord.H1[0] = coord.H1[0] - dh;
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -