坐标转换（c#） 原 荐

坐标转换  

坐标转换，简言之就是讲一个坐标系的坐标通过特定模型转换到另一坐标下；转换中，需要用到不同的转换模型，例如三维空间下的布尔莎七参数、莫洛金斯模型；二维空间下的平面三参、四参、多项式拟合等模型。本文中，主要使用布尔莎七参数模型进行坐标转换。

坐标转换关系

布尔莎模型（七参数模型）

详细见百科： https://baike.baidu.com/item/%E4%B8%83%E5%8F%82%E6%95%B0

程序实现

1.编写数据转换类

把各种格式表示的度分秒转换为弧度形式进行计算，最后一个方法把弧度转为度分秒连写的格式。

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
using System.Text.RegularExpressions;
namespace CoordTF
{ 
    class DataCnversion
    {
        #region 
        /// <summary>
        /// 数据格式化，全转为度分秒（x.xxxx）连写格式
        /// </summary>
        /// <param name="oldCoord">大地坐标</param>
        /// <param name="dataFormat">数据格式</param>
        /// <returns></returns>
        public double DataFormatting(string oldCoord,string dataFormat)
        {
            double newCoord = 0.00;
            switch (dataFormat)
            {
                case "度°分'秒″":
                    newCoord = DMSA_RAD(oldCoord);
                    break;
                case "度:分:秒":
                    newCoord = DMSB_RAD(oldCoord);
                    break;
                case "度:分":

                    break;
                case "度":
                    newCoord = D_RAD(oldCoord); 
                    break;
                case "度分秒连写":
                    newCoord = DMSLX_RAD(oldCoord);
                    break;
            }
            return newCoord;
        }
        #endregion
        #region 各形式角度转弧度
        double d, m, s;     //度分秒
        /// <summary>
        /// 度°分'秒"格式转弧度
        /// </summary>
        /// <param name="dms">格式为XX ° XX ' XX.XX " </param>
        /// <returns></returns>
        public double DMSA_RAD(string dms)
        {
            double rad;
            var newD = dms.Substring(0, dms.IndexOf("°"));
            var newM = dms.Substring(dms.IndexOf("°")+1,dms.IndexOf("′")-3);
            var newS = dms.Substring(dms.IndexOf("′")+1,dms.Length-dms.IndexOf("′")-2);
            rad = (double.Parse(newD) + double.Parse(newM)/ 60 + double.Parse(newS) / 3600) * 
            Math.PI / 180;
            return rad;
        }
        /// <summary>
        /// 度:分:秒 格式转弧度
        /// </summary>
        /// <param name="dms">格式为度:分:秒</param>
        /// <returns></returns>
        public double DMSB_RAD(string dms)
        {
            double rad;
            string[] newDms = dms.Split(':');
            rad = (double.Parse(newDms[0].ToString()) + double.Parse(newDms[1].ToString()) / 60 + 
            double.Parse(newDms[2].ToString()) / 3600) * Math.PI / 180;
            return rad;
        }
        /// <summary>
        /// 度分秒连写（x.xxxx）转弧度
        /// </summary>
        /// <param name="dms">格式为x.xxxx</param>
        /// <returns></returns>
        public double DMSLX_RAD(string dms)
        {
            double rad;
            var newDms = double.Parse(dms);
            d = Math.Floor(newDms); //度
            var x = (newDms - d) * 100;
            x = double.Parse(x.ToString("N8"));
            m = Math.Floor(x);  //分
            s = (x - m) * 100;
            m /= 60; //秒
            s /= 3600;
            rad = (d + m + s) / 180 * (Math.PI);    //弧度
            return rad;
        }
        /// <summary>
        /// xx.xxxx° 格式转弧度
        /// </summary>
        /// <param name="dms">角度，格式为XX.XXXX°</param>
        /// <returns></returns>
        public double D_RAD(string dms)
        {
            double rad;
            rad = (double.Parse(dms) * Math.PI / 180);
            return rad;
        }
        #endregion
        /// <summary>
        /// 弧度转角度（结果为xx.xxxxxx连写格式）
        /// </summary>
        /// <param name="rad">弧度值</param>
        /// <returns></returns>
        public double RAD_DMS(string rad)
        {
            double dms;
            var newRad = double.Parse(rad);
            var d = newRad / Math.PI * 180; //度
            var dStr = d.ToString("N8");
            d = double.Parse(dStr);
            var x = Math.Floor(d);
            var m = (d - x) * 60;   //分
            var mStr = m.ToString("N8");
            m = double.Parse(mStr);
            var f = Math.Floor(m);
            var s = (m - f) * 60;   //秒
            var sStr = s.ToString("N8");
            s = double.Parse(sStr);
            dms = x + f / 100 + s / 10000;
            return dms;
        }
    }
}

