本设计实现了一个大地线长度计算工具,用户可以输入两个点的经纬度坐标,然后点击计算按钮,程序会根据输入的经纬度坐标计算出这两个点之间的大地线长度,并将结果显示在界面上。如下图所示。
附录完整代码如下:
# 导入所需的库
import sys
import math
from PyQt5.QtWidgets import QApplication, QWidget, QVBoxLayout, QHeaderView, QPushButton, QTableWidget, QTableWidgetItem# 定义一个名为GeodeticLine的类,初始化方法接受一个参数e2(椭球体第二偏心率平方)
class GeodeticLine:def __init__(self, e2):self.e2 = e2# 定义名为bessel_inverse的方法,接受输入数据input_data和一个整型变量geodesy。def bessel_inverse(self, input_data, geodesy):# 创建一个名为output的空字典,用于存储计算结果。output = {}# 调用calculate_geodetic_parameters方法计算输入点的辅助参数sin_u和cos_u。sin_u1, cos_u1 = self.calculate_geodetic_parameters(input_data['B1'])sin_u2, cos_u2 = self.calculate_geodetic_parameters(input_data['B2'])# 计算两点经度差L,并将lamda初始值设为L,初始化delta0为0。L = input_data['L2'] - input_data['L1']lamda = Ldelta0 = 0# 在一个无限循环中(后面会有跳出条件),计算p和q的值,然后计算A12(A至B的方位角)while True:p = cos_u2 * math.sin(lamda)q = cos_u1 * sin_u2 - sin_u1 * cos_u2 * math.cos(lamda)A12 = math.atan2(p, q)# 根据p和q的符号调整A12的值。if p > 0 and q > 0:A12 = abs(A12)elif p > 0 and q < 0:A12 = math.pi - abs(A12)elif p < 0 and q < 0:A12 = abs(A12) + math.pielse:A12 = 2 * math.pi - abs(A12)# 计算sin_sigma,cos_sigma和sigma。sin_sigma = p * math.sin(A12) + q * math.cos(A12)cos_sigma = sin_u1 * sin_u2 + cos_u1 * cos_u2 * math.cos(lamda)sigma = math.atan2(sin_sigma, cos_sigma)# 计算sin_A0和cos_sq_A0。sin_A0 = cos_u1 * math.sin(A12)cos_sq_A0 = 1 - sin_A0 ** 2# 计算x、alpha和beta。x = 2 * sin_u1 * sin_u2 - cos_sq_A0 * cos_sigmaalpha = (33523299 - (28189 - 70 * cos_sq_A0) * cos_sq_A0) * 1e-10beta = (28189 - 94 * cos_sq_A0) * 1e-10# 计算delta并更新lamda的值。delta = (alpha * sigma - beta * x * sin_sigma) * sin_A0lamda = L + delta# 判断delta与delta0之间的差值是否小于1e - 12,如果是,则跳出循环;否则,更新delta0的值。if abs(delta - delta0) < 1e-12:breakdelta0 = delta# 对A21(从点B到点A的方位角)进行计算和调整。A21 = A12 + math.piif A21 < 0:A21 += 2 * math.piif A21 > 2 * math.pi:A21 -= 2 * math.piif A12 < math.pi and A21 < math.pi:A21 += math.piif A12 > math.pi and A21 > math.pi:A21 -= math.pi# 根据geodesy的值计算A、B和C。if geodesy == 1:A = 6356863.020 + (10708.949 - 13.474 * cos_sq_A0) * cos_sq_A0B = 10708.938 - 17.956 * cos_sq_A0C = 4.487elif geodesy == 2:A = 6356755.288 + (10710.341 - 13.534 * cos_sq_A0) * cos_sq_A0B = 10710.342 - 18.046 * cos_sq_A0C = 4.512elif geodesy == 3:A = 6356755.288 + (10710.341 - 13.534 * cos_sq_A0) * cos_sq_A0B = 10710.342 - 18.046 * cos_sq_A0C = 4.512# 计算y和S(两点间的大地线长度)。y = (cos_sq_A0 ** 2 - 2 * x * x) * cos_sigmaS = A * sigma + (B * x + C * y) * sin_sigma# 将计算结果存入output字典并返回。output['S'] = Soutput['A12'] = A12output['A21'] = A21return output# 定义一个名为calculate_geodetic_parameters的辅助方法,接受一个参数B(纬度值),计算并返回sin_u和cos_u。def calculate_geodetic_parameters(self, B):W = math.sqrt(1 - self.e2 * math.pow(math.sin(B), 2))sin_u = math.sin(B) * math.sqrt(1 - self.e2) / Wcos_u = math.cos(B) / Wreturn sin_u, cos_uclass MainWindow(QWidget):def __init__(self):super().__init__()self.initUI()def initUI(self):self.setWindowTitle('大地线长度计算')self.setGeometry(300, 300, 600, 300)layout = QVBoxLayout()self.table = QTableWidget(2, 2)self.table.setHorizontalHeaderLabels(['L(dd.mmss)', 'B(dd.mmss)'])self.table.setVerticalHeaderLabels(['点P1', '点P2'])header = self.table.horizontalHeader()header.setSectionResizeMode(QHeaderView.Stretch)layout.addWidget(self.table)self.calc_button = QPushButton('计算')self.calc_button.clicked.connect(self.calculate_distance)layout.addWidget(self.calc_button)self.result_table = QTableWidget(1, 1)self.result_table.setHorizontalHeaderLabels(['大地线长度S (m)'])self.result_table.horizontalHeader().setStretchLastSection(True)layout.addWidget(self.result_table)self.setLayout(layout)def calculate_distance(self):try:lat1 = float(self.table.item(0, 0).text())lon1 = float(self.table.item(0, 1).text())lat2 = float(self.table.item(1, 0).text())lon2 = float(self.table.item(1, 1).text())input_data = {'e2': 0.00669437999013,'L1': math.radians(lon1),'B1': math.radians(lat1),'L2': math.radians(lon2),'B2': math.radians(lat2)}geodetic_line = GeodeticLine(input_data['e2'])output = geodetic_line.bessel_inverse(input_data, 1)self.result_table.setItem(0, 0, QTableWidgetItem(f"{output['S']:.2f}"))except Exception as e:print(f'Error: {e}')if __name__ == '__main__':app = QApplication(sys.argv)mainWindow = MainWindow()mainWindow.show()sys.exit(app.exec_())