同一基准下坐标形式转换
测绘中各种坐标系：
1.大地坐标系：以参考椭球为基准面建立起来的坐标系。
将某点投影到椭球面上的位置用大地经纬度L和大地纬度表示B，（B,L,H）
2.空间直角坐标系：以椭球体中心为原点，起始子午面与赤道面为X轴，在赤道面上与X轴正交的方向为Y轴，椭球体的旋转轴为Z轴
某点位置用(X,Y,Z)表示
3.平面直角坐标系：通过高斯投影到平面上，以x轴为纵轴，y轴为横轴的平面直角坐标系表示(x,y)
大地坐标系---->空间直角坐标系：
/// <summary>
/// 大地坐标转空间直角坐标
/// </summary>
/// <param name="GeodeticPoint">大地坐标(度)</param>
/// <returns>空间直角坐标(米)</returns>
public SpacePoint BLH_XYZ(GeodeticPoint GeodeticPoint)
{
SpacePoint spacePoint = new SpacePoint();
double B = GeodeticPoint.B;
double L = GeodeticPoint.L;
double H = GeodeticPoint.H;
double a = this._ellipsoidParamets.a;
double e = Math.Sqrt(this._ellipsoidParamets.e2);
double Bhd = B * Math.PI / 180;
double Lhd = L * Math.PI / 180;
double BLH_XYZ_N = a / Math.Sqrt(1 - Math.Pow(e, 2) * Math.Pow(Math.Sin(Bhd), 2));
spacePoint.X = (BLH_XYZ_N + H) * Math.Cos(Bhd) * Math.Cos(Lhd);
spacePoint.Y = (BLH_XYZ_N + H) * Math.Cos(Bhd) * Math.Sin(Lhd);
spacePoint.Z = (BLH_XYZ_N * (1 - Math.Pow(e, 2)) + H) * Math.Sin(Bhd);
spacePoint.Name = GeodeticPoint.Name;
spacePoint.Property = GeodeticPoint.Property;
return spacePoint;
}
空间直角坐标----->大地坐标系：
/// <summary>
/// 空间直角坐标转大地坐标
/// </summary>
/// <param name="spacePoint">空间直角坐标(米)</param>
/// <returns>大地坐标(度)</returns>
public GeodeticPoint XYZ_BLH(SpacePoint spacePoint)
{
GeodeticPoint GeodeticPoint = new GeodeticPoint();
double X = spacePoint.X;
double Y = spacePoint.Y;
double Z = spacePoint.Z;
double a = this._ellipsoidParamets.a;
double e = Math.Sqrt(this._ellipsoidParamets.e2);
double R = Math.Sqrt(Math.Pow(X, 2) + Math.Pow(Y, 2));
double[] b = new double[20000];
b[0] = Math.Atan(Z / R);
int j = 1;
double N;
do
{
N = a / Math.Sqrt(1 - Math.Pow(e, 2) * Math.Pow(Math.Sin(b[j - 1]), 2));
b[j] = Math.Atan((Z + N * Math.Pow(e, 2) * Math.Sin(b[j - 1])) / R);
j++;
} while (Math.Abs(b[j - 1] - b[j - 2]) > 0.0000000001);
double[] BLH = new double[3];
BLH[0] = b[j - 1];
BLH[1] = Math.Atan(Y / X);
if (X < 0)
{
BLH[1] = BLH[1] + Math.PI;
}
BLH[2] = R / Math.Cos(BLH[0]) - N;
BLH[0] = BLH[0] * 180 / Math.PI;
BLH[1] = BLH[1] * 180 / Math.PI;
GeodeticPoint.B = BLH[0];
GeodeticPoint.L = BLH[1];
GeodeticPoint.H = BLH[2];
GeodeticPoint.Name = spacePoint.Name;
GeodeticPoint.Property = spacePoint.Property;
return GeodeticPoint;
}
大地坐标--->平面直角坐标(高斯正算):
1 /// <summary> 2 /// 大地坐标BLH转平面坐标xyh 3 /// </summary> 4 /// <param name="_L0">中央子午线(度)</param> 5 /// <param name="GeodeticPoint">大地坐标(度)</param> 6 /// <returns>平面坐标(米 y=y+500000)</returns> 7 public PlanePoint BLH_xyh(double _L0, GeodeticPoint GeodeticPoint) 8 { 9 PlanePoint planePoint = new PlanePoint(); 10 double B = GeodeticPoint.B; 11 double L = GeodeticPoint.L; 12 B = B * Math.PI / 180.0; 13 L = L * Math.PI / 180.0; 14 15 double e = Math.Sqrt(this._ellipsoidParamets.e2); 16 17 double AA, BB, CC, DD, EE, FF, GG, p1; 18 p1 = 3600.0 / Math.PI * 180.0; 19 AA = 1.0 + 3.0 / 4.0 * Math.Pow(e, 2) + 45.0 / 64.0 * Math.Pow(e, 4) + 175.0 / 256.0 * Math.Pow(e, 6) + 11025.0 / 16384.0 * Math.Pow(e, 8) + 43659.0 / 65536.0 * Math.Pow(e, 10) + 693693.0 / 1048576.0 * Math.Pow(e, 12); 20 BB = 3.0 / 8 * Math.Pow(e, 2) + 15.0 / 32.0 * Math.Pow(e, 4) + 525.0 / 1024.0 * Math.Pow(e, 6) + 2205.0 / 4096.0 * Math.Pow(e, 8) + 72765.0 / 131072.0 * Math.Pow(e, 10) + 297297.0 / 524288.0 * Math.Pow(e, 12); 21 CC = 15.0 / 256.0 * Math.Pow(e, 4) + 105.0 / 1024.0 * Math.Pow(e, 6) + 2205.0 / 16384.0 * Math.Pow(e, 8) + 10395.0 / 65536.0 * Math.Pow(e, 10) + 1486485.0 / 8388608.0 * Math.Pow(e, 12); 22 DD = 35.0 / 3072.0 * Math.Pow(e, 6) + 105.0 / 4096.0 * Math.Pow(e, 8) + 10395.0 / 262144.0 * Math.Pow(e, 10) + 55055.0 / 1048576.0 * Math.Pow(e, 12); 23 EE = 315.0 / 131072.0 * Math.Pow(e, 8) + 3465.0 / 524288.0 * Math.Pow(e, 10) + 99099.0 / 8388608.0 * Math.Pow(e, 12); 24 FF = 693.0 / 1310720.0 * Math.Pow(e, 10) + 9009.0 / 5242880.0 * Math.Pow(e, 12); 25 GG = 1001.0 / 8388608.0 * Math.Pow(e, 12); 26 27 double ee = Math.Sqrt(this._ellipsoidParamets.ee2); 28 double a = this._ellipsoidParamets.a; 29 double H, T, N, XX; 30 H = ee * Math.Cos(B); 31 T = Math.Tan(B); 32 N = a / Math.Sqrt(1 - Math.Pow(e, 2) * Math.Pow(Math.Sin(B), 2)); 33 XX = a * (1 - Math.Pow(e, 2)) * (AA * B - BB * Math.Sin(2.0 * B) + CC * Math.Sin(4.0 * B) - DD * Math.Sin(6.0 * B) + EE * Math.Sin(8.0 * B) - FF * Math.Sin(10.0 * B) + GG * Math.Sin(12.0 * B)); 34 35 double xl; 36 xl = L - _L0 * Math.PI / 180.0; 37 xl = xl * 180.0 / Math.PI * 3600.0; 38 39 double zjbl10; 40 zjbl10 = (0.5 + 1.0 / 24.0 * (5.0 - Math.Pow(T, 2) + 9.0 * Math.Pow(H, 2) + 4.0 * Math.Pow(H, 4)) * Math.Pow(Math.Cos(B), 2) * Math.Pow(xl, 2) / Math.Pow(p1, 2) + 1.0 / 720.0 * (61.0 - 58.0 * Math.Pow(T, 2) + Math.Pow(T, 4)) * Math.Pow(Math.Cos(B), 4) * Math.Pow(xl, 4) / Math.Pow(p1, 4)); 41 42 double x, y; 43 x = XX + N * T * Math.Pow(Math.Cos(B), 2) * Math.Pow(xl, 2) / Math.Pow(p1, 2) * zjbl10; 44 zjbl10 = (1.0 + 1.0 / 6.0 * (1.0 - Math.Pow(T, 2) + Math.Pow(H, 2)) * Math.Pow(Math.Cos(B), 2) * Math.Pow(xl, 2) / Math.Pow(p1, 2) + 1.0 / 120.0 * (5.0 - 18.0 * Math.Pow(T, 2) + Math.Pow(T, 4) + 14.0 * Math.Pow(H, 2) - 58.0 * Math.Pow(H, 2) * Math.Pow(T, 2)) * Math.Pow(Math.Cos(B), 4) * Math.Pow(xl, 4) / Math.Pow(p1, 4)); 45 y = N * Math.Cos(B) * xl / p1 * zjbl10; 46 y = y + 500000.0; 47 48 planePoint.x = x; 49 planePoint.y = y; 50 planePoint.h = GeodeticPoint.H; 51 planePoint.Name = GeodeticPoint.Name; 52 planePoint.Property = GeodeticPoint.Property; 53 54 return planePoint; 55 }
平面直角坐标---->大地坐标(高斯反算)
/// <summary> /// 平面坐标xyh(米 y=y+500000)转大地坐标BLH（度）----同一基准 /// </summary> /// <param name="_L0">中央子午线(度)</param> /// <param name="planePoint">平面坐标点(米 y=y+500000)</param> /// <returns>大地坐标(度)</returns> public GeodeticPoint xyh_BLH(double _L0, PlanePoint planePoint) { GeodeticPoint elloipsoidPoint = new GeodeticPoint(); double pi = Math.PI; double C = this._ellipsoidParamets.a / Math.Sqrt(1 - this._ellipsoidParamets.e2); double x2 = this._ellipsoidParamets.e2; double x4 = x2 * x2; double x6 = x2 * x4; double x8 = x4 * x4; double x10 = x2 * x8; double AA = 1.0 + 3.0 * x2 / 4.0 + 45.0 * x4 / 64.0 + 175.0 * x6 / 256.0; AA = AA + 11025.0 * x8 / 16384.0 + 43659.0 * x10 / 65536.0; double BB = 3.0 * x2 / 4.0 + 15.0 * x4 / 16.0 + 525.0 * x6 / 512.0; BB = BB + 2205.0 * x8 / 2048.0 + 72765.0 * x10 / 65536.0; double CC = 15.0 * x4 / 64.0 + 105.0 * x6 / 256.0; CC = CC + 2205.0 * x8 / 4096.0 + 10395.0 * x10 / 16384.0; double DD = 35.0 * x6 / 512.0 + 315.0 * x8 / 2048.0 + 31185.0 * x10 / 13072.0; double a1 = AA * this._ellipsoidParamets.a * (1 - this._ellipsoidParamets.e2); double a2 = -BB * this._ellipsoidParamets.a * (1 - this._ellipsoidParamets.e2) / 2.0; double a3 = CC * this._ellipsoidParamets.a * (1 - this._ellipsoidParamets.e2) / 4.0; double a4 = -DD * this._ellipsoidParamets.a * (1 - this._ellipsoidParamets.e2) / 6.0; double[] b11 = new double[6]; double[] r11 = new double[6]; double[] d11 = new double[6]; b11[1] = -a2 / a1; r11[1] = -a3 / a1; d11[1] = -a4 / a1; for (int i = 1; i <= 4; i++) { b11[i + 1] = b11[1] + b11[1] * r11[i] - 2.0 * r11[1] * b11[i] - 3.0 * b11[1] * b11[i] * b11[i] / 2.0; r11[i + 1] = r11[1] + b11[1] * b11[i]; d11[i + 1] = d11[1] + b11[1] * r11[i] + 2.0 * r11[1] * b11[i] + b11[1] * b11[i] * b11[i] / 2.0; } double K1, K2, K3; K1 = 2.0 * b11[5] + 4.0 * r11[5] + 6.0 * d11[5]; K2 = -8.0 * r11[5] - 32.0 * d11[5]; K3 = 32.0 * d11[5]; double y = planePoint.y - 500000.0; double x = planePoint.x; double Bf, Tf, y2, Vf2, Nf, M1, Mb; Bf = x / a1; Bf = Bf + Math.Cos(Bf) * Math.Sin(Bf) * (K1 + Math.Sin(Bf) * Math.Sin(Bf) * (K2 + Math.Sin(Bf) * Math.Sin(Bf) * K3)); Tf = Math.Tan(Bf); y2 = x2 / (1 - x2) * Math.Cos(Bf) * Math.Cos(Bf); Vf2 = 1.0 + x2 / (1.0 - x2) * Math.Cos(Bf) * Math.Cos(Bf); Nf = C / Math.Sqrt(Vf2); M1 = y / Nf; Mb = (y / Nf) * (y / Nf); double B, L; B = Bf - Mb * Vf2 * Tf / 2.0 + Mb * Mb * Vf2 * Tf / 24.0 * (5.0 + 3.0 * Tf * Tf + y2 - 9.0 * y2 * Tf * Tf) - Mb * Mb * Mb * Vf2 * Tf / 720.0 * (61.0 + (90.0 + 45.0 * Tf * Tf) * Tf * Tf); L = M1 / Math.Cos(Bf) - M1 * M1 * M1 / 6.0 / Math.Cos(Bf) * (1.0 + y2 + 2.0 * Tf * Tf) + M1 * Mb * Mb / 120.0 / Math.Cos(Bf) * (5.0 + (6.0 + 8.0 * Tf * Tf) * y2 + (28.0 + 24.0 * Tf * Tf) * Tf * Tf); L = L + _L0 * pi / 180.0; elloipsoidPoint.B = B * 180.0 / pi; elloipsoidPoint.L = L * 180.0 / pi; elloipsoidPoint.H = planePoint.h; elloipsoidPoint.Name = planePoint.Name; elloipsoidPoint.Property = planePoint.Property; return elloipsoidPoint; }
