坐标转换1------同基准坐标换算
测绘中各种坐标系:
1.大地坐标系:以参考椭球为基准面建立起来的坐标系。
将某点投影到椭球面上的位置用大地经纬度L和大地纬度表示B,(B,L,H)
2.空间直角坐标系:以椭球体中心为原点,起始子午面与赤道面为X轴,在赤道面上与X轴正交的方向为Y轴,椭球体的旋转轴为Z轴
某点位置用(X,Y,Z)表示
3.平面直角坐标系:通过高斯投影到平面上,以x轴为纵轴,y轴为横轴的平面直角坐标系表示(x,y)
大地坐标系---->空间直角坐标系:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
|
/// <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; } |
空间直角坐标----->大地坐标系:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
|
/// <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; } |
大地坐标--->平面直角坐标(高斯正算):
/// <summary> /// 大地坐标BLH转平面坐标xyh /// </summary> /// <param name="_L0">中央子午线(度)</param> /// <param name="GeodeticPoint">大地坐标(度)</param> /// <returns>平面坐标(米 y=y+500000)</returns> public PlanePoint BLH_xyh(double _L0, GeodeticPoint GeodeticPoint) { PlanePoint planePoint = new PlanePoint(); double B = GeodeticPoint.B; double L = GeodeticPoint.L; B = B * Math.PI / 180.0; L = L * Math.PI / 180.0; double e = Math.Sqrt(this._ellipsoidParamets.e2); double AA, BB, CC, DD, EE, FF, GG, p1; p1 = 3600.0 / Math.PI * 180.0; 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); 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); 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); 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); 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); FF = 693.0 / 1310720.0 * Math.Pow(e, 10) + 9009.0 / 5242880.0 * Math.Pow(e, 12); GG = 1001.0 / 8388608.0 * Math.Pow(e, 12); double ee = Math.Sqrt(this._ellipsoidParamets.ee2); double a = this._ellipsoidParamets.a; double H, T, N, XX; H = ee * Math.Cos(B); T = Math.Tan(B); N = a / Math.Sqrt(1 - Math.Pow(e, 2) * Math.Pow(Math.Sin(B), 2)); 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)); double xl; xl = L - _L0 * Math.PI / 180.0; xl = xl * 180.0 / Math.PI * 3600.0; double zjbl10; 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)); double x, y; x = XX + N * T * Math.Pow(Math.Cos(B), 2) * Math.Pow(xl, 2) / Math.Pow(p1, 2) * zjbl10; 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)); y = N * Math.Cos(B) * xl / p1 * zjbl10; y = y + 500000.0; planePoint.x = x; planePoint.y = y; planePoint.h = GeodeticPoint.H; planePoint.Name = GeodeticPoint.Name; planePoint.Property = GeodeticPoint.Property; return planePoint; }
平面直角坐标---->大地坐标(高斯反算)
/// <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; }