5
$\begingroup$

Im trying to get the historical Lon/Lat/Alt of a starlink satellite. I'm using TLE's from space-track and Skyfield to calculate latitude, longitude and altitude of satellites for a given date range. The methodology I'm using:

  1. Get TLE for date range 06/01 - 06/07 (say you get 30 data points). EDIT: got 19 historical TLE's.
  2. Use Skyfield to generate the LLA (lat,lon,alt) for those 30 data points + timestamps between those 30 TLE's.
  3. Plot the altitude as predicted using Skyfield (that uses SGP4 under the hood?) over time.

Here's what I got for STARLINK-55775 (y-axis in kilometres).

satellite bouncing

I'm not familiar with astrodynamics, just trying to learn something new.

--EDITS--

Code:

from sgp4.api import Satrec
from sgp4.api import SGP4_ERRORS
from skyfield.api import EarthSatellite
import requests
from skyfield.api import load, wgs84
import pandas as pd
import datetime
from dateutil import parser, tz
import pytz
import plotly.express as px

def create_session():
    user_name = 'username'
    password = 'password'
    base_url = 'https://www.space-track.org/'
    request_login = "/ajaxauth/login"
    request_cmd_action = "/basicspacedata/query" 
    request_find_starlinks = "/class/tle_latest/NORAD_CAT_ID/53871/ORDINAL/1/OBJECT_NAME/STARLINK~~/format/json/orderby/NORAD_CAT_ID%20asc"
    session = requests.Session() 
    # need to log in first. note that we get a 200 to say the web site got the data, not that we are logged in
    resp = session.post(base_url + request_login, data={'identity': user_name, 'password': password})
    if resp.status_code != 200:
        raise MyError(resp, "POST fail on login")

    # this query picks up all Starlink satellites from the catalog. Note - a 401 failure shows you have bad credentials 
    resp = session.get(base_url + request_cmd_action + request_find_starlinks)
    if resp.status_code != 200:
        print(resp)
        raise MyError(resp, "GET fail on request for Starlink satellites")
    else:
        print("########### SUCCESSFULLY SET UP SESSION ###############")
        return session

def get_tle(start_dt, end_dt, session, id=53871, ):
    url = f'https://www.space-track.org/basicspacedata/query/class/gp_history/NORAD_CAT_ID/{id}/orderby/EPOCH ASC/EPOCH/{start_dt}--{end_dt}/format/json'
    result = session.get(url) #logs into space-track to get data.
    if result.status_code == 200: 
        docs = result.json()
        return_doc = []
        for doc in docs:
            _t = {}
            _t['timestamp'] = doc['EPOCH']
            _t['TLE_LINE1'] = doc['TLE_LINE1']
            _t['TLE_LINE2'] = doc['TLE_LINE2']
            return_doc.append(_t)
        return return_doc
    else:
        print(result.text)

def get_satellite_lla(tle_1, tle_2, timestamp, **kwargs):
    ts = load.timescale()
    if type(timestamp) == str:
        t = ts.utc(parser.parse(timestamp))
    else:
        t = ts.utc(timestamp)
    satellite = EarthSatellite(tle_1, tle_2, 'StarLink')
    geocentric = satellite.at(t)
    lat, lon = wgs84.latlon_of(geocentric)
    alt = wgs84.height_of(geocentric).km
    session.close()
    return pd.Series([lat.degrees, lon.degrees, alt], index=['lat', 'lon', 'alt'])

session = create_session()
id = 56529
tle_data = get_tle(id=id,start_dt='2023-05-31', end_dt='2023-06-07', session=session)
df = pd.DataFrame.from_dict(tle_data)
df.drop_duplicates(subset='timestamp', inplace=True)
resampled = df.copy()
resampled.timestamp = resampled.timestamp.apply(lambda x: parser.parse(x + ' UTC'))
resampled = resampled.set_index('timestamp').resample('5T').ffill().reset_index().dropna() # forward filling the TLE data for newly generated timestamps.
resampled[['lat','lon', 'alt']] = resampled.apply(lambda x: get_satellite_lla(tle_1=x.TLE_LINE1, tle_2=x.TLE_LINE2, timestamp=x.timestamp), axis=1)
resampled.to_csv(f"Starlink_{id}.csv", index=False)
fix = px.line(resampled, x='timestamp', y='alt')
fix.show()
session.close()

Chose a different satellite, which isn't being raised (but hasn't been put into station yet) and changed the upsampling from 30minutes to 5 minutes. enter image description here

The satellite still seems to swing up and down by close to 10km!

