NDBGJK5.0/Kingo.Plugins/Kingo.Plugin.ProofManager/Common/CoordinateTrans.cs

using ESRI.ArcGIS.Client.Geometry;
using System;
using System.Collections.Generic;
using System.ComponentModel;
using System.Data;
using System.IO;
using System.Runtime.Serialization.Json;
using System.Text;

namespace Kingo.Plugin.ProofManager.Common
{
    public class CoordinateTrans
    {


        #region 四参数集合操作

        /// <summary>
        /// 初始化四参数
        /// </summary>
        /// <param name="dbHelper"></param>
        public static void InitParams()
        {
            try
            {
                string selectSql = "SELECT * FROM TB_CSSZ";
                DataTable dt = null;// Kingo.SanDiaoCheck.Common.CommonAPI.GetSystemTable(selectSql);
                if (dt != null && dt.Rows.Count > 0)
                {
                    ParamsEntity hcryEntity = null;
                    foreach (DataRow item in dt.Rows)
                    {
                        hcryEntity = new ParamsEntity();
                        hcryEntity.BSM = ObjectConvertOperate.ObjectToInt(item["BSM"]);
                        hcryEntity.XZQDM = ObjectConvertOperate.ObjectToString(item["XZQDM"]);
                        hcryEntity.X80OFFSET = ObjectConvertOperate.ObjectToDouble(item["X80OFFSET"]);
                        hcryEntity.X80ROTATE = ObjectConvertOperate.ObjectToDouble(item["X80ROTATE"]);
                        hcryEntity.X84OFFSET = ObjectConvertOperate.ObjectToDouble(item["X84OFFSET"]);
                        hcryEntity.X84ROTATE = ObjectConvertOperate.ObjectToDouble(item["X84ROTATE"]);
                        hcryEntity.Y80OFFSET = ObjectConvertOperate.ObjectToDouble(item["Y80OFFSET"]);
                        hcryEntity.Y80ROTATE = ObjectConvertOperate.ObjectToDouble(item["Y80ROTATE"]);
                        hcryEntity.Y84OFFSET = ObjectConvertOperate.ObjectToDouble(item["Y84OFFSET"]);
                        hcryEntity.Y84ROTATE = ObjectConvertOperate.ObjectToDouble(item["Y84ROTATE"]);
                        hcryEntity.DH = ObjectConvertOperate.ObjectToString(item["DH"]);
                        AddOrUpdateParams(hcryEntity);
                    }
                }
            }
            catch (Exception ex)
            {
            }
        }

        /// <summary>
        /// 四参数集合
        /// </summary>
        private static Dictionary<string, ParamsEntity> paramsEntityList = new Dictionary<string, ParamsEntity>();
        /// <summary>
        /// 增加或修改四参数
        /// </summary>
        /// <param name="paramEntity"></param>
        /// <returns></returns>
        public static bool AddOrUpdateParams(ParamsEntity paramEntity)
        {
            lock (paramsEntityList)
            {
                if (paramsEntityList.ContainsKey(paramEntity.XZQDM))
                    paramsEntityList[paramEntity.XZQDM] = paramEntity;
                else
                    paramsEntityList.Add(paramEntity.XZQDM, paramEntity);
            }
            return true;
        }

        /// <summary>
        /// 获取四参数根据行政区
        /// </summary>
        /// <param name="xzq"></param>
        /// <returns></returns>
        public static ParamsEntity GetParamsByXzq(string xzq)
        {
            if (string.IsNullOrWhiteSpace(xzq)) return null;
            if (paramsEntityList.ContainsKey(xzq))
                return paramsEntityList[xzq];
            return null;
        }
        public static Param4 Get80To84(string xzq)
        {
            ParamsEntity paramsEntity = GetParamsByXzq(xzq);
            Param4 param4 = null;
            if (paramsEntity != null)
            {
                if (paramsEntity.X80OFFSET == 0 && paramsEntity.X80ROTATE == 0 && paramsEntity.Y80OFFSET == 0 && paramsEntity.Y80ROTATE == 0)
                    return null;
                param4 = new Param4();
                param4.PX = paramsEntity.X80OFFSET;
                param4.PY = paramsEntity.Y80OFFSET;
                param4.SX = paramsEntity.X80ROTATE;
                param4.SY = paramsEntity.Y80ROTATE;
            }
            return param4;
        }

        public static Param4 Get84To80(string xzq)
        {
            ParamsEntity paramsEntity = GetParamsByXzq(xzq);
            Param4 param4 = null;
            if (paramsEntity != null)
            {
                if (paramsEntity.X84OFFSET == 0 && paramsEntity.X84ROTATE == 0 && paramsEntity.Y84OFFSET == 0 && paramsEntity.Y84ROTATE == 0)
                    return null;
                param4 = new Param4();
                param4.PX = paramsEntity.X84OFFSET;
                param4.PY = paramsEntity.Y84OFFSET;
                param4.SX = paramsEntity.X84ROTATE;
                param4.SY = paramsEntity.Y84ROTATE;
            }
            return param4;
        }

        #endregion

