C#计算wgs84大地坐标转换为空间直角坐标

名词解释：

        大地坐标表示：BLH(大地纬度B、大地经度L、大地高H)

        空间坐标表示：XYZ

正常流程：

        将BLH转换为XYZ,然后将XYZ通过四(七)参数的办法转换为xyz

WGS84坐标系的参数

长半轴：a = 6378137

WGS84椭圆扁率：f = 1/298.257223563(f = (a-b)/a)

短半轴：b = a * (1 - f);

BLH转XYZ的公式

其中公式中e为椭圆的第一偏心率，N为椭圆为曲率半径,计算参数为第一辅助系数W

核心代码（C#实现）

double a = 6378137;
double f = 1 / 298.257223563;
double b = a * (1 - f);
double e = Math.Sqrt(a * a - b * b) / a;
double N = a / Math.Sqrt(1 - e * e * Math.Sin(lat * Math.PI / 180) * Math.Sin(lat * Math.PI / 180));
double WGS84_X = (N + H1) * Math.Cos(lat * Math.PI / 180) * Math.Cos(lon * Math.PI / 180);
double WGS84_Y = (N + H1) * Math.Cos(lat * Math.PI / 180) * Math.Sin(lon * Math.PI / 180);
double WGS84_Z = (N * (1 - (e * e)) + H1) * Math.Sin(lat * Math.PI / 180);

XYZ => xyz (xyz是经过投影转换的后的坐标)

 //WGS84坐标系下6度带投影坐标
double iPI = Math.PI / 180.0; 3.1415926535898/180.0;   
double ZoneWide = 6; 6度带宽   
double a1 = 6378245.0; f1 = 1.0 / 298.257223563; //WGS84坐标系参数   
 //a1 = 6378245.0; f1 = 1.0 / 298.3; //54年北京坐标系参数
 //a2=6378140.0; f2=1/298.257; //80年西安坐标系参数   
double ProjNo = (int)(lon / ZoneWide);
double longitude0 = ProjNo * ZoneWide + ZoneWide / 2;
double longitude0 = longitude0 * iPI;
double latitude0 = 0;
double longitude1 = lon * iPI; //经度转换为弧度  
double latitude1 = lat * iPI; //纬度转换为弧度  
double e2 = 2 * f - f * f;
double ee = e2 * (1.0 - e2);
double NN = a / Math.Sqrt(1.0 - e2 * Math.Sin(latitude1) * Math.Sin(latitude1));
double T = Math.Tan(latitude1) * Math.Tan(latitude1);
double C = ee * Math.Cos(latitude1) * Math.Cos(latitude1);
double A = (longitude1 - longitude0) * Math.Cos(latitude1);
 //参数转换
double  M = a * ((1 - e2 / 4 - 3 * e2 * e2 / 64 - 5 * e2 * e2 * e2 / 256) * latitude1 - (3 * e2 / 8 + 3 * e2 * e2 / 32 + 45 * e2 * e2 * e2 / 1024) * Math.Sin(2 * latitude1) + (15 * e2 * e2 / 256 + 45 * e2 * e2 * e2 / 1024) * Math.Sin(4 * latitude1) - (35 * e2 * e2 * e2 / 3072) * Math.Sin(6 * latitude1));
double xval = NN * (A + (1 - T + C) * A * A * A / 6 + (5 - 18 * T + T * T + 72 * C - 58 * ee) * A * A * A * A * A / 120);
double yval = M + NN * Math.Tan(latitude1) * (A * A / 2 + (5 - T + 9 * C + 4 * C * C) * A * A * A * A / 24 + (61 - 58 * T + T * T + 600 * C - 330 * ee) * A * A * A * A * A * A / 720);
double X0 = 1000000L * (ProjNo + 1) + 500000L;
double Y0 = 0;
xval = xval + X0;
yval = yval + Y0;
double yval1 = Convert.ToDouble(xval.ToString().Substring(2, xval.ToString().Length - 3));
double xval1 = Convert.ToDouble(yval.ToString().Substring(0, yval.ToString().Length - 5));