基于C#实现的高斯正反算公式实现的国家大地2000坐标转换功能

基于C#实现的高斯正反算公式实现的国家大地2000坐标转换功能

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

namespace ConsoleApp_CoordConvert2000
{
    public interface IProjectionConversion2
    {
        void GetBLFromXY(decimal x, decimal y, ref decimal B, ref decimal L);
        void GetXYFromBL(decimal B, decimal L, ref decimal x, ref decimal y);
    }
    public enum EnumStrip
    {
        Strip3=3,
        Strip6=6,
    }
    public abstract class ProjectionConversion2 : IProjectionConversion2
    {
        protected decimal a = 6378137M;           //2000椭球长半轴
        protected decimal b = 6356752.3141M;      //2000椭球短半轴
        protected decimal f = 1 / 298.257222101M;     //椭球扁率f
        protected decimal e = 0.0818191910428M;       //第一偏心率
        protected decimal e2 = 6.693421622966E-03M;   //椭球第一偏心率平方 e2
        protected decimal e12 = 6.738525414683E-03M;  //椭球第二偏心率平方 e12
        protected decimal c = 6399593.62586M;         //极点子午圈曲率半径 c
        protected decimal p = 206264.806252992M;      //弧度秒=180*3600/pi
        protected decimal pi = 3.1415926535M;

        //需要三个配置参数
        private bool m_IsBigNumber = true;
        private EnumStrip m_Strip = EnumStrip.Strip3;
        private decimal m_L0 = 105;

        public bool IsBigNumber
        {
            get
            {
                return m_IsBigNumber;
            }
            set
            {
                m_IsBigNumber = value;
            }
        }
        public EnumStrip Strip
        {
            get
            {
                return m_Strip;
            }
            set
            {
                m_Strip = value;
            }
        }
        public decimal L0
        {
            get
            {
                return m_L0;
            }
            set
            {
                m_L0 = value;
            }
        }
        protected double toNum(decimal num)
        {
            //return double.Parse(num.ToString());
            return (double)num;
        }
        protected decimal toDec(double num)
        {
            //return decimal.Parse(num.ToString());
            return (decimal)num;
        }
        protected decimal sin(decimal num)
        {
            return this.toDec(Math.Sin(this.toNum(num)));
        }
        protected decimal cos(decimal num)
        {
            return this.toDec(Math.Cos(this.toNum(num)));
        }
        protected decimal sqrt(decimal num)
        {
            return this.toDec(Math.Sqrt(this.toNum(num)));
        }
        protected decimal tan(decimal num)
        {
            return this.toDec(Math.Tan(this.toNum(num)));
        }
        public  decimal pow(decimal x, decimal y)
        {
            return toDec(Math.Pow(toNum(x), toNum(y)));
        }