        /// <summary>
        /// 计算四参数
        /// </summary>
        /// <param name="N">重合点个数</param>
        /// <param name="oldp">重合点在旧坐标系中的坐标</param>
        /// <param name="newp1">重合点的在新坐标系中的坐标</param>
        /// <returns></returns>
        private static double[][] Calculate4Param(int N, double[][] oldp, double[][] newp1)
        {
            int i, j;
            double[][] coneB;                      //系数矩阵B   
            coneB = new double[2 * N][];
            for (i = 0; i < 2 * N; i++)
                coneB[i] = new double[4];
            double[][] TconeB;                      //系数矩阵B的转置   
            TconeB = new double[4][];
            for (i = 0; i < 4; i++)
                TconeB[i] = new double[2 * N];
            double[][] consL;                  //常数项矩阵L   
            consL = new double[2 * N][];
            for (i = 0; i < 2 * N; i++)
                consL[i] = new double[1];
            double[][] Nbb;                     //法方程系数矩阵Nbb   
            Nbb = new double[4][];
            for (i = 0; i < 4; i++)
                Nbb[i] = new double[4];
            double[][] DNbb;                 //法方程系数矩阵Nbb的逆矩阵   
            DNbb = new double[4][];
            for (i = 0; i < 4; i++)
                DNbb[i] = new double[4];
            double[][] ENbb;                 //法方程系数矩阵Nbb的扩展矩阵   
            ENbb = new double[4][];
            for (i = 0; i < 4; i++)
                ENbb[i] = new double[8];
            double[][] W;               //用于表示Nbb*TconeB   
            W = new double[4][];
            for (i = 0; i < 4; i++)
                W[i] = new double[2 * N];
            double[][] parameter;              //平差参数   
            parameter = new double[4][];
            for (i = 0; i < 4; i++)
                parameter[i] = new double[1];
            double[][] V;                      //平差改正数   
            V = new double[2 * N][];
            for (i = 0; i < 2 * N; i++)
                V[i] = new double[1];
            double[][] result2;                //用于存储系数B与平差参数parameter的乘积   
            result2 = new double[2 * N][];
            for (i = 0; i < 2 * N; i++)
                result2[i] = new double[1];

            for (i = 0; i < 2 * N; i++)                        //写出方程系数矩阵B   
            {
                for (j = 0; j < 4; j++)
                {
                    coneB[i][0] = (0.5 + Math.Pow(-1, i) * 0.5);
                    coneB[i][1] = (0.5 - Math.Pow(-1, i) * 0.5);
                    coneB[i][2] = oldp[i / 2][0] * (0.5 + Math.Pow(-1, i) * 0.5) + oldp[i / 2][1] * (0.5 - Math.Pow(-1, i) * 0.5);
                    coneB[i][3] = oldp[i / 2][0] * (0.5 - Math.Pow(-1, i) * 0.5) + oldp[i / 2][1] * (-0.5 - Math.Pow(-1, i) * 0.5);
                }
            }
            for (i = 0; i < 2 * N; i++)              //写出常数项矩阵L   
                for (j = 0; j < 1; j++)
                    consL[i][j] = newp1[i / 2][0] * (0.5 + Math.Pow(-1, i) * 0.5) + newp1[i / 2][1] * (0.5 - Math.Pow(-1, i) * 0.5);
            transform(coneB, TconeB, 2 * N, 4);    //调用转置函数   
            multiplication(TconeB, coneB, Nbb, 4, 2 * N, 4);    //调用乘法函数计算法方程系数   
            det(Nbb, ENbb, DNbb, 4);           //调用求逆函数，计算法方程系数的逆阵   
            multiplication(TconeB, consL, W, 4, 2 * N, 1);//调用乘法函数，计算常数项W   
            multiplication(DNbb, W, parameter, 4, 4, 1);    //调用乘法函数,解算转换系数a,b,c,d  
            return parameter;
        }
        /// <summary>
        /// 四参数转换
        /// </summary>
        /// <param name="oldpoint"></param>
        /// <param name="parameter"></param>
        /// <param name="M"></param>
        /// <returns></returns>
        private static double[][] Tras4Param(double[][] oldpoint, double[][] parameter, int M)
        {
            int i;
            double[][] newpoint;                    //用于存储旧坐标经转换后的新坐标   
            newpoint = new double[M][];
            for (i = 0; i < M; i++)
            {
                newpoint[i] = new double[2];
            }
            for (i = 0; i < M; i++)
            {
                newpoint[i][0] = parameter[0][0] + oldpoint[i][0] * parameter[2][0] - oldpoint[i][1] * parameter[3][0];
                newpoint[i][1] = parameter[1][0] + oldpoint[i][1] * parameter[2][0] + oldpoint[i][0] * parameter[3][0];
            }
            return newpoint;
        }

        //求逆函数   
        private static void det(double[][] p, double[][] b, double[][] c, int n)
        {
            int i, j, k;
            double m;
            //扩展为增广矩阵   
            for (i = 0; i < n; i++)
            {
                for (j = 0; j < n; j++)
                {
                    b[i][j] = p[i][j];
                }
                for (j = n; j < 2 * n; j++)
                {
                    if (i + n == j)
                        b[i][j] = 1;
                    else
                        b[i][j] = 0;
                }
            }
            //将左半部分化为单位矩阵   
            for (i = 0; i < n; i++)
            {
                m = b[i][i];
                for (j = i; j < n * 2; j++)
                    b[i][j] /= m;
                for (j = i + 1; j < n; j++)
                {
                    m = b[j][i];
                    for (k = i; k < n * 2; k++)
                        b[j][k] = b[j][k] - m * b[i][k];
                }
            }
            for (i = n - 1; i >= 0; i--)
            {
                m = b[i][i];
                for (j = i - 1; j >= 0; j--)
                {
                    m = b[j][i];
                    for (k = i; k < n * 2; k++)
                        b[j][k] = b[j][k] - m * b[i][k];
                }
            }
            for (i = 0; i < n; i++)
            {
                for (j = n; j < 2 * n; j++)
                {
                    c[i][j - n] = b[i][j];
                }
            }
        }
        //转置函数   
        private static void transform(double[][] a, double[][] b, int m, int n)
        {
            int i, j;
            for (i = 0; i < m; i++)
                for (j = 0; j < n; j++)
                {
                    b[j][i] = a[i][j];
                }
        }
        //乘法函数   
        public static void multiplication(double[][] a, double[][] b, double[][] c, int m, int s, int n)
        {
            int i, j, k;
            double sum;
            for (i = 0; i < m; i++)
                for (j = 0; j < n; j++)
                {
                    sum = 0.0;
                    for (k = 0; k < s; k++)
                        sum += a[i][k] * b[k][j];
                    c[i][j] = sum;
                }
        }
        //矩阵加法函数   
        private static void add(double[][] a, double[][] b, double[][] ab, int m, int n, int sign)
        {
            for (int i = 0; i < m; i++)
                for (int j = 0; j < n; j++)
                    ab[i][j] = a[i][j] + Math.Pow(-1, sign) * b[i][j];     //sign为偶数是表示矩阵相加，奇数表示相减   
        }
        /// <summary>
        /// 根据同位点计算平面转换四参数
        /// </summary>
        /// <param name="sourcePoints"></param>
        /// <param name="targetPoints"></param>
        /// <returns></returns>
        public static Param4 Calculate4Param(MapPoint[] sourcePoints, MapPoint[] targetPoints)
        {
            if (sourcePoints == null || targetPoints == null || sourcePoints.Length != targetPoints.Length || sourcePoints.Length == 0)
                throw new ArgumentException("参数为空或者个数不一致", "sourcePoints,targetPoints");
            int n = sourcePoints.Length;
            double[][] oldp1 = new double[n][];
            double[][] newp1 = new double[n][];
            for (int i = 0; i < n; i++)
            {
                oldp1[i] = new double[] { sourcePoints[i].X, sourcePoints[i].Y };
                newp1[i] = new double[] { targetPoints[i].X, targetPoints[i].Y };
            }
            double[][] param4 = Calculate4Param(n, oldp1, newp1);
            return new Param4()
            {
                PX = param4[0][0],
                PY = param4[1][0],
                SX = param4[2][0],
                SY = param4[3][0]
            };
        }

        public static MapPoint[] Tras4Param(MapPoint[] sourcePoints, Param4 param4)
        {
            if (sourcePoints == null || sourcePoints.Length == 0)
                throw new ArgumentNullException("sourcePoints");
            if (param4 == null)
                throw new ArgumentNullException("param4");
            double[][] parameter = new double[4][];
            parameter[0] = new double[] { param4.PX };
            parameter[1] = new double[] { param4.PY };
            parameter[2] = new double[] { param4.SX };
            parameter[3] = new double[] { param4.SY };
            int m = sourcePoints.Length;
            double[][] oldp = new double[m][];
            for (int i = 0; i < m; i++)
            {
                oldp[i] = new double[] { sourcePoints[i].X, sourcePoints[i].Y };
            }
            double[][] newp = Tras4Param(oldp, parameter, m);
            MapPoint[] targetPoints = new MapPoint[m];
            for (int i = 0; i < m; i++)
            {
                targetPoints[i] = new MapPoint(newp[i][0], newp[i][1]);
            }
            return targetPoints;
        }
        /// <summary>
        /// 平面四参数转换方法
        /// </summary>
        /// <param name="sourcePoint"></param>
        /// <param name="param4"></param>
        /// <returns></returns>
        public static MapPoint Tras4Param(MapPoint sourcePoint, Param4 param4)
        {
            return Tras4Param(new MapPoint[] { sourcePoint }, param4)[0];
        }

        //高斯投影正、反算 
        //////3度带宽  
        //高斯投影由经纬度(Unit:DD)反算大地坐标(含带号，Unit:Metres) 
        private static void GaussProjCal(double longitude, double latitude, double a, double f, out double X, out double Y, int ProjNo = 0)
        {
            //int ProjNo = 0; 
            int ZoneWide;   ////带宽 
            double longitude1, latitude1, longitude0, latitude0, X0, Y0, xval, yval;
            double e2, ee, NN, T, C, A, M, iPI;
            iPI = Math.PI / 180;   ////3.1415926535898/180.0; 
            ZoneWide = 3;     ////3度带宽 
            //a = 6378245.0; f = 1.0 / 298.3;      //54年北京坐标系参数 
            //a=6378140.0;   f=1/298.257;      //80年西安坐标系参数 
            //a=6378137m;f=1/298.257223563;//WGS-84坐标系
            if (ProjNo == 0)
            {
                ProjNo = (int)Math.Round(longitude / ZoneWide);
            }
            longitude0 = ProjNo * ZoneWide;
            longitude0 = longitude0 * iPI;
            latitude0 = 0;
            longitude1 = longitude * iPI;   //经度转换为弧度
            latitude1 = latitude * iPI;     //纬度转换为弧度
            e2 = 2 * f - f * f;
            ee = e2 * (1.0 - e2);
            NN = a / Math.Sqrt(1.0 - e2 * Math.Sin(latitude1) * Math.Sin(latitude1));
            T = Math.Tan(latitude1) * Math.Tan(latitude1);
            C = ee * Math.Cos(latitude1) * Math.Cos(latitude1);
            A = (longitude1 - longitude0) * Math.Cos(latitude1);
            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));
            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);
            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);
            X0 = 1000000L * (ProjNo) + 500000L;
            Y0 = 0;
            xval = xval + X0; yval = yval + Y0;
            X = xval;
            Y = yval;
        }

        /// <summary>
        /// 高斯投影正算
        /// </summary>
        /// <param name="sourcePoint"></param>
        /// <param name="param"></param>
        /// <returns></returns>
        public static MapPoint GaussProjCal(MapPoint sourcePoint, EarthParam param, int ProjNo = 0)
        {
            double x, y;
            GaussProjCal(sourcePoint.X, sourcePoint.Y, param.A, param.F, out x, out y, ProjNo);
            return new MapPoint(x, y);
        }

        //高斯投影由大地坐标(Unit:Metres)反算经纬度(Unit:DD)
        private static void GaussProjInvCal(double X, double Y, double a, double f, out double longitude, out double latitude, int ProjNo = 0)
        {
            //int ProjNo; 
            int ZoneWide;   ////带宽 
            double longitude1, latitude1, longitude0, latitude0, X0, Y0, xval, yval;
            double e1, e2, ee, NN, T, C, M, D, R, u, fai, iPI;
            iPI = Math.PI / 180;   ////3.1415926535898/180.0; 
            //a = 6378245.0; f = 1.0 / 298.3;    //54年北京坐标系参数 
            //a=6378140.0;   f=1/298.257;      //80年西安坐标系参数 
            //a=6378137m;f=1/298.257223563;//WGS-84坐标系
            ZoneWide = 3;     ////3度带宽 
            if (ProjNo == 0)
            {
                ProjNo = (int)(X / 1000000L);    //查找带号
            }
            longitude0 = ProjNo * ZoneWide;
            longitude0 = longitude0 * iPI;     //中央经线
            X0 = ProjNo * 1000000L + 500000L;
            Y0 = 0;
            xval = X - X0; yval = Y - Y0;      //带内大地坐标
            e2 = 2 * f - f * f;
            e1 = (1.0 - Math.Sqrt(1 - e2)) / (1.0 + Math.Sqrt(1 - e2));
            ee = e2 / (1 - e2);
            M = yval;
            u = M / (a * (1 - e2 / 4 - 3 * e2 * e2 / 64 - 5 * e2 * e2 * e2 / 256));
            fai = u + (3 * e1 / 2 - 27 * e1 * e1 * e1 / 32) * Math.Sin(2 * u) + (21 * e1 * e1 / 16 - 55 * e1 * e1 * e1 * e1 / 32) * Math.Sin(4 * u)
             + (151 * e1 * e1 * e1 / 96) * Math.Sin(6 * u) + (1097 * e1 * e1 * e1 * e1 / 512) * Math.Sin(8 * u);
            C = ee * Math.Cos(fai) * Math.Cos(fai);
            T = Math.Tan(fai) * Math.Tan(fai);
            NN = a / Math.Sqrt(1.0 - e2 * Math.Sin(fai) * Math.Sin(fai));
            R = a * (1 - e2) / Math.Sqrt((1 - e2 * Math.Sin(fai) * Math.Sin(fai)) * (1 - e2 * Math.Sin(fai) * Math.Sin(fai)) * (1 - e2 * Math.Sin(fai) * Math.Sin(fai)));
            D = xval / NN;
            //计算经度(Longitude) 纬度(Latitude)
            longitude1 = longitude0 + (D - (1 + 2 * T + C) * D * D * D / 6 + (5 - 2 * C + 28 * T - 3 * C * C + 8 * ee + 24 * T * T) * D * D * D * D * D / 120) / Math.Cos(fai);
            latitude1 = fai - (NN * Math.Tan(fai) / R) * (D * D / 2 - (5 + 3 * T + 10 * C - 4 * C * C - 9 * ee) * D * D * D * D / 24
             + (61 + 90 * T + 298 * C + 45 * T * T - 256 * ee - 3 * C * C) * D * D * D * D * D * D / 720);
            //转换为度 DD
            longitude = longitude1 / iPI;
            latitude = latitude1 / iPI;
        }
        /// <summary>
        /// 高斯投影反算
        /// </summary>
        /// <param name="sourcePoint"></param>
        /// <param name="param"></param>
        /// <returns></returns>
        public static MapPoint GaussProjInvCal(MapPoint sourcePoint, EarthParam param, int ProjNo = 0)
        {
            double x, y;
            GaussProjInvCal(sourcePoint.X, sourcePoint.Y, param.A, param.F, out x, out y, ProjNo);
            return new MapPoint(x, y);
        }
        /// <summary>
        /// 西安80平面坐标转84经纬度坐标
        /// </summary>
        /// <param name="sourcePoint"></param>
        /// <param name="param"></param>
        /// <returns></returns>
        public static MapPoint Xian80ToWgs84(MapPoint sourcePoint, Param4 param = null, int ProjNo = 0)
        {
            if (sourcePoint == null)
                throw new ArgumentNullException("sourcePoint");
            if (param != null && param.SX != 0 && param.PX != 0)
            {
                MapPoint targetPoint = Tras4Param(sourcePoint, param);
                return GaussProjInvCal(targetPoint, EarthParam.WGS84, ProjNo);
            }
            else
            {
                return GaussProjInvCal(sourcePoint, EarthParam.WGS84, ProjNo);
            }
        }
        /// <summary>
        /// 84经纬度坐标转西安80平面坐标
        /// </summary>
        /// <param name="sourcePoint"></param>
        /// <param name="param"></param>
        /// <returns></returns>
        public static MapPoint Wgs84ToXian80(MapPoint sourcePoint, Param4 param = null, int ProjNo = 0)
        {
            if (sourcePoint == null)
            {
                throw new ArgumentNullException("sourcePoint");
            }
            MapPoint targetPoint = GaussProjCal(sourcePoint, EarthParam.WGS84, ProjNo);
            if (param != null && param.PX != 0 && param.SX != 0)
            {
                return Tras4Param(targetPoint, param);
            }
            else
            {
                return targetPoint;
            }
        }
        public static MapPoint Wgs84ToXian80(MapPoint sourcePoint, string xzqdm, int ProjNo = 0)
        {
            return Wgs84ToXian80(sourcePoint, Get84To80(xzqdm), ProjNo);
        }

        public static Envelope Xian80ToWgs84(Envelope source, Param4 param = null, int ProjNo = 0)
        {
            if (source == null)
            {
                throw new ArgumentNullException("source");
            }
            if (ProjNo == 0)
            {
                MapPoint pt = source.GetCenter();
                ProjNo = (int)(pt.X / 1000000L);
            }
            MapPoint ptLT = Xian80ToWgs84(new MapPoint(source.XMin, source.YMax), param, ProjNo);
            MapPoint ptRB = Xian80ToWgs84(new MapPoint(source.XMax, source.YMin), param, ProjNo);
            return new Envelope(ptLT.X, ptRB.Y, ptRB.X, ptLT.Y);
        }

        public static Envelope Wgs84ToXian80(Envelope source, Param4 param = null, int ProjNo = 0)
        {
            if (source == null)
                throw new ArgumentNullException("source");
            if (ProjNo == 0)
            {
                MapPoint pt = source.GetCenter();
                ProjNo = (int)Math.Round(pt.X / 3);
            }
            MapPoint ptLT = Wgs84ToXian80(new MapPoint(source.XMin, source.YMax), param, ProjNo);
            MapPoint ptRB = Wgs84ToXian80(new MapPoint(source.XMax, source.YMin), param, ProjNo);
            return new Envelope(ptLT.X, ptRB.Y, ptRB.X, ptLT.Y);
        }

        /// <summary>
        /// 度分秒转成度
        /// </summary>
        /// <param name="point"></param>
        /// <returns></returns>
        public static MapPoint DFM2D(MapPoint point)
        {
            return new MapPoint(DFM2D(point.X), DFM2D(point.Y));
        }
        /// <summary>
        /// 度转成度分秒
        /// </summary>
        /// <param name="point"></param>
        /// <returns></returns>
        public static MapPoint D2DFM(MapPoint point)
        {
            return new MapPoint(D2DFM(point.X), D2DFM(point.Y));
        }

        public static double DFM2D(double dfm)
        {
            double d = Math.Floor(dfm);
            double f = Math.Floor((dfm - d) * 100);
            double m = (dfm - d - f / 100) * 10000;
            return d + f / 60 + m / 3600;
        }

        public static double D2DFM(double degree)
        {
            double d = Math.Floor(degree);
            double f = Math.Floor((degree - d) * 60);
            double m = (degree - d - f / 60) * 3600;
            return d + f / 100 + m / 10000;
        }

        public static int Get80PrjNo(double x)
        {
            if (x > 1000000 && x < 70000000)
            {
                return (int)x / 1000000;
            }
            return 0;
        }
    }

    public class Param4 : INotifyPropertyChanged
    {
        public event PropertyChangedEventHandler PropertyChanged;

        protected void PropertyChangedNotify(string propName)
        {
            if (this.PropertyChanged != null)
                PropertyChanged(this, new PropertyChangedEventArgs(propName));
        }

        private double _PX;

        public double PX
        {
            get { return _PX; }
            set
            {
                _PX = value;
                PropertyChangedNotify("PX");
            }
        }
        private double _PY;

        public double PY
        {
            get { return _PY; }
            set { _PY = value; PropertyChangedNotify("PY"); }
        }
        private double _SX;

        public double SX
        {
            get { return _SX; }
            set { _SX = value; PropertyChangedNotify("SX"); }
        }
        private double _SY;

        public double SY
        {
            get { return _SY; }
            set { _SY = value; PropertyChangedNotify("SY"); }
        }

        public Param4 Reverse()
        {
            return new Param4()
            {
                PX = -this.PX,
                PY = -this.PY,
                SX = -this.SX,
                SY = -this.SY
            };
        }

        public bool IsEmpty
        {
            get
            {
                return PX == 0 && PY == 0 && SX == 0 && SY == 0;
            }
        }
    }

    public class EarthParam
    {
        /// <summary>
        /// 长轴 eg:6378137
        /// </summary>
        public double A { get; set; }
        /// <summary>
        /// 扁心率 eg:1 / 298.257223563
        /// </summary>
        public double F { get; set; }

        public static EarthParam WGS84 { get; private set; }
        public static EarthParam XIAN80 { get; private set; }
        public static EarthParam BJ54 { get; private set; }
        static EarthParam()
        {
            WGS84 = new EarthParam() { A = 6378137, F = 1 / 298.257223563 };
            XIAN80 = new EarthParam() { A = 6378140, F = 1.0 / 298.3 };
            BJ54 = new EarthParam() { A = 6378245.0, F = 1 / 298.257 };
        }
    }

    public class YPRMethod
    {
        public static double GetAngle(double yaw, double pitch, double roll)
        {
            yaw = yaw * Math.PI / 180;
            pitch = pitch * Math.PI / 180;
            roll = roll * Math.PI / 180;
            //MapPoint pt1 = new MapPoint(0, 0, 0, 0);
            //MapPoint pt2 = new MapPoint(0, 1, 0, 0);
            //MapPoint pt22 = new MapPoint(0, 1, 0, 0);
            //MapPoint pt3 = new MapPoint(1, 0, 0, 0);
            //rotate(pt1, yaw, pitch, roll);
            //rotate(pt2, yaw, pitch, roll);
            //rotate2(pt22, yaw, pitch, roll);
            //rotate(pt3, yaw, pitch, roll);
            //MapPoint vec = get_Normal(pt1, pt2, pt3);
            //if (vec.X == 0 && vec.Y == 0)
            //{
            //    return (360 + yaw / Math.PI * 180) % 360;
            //}
            //double angle = Math.Atan2(vec.Y, vec.X);
            //if (double.IsNaN(angle))
            //{
            //    return angle;
            //}
            //else
            //{
            //    return (360 + angle * 180 / Math.PI) % 360;
            //}
            MapPoint pt = RotateOL(yaw, pitch, roll);
            double angle = Math.Atan2(pt.X, pt.Z);
            if (double.IsNaN(angle) || angle == 0)
            {
                return yaw * 180 / Math.PI;
            }
            else
            {
                return -angle * 180 / Math.PI;
            }
            //return 0;
        }

        private static void rotate(MapPoint pt, double yaw, double pitch, double roll)
        {
            double x1 = pt.X * Math.Cos(yaw) - pt.Y * Math.Sin(yaw);
            double y1 = pt.X * Math.Sin(yaw) + pt.Y * Math.Cos(yaw);
            double z1 = pt.Z;
            double y2 = y1 * Math.Cos(roll) - z1 * Math.Sin(roll);
            double z2 = y1 * Math.Sin(roll) + z1 * Math.Cos(roll);
            double x2 = x1;
            pt.Z = z2 * Math.Cos(pitch) - x2 * Math.Sin(pitch);
            pt.X = z2 * Math.Sin(pitch) + x2 * Math.Cos(pitch);
            pt.Y = y2;
        }

        private static void rotate2(MapPoint pt, double yaw, double pitch, double roll)
        {
            double y1 = pt.Y * Math.Cos(roll) - pt.Z * Math.Sin(roll);
            double z1 = pt.Y * Math.Sin(roll) + pt.Z * Math.Cos(roll);
            double x1 = pt.X;

            double x2 = x1 * Math.Cos(yaw) - y1 * Math.Sin(yaw);
            double y2 = x1 * Math.Sin(yaw) + y1 * Math.Cos(yaw);
            double z2 = z1;

            pt.Z = z2 * Math.Cos(pitch) - x2 * Math.Sin(pitch);
            pt.X = z2 * Math.Sin(pitch) + x2 * Math.Cos(pitch);
            pt.Y = y2;
        }

        private static MapPoint get_Normal(MapPoint p1, MapPoint p2, MapPoint p3)
        {

            double a = ((p2.Y - p1.Y) * (p3.Z - p1.Z) - (p2.Z - p1.Z) * (p3.Y - p1.Y));

            double b = ((p2.Z - p1.Z) * (p3.X - p1.X) - (p2.X - p1.X) * (p3.Z - p1.Z));

            double c = ((p2.X - p1.X) * (p3.Y - p1.Y) - (p2.Y - p1.Y) * (p3.X - p1.X));

            return new MapPoint(a, b, c, 0);
        }

        private static MapPoint RotateOL(double yaw, double pitch, double roll)
        {
            //double[][] mb = GetMaxtrixB(roll);
            //double[][] mc = GetMaxtrixC(pitch);
            //double[][] md = GetMaxtrixD(yaw);
            double[][] ma = GetMaxtrixA(pitch, yaw, roll);

            double[][] xyz = new double[3][];
            xyz[0] = new double[1] { 0 };
            xyz[1] = new double[1] { 1 };
            xyz[2] = new double[1] { 0 };
            //double[][] bc = new double[3][];
            //for (int i = 0; i < 3; i++)
            //{
            //    bc[i] = new double[3];
            //}
            //CoordinateTrans.multiplication(mb, mc, bc, 3, 3, 3);
            //double[][] bcd = new double[3][];
            //for (int i = 0; i < 3; i++)
            //{
            //    bcd[i] = new double[3];
            //}
            //CoordinateTrans.multiplication(bc, md, bcd, 3, 3, 3);

            double[][] xyz1 = new double[3][];
            for (int i = 0; i < 3; i++)
            {
                xyz1[i] = new double[1];
            }
            CoordinateTrans.multiplication(ma, xyz, xyz1, 3, 3, 1);
            return new MapPoint(xyz1[0][0], xyz1[1][0], xyz1[2][0], 0);
        }

        private static double[][] GetMaxtrixB(double d)
        {
            double[][] mxtrix = new double[3][];
            mxtrix[0] = new double[] { Math.Cos(d), Math.Sin(d), 0 };
            mxtrix[1] = new double[] { -Math.Sin(d), Math.Cos(d), 0 };
            mxtrix[2] = new double[] { 0, 0, 1 };
            return mxtrix;
        }

        private static double[][] GetMaxtrixC(double d)
        {
            double[][] mxtrix = new double[3][];
            mxtrix[0] = new double[] { 1, 0, 0 };
            mxtrix[1] = new double[] { 0, Math.Cos(d), Math.Sin(d) };
            mxtrix[2] = new double[] { 0, -Math.Sin(d), Math.Cos(d) };
            return mxtrix;
        }

        private static double[][] GetMaxtrixD(double d)
        {
            double[][] mxtrix = new double[3][];
            mxtrix[0] = new double[] { Math.Cos(d), Math.Sin(d), 0 };
            mxtrix[1] = new double[] { -Math.Sin(d), Math.Cos(d), 0 };
            mxtrix[2] = new double[] { 0, 0, 1 };
            return mxtrix;
        }
        private static double[][] GetMaxtrixA(double x, double y, double z)
        {
            double[][] mxtrix = new double[3][];
            mxtrix[0] = new double[] {
                Math.Cos(y)*Math.Cos(z)-Math.Sin(x)*Math.Sin(y)*Math.Sin(z),
                Math.Cos(y)*Math.Sin(z)+Math.Sin(x)*Math.Sin(y)*Math.Cos(z),
                -Math.Sin(y)*Math.Cos(x) };
            mxtrix[1] = new double[] {
                -Math.Cos(x)*Math.Sin(z),
                Math.Cos(x)*Math.Cos(z),
                Math.Sin(x) };
            mxtrix[2] = new double[] {
                Math.Sin(y)*Math.Cos(z)+Math.Sin(x)*Math.Cos(y)*Math.Sin(z),
                Math.Sin(y)*Math.Sin(z)+Math.Sin(x)*Math.Cos(y)*Math.Cos(z),
                Math.Cos(x)*Math.Cos(y) };
            return mxtrix;
        }
    }

    /// <summary>
    /// 墨卡托与经纬度转换操作
    /// </summary>
    public class MercatorToLonLat
    {
        /// <summary>
        /// 经纬度转墨卡托坐标
        /// </summary>
        /// <param name="lonLat"></param>
        /// <returns></returns>
        public static ESRI.ArcGIS.Client.Geometry.MapPoint lonLat2Mercator(ESRI.ArcGIS.Client.Geometry.MapPoint lonLat)
        {
            double x1 = 0, y1 = 0;
            transform(lonLat.Y, lonLat.X, out y1, out x1);
            lonLat = new MapPoint(x1, y1);
            ESRI.ArcGIS.Client.Geometry.MapPoint mercator = new ESRI.ArcGIS.Client.Geometry.MapPoint();
            double x = lonLat.X * 20037508.34 / 180;
            double y = Math.Log(Math.Tan((90 + lonLat.Y) * Math.PI / 360)) / (Math.PI / 180);
            y = y * 20037508.34 / 180;
            mercator.X = x;
            mercator.Y = y;
            return mercator;
        }
        /// <summary>
        /// 墨卡托坐标转经纬度
        /// </summary>
        /// <param name="mercator"></param>
        /// <returns></returns>
        public static ESRI.ArcGIS.Client.Geometry.MapPoint Mercator2lonLat(ESRI.ArcGIS.Client.Geometry.MapPoint mercator)
        {
            ESRI.ArcGIS.Client.Geometry.MapPoint lonLat = new ESRI.ArcGIS.Client.Geometry.MapPoint();
            double x = mercator.X / 20037508.34 * 180;
            double y = mercator.Y / 20037508.34 * 180;
            y = 180 / Math.PI * (2 * Math.Atan(Math.Exp(y * Math.PI / 180)) - Math.PI / 2);
            lonLat.X = x;
            lonLat.Y = y;
            return lonLat;
        }

        public static ESRI.ArcGIS.Client.Geometry.Envelope Mercator2lonLat(ESRI.ArcGIS.Client.Geometry.Envelope mercator)
        {
            double xMax = MercatorX2Lon(mercator.XMax);
            double xMin = MercatorX2Lon(mercator.XMin);
            double yMax = MercatorY2Lat(mercator.YMax);
            double yMin = MercatorY2Lat(mercator.YMin);
            return new Envelope()
            {
                XMax = xMax,
                XMin = xMin,
                YMax = yMax,
                YMin = yMin
            };
        }

        public static ESRI.ArcGIS.Client.Geometry.Envelope LonLat2Mercator(ESRI.ArcGIS.Client.Geometry.Envelope lonlat)
        {
            double xMax = Lon2Mercator(lonlat.XMax);
            double xMin = Lon2Mercator(lonlat.XMin);
            double yMax = Lat2Mercator(lonlat.YMax);
            double yMin = Lat2Mercator(lonlat.YMin);
            return new Envelope()
            {
                XMax = xMax,
                XMin = xMin,
                YMax = yMax,
                YMin = yMin
            };
        }

        private static double MercatorX2Lon(double x)
        {
            return x / 20037508.34 * 180;
        }

        private static double MercatorY2Lat(double y)
        {
            double lat = y / 20037508.34 * 180;
            return 180 / Math.PI * (2 * Math.Atan(Math.Exp(lat * Math.PI / 180)) - Math.PI / 2);
        }

        private static double Lon2Mercator(double x)
        {
            return x * 20037508.34 / 180;
        }

        private static double Lat2Mercator(double y)
        {
            double mercatorY = Math.Log(Math.Tan((90 + y) * Math.PI / 360)) / (Math.PI / 180);
            return mercatorY * 20037508.34 / 180;
        }


        #region  经纬度转火星坐标
        const double pi = 3.14159265358979324;

        //
        // Krasovsky 1940
        //
        // a = 6378245.0, 1/f = 298.3
        // b = a * (1 - f)
        // ee = (a^2 - b^2) / a^2;
        const double a = 6378245.0;
        const double ee = 0.00669342162296594323;

        //
        // World Geodetic System ==> Mars Geodetic System
        public static void transform(double wgLat, double wgLon, out double mgLat, out double mgLon)
        {
            if (outOfChina(wgLat, wgLon))
            {
                mgLat = wgLat;
                mgLon = wgLon;
                return;
            }
            double dLat = transformLat(wgLon - 105.0, wgLat - 35.0);
            double dLon = transformLon(wgLon - 105.0, wgLat - 35.0);
            double radLat = wgLat / 180.0 * pi;
            double magic = Math.Sin(radLat);
            magic = 1 - ee * magic * magic;
            double sqrtMagic = Math.Sqrt(magic);
            dLat = (dLat * 180.0) / ((a * (1 - ee)) / (magic * sqrtMagic) * pi);
            dLon = (dLon * 180.0) / (a / sqrtMagic * Math.Cos(radLat) * pi);
            mgLat = wgLat + dLat;
            mgLon = wgLon + dLon;
        }

        static bool outOfChina(double lat, double lon)
        {
            if (lon < 72.004 || lon > 137.8347)
                return true;
            if (lat < 0.8293 || lat > 55.8271)
                return true;
            return false;
        }

        static double transformLat(double x, double y)
        {
            double ret = -100.0 + 2.0 * x + 3.0 * y + 0.2 * y * y + 0.1 * x * y + 0.2 * Math.Sqrt(Math.Abs(x));
            ret += (20.0 * Math.Sin(6.0 * x * pi) + 20.0 * Math.Sin(2.0 * x * pi)) * 2.0 / 3.0;
            ret += (20.0 * Math.Sin(y * pi) + 40.0 * Math.Sin(y / 3.0 * pi)) * 2.0 / 3.0;
            ret += (160.0 * Math.Sin(y / 12.0 * pi) + 320 * Math.Sin(y * pi / 30.0)) * 2.0 / 3.0;
            return ret;
        }

        static double transformLon(double x, double y)
        {
            double ret = 300.0 + x + 2.0 * y + 0.1 * x * x + 0.1 * x * y + 0.1 * Math.Sqrt(Math.Abs(x));
            ret += (20.0 * Math.Sin(6.0 * x * pi) + 20.0 * Math.Sin(2.0 * x * pi)) * 2.0 / 3.0;
            ret += (20.0 * Math.Sin(x * pi) + 40.0 * Math.Sin(x / 3.0 * pi)) * 2.0 / 3.0;
            ret += (150.0 * Math.Sin(x / 12.0 * pi) + 300.0 * Math.Sin(x / 30.0 * pi)) * 2.0 / 3.0;
            return ret;
        }
        #endregion

        #region 百度坐标转火星坐标

        const double x_pi = 3.14159265358979324 * 3000.0 / 180.0;
        static MapPoint LonLatToBaidu(MapPoint mapPointLon)
        {
            MapPoint rst = new MapPoint();
            double bd_lat = 0;
            double bd_lon = 0;
            LonLatToBaidu(mapPointLon.X, mapPointLon.Y, ref bd_lat, ref bd_lon);
            rst.X = bd_lon;
            rst.Y = bd_lat;
            return rst;
        }

        static void LonLatToBaidu(double lon, double lat, ref double bd_lat, ref double bd_lon)
        {
            double x = lon, y = lat;
            double z = Math.Sqrt(x * x + y * y) + 0.00002 * Math.Sin(y * x_pi);
            double theta = Math.Atan2(y, x) + 0.000003 * Math.Cos(x * x_pi);
            bd_lon = z * Math.Cos(theta) + 0.0065;
            bd_lat = z * Math.Sin(theta) + 0.006;
        }

        static void BaiduToLonLat(double bd_lat, double bd_lon, ref double lon, ref double lat)
        {
            double x = bd_lon - 0.0065, y = bd_lat - 0.006;
            double z = Math.Sqrt(x * x + y * y) - 0.00002 * Math.Sin(y * x_pi);
            double theta = Math.Atan2(y, x) - 0.000003 * Math.Cos(x * x_pi);
            lon = z * Math.Cos(theta);
            lat = z * Math.Sin(theta);
        }
        #endregion
    }



    /// <summary>
    /// 四参数实体
    /// </summary>
    public class ParamsEntity
    {
        /// <summary>
        /// 标识码
        /// </summary>
        public int BSM { get; set; }
        /// <summary>
        /// 行政区代码
        /// </summary>
        public string XZQDM { get; set; }
        /// <summary>
        /// X平移84转80
        /// </summary>
        public double X84OFFSET { get; set; }
        /// <summary>
        /// Y平移84转80
        /// </summary>
        public double Y84OFFSET { get; set; }
        /// <summary>
        /// X旋转84转80
        /// </summary>
        public double X84ROTATE { get; set; }
        /// <summary>
        /// Y旋转84转80
        /// </summary>
        public double Y84ROTATE { get; set; }
        /// <summary>
        /// X平移80转84
        /// </summary>
        public double X80OFFSET { get; set; }
        /// <summary>
        /// Y平移80转84
        /// </summary>
        public double Y80OFFSET { get; set; }
        /// <summary>
        /// X旋转80转84
        /// </summary>
        public double X80ROTATE { get; set; }
        /// <summary>
        /// Y旋转80转84
        /// </summary>
        public double Y80ROTATE { get; set; }
        /// <summary>
        ///代号
        /// </summary>
        public string DH { get; set; }



    }


    public class ObjectConvertOperate
    {


        /// <summary>
        /// Object转string
        /// </summary>
        /// <param name="obj"></param>
        /// <returns></returns>
        public static string ObjectToString(object obj)
        {
            if (obj == null)
                return "";
            return obj.ToString();
        }

        /// <summary>
        /// Object转double
        /// </summary>
        /// <param name="obj"></param>
        /// <returns></returns>
        public static double ObjectToDouble(object obj)
        {
            if (obj == null)
                return 0.0;
            double rst = 0.0;
            double.TryParse(obj.ToString(), out rst);
            return rst;
        }

        /// <summary>
        /// Object转double
        /// </summary>
        /// <param name="obj"></param>
        /// <returns></returns>
        public static int ObjectToInt(object obj)
        {
            if (obj == null)
                return 0;
            int rst = 0;
            int.TryParse(obj.ToString(), out rst);
            return rst;
        }

        /// <summary>
        /// Object转bool
        /// </summary>
        /// <param name="obj"></param>
        /// <returns></returns>
        public static bool ObjectToBool(object obj)
        {
            if (obj == null)
                return false;
            bool rst = false;
            bool.TryParse(obj.ToString(), out rst);
            return rst;
        }

        /// <summary>
        /// Object转datetime
        /// </summary>
        /// <param name="obj"></param>
        /// <returns></returns>
        public static DateTime ObjectToDate(object obj)
        {
            if (obj == null)
                return new DateTime();
            DateTime dt = new DateTime();
            DateTime.TryParse(obj.ToString(), out dt);
            return dt;
        }


    }

    public class JsonRings
    {
        public List<List<double[]>> rings;

        public static JsonRings GetTemp()
        {
            JsonRings jr = new JsonRings() { rings = new List<List<double[]>>() };
            for (int i = 0; i < 10; i++)
            {
                List<double[]> ring = new List<double[]>();
                for (int j = 0; j < 10; j++)
                {
                    ring.Add(new double[] { 1111, 1111 });
                }
                jr.rings.Add(ring);
            }
            return jr;
        }

        public string GetJsonString()
        {
            var serializer = new DataContractJsonSerializer(this.GetType());
            var stream = new MemoryStream();
            serializer.WriteObject(stream, this);
            byte[] dataBytes = new byte[stream.Length];
            stream.Position = 0;
            stream.Read(dataBytes, 0, (int)stream.Length);
            return Encoding.UTF8.GetString(dataBytes);
        }

        public static JsonRings FromJson(string json)
        {
            try
            {
                if (string.IsNullOrWhiteSpace(json)) return null;
                var serializer = new DataContractJsonSerializer(typeof(JsonRings));
                var mStream = new MemoryStream(Encoding.Default.GetBytes(json));
                return (JsonRings)serializer.ReadObject(mStream);
            }
            catch
            {
                return null;
            }

        }
    }
}