首页 星云 工具 资源 星选 资讯 热门工具
:

PDF转图片 完全免费 小红书视频下载 无水印 抖音视频下载 无水印 数字星空

使用广播星历计算卫星坐标(Python)

编程知识
2024年08月31日 15:06

前言

  1. 本代码为GNSS课程设计代码,仅供参考,使用的计算方法与公式均来源于王坚主编的《卫星定位原理与应用(第二版)》。
  2. 本代码计算结果可以通过下载精密星历进行比照,误差在1-10m左右。
  3. 实现功能:读取卫星广播星历,并将其计算为WGS-84坐标系下的坐标,每颗卫星,每15分钟输出一次。

广播星历下载

  1. 有多重方法进行下载,由于网络原因以及使用的便捷程度,建议使用武汉大学的IGS网站进行下载。(http://www.igs.gnsswhu.cn/index.php)
  2. 文件名不需要填写,在选择时间时注意,即使只选择一天的数据,设置结束时间也要到第二天,否则会显示错误,检索结果中的数字代表该日期是所选年份的第几天。

Python函数库

本代码使用numpy,pandas,gnss_lib_py,matplotlib四个函数库,请提前安装。

安装代码

#两条命令根据使用环境进行选择
pip install gnss-lib-py pandas matplotlib #Python环境安装代码
conda install gnss-lib-py pandas matplotlib -c conda-forge #conda环境安装代码

gnss_lib_py库

GitHub主页:https://github.com/Stanford-NavLab/gnss_lib_py?tab=readme-ov-file
文档主页:https://gnss-lib-py.readthedocs.io/en/latest/index.html
本文主要使用该库的读取以及转化为DataFrame功能,其中参数的命名规则以及时间转换规则可以在文档中找到。

具体代码

建议使用jupyter进行执行,下方代码为分单元块格式,如果没有jupyter环境可以直接粘贴到一个python文件进行运行。

代码块一:

import gnss_lib_py as glp
import datetime
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

代码块二:

# 导入23n文件
file_path = 'brdc2550.23n'
data = glp.RinexNav(file_path)
data_df = data.pandas_df()

代码块三:

# 寻找最小差值的参考时刻
def find(inweekmilli, refers, insv):
    filter_sv = refers[refers['gnss_sv_id'] == insv]
    defference = np.abs(inweekmilli - filter_sv['t_oe'])
    return defference.idxmin()


times = np.array([None] * 24 * 4)
gpsmillis = np.array([None] * 24 * 4)

n = 0
for hour in range(0, 24):
    minut = 0
    while minut < 60:
        times[n] = datetime.datetime(2023, 9, 12, hour, minut, 0, tzinfo=datetime.timezone.utc)
        gpsmillis[n] = glp.datetime_to_gps_millis(times[n])
        minut += 15
        n += 1
gpsmillis = np.array(gpsmillis)

GM=3.986005E+14
sqrtGM = np.sqrt(GM)
sv_list = [f'G{str(i).zfill(2)}' for i in range(1, 33)]
outdata = pd.DataFrame(columns=['data', 'gnss_sv_id', 'X', 'Y', 'Z'], index=range(24 * 4 * 32))
orbit = pd.DataFrame(columns=['data', 'gnss_sv_id', 'x', 'y'], index=range(24 * 4 * 32))
m = 0
j = 0
print("正在计算,请稍候。")
for gpsmilli in gpsmillis:
    for sv in sv_list:
        week, milli_week = glp.gps_millis_to_tow(gpsmilli)
        milli_week=milli_week-18
        index = find(milli_week, data_df, sv)
        print(f'sv:{sv},time:{times[j]},index:{index}')
        print(milli_week,data_df.iloc[index]['t_oe'])
        a = np.power(data_df.iloc[index]['sqrtA'], 2)
        n0 = sqrtGM / np.power(a, 3 / 2)
        n = n0 + data_df.iloc[index]['deltaN']
        tk = milli_week - data_df.iloc[index]['t_oe']
        M = data_df.iloc[index]['M_0'] + n * tk
        e = data_df.iloc[index]['e']
        # 打印中间结果M和e
        print(f"M: {M}, e: {e}")

        # 解开普勒方程
        E = M
        for _ in range(50):  # 使用迭代方法求解E
            E = M + e * np.sin(E)
        
        # 打印中间结果E
        print(f"E: {E}")
        f = np.arctan((np.sqrt(1 - e**2) * np.sin(E)) / (np.cos(E) - e))
        if E > np.pi*0.5:
            f=f+np.pi
        if E < -np.pi*0.5:
            f=f-np.pi
        if np.pi*0.5 > E > 0 > f:
            f=f+np.pi
        if -np.pi*0.5 < E < 0 < f:
            f=f-np.pi
        print(f"arctan({(np.sqrt(1 - e**2) * np.sin(E)) / (np.cos(E) - e)}),f:{f}")
        u_pie = data_df.iloc[index]['omega'] + f
        r_pie = a * (1 - e * np.cos(E))
        C_uc = data_df.iloc[index]['C_uc']
        C_us = data_df.iloc[index]['C_us']
        C_rc = data_df.iloc[index]['C_rc']
        C_rs = data_df.iloc[index]['C_rs']
        C_ic = data_df.iloc[index]['C_ic']
        C_is = data_df.iloc[index]['C_is']
        delta_u = C_uc * np.cos(2 * u_pie) + C_us * np.sin(2 * u_pie)
        delta_r = C_rc * np.cos(2 * u_pie) + C_rs * np.sin(2 * u_pie)
        delta_i = C_ic * np.cos(2 * u_pie) + C_is * np.sin(2 * u_pie)
        u = u_pie + delta_u
        r = r_pie + delta_r
        i = data_df.iloc[index]['i_0'] + delta_i + data_df.iloc[index]['IDOT'] * tk
        print(f'u:{u}')
        x = r * np.cos(u)
        y = r * np.sin(u)
        w_e = 7.292115E-5
        L = data_df.iloc[index]['Omega_0'] + (data_df.iloc[index]['OmegaDot']- w_e )* milli_week - data_df.iloc[index]['OmegaDot']*data_df.iloc[index]['t_oe']
        X = x * np.cos(L) - y * np.cos(i) * np.sin(L)
        Y = x * np.sin(L) + y * np.cos(i) * np.cos(L)
        Z = y * np.sin(i)
        orbit.iloc[m,:] = [times[j],sv,x,y]
        outdata.iloc[m, :] = [times[j], sv, X, Y, Z]
        m += 1
    j += 1
print("由于结果较长,请到Excel中查看,文件位于代码同级目录下outdata.csv。")
outdata.to_csv('outdata.csv')
print("导出成功。")

代码块四:

# 三维坐标可视化显示
out_sv = 'G20'
fig = plt.figure()
ax = plt.axes(projection='3d')
X = outdata[outdata['gnss_sv_id'] == out_sv]['X']
Y = outdata[outdata['gnss_sv_id'] == out_sv]['Y']
Z = outdata[outdata['gnss_sv_id'] == out_sv]['Z']
ax.plot(X, Y, Z, label=out_sv)
ax.legend()
plt.show()
From:https://www.cnblogs.com/BIGWangqz/p/18390417
本文地址: http://www.shuzixingkong.net/article/1613
0评论
提交 加载更多评论
其他文章 使用 Quickwit 的搜索流功能为 ClickHouse 添加全文搜索
本指南将帮助您使用 Quickwit 的搜索流功能为知名的 OLAP 数据库 ClickHouse 添加全文搜索。Quickwit 暴露了一个 REST 端点,可以极快地(每秒最多 5000 万条)流式传输匹配搜索查询的 ID 或其他属性,ClickHouse 可以轻松地使用它们进行连接查询。 我们
使用 Quickwit 的搜索流功能为 ClickHouse 添加全文搜索
探索一下 Enum 优化
探索一下 Enum 优化 SV.Enums主要是探索如何让 enum 更高效 其中涉及的优化手段并非完全自创 很多内容参考于以下项目 NetEscapades.EnumGenerators FastEnum runtime 主要优化手段 其实主要全是 空间换时间,大量缓存 封装入口方法以及 sour
使用 nuxi build-module 命令构建 Nuxt 模块
title: 使用 nuxi build-module 命令构建 Nuxt 模块 date: 2024/8/31 updated: 2024/8/31 author: cmdragon excerpt: nuxi build-module 命令是构建 Nuxt 模块的核心工具,它将你的模块打包成适合
使用 nuxi build-module 命令构建 Nuxt 模块 使用 nuxi build-module 命令构建 Nuxt 模块
手把手在STM32F103C8T6上构建可扩展可移植的DHT11驱动
前言 如何驱动一个你陌生的传感器呢?别看我,也别在网上死马当活马医!你需要做的,首先是明确你的传感器的名称,在这里,我们想要使用的是DHT11温湿度传感器 可能需要的前置知识 简单的OLED驱动原理 简单的IIC通信知识 基础的查手册能力 相对稳固的C语言基础 不会没关系,我会详细说明的! 一种可能
手把手在STM32F103C8T6上构建可扩展可移植的DHT11驱动 手把手在STM32F103C8T6上构建可扩展可移植的DHT11驱动 手把手在STM32F103C8T6上构建可扩展可移植的DHT11驱动
设计模式之迭代器模式
迭代器模式很多人都熟悉,但是什么是迭代器,为什么要用迭代器?这个很多人就很难做出具体回答了,只是知道如果有了迭代器,那么我们就能foreach遍历了,方便循环处理。这只是对迭代器的用途,进行了回答,foreach语法是java1.5时加入的语法糖,那么在这之前呢,之前是怎么做的?要知道并不是所有容器
设计模式之迭代器模式
2024 NepCTF
NepCTF NepMagic —— CheckIn 直接玩游戏就能出 注意有一关要把隐藏的方块全找到 NepCamera 先使用tshark读取数据 结果文件中发现大量jpeg头ffd8ffe0。 猜测是多个图片,编写脚本,提取出来。 脚本: import re inputFilePath=&qu
2024 NepCTF 2024 NepCTF 2024 NepCTF
Go plan9 汇编: 打通应用到底层的任督二脉
原创文章,欢迎转载,转载请注明出处,谢谢。 0. 前言 作为一个严肃的 Gopher,了解汇编是必须的。本汇编系列文章会围绕基本的 Go 程序介绍汇编的基础知识。 1. Go 程序到汇编 首先看一个简单到令人发指的示例: package main func main() { a := 1 print
Go plan9 汇编: 打通应用到底层的任督二脉 Go plan9 汇编: 打通应用到底层的任督二脉
[postgres]使用pgbench进行基准测试
前言 pgbench是一种在postgres上进行基准测试的简单程序,一般安装后就会自带。pgbench可以再并发的数据库绘画中一遍遍地进行相同序列的SQL语句,并且计算平均事务率。 测试准备 既然要测postgres,肯定要先有个postgres。安装过程略过。 一些环境信息: postgres版