        //高斯反算方法  OK   (x,y)=>(B,L)
        public virtual  void GetBLFromXY(decimal x, decimal y, ref decimal B, ref decimal L)
        {
            //去掉大数和东移500公里
            decimal y1 = y - 500000.0M;
            if (this.IsBigNumber == true)
            {
                y1 = y1 - (this.L0 / (int)this.Strip) * 1000000M;
            }
            y = y1;
            decimal l0 = this.L0 * 3600;  //中央子午线转为秒值 如=105*3600
            //计算临时值
            //decimal e4 = e2 * e2;
            //decimal e6 = e4 * e2;
            //decimal e8 = e6 * e2;
            //decimal e10 = e8 * e2;
            //decimal e_12 = e10 * e2;
            decimal e4 = pow(e2, 2); // e2 * e2;
            decimal e6 = pow(e2, 3); //e4 * e2;
            decimal e8 = pow(e2, 4); //e6 * e2;
            decimal e10 = pow(e2, 5); //e8 * e2;
            decimal e_12 = pow(e2, 6); //e10 * e2;
            //
            decimal A1 = 1 + 3 * e2 / 4 + 45 * e4 / 64 + 175 * e6 / 256 + 11025 * e8 / 16384 + 43659 * e10 / 65536 + 693693 * e_12 / 1048576;
            decimal B1 = 3 * e2 / 8 + 15 * e4 / 32 + 525 * e6 / 1024 + 2205 * e8 / 4096 + 72765 * e10 / 131072 + 297297 * e_12 / 524288;
            decimal C1 = 15 * e4 / 256 + 105 * e6 / 1024 + 2205 * e8 / 16384 + 10395 * e10 / 65536 + 1486485 * e_12 / 8388608;
            decimal D1 = 35 * e6 / 3072 + 105 * e8 / 4096 + 10395 * e10 / 262144 + 55055 * e_12 / 1048576;
            decimal E1 = 315 * e8 / 131072 + 3465 * e10 / 524288 + 99099 * e_12 / 8388608;
            decimal F1 = 693 * e10 / 1310720 + 9009 * e_12 / 5242880;
            decimal G1 = 1001 * e_12 / 8388608;
            //求底点纬度值Bf
            decimal B0 = x / (a * (1 - e2) * A1);
            decimal Bf = 0.0M;
            decimal FB = 0.0M;
            decimal FB1 = 0.0M;
            decimal a0 = a * (1 - e2);
            decimal delta = Math.Abs(Bf - B0);
            while (delta > 4.8E-11M)   //0.000000000048M
            {
                Bf = B0;
                FB = a0 * (A1 * Bf - B1 * sin(2 * Bf) + C1 * sin(4 * Bf) - D1 * sin(6 * Bf) + E1 * sin(8 * Bf) - F1 * sin(10 * Bf) + G1 * sin(12 * Bf));
                FB1 = a0 * (A1 - 2 * B1 * cos(2 * Bf) + 4 * C1 * cos(4 * Bf) - 6 * D1 * cos(6 * Bf) + 8 * E1 * cos(8 * Bf) - 10 * F1 * cos(10 * Bf) + 12 * G1 * cos(12 * Bf));
                B0 = Bf + (x - FB) / FB1;
                //
                delta = Math.Abs(Bf - B0);
            }
            //
            decimal sinBf = sin(Bf);
            decimal sinBf2 = sinBf * sinBf;
            decimal W = sqrt(1 - e2 * sinBf2);
            decimal W3 = W * W * W;
            decimal N = a / W;
            decimal t = tan(Bf);
            decimal t2 = t * t;
            decimal t4 = t2 * t2;
            decimal cosBf = cos(Bf);
            decimal cosBf2 = cosBf * cosBf;
            decimal yy = e12 * cosBf2;   //η2
            decimal Mf = a0 / W3;
            //
            decimal y_N = y / N;
            decimal y_N2 = y_N * y_N;
            decimal y_N4 = y_N2 * y_N2;
            //
            //计算经伟度值B,L
            decimal t_B = Bf*p  - (p * t / (2 * Mf) * y * y_N) * (1 - (5 + 3 * t2 + yy - 9 * yy * t2) * y_N2 + (61 + 90 * t2 + 45 * t4) * y_N4 / 360);
            decimal t_L = (p / cosBf) * y_N * (1 - (1 + 2 * t2 + yy) * y_N2 / 6 + (5 + 28 * t2 + 24 * t4 + 6 * yy + 8 * yy * t2) * y_N4 / 120);
            //
            L = t_L + l0;
            //
            B = t_B / 3600;   //转为度
            L = L / 3600;   //转为度
            //--the--end--
        }

        //高斯正算方法 (B,L)=>(x,y)
        public virtual void GetXYFromBL(decimal B, decimal L, ref decimal x, ref decimal y)
        {
            //计算临时值
            decimal e4 = pow(e2, 2); // e2 * e2;
            decimal e6 = pow(e2, 3); //e4 * e2;
            decimal e8 = pow(e2, 4); //e6 * e2;
            decimal e10 = pow(e2, 5); //e8 * e2;
            decimal e_12 = pow(e2, 6); //e10 * e2;
            //
            decimal A1 = 1 + 3 * e2 / 4 + 45 * e4 / 64 + 175 * e6 / 256 + 11025 * e8 / 16384 + 43659 * e10 / 65536 + 693693 * e_12 / 1048576;
            decimal B1 = 3 * e2 / 8 + 15 * e4 / 32 + 525 * e6 / 1024 + 2205 * e8 / 4096 + 72765 * e10 / 131072 + 297297 * e_12 / 524288;
            decimal C1 = 15 * e4 / 256 + 105 * e6 / 1024 + 2205 * e8 / 16384 + 10395 * e10 / 65536 + 1486485 * e_12 / 8388608;
            decimal D1 = 35 * e6 / 3072 + 105 * e8 / 4096 + 10395 * e10 / 262144 + 55055 * e_12 / 1048576;
            decimal E1 = 315 * e8 / 131072 + 3465 * e10 / 524288 + 99099 * e_12 / 8388608;
            decimal F1 = 693 * e10 / 1310720 + 9009 * e_12 / 5242880;
            decimal G1 = 1001 * e_12 / 8388608;
            //
            decimal l0 = this.L0 * 3600;  //中央子午线 度转为秒值 如=105*3600
            decimal LL = L * 3600;                   //转为秒值
            //
            decimal t_B = B * this.pi / 180;     //转为弧度值  b
            decimal t_L = (LL - l0) / p;          //转为秒值    l
            decimal L2 = pow(t_L, 2);// t_L * t_L;
            decimal L4 = pow(t_L, 4);// L2 * L2;
            //
            decimal sinB = sin(t_B);
            decimal sinB2 = sinB * sinB;
            decimal W = sqrt(1 - e2 * sinB2);
            //decimal W3 = pow(W, 3);// W * W * W;
            decimal N = a / W;
            decimal t = tan(t_B);
            decimal t2 = t * t;
            decimal t4 = t2 * t2;
            decimal cosB = cos(t_B);
            decimal cosB2 = cosB * cosB;
            decimal cosB4 = cosB2 * cosB2;
            decimal y2 = e12 * cosB2;   //η2
            decimal y4 = y2 * y2;
            //
            decimal l_p = t_L;  //t_L/p;  //上面t_L已经除了p值,这里就不再除p值
            decimal l2_p2 = L2;   //l2/(p*p);
            decimal l4_p4 = L4;   //l4/(p*p*p*p);
            //
            decimal a0 = a * (1 - e2);
            //计算子午弧长公式xx
            decimal xx = a0 * (A1 * t_B - B1 * sin(2 * t_B) + C1 * sin(4 * t_B) - D1 * sin(6 * t_B) + E1 * sin(8 * t_B) - F1 * sin(10 * t_B) + G1 * sin(12 * t_B));
            //计度平面坐标值x,y
            x = xx + N * t * cosB2 * l2_p2 * (0.5M + (5 - t2 + 9 * y2 + 4 * y4) * cosB2 * l2_p2 / 24 + (61 - 58 * t2 + t4) * cosB4 * l4_p4 / 720);
            y = N * cosB * l_p * (1 + (1 - t2 + y2) * cosB2 * l2_p2 / 6 + (5 - 18 * t2 + t4 + 14 * y2 - 58 * y2 * t2) * cosB4 * l4_p4 / 120);
            //
            if (IsBigNumber == true)  //转为高斯投影是大数投影吗?即Zone 35带数  (小数投影为CM_105E)
            {
                y = y + (this.L0 / (int)this.Strip) * 1000000M;
            }
            y = y + 500000.0M;
            //--the--end--
        }

    }
}

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;


namespace ConsoleApp_CoordConvert2000
{
    /// <summary>
    /// 2000坐标系
    /// vp:hsg
    /// create date:2015-12
    /// </summary>
    public class GSCoordConvertionClass_2000 : ProjectionConversion2, IProjectionConversion2
    {
        public GSCoordConvertionClass_2000()
        {
            this.a = 6378137M;           //2000椭球长半轴
            this.b = 6356752.3141M;      //2000椭球短半轴
            this.f = 1 / 298.257222101M;     //椭球扁率f
            this.e = 0.0818191910428M;       //第一偏心率
            this.e2 = 6.693421622966E-03M;   //椭球第一偏心率平方 e2
            this.e12 = 6.738525414683E-03M;  //椭球第二偏心率平方 e12
            this.c = 6399593.62586M;         //极点子午圈曲率半径 c
            this.p = 206264.806252992M;      //弧度秒=180*3600/pi
            this.pi = 3.1415926535M;
        }
        
    }
    /// <summary>
    /// 西安80坐标
    /// vp:hsg
    /// create date:2015-12
    /// </summary>
    public class GSCoordConvertionClass_Xian80 : ProjectionConversion2, IProjectionConversion2
    {
        public GSCoordConvertionClass_Xian80()
        {
            this.a = 6378140M;           //西安80椭球长半轴
            this.b = 6356755.29M;        //西安80椭球短半轴
            this.f = 1 / 298.257M;     //椭球扁率f
            this.e = 0.081819221455523M;       //第一偏心率
            this.e2 = 6.69438499958795E-03M;   //椭球第一偏心率平方 e2
            this.e12 = 6.73950181947292E-03M;  //椭球第二偏心率平方 e12
            this.c = 6399596.65198801M;         //极点子午圈曲率半径 c
            this.p = 206264.806252992M;      //弧度秒=180*3600/pi
            this.pi = 3.1415926535M;
        }
    }
}

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

namespace ConsoleApp_CoordConvert2000
{
    class Program
    {
        static void Main(string[] args)
        {
            decimal x = 0;
            decimal y = 0;
            decimal B = 0;
            decimal L = 0;

            GSCoordConvertionClass_2000 cc = new GSCoordConvertionClass_2000();
            cc.IsBigNumber = true;
            cc.Strip = EnumStrip.Strip3;
            cc.L0 = 105;
            //
            //---------------------------------
            //反算 OK
            x = 4016159.7706M;
            y = 35358852.807M;
            cc.GetBLFromXY(x, y, ref B, ref L);
            //
            B = Math.Round(B, 8);
            L = Math.Round(L, 8);
            System.Console.WriteLine("X, Y = (" + x + ", " + y + "=>B,L=(" + B + "," + L + ")");
            //
            //正算
            //B = 36.155619734M;
            //L = 105.254854607M;
            cc.GetXYFromBL(B, L, ref x, ref y);
            System.Console.WriteLine("B,L=(" + B + "," + L + ")=>X,Y=(" + x + "," + y + "");
            //
            System.Console.Read();

        }
    }
}


已标记关键词 清除标记
相关推荐
©️2020 CSDN 皮肤主题: 大白 设计师:blogdevteam 返回首页