C#程式設計練習(02):大地座標系(LBH)向空間直角座標系(XYZ)的轉換及其逆轉換
阿新 • • 發佈:2018-12-19
需求說明:以WGS-84軟體為例,實現大地座標系(LBH)向空間直角座標系(XYZ)的轉換及其逆轉換
原理說明:
程式原始碼:
using System; using System.Collections.Generic; using System.Linq; using System.Text; namespace XYZ2BLH { class Program { //輸出為度格式的字串 (度分秒.秒的小數部分) //例如:35度23分32.9秒 表述為352332.9 static string[] XYZ2LBH(double x,double y,double z) { double[] blh = {0.0, 0.0, 0.0}; double a = 6378137; double b = 6356752.3142; double e2 = 0.0066943799013; double ePie2 = 0.00673949674227; double xita = Math.Atan(z * a / (Math.Sqrt(x * x + y * y) * b)); double xitaSin3 = Math.Sin(xita) * Math.Sin(xita) * Math.Sin(xita); double xitaCos3 = Math.Cos(xita) * Math.Cos(xita) * Math.Cos(xita); double B = Math.Atan((z + ePie2 * b * xitaSin3) / ((Math.Sqrt(x * x + y * y) - e2 * a * xitaCos3))); double L = Math.Atan(y / x); double N = a / Math.Sqrt(1 - e2 * Math.Sin(B) * Math.Sin(B)); ; double H = Math.Sqrt(x * x + y * y) / Math.Cos(B) - N; //轉化為角度 B = 180 * B / Math.PI; if (B < 0.0) { B = B + 90; } L = 180 * L / Math.PI; if (L < 0.0) { L = L + 180; } return new string[] {rad2dms(B), rad2dms(L), H.ToString()}; } //獲取一個小數的整數部分和小數部分 static double[] modf(double t) { t += 1.0e-14;//避免角度轉換過程產生無效近似 double t_integer = Math.Floor(t); return new double[] { t_integer, t - t_integer }; } //角度轉弧度 static double turn1(double t) { double[] t1 = modf(t); double[] t2 = modf(t1[1] * 100); double x1 = t1[0]; double x2 = t2[0]; double x3 = t2[1] * 100; //Console.Write("{0:N} \t {1:N}\t {2:N}\r\n", x1, x2, x3); return (x1 + x2 / 60 + x3 / 3600) / 180 * Math.PI; } //弧度轉角度 static string rad2dms(double rad) { int du, fen; double miao; du = (int)rad; fen = (int)((rad - (double)du) * 60); miao = ((rad - (double)du) * 60 - (double)fen) * 60; return du.ToString() + fen.ToString().PadLeft(2,'0') + miao.ToString(); } //經度、緯度和大地高轉向空間直角座標 static double[] LBH2XYZ(double[] lbh) { double B = turn1(lbh[1]); double L = turn1(lbh[0]); double H = lbh[2]; double A = 1 / 298.257223563; double a = 6378137.0000000000; double E = 2 * A - A * A; double W = Math.Sqrt(1 - E * Math.Sin(B) * Math.Sin(B)); double N = a / W; double[] xyz = { 0.0, 0.0, 0.0 }; xyz[0] = (N + H) * Math.Cos(B) * Math.Cos(L); xyz[1] = (N + H) * Math.Cos(B) * Math.Sin(L); xyz[2] = (N * (1 - E) + H) * Math.Sin(B); return xyz; } static void Main(string[] args) { double[] lbh = { 108.5743, 32.2352, 100.59 }; double[] xyz = LBH2XYZ(lbh); Console.WriteLine("{0:N6} {1:N6} {2:N6}", xyz[0], xyz[1], xyz[2]); string[] lbhTmp = XYZ2LBH(xyz[0], xyz[1], xyz[2]); Console.WriteLine("{0:N6} {1:N6} {2:N6}", double.Parse(lbhTmp[0]), double.Parse(lbhTmp[1]), double.Parse(lbhTmp[2])); } } }
測試結果: