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

最新推荐文章于 2024-08-08 22:15:04 发布
最新推荐文章于 2024-08-08 22:15:04 发布
阅读量1w 收藏 50
点赞数 2
分类专栏： 国土资源 算法
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/hsg77/article/details/50200585
版权

基于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();

        }
    }
}


hsg77
关注
专栏目录
大地坐标与经纬度坐标转换
03-20
大地坐标与经纬度坐标转换
高斯正反算（C#实现）
10-25
该程序提供了在不同坐标系统下进行高斯正算和反算的功能，用C#编码，代码简洁，界面漂亮。
矢量文件坐标转换：2000坐标系转换为wgs84坐标系，具体代码实现
最新发布
wavemenu的博客
08-08 693
在上述代码的主程序部分，我们指定了一个Shapefile文件路径和保存目录，并调用了 VectorTranslate 函数将其转换为ESRI Shapefile格式，并重新投影到WGS84坐标系。这些步骤共同构成了一个完整的矢量数据转换流程，确保了数据的准确性和实用性。最近在处理矢量样本的时候，遇到一些shp文件的坐标系为2000坐标，需要统一地把非WGS84坐标系的矢量转换为WGS84坐标系。首先，在转换之前，通过检查并创建目标保存目录，保证了正确的目录结构，这对于避免转换失败或产生错误输出至关重要。
使用c#设计winform窗体实现高斯坐标正反算
Cyy2511250589的博客
06-09 1402
但是按照他们的思路正算x坐标是对的，y是错的。害，计算不出来正确的，搞得花里胡哨的希望搞个高分。注意：反算算不出来啊！主要就是公式，这是老师发的。
最新坐标转换含（2000坐标系转换）
04-16
最新2000坐标系转换，转换功能强大，坐标系转换，54坐标系、80坐标系、2000坐标系,相互转换，免费
利用C#编写一个高斯正反算程序
Zj1638的博客
07-11 6422
利用C#编写一个高斯正反算程序
GPS经纬度转工程坐标，包含七参数转换，c#实现代码，超详细
qq_41816368的博客
08-12 3332
完整代码 using System; using System.Collections.Generic; using System.Linq; using System.Text; using MathNet.Numerics.LinearAlgebra; using MathNet.Numerics.LinearAlgebra.Double; namespace CoordTrans { public class MyPoint { public double L =
C#坐标转换（大地坐标与空间坐标互换、高斯正反算）.zip
05-28
本文将详细介绍C#编程语言中关于大地坐标与空间直角坐标系之间的转换，以及如何进行高斯正反算。这些知识对于理解和处理地理定位数据至关重要。 首先，我们要理解两个核心概念：大地坐标和空间直角坐标。 1. **...
C# -大地测量学-高斯投影正反算公式程序
06-20
本项目是一个用C#编程语言实现的高斯投影正反算公式程序，可以帮助我们理解和应用这个关键的大地测量技术。 首先，我们要理解高斯投影的基本原理。高斯投影将地球表面的经纬度转换为平面坐标，通过一个椭圆柱面与...
C#实现高斯正反算
09-13
这个“C#实现高斯正反算”项目显然是一个利用C#编程语言来实现高斯投影的计算过程，包括正算和反算两个部分。下面我们将深入探讨这两个概念以及它们在实际应用中的重要性。 **高斯投影正算**是将地理坐标（经度L和...
C# 高斯正反算，坐标转换源代码，（WGS84，西安80，北京54）
02-26
这里我们将详细探讨C#语言实现高斯正反算以及坐标转换的相关知识点。 1. **高斯投影**：高斯-克吕格投影是地图投影的一种，由德国数学家弗里德里希·高斯提出，用于将球面地理坐标转换为平面坐标。这种投影方式保持...
利用C#编写一个大地坐标与空间直角坐标相互转换程序
06-21
利用c#编写的一个window窗体应用，可以实现大地坐标与空间直角坐标相互转换，转换结果利用https://blog.csdn.net/starmings/article/details/123034189的程序进行验证，结果一致。转换公式教材和csdn上都有，但是需要主要角度与弧度的转换。
C#经纬度坐标与大地坐标相互转换工具
08-09
经纬度转换工具，对手持设备获取的经纬度坐标转换为真实的大地坐标。
高斯投影正反算及邻带换算
03-26
话不多说，直接上代码 using System; using System.Collections.Generic; using System.ComponentModel; using System.Data; using System.Drawing; using System.Linq; using System.Text; using System.Threading.Tasks; using System.Windows.Forms; namespace _高斯投影 { public partial class Form2 : Form { public Form2() { InitializeComponent(); } double DD2RAD(double n) { double DD; double MM; double SS; DD = Math.Floor(n); MM = Math.Floor((n - DD) * 100); SS = ((n - DD) * 100 - MM) * 100; n = (DD + MM / 60.0 + SS / 3600.0) * Math.PI / 180.0; return n; } private void Form2_Load(object sender, EventArgs e) { } private void button1_Click(object sender, EventArgs e) { double B, L; B = double.Parse(textBox1.Text); B = DD2RAD(B); L = double.Parse(textBox2.Text); L = DD2RAD(L); double L0 = double.Parse(textBox3.Text); L0 = DD2RAD(L0); double a = double.Parse(textBoxa.Text); double e2 = double.Parse(textBoxe2.Text); double A = 1.0 + 3.0 * e2 / 4 + 45.0 * e2 * e2 / 64 + 175.0 * Math.Pow(e2, 3) / 256 + 11025.0 * Math.Pow(e2, 4) / 16384 + 43659.0 * Math.Pow(e2, 5) / 65536; double B0 = 3.0 * e2 / 4 + 15.0 * e2 * e2 / 16 + 525.0 * Math.Pow(e2, 3) / 512 + 2205.0 * Math.Pow(e2, 4) / 2048 + 72765.0 * Math.Pow(e2, 5) / 65536; double C = 15.0 * e2 * e2 / 64 + 105.0 * Math.Pow(e2, 3) / 256 + 2205.0 * Math.Pow(e2, 4) / 4096 + 10395.0 * Math.Pow(e2, 5) / 16384; double D = 35.0 * Math.Pow(e2, 3) / 512 + 315.0 * Math.Pow(e2, 4) / 2048 + 31185.0 * Math.Pow(e2, 5) / 131072; double α = A * a * (1 - e2);//α double β = -B0 * a * (1 - e2) / 2.0;//β double γ = C * a * (1 - e2) / 4.0; double σ = -D * a * (1 - e2) / 6.0; double X
大地经纬度坐标转换
04-06
转换大地经纬度坐标使用。可用于西安80北京54和经纬度坐标之间的转换
工程坐标转换方法C#代码实现
Jimmysaid的博客
12-04 1529
在前面的文章中系统的阐述了工程坐标的转换类别和转换的方法。关于转换代码实现，有很多的类库。这里针对GPS接收的WGS84椭球的经纬度转换为地方坐标系的问题，利用C#，对工程坐标转换方法和步骤做出详细的解答。不基于任何类库和函数库，也未使用矩阵库，可以便利的将代码移植到任何语言。
Autolisp高斯反算，CGCS2000大地坐标系平面坐标转换为经纬度坐标
泥鳅也是鱼的博客
05-10 6440
CAD制图中，经常要用到根据图上某点坐标标注经纬度的功能，之前在网上找了一个西安80坐标系下标注的，现在用2000的比较多，网上没有找到合适的，所以自己照着三调的培训资料写了一下，三调的资料内有高斯投影的计算公式以及图幅分幅和随球面积等计算内容。网上有很多.net下写的2000坐标系下高斯反算的代码，但是我没看懂，没理解需要迭代循环的一部分，如果有懂的，可以给科普一下。谢谢。 五个常数自行翻阅三调培训资料，我不确定是不是涉密。 (defun xytobl2000 (y x cm / a bb f e2 e1
android 获取经纬度_XOMAP-爬虫[4]-零代码实现将"地址"转换为“经纬度”坐标及“结构化地址”...
weixin_39790686的博客
10-30 131
小Ｏ地图提供将地址转换为经纬度坐标的功能，例如：将 “北京市西城区金融街” 转换为　“116.359861, 39.917225”　经纬度坐标，并将地址分解为，北京市｜西城区｜金融街等结构。文末有小O地图介绍及下载说明。【１】新建任务​【２】导入数据导入要处理的地址数据，支持导入CSV和EXCEL格式的数据源。​选择导入数据文件后，然后进行导入设置，如下图，完成后单击确定导入，提示成功后，关闭设置...
C#实现的高斯坐标转换工具
本文将详细介绍如何使用C#语言实现高斯投影坐标转换，包括正算（地理坐标到平面坐标）和反算（平面坐标到地理坐标）的过程。在地理信息系统（GIS）中，高斯投影是一种广泛使用的地图投影方法，尤其在中国的国家统一...
评论 7
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值