坐标转换(c#)

C#

浏览数:300

2019-7-4

AD:资源代下载服务

坐标转换  

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

坐标转换关系

布尔莎模型(七参数模型)

详细见百科: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>
        /// 数据格式化,全转为弧度格式
        /// </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 = {X,Y,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 = {B,L,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 = {x,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 = {B,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);
            double[] BL2xy = BL_xy(xy2BL[0],xy2BL[1], l2, a, f1, xCon, yCon);
            return BL2xy;
        }
    }
}

   

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

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

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

作者:干脆面两毛五