坐标转换（空间直角坐标系与大地坐标系）

最新推荐文章于 2024-01-03 08:41:29 发布

MqtGhj

最新推荐文章于 2024-01-03 08:41:29 发布

阅读量1.8k

点赞数 6

分类专栏：测绘专业相关编程文章标签： python pycharm 职场和发展经验分享学习

本文链接：https://blog.csdn.net/weixin_64338372/article/details/127650193

版权

测绘专业相关编程专栏收录该内容

17 篇文章 26 订阅

订阅专栏

代码中有详细的注释

根据《大地测量学（第三版）》武汉大学出版社的教材进行编写

import math as m
import sys as s
import os


class Earth():
    # 建立椭球类，采用CGCS2000中国大地坐标系数据
    # 《大地测量学基础》P.106的表4-1
    def __init__(self):  # 给定椭球数据
        # 长半轴
        self.a = 6378137.0
        # 短半轴
        self.b = 6356752.31414
        # 扁率
        self.f = 1 / 298.257222101
        # 第一偏心率的平方
        self.e1s = m.pow(0.0818191910428, 2)
        # 第二偏心率的平方
        self.e2s = self.e1s / (1 - self.e1s)


class Transform():
    # 建立转换类
    def __init__(self):  # 导入椭球数据
        self.earth = Earth()


    def BLH_XYZ(self, B, L, H):  # 大地坐标转换为空间直角坐标
        # 将角度转换为弧度
        self.b = m.radians(B)
        self.l = m.radians(L)
        # 这里用个python的包就很简便的实现了角度到弧度间的转换

        # 计算辅助函数
        self.W = m.sqrt(1 - self.earth.e1s * m.pow(m.sin(self.b), 2))
        #《大地测量学基础》P.106的（4-5）

        # 转换为空间直角坐标
        # 《大地测量学基础》P.111的（4-30）
        self.N = self.earth.a / self.W
        self.X = (self.N + H) * m.cos(self.b) * m.cos(self.l)
        self.Y = (self.N + H) * m.cos(self.b) * m.sin(self.l)
        self.Z = (self.N * (1 - self.earth.e1s) + H) * m.sin(self.b)

    def XYZ_BLH(self, X, Y, Z):  # 空间直角坐标转换为大地坐标
        # 求出大地经度
        self.l = m.atan(Y / X)  #《大地测量学基础》P.111的（4-31）

        # 求出大地纬度
        # 《大地测量学基础》P.112
        self.r = m.sqrt(X * X + Y * Y)
        self.tb1 = Z / self.r    #《大地测量学基础》P.112的（4-32）中的sinB的B取0得到  self.e1s是第一偏心率的平方
        while True:
            '''
            self.tb2 = 1 / self.r * (
                    Z + self.earth.a * self.earth.e1s * self.tb1 /
                    m.sqrt(1 + self.tb1 * self.tb1 * (1 - self.earth.e1s)))  #《大地测量学基础》P.112的（4-38）的变式
            '''
# 这里的俩个方法的tb2的结果是一样，本质是公式的转变，下面这个是课本的
            self.tb2 = 1 / self.r * (
                    Z + self.earth.a * self.earth.a * self.earth.e1s * self.tb1 /
                    (self.earth.b * m.sqrt(1 + self.earth.e2s + self.tb1 * self.tb1)))# 《大地测量学基础》P.112的（4-38）

            if abs(self.tb2 - self.tb1) <= 5e-10:  # 对于这个范围特地查找了下资料
                break
            self.tb1 = self.tb2
        self.b = m.atan(self.tb2)   # 利用迭代的方法求出B的值，直到符合要求时返回

        # 求出大地高
        self.W = m.sqrt(1 - self.earth.e1s * m.pow(m.sin(self.b), 2))
        self.N = self.earth.a / self.W
        self.H = self.r / m.cos(self.b) - self.N  #《大地测量学基础》P.112的（4-35）


        # 将弧度转换为角度
        self.B = m.degrees(self.b)
        self.L = m.degrees(self.l)
        # 这里用个python的包就很简便的实现了弧度到角度间的转换


class Point():
    # 建立点类
    def BLH(self, B, L, H):  # 建立大地坐标点类
        self.B = B
        self.L = L
        self.H = H
        self.P = Transform()
        self.P.BLH_XYZ(B, L, H)
        self.X = self.P.X
        self.Y = self.P.Y
        self.Z = self.P.Z

    def XYZ(self, X, Y, Z):  # 建立直角坐标点类
        self.X = X
        self.Y = Y
        self.Z = Z
        self.P = Transform()
        self.P.XYZ_BLH(X, Y, Z)
        self.B = self.P.B
        self.L = self.P.L
        self.H = self.P.H

SolveType = input('请输入解算类型：（键入A为大地坐标->直角坐标，B为直角坐标->大地坐标）\n')
while True:
    if SolveType == 'A':
        # 当选择'大地坐标->直角坐标'模式时
        # 输入大地坐标
        I = input('请输入大地坐标B、L、H，中间用空格隔开：\n')
        while True:
            if I:
                pass
            else:
                s.exit(0)
            try:
                [B, L, H] = [float(n) for n in I.split()]
                break
            except:
                I = input('输入错误，请重新输入大地坐标B、L、H，中间用空格隔开：\n')
            else:
                pass

        # 开始转换
        P = Point()  # 创建一个实例对象
        P.BLH(B, L, H)
        B = str(B)
        L = str(L)
        H = str(H)
        X = str(P.X)
        Y = str(P.Y)
        Z = str(P.Z)
        massege = '\n您输入的大地坐标为：（' + B + ', ' + L + ', ' + H + '）\n转换得到直角坐标为：（' + X + ', ' + Y + ', ' + Z + '）'
        print(massege)
        break

    elif SolveType == 'B':
        # 当选择'直角坐标->大地坐标'模式时
        # 输入直角坐标
        I = input('请输入直角坐标X、Y、Z，中间用空格隔开：\n')
        while True:
            if I:
                pass
            else:
                s.exit(0)
            try:
                [X, Y, Z] = [float(n) for n in I.split()]
                1/X
                break
            except:
                I = input('输入错误，请重新输入直角坐标X、Y、Z，中间用空格隔开：\n')
            else:
                pass

        # 开始转换
        P = Point()   # 创建一个实例对象
        P.XYZ(X, Y, Z)
        X = str(X)
        Y = str(Y)
        Z = str(Z)
        B = str(P.B)
        L = str(P.L)
        H = str(P.H)
        massege = '\n您输入的直角坐标为：（' + X + ', ' + Y + ', ' + Z + '）\n转换得到大地坐标为：（' + B + ', ' + L + ', ' + H + '）'
        print(massege)
        break

    else:
        SolveType = input('输入错误，请重新输入解算类型：（键入C为正算，D为反算）\n')
os.system('pause')  # 暂停命令

运行结果示例：