2.编写坐标转换类

          包含空间与大地互转、高斯正反算、同一坐标系下换带计算、布尔莎七参数转换计算。

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
using System.Windows.Forms;
using MathNet.Numerics.LinearAlgebra.Double;
namespace CoordTF
{
    class CoordTransform
    {
        /// <summary>
        /// 大地坐标转空间直角
        /// </summary>
        /// <param name="B">纬度</param>
        /// <param name="L">经度</param>
        /// <param name="H">高度</param>
        /// <param name="a">椭球长半轴</param>
        /// <param name="fl">扁率倒数</param>
        /// <returns></returns>
        public double[] BLH_XYZ(double B, double L, double H, double a, double fl,string 
        dataFormat)
        {
            var f = 1 / fl;
            var e2 = 2 * f - f * f;
            DataCnversion dataCnversion = new DataCnversion();
            B = dataCnversion.DataFormatting(B.ToString(), dataFormat);
            L = dataCnversion.DataFormatting(L.ToString(), dataFormat);
            var W = Math.Sqrt(1 - e2 * (Math.Sin(B) * Math.Sin(B)));
            var N = a / W;
            var X = (N + H) * Math.Cos(B) * Math.Cos(L);
            var Y = (N + H) * Math.Cos(B) * Math.Sin(L);
            var Z = (N * (1 - e2) + H) * Math.Sin(B);
            double[] XYZ = new double[3];
            XYZ[0] = X;
            XYZ[1] = Y;
            XYZ[2] = Z;
            return XYZ;
        }
        /// <summary>
        /// 空间直角转大地坐标
        /// </summary>
        /// <param name="X">空间X</param>
        /// <param name="Y">空间Y</param>
        /// <param name="Z">空间Z</param>
        /// <param name="a">椭球长半轴</param>
        /// <param name="fl">扁率倒数</param>
        /// <returns></returns>
        public double[] XYZ_BLH(double X, double Y, double Z, double a, double fl)
        {
            var f = 1 / fl;
            var e2 = 2 * f - f * f;
            var S = Math.Sqrt(X * X + Y * Y);
            var L = Math.Acos(X / S);
            if (Y < 0)
                L = -L;
            var Rho = 206264.806247096363;
            var Error = 0.00001;
            var B1 = Math.Atan(Z / S);
            double B2;
            while (true)
            {
                var W1 = Math.Sqrt(1 - e2 * (Math.Sin(B1) * Math.Sin(B1)));
                var N1 = a / W1;
                B2 = Math.Atan((Z + N1 * e2 * Math.Sin(B1)) / S);
                if (Math.Abs(B2 - B1) * Rho <= Error)
                    break;
                B1 = B2;
            }
            var B = B2;
            var W = Math.Sqrt(1 - e2 * (Math.Sin(B) * Math.Sin(B)));
            var N = a / W;
            var H = S / Math.Cos(B) - N;
            DataCnversion dataCnversion = new DataCnversion();
            B = dataCnversion.RAD_DMS(B.ToString("N8"));
            L = dataCnversion.RAD_DMS(L.ToString("N8"));
            double[] BLH = new double[3];
            BLH[0] = B;
            BLH[1] = L;
            BLH[2] = H;
            return BLH;
        }
        /// <summary>
        /// 布尔莎七参转换
        /// </summary>
        /// <param name="X">空间直角X</param>
        /// <param name="Y">空间直角Y</param>
        /// <param name="Z">空间直角Z</param>
        /// <param name="dx">平移参数x</param>
        /// <param name="dy">平移参数y</param>
        /// <param name="dz">平移参数z</param>
        /// <param name="rpx">旋转参数x</param>
        /// <param name="rpy">旋转参数y</param>
        /// <param name="rpz">旋转参数z</param>
        /// <param name="k">尺度参数k</param>
        /// <returns></returns>
        public double[] BursaTF(double X, double Y, double Z, double dx, double dy, double dz, 
        double rpx, double rpy, double rpz, double k)
        {
            var Rho = 206264.806247096355;//常数          
            double A15, A16, A17, A24, A26, A27, A34, A35, A37;
            A15 = Z / Rho;  A16 = Y / Rho;  A17 = X / 1000000;
            A24 = Z / Rho;  A26 = X / Rho;  A27 = Y/ 1000000;
            A34 = Y / Rho;  A35 = X / Rho;  A37 = Z / 1000000;
            var newX = X + dx - A15 * rpy + A16 * rpz + A17 * k;
            var newY = Y + dy + A24 * rpx - A26 * rpz + A27 * k;
            var newZ = Z + dz - A34 * rpx + A35 * rpy + A37 * k;
            double[] result = { newX, newY, newZ };

            return result;
        }
        /// <summary>
        /// 高斯正算
        /// </summary>
        /// <param name="B">大地纬度</param>
        /// <param name="L">大地经度</param>
        /// <param name="L0">中央子午线经度(度°分′秒″格式)</param>
        /// <param name="a">椭球长半径</param>
        /// <param name="f1">扁率倒数</param>
        /// <param name="xCon">x常参</param>
        /// <param name="yCon">y常参</param>
        /// <returns></returns>
        public double[] BL_xy(double B, double L, double L0, double a, double f1, double xCon, 
        double yCon)
        {
            DataCnversion dataCnversion = new DataCnversion();
            var f = 1 / f1;
            var e2 = 2 * f - f * f;
            var e12 = e2 / (1 - e2);
            B = dataCnversion.DMSLX_RAD(B.ToString("N8"));
            L = dataCnversion.DMSLX_RAD(L.ToString("N8"));
            L0 = dataCnversion.DMSLX_RAD(L0.ToString("N8"));
            var l = L - L0;
            var A0 = 1 + 3 * e2 / 4 + 45 * Math.Pow(e2, 2) / 64 + 350 * Math.Pow(e2, 3) / 512 + 
            11025 * Math.Pow(e2, 4) / 16384;
            var A2 = -(3 * e2 / 4 + 60 * Math.Pow(e2, 2) / 64 + 525 * Math.Pow(e2, 3) / 512 + 17640 
            * Math.Pow(e2, 4) / 16384) / 2;
            var A4 = (15 * Math.Pow(e2, 2) / 64 + 210 * Math.Pow(e2, 3) / 512 + 8820 * Math.Pow(e2, 
            4) / 16384) / 4;
            var A6 = -(35 * Math.Pow(e2, 3) / 512 + 2520 * Math.Pow(e2, 4) / 16384) / 6;
            var A8 = (315 * Math.Pow(e2, 4) / 16384) / 8;
            var X = a * (1 - e2) * (A0 * B + A2 * Math.Sin(2 * B) + A4 * Math.Sin(4 * B) + A6 * 
            Math.Sin(6 * B) + A8 * Math.Sin(8 * B));
            var m0 = l * Math.Cos(B);
            var t = Math.Tan(B);
            var n2 = e12 * Math.Pow(Math.Cos(B), 2);
            var W = Math.Sqrt(1 - e2 * Math.Pow(Math.Sin(B), 2));
            var N = a / W;
            var Nt = N * t;
            double x = X + (Nt * Math.Pow(m0, 2) / 2) + (Nt * (5 - Math.Pow(t, 2) + (9 * n2) + (4 * 
            Math.Pow(n2, 2))) * Math.Pow(m0, 4) / 24) + (Nt * (61 - (58 * Math.Pow(t, 2)) + 
            Math.Pow(t, 4) + (270 * n2) - (330 * n2 * Math.Pow(t, 2))) * Math.Pow(m0, 6) / 720);
            double y = (N * m0) + (N * (1 - Math.Pow(t, 2) + n2) * Math.Pow(m0, 3) / 6) + (N * (5 - 
            (18 * Math.Pow(t, 2)) + Math.Pow(t, 4) + (14 * n2) - (58 * n2 * Math.Pow(t, 2))) * 
            Math.Pow(m0, 5) / 120);
            x += xCon;
            y += yCon;
            double[] point = new double[2];
            point[0] = x;
            point[1] = y;
            return point;
        }
        /// <summary>
        /// 高斯反算
        /// </summary>
        /// <param name="x">平面坐标x</param>
        /// <param name="y">平面坐标y</param>
        /// <param name="L0">中央子午线经度</param>
        /// <param name="a">椭球长半轴</param>
        /// <param name="f1">扁率倒数</param>
        /// <param name="xCon">x常参数</param>
        /// <param name="yCon">y常参数</param>
        /// <returns></returns>
        public double[] Xy_BL(double x, double y, double L0, double a, double f1, double xCon, 
        double yCon)
        {
            var k = 1000000;
            x -= xCon;
            if (Math.Abs(y) >= k)
            {
                var zoneNumber = Math.Floor(y / k);
                y -= zoneNumber;
            }
            y -= yCon;
            var f = 1 / f1;
            var e2 = 2 * f - f * f;
            var e12 = e2 / (1 - e2);
            DataCnversion dataCnversion = new DataCnversion();
            L0 = dataCnversion.DMSLX_RAD(L0.ToString("N8"));
            var A0 = 1 + 3 * e2 / 4 + 45 * Math.Pow(e2, 2) / 64 + 350 * Math.Pow(e2, 3) / 512 + 
            11025 * Math.Pow(e2, 4) / 16384;
            var B0 = x / (a * (1 - e2) * A0);
            var K0 = (3 * e2 / 4 + 45 * Math.Pow(e2, 2) / 64 + 350 * Math.Pow(e2, 3) / 512 + 11025 
            * Math.Pow(e2, 4) / 16384) / 2;
            var K2 = -(63 * Math.Pow(e2, 2) / 64 + 1108 * Math.Pow(e2, 3) / 512 + 58239 * 
            Math.Pow(e2, 4) / 16384) / 3;
            var K4 = (604 * Math.Pow(e2, 3) / 512 + 68484 * Math.Pow(e2, 4) / 16384) / 3;
            var K6 = -(26328 * Math.Pow(e2, 4) / 16384) / 3;
            var sB0 = Math.Pow(Math.Sin(B0), 2);
            var Bf = B0 + Math.Sin(2 * B0) * (K0 + sB0 * (K2 + sB0 * (K4 + K6 * sB0)));
            var nf2 = e12 * Math.Pow(Math.Cos(Bf), 2);
            var tf = Math.Tan(Bf);
            var tf2 = tf * tf;
            var tf4 = Math.Pow(tf, 4);
            var Wf = Math.Sqrt(1 - e2 * (Math.Pow(Math.Sin(Bf), 2)));
            var Nf = a / Wf;
            var yNf = y / Nf;
            var Vf2 = 1 + nf2;
            var B = Bf - 0.5 * Vf2 * tf * (Math.Pow(yNf, 2) - (5 + 3 * tf2 + nf2 - 9 * nf2 * tf2) * 
            Math.Pow(yNf, 4) / 12 + (61 + 90 * tf2 + 45 * tf4) * Math.Pow(yNf, 6) / 360);
            var l = (yNf - (1 + 2 * tf2 + nf2) * Math.Pow(yNf, 3) / 6 + (5 + 28 * tf2 + 24 * tf4 + 
            6 * nf2 + 8 * nf2 * tf2) * Math.Pow(yNf, 5) / 120) / Math.Cos(Bf);
            var L = L0 + l;
            B = dataCnversion.RAD_DMS(B.ToString("N8"));
            L = dataCnversion.RAD_DMS(L.ToString("N8"));
            double[] point = new double[2];
            point[0] = B;
            point[1] = L;
            return point;
        }
        /// <summary>
        /// 同一坐标系下坐标换带
        /// </summary>
        /// <param name="x"></param>
        /// <param name="y"></param>
        /// <param name="l1">换带前经度</param>
        /// <param name="l2">换带后经度</param>
        /// <param name="a"></param>
        /// <param name="f1"></param>
        /// <param name="xCon"></param>
        /// <param name="yCon"></param>
        /// <returns></returns>
        public double[] CoordChangeBelt(double x, double y, double l1, double l2, double a, double 
        f1, double xCon, double yCon)
        {
            double[] xy2BL = Xy_BL(x, y, l1, a, f1, xCon, yCon);
            var B = xy2BL[0];
            var L = xy2BL[1];
            double[] BL2xy = BL_xy(B, L, l2, a, f1, xCon, yCon);
            return BL2xy;
        }
    }
}

到此，所有的基础转换类完成，使用时依据坐标转换关系，依次调用类中的方法即可。

本人仅与中海达 HGO 数据处理软件包中CoordTool工具，做过七参数转换对比测试，同样参数下转换后坐标差值如下图（数据就不贴出来了）。