$\endgroup$
7
  • 1
    $\begingroup$ That looks reasonable to me. SGP4 is an analytic theory, expressed as polynomials of trig functions, which means it usually looks like several sine waves at different frequencies. $\endgroup$
    – Ryan C
    Jun 16, 2023 at 2:46
  • 3
    $\begingroup$ I don't really understand your process "generate the LLA (lat,lon,alt) for those 30 data points + timestamps between those 30 TLE's." Are you really using 30 different Two line element sets for this six day period? Also, I wonder if you are just seeing aliasing due to undersampling - the orbit is elliptical but your sparse sample points sometimes catch the max/min in altitude and sometimes miss them? Can you calculate a lot more closely spaced points and see what it looks like? $\endgroup$
    – uhoh
    Jun 16, 2023 at 5:04
  • 1
    $\begingroup$ You are perfectly welcome to post your code or script if that helps to explain how you're handling 30 TLEs. Paste a code block all indented by four spaces and it will format nicely. I mean if you have 30 TLEs equally spaced and 30 data points for each TLE, that's 900 points or one point every 10 minutes for six days. But I don't think you are showing that. I wonder if you are OVERLAPPING data from different TLEs such that adjacent points are coming from TLEs with very different epochs? I think you have to show your code/script to sort this out. That beating is to large and fast for LEO I think $\endgroup$
    – uhoh
    Jun 16, 2023 at 5:05
  • 3
    $\begingroup$ I agree with @uhoh this is mostly due to aliasing. There's a slight eccentricity of the orbit of about 20km. And the orbit is raised during the last two days. The rest is artifacts from your plotting technique. Try rerunning the simulation generating more points in step 3. $\endgroup$
    – asdfex
    Jun 16, 2023 at 11:15
  • 1
    $\begingroup$ So the question is about the larger 26-ish hour waveform you're getting? I agree that's likely an artifact from your low sample rate. The starlinks orbit every 90 minutes, so taking only a few readings per day is likely to introduce artifacts -- you're basically superposing the orbit frequency with the sampling frequency and getting a new waveform out of the interference between them. If you're into music, this is pretty much the same effect as the 'beating' effect you hear when you play a two pure tones together with one slightly out of tune. $\endgroup$ Jun 16, 2023 at 14:22

1 Answer 1

4
$\begingroup$

Edit: Responding to the edit/reframing of your question - you are calculating things correctly. The behavior you are seeing is real, although there are some caveats. These are described in my initial response below.

The primary cause of the plotted variance in altitude is the oblateness of the earth. The polar radius of the earth is 6357 km vs the equatorial radius which is 6378 km. A difference of 21 km.

Plotting true altitude at the instantaneous location of the object will hence result in changes due to the ground below moving further away from the satellite. A satellite in a pure 90 degree inclination polar orbit, would see a perfect sinusoid with 21 km amplitude in its altitude. Object 55775 has a fairly high inclination of 70 degrees, so you are seeing almost the full 21 km swing (361 to 382 km during June 2). Note that some software will report "altitude" not as the instantaneous value but rather referenced against a sphere with radius equal to the equatorial radius. This is usually also the way we talk about altitudes. When we say the ISS is at 400 km orbit, it is implied that this is altitude over the equator, even though altitude at northern most latitudes could be 15ish km more. To avoid confusion, you should see if you can plot the satellite's radius (measured from the center of the earth rather than altitude). This is why altitude is evil when dealing with orbit plots.

A second effect is the eccentricity of the orbit. Object 55775 seems to have around 20 km difference between apogee and perigee. Why the orbit has this eccentricity is a question for SpaceX and their concept of operations. However, you will also see that variation in the altitude plots. Because in the time frame of a few days perigee and apogee will be at constant latitude, the variation due to eccentricity and the variation due to earth oblateness will have the same exact frequency (once per orbit), but may be out of phase. That is why you see the double dips second plot.

Both of these are real and represent the true behavior of the satellite.

An artificial effect, as some pointed out in the comments, is likely due to your sampling (not plotting enough data points). That leads to the lower frequency oscillations in your initial plot. This is a process known as aliasing, where you are seeing the higher frequency oscillation due to altitude "aliased" to a lower frequency oscillation because you don't have enough data points. This does not appear in your second plot, because you changed the sampling to 5 minutes. Try changing the sampling from 30 minutes to 45 minutes, and things will get really weird as now you are sampling only two times per orbit. You are now at the "Nyquist limit" sampling at 2x times the signal frequency. That means that you will start to lose data and you may not even see the orbital period, or the true amplitude of the signal (altitude variation).

Other things you may notice, especially if you plot radius or semi major axis instead of altitude:

  • TLE's are not very accurate. Every time SGP4 changes to a new TLE, you will see a "jump" in the orbital parameters. 1 km change in semi major axis is not uncommon when switching from one TLE to another, especially if the epoch difference between TLEs is several days
  • The satellite may be doing propulsive maneuvers
  • Around June 4 this satellite started the process of increasing its altitude
  • If you were to use a higher accuracy propagator (not SGP4) that takes into account higher order gravity effects as well as 3rd body perturbations (moon, sun) you would see even more "jumping around"
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.