Location>code7788 >text

Calculation of satellite coordinates using broadcast ephemeris (Python)

Popularity：877 ℃/2024-08-31 16:11:17

preamble

  1. This code is a GNSS course design code for reference only. The calculation methods and formulas used are from "Satellite Positioning Principles and Applications (Second Edition)" edited by Wang Jian.
  2. The results of this code calculation can be compared by downloading a precision ephemeris with an error of about 1-10m.
  3. Function to be realized: Read the satellite broadcast ephemeris and calculate it into coordinates under the WGS-84 coordinate system, per satellite, and output it every 15 minutes.

Broadcast Ephemeris Download

  1. There are multiple ways to download, and due to network reasons and ease of use, it is recommended to use Wuhan University's IGS website to download. (/）
  2. The file name does not need to be filled in. When selecting the time, pay attention to the fact that even if only one day's data is selected, the end time of the setup should be up to the next day, otherwise it will be displayed incorrectly, and the number in the retrieval result represents that the date is the first day of the selected year.

Python Function Library

This code uses numpy, pandas, gnss_lib_py,matplotlib four libraries, please install in advance.

Installation Code

#The two commands are chosen according to the environment in which they will be used
pip install gnss-lib-py pandas matplotlib #Python environment install code
conda install gnss-lib-py pandas matplotlib -c conda-forge #conda environment installation code

gnss_lib_py library

GitHub homepage:/Stanford-NavLab/gnss_lib_py?tab=readme-ov-file
Documentation home page:/en/latest/
This article mainly uses the library's read and transformed into a DataFrame function, which parameter naming rules and time conversion rules can be found in the documentation.

specific code

It is recommended to use jupyter for execution, the code below is in sub-cell block format, if you do not have a jupyter environment can be directly pasted into a python file to run.

Code block one:

import gnss_lib_py as glp
import datetime
import pandas as pd
import numpy as np
import  as plt

Code block two:

# Import the 23n file
file_path = 'brdc2550.23n'
data = (file_path)
data_df = data.pandas_df()

Code block three:

# Finding the reference moment for the minimum difference
def find(inweekmilli, refers, insv):
    filter_sv = refers[refers['gnss_sv_id'] == insv]
    defference = (inweekmilli - filter_sv['t_oe'])
    return ()


times = ([none] * 24 * 4)
gpsmillis = ([None] * 24 * 4)

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

GM=3.986005E+14
sqrtGM = (GM)
sv_list = [f'G{str(i).zfill(2)}' for i in range(1, 33)]
outdata = (columns=['data', 'gnss_sv_id', 'X', 'Y', 'Z'], index=range(24 * 4 * 32))
orbit = (columns=['data', 'gnss_sv_id', 'x', 'y'], index=range(24 * 4 * 32))
m = 0
j = 0
print("computation in progress，Please wait.。")
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 = (data_df.iloc[index]['sqrtA'], 2)
        n0 = sqrtGM / (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']
        # Print intermediate resultsMcap (a poem)e
        print(f "M: {M}, e: {e}")

        # solve Kepler's equations
        E = M
        for _ in range(50):  # Solve using iterative methodsE
            E = M + e * (E)
        
        # Print intermediate resultsE
        print(f "E: {E}")
        f = (((1 - e**2) * (E)) / ((E) - e))
        if E > *0.5:
            f=f+
        if E < -*0.5:
            f=
        if *0.5 > E > 0 > f:
            f=f+
        if -*0.5 < E < 0 < f:
            f=
        print(f "arctan({((1 - e**2) * (E)) / ((E) - e)}),f:{f}")
        u_pie = data_df.iloc[index]['omega'] + f
        r_pie = a * (1 - e * (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 * (2 * u_pie) + C_us * (2 * u_pie)
        delta_r = C_rc * (2 * u_pie) + C_rs * (2 * u_pie)
        delta_i = C_ic * (2 * u_pie) + C_is * (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 * (u)
        y = r * (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 * (L) - y * (i) * (L)
        Y = x * (L) + y * (i) * (L)
        Z = y * (i)
        [m, :] = [times[j],sv,x,y]
        [m, :] = [times[j], sv, X, Y, Z]
        m += 1
    j += 1
print("Due to longer results，please go toExcellook sth. up，The file is located in the same level directory as the code。")
outdata.to_csv('')
print("Export Successful。")

Code block four:

# Visualization of 3D coordinates
out_sv = 'G20'
fig = ()
ax = (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']
(X, Y, Z, label=out_sv)
()
()
<Previous Next>