需求说明:以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]));
}
}
}
测试结果: