测绘程序设计学习—Day3-大地坐标系(BLH)-空间直角坐标系(XYZ)转换

于 2025-07-09 17:46:11 发布
阅读量921 收藏 6
点赞数 22
CC 4.0 BY-SA版权
文章标签： 学习
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/chatcathat/article/details/149225732

今天做到这，发现这个测绘程序设计就是把用编程把公式表达出来。

首先看题目：这题是把大地坐标系(BLH)转换到空间直角坐标系(XYZ)。

你不懂公式？没关系，不懂公式也可以写！ （看得懂数学符号就行）

 首先明确一哈输入输出，中间就是公式了。

  • 输入：纬度(B)、经度(L)、高程(H)

  • 输出：空间直角坐标(X, Y, Z)

1.大地坐标系(BLH)转换到空间直角坐标系(XYZ)

类成员定义 

这两个参数永远不会变吧，是固定值，为了避免运行时被修，我们用const定义。

 const double a = 6378137;
 const double f = 1 / 298.257222101;

辅助方法：

1.我们注意计算公式里面有三角函数，三角函数需弧度输入。我们需要两个方法使得弧度的度可以互转。

        //Convert Degree to rad.
        private double  Deg2Rad(double deg)//输入角度
        {
            double rad = deg * Math.PI / 180.0;// Math.PI =π
            return rad;//输出弧度
        }
        //Convert rad to Degree.
        private double  Rad2Deg(double rad)
        {
            double deg = rad * 180 / Math.PI;
            return deg;
        }

2.为了方便我们可以把e^{2}计算出来：

  private double  ComPuteE2()
        {
            double e2 = 2 * f - f * f;
            return e2;//调用这个就返回e²
        }

3.计算N ：计算时需要输入B，因为是sinB，所以把度转化为弧度，这就用到上面的方法了。

  private double  computeN(double B)
        {
            double e2 = ComPuteE2();// 计算e²
            double b = Deg2Rad(B);//b在sin里面，转弧度
            double W = Math.Sqrt(1 - e2 * Math.Sin(b) * Math.Sin(b));//公式照抄
            //double b=Math.Sin(Deg2Rad(B));
            //double W =Math.Sqrt(1-e2*b*b);
            double N = a / W;
            return N;
        }

核心方法：

输入BLH，outXYZ。调用之前的方法，表达公式。

public void BLH2XYZ(double B, double L, double H, out double X, out double Y, out double Z)
        {
            //sin里面是弧度制 谁在sin里面谁转
            double b = Deg2Rad(B);
            double l = Deg2Rad(L);

            double e2 = ComPuteE2();
            double N = computeN(B);
            
            X = (N + H) * Math.Cos(b) * Math.Cos(l);
            Y = (N + H) * Math.Cos(b) * Math.Sin(l);
            Z = (N * (1 - e2) + H) * Math.Sin(b);
        }

运行

main里面调用一下方法：

 class Program
    {
        static void Main(string[] args)
        {
            double B = 30.30, L = 114.45, H = 50;
            CoorTrans trans = new CoorTrans();
            Console.WriteLine("原始BLH坐标:");
            Console.WriteLine($"纬度 B: {B}");
            Console.WriteLine($"经度 L: {L}");
            Console.WriteLine($"高程 H: {H}");
            Console.WriteLine();
            double X, Y, Z;
            trans.BLH2XYZ(B, L, H, out X, out Y, out Z);
            Console.WriteLine("转换后的XYZ坐标:");
            Console.WriteLine($"X: {X:F3} 米");
            Console.WriteLine($"Y: {Y:F3} 米");
            Console.WriteLine($"Z: {Z:F3} 米");
            Console.WriteLine();
          
        }
    }

2.空间直角坐标系(XYZ) 转换到大地坐标系(BLH)

在看这到题目：空间直角坐标系(XYZ)转换到大地坐标系(BLH)

计算：（说实话各位，这个L的判断处理我也不太懂，以后懂了在做解释）

  private double computeL(double X,double Y)
        {
            //特殊情况处理
            double L = 0;
            if (Math.Abs(X) < 1e-15 && Y > 0)
                L = 0;
            else if (Math.Abs(X) < 1e-15 && Y < 0)
                L = 180;
            //一般情况处理
            else
            {
                L = Math.Atan(Y / X);
                L = Rad2Deg(L);
                if (Y < 0 && X > 0)
                    L = L + 360;
                if (Y > 0 && X < 0)
                    L = 180 + L;
                if (X < 0 && Y < 0)
                    L = 180 + L;
            }
            return L;
        }
 //B can not be 90 deg
        private double computeH(double X,double Y,double B)
        {
            double b = Deg2Rad(B);
            double N = computeN(B);
            double H = Math.Sqrt(X * X + Y * Y) / Math.Cos(b) - N;
            return H;
        }
        private double ComputeB(double Z,double H,double Y,double B0,double X)
        {
            double N = computeN(B0);
            double up = Z * (N + H);
            double under = Math.Sqrt(X * X + Y * Y) * (N*(1 - ComPuteE2() * ComPuteE2()) + H);
            double B = Math.Atan(up / under);
            B = Rad2Deg(B);
            return B;
        }

方法：

public  void XYZ2BLH(double X,double Y,double Z,out  double B,out double L,out double H)
        {
            L = computeL(X,Y);
            double B0 = 0;
            double detaB = 90;
            //纬度 B 和高程 H 的公式互相依赖（H 的计算需要 B，B 的计算又需要 H），无法直接求解，必须通过迭代逼近。
            do
            {
                H = computeH(X, Y, B0);
                B = ComputeB(Z, H, Y, B0, X);
                detaB = B - B0;
                B0 = B;
            } while (Math.Abs(detaB ) < 1e-10);

        }

3.所有代码

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

namespace Day1_4
{
    class CoorTrans
    {
        const double a = 6378137;
        const double f = 1 / 298.257222101;

        public void BLH2XYZ(double B, double L, double H, out double X, out double Y, out double Z)
        {
            //sin里面是弧度制
            double b = Deg2Rad(B);
            double l = Deg2Rad(L);

            double e2 = ComPuteE2();
            double N = computeN(B);
            X = (N + H) * Math.Cos(b) * Math.Cos(l);
            Y = (N + H) * Math.Cos(b) * Math.Sin(l);
            Z = (N * (1 - e2) + H) * Math.Sin(b);
        }
        public  void XYZ2BLH(double X,double Y,double Z,out  double B,out double L,out double H)
        {
            L = computeL(X,Y);
            double B0 = 0;
            double detaB = 90;
            do
            {
                H = computeH(X, Y, B0);
                B = ComputeB(Z, H, Y, B0, X);
                detaB = B - B0;
                B0 = B;
            } while (Math.Abs(detaB ) < 1e-10);

        }
        private double computeL(double X,double Y)
        {
            //特殊情况处理
            double L = 0;
            if (Math.Abs(X) < 1e-15 && Y > 0)
                L = 0;
            else if (Math.Abs(X) < 1e-15 && Y < 0)
                L = 180;
            //一般情况处理
            else
            {
                L = Math.Atan(Y / X);
                L = Rad2Deg(L);
                if (Y < 0 && X > 0)
                    L = L + 360;
                if (Y > 0 && X < 0)
                    L = 180 + L;
                if (X < 0 && Y < 0)
                    L = 180 + L;
            }
            return L;
        }
        //B can not be 90 deg
        private double computeH(double X,double Y,double B)
        {
            double N = computeN(B);
            double H = Math.Sqrt(X * X + Y * Y) / Math.Cos(Deg2Rad(B)) - N;
            return H;
        }
        private double ComputeB(double Z,double H,double Y,double B0,double X)
        {
            double N = computeN(B0);
            double up = Z * (N + H);
            double under = Math.Sqrt(X * X + Y * Y) * (N*(1 - ComPuteE2() * ComPuteE2()) + H);
            double B = Math.Atan(up / under);
            B = Rad2Deg(B);
            return B;
        }
        private double  ComPuteE2()
        {
            double e2 = 2 * f - f * f;
            return e2;
        }



        private double  computeN(double B)
        {
            double e2 = ComPuteE2();
            double b = Deg2Rad(B);
            double W = Math.Sqrt(1 - e2 * Math.Sin(b) * Math.Sin(b));
            //double b=Math.Sin(Deg2Rad(B));
            //double W =Math.Sqrt(1-e2*b*b);
            double N = a / W;
            return N;
        }
        
        
        //Convert Degree to rad.
        private double  Deg2Rad(double deg)
        {
            double rad = deg * Math.PI / 180.0;
            return rad;
        }
        //Convert rad to Degree.
        private double  Rad2Deg(double rad)
        {
            double deg = rad * 180 / Math.PI;
            return deg;
        }

    }
}

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

namespace Day1_4
{
    class Program
    {
        static void Main(string[] args)
        {
            double B = 30.30, L = 114.45, H = 50;
            CoorTrans trans = new CoorTrans();
            Console.WriteLine("原始BLH坐标:");
            Console.WriteLine($"纬度 B: {B}");
            Console.WriteLine($"经度 L: {L}");
            Console.WriteLine($"高程 H: {H}");
            Console.WriteLine();
            double X, Y, Z;
            trans.BLH2XYZ(B, L, H, out X, out Y, out Z);
            Console.WriteLine("转换后的XYZ坐标:");
            Console.WriteLine($"X: {X:F3} 米");
            Console.WriteLine($"Y: {Y:F3} 米");
            Console.WriteLine($"Z: {Z:F3} 米");
            Console.WriteLine();
            // XYZ转BLH
            trans.XYZ2BLH(X, Y, Z, out B, out L, out H);

            Console.WriteLine("转换回的BLH坐标:");
            Console.WriteLine($"纬度 B: {B:F9}°");
            Console.WriteLine($"经度 L: {L:F9}°");
            Console.WriteLine($"高程 H: {H:F3}米");

            Console.ReadKey();
        }
    }
}

注意 XYZ转BLH中H算出来有错误，我也不知道怎么修改，整体代码有待优化。

chatcathat
关注
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

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

抵扣说明：

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

余额充值