Skyfield and JPL disagreement on RA Dec for EarthSatellite #1093
-
|
Dear Brandon, first of all let me thank you for writing Skyfield. It is extremely useful and very fast! That is amazing! My background is research in solar system astronomy, mainly comets, so I know JPL Horizons and the MPC very well. Skyfield is a new tool for me but I love it already. I read the comment of the person who got their PhD in great part thanks to Skyfield. Now on to my problem. I am not being able to get JPL Horizons and Skyfield to agree on the RA and Dec of an earth satellite, given the same TLE. I tried to find a similar problem and the associated solution, but failed. Here's the output of the code below: The reason I need the precision in the RA and Dec is so that I can compare with, e.g. stars in the Gaia catalog, comets, or asteroids in the background. I appreciate any help anyone can give me in pointing out what I am doing wrong. Best regards, Here is the code. The ugly function is just to get the JPL Horizons coordinates using their API. The TLE and date are hardcoded so you can run it anytime. import json
import sys
from urllib.parse import quote
import pandas as pd
import requests
from skyfield.api import EarthSatellite, load, wgs84
def query_jpl_for_tle_and_time(t0, t1, tle) -> pd.DataFrame:
"""
Query the JPL Horizons API for the positional data of a satellite in a given time range.
Parameters
----------
t0 : skyfield.api.Time
The start time of the time range.
t1 : skyfield.api.Time
The end time of the time range.
tle : dict
The Two Line Element (TLE) data for the satellite.
Returns
-------
dict
The positional data of the satellite in a dictionary.
"""
# Convert times to the required UTC format: 'YYYY-MM-DD HH:MM:SS'
t0 = t0.utc_strftime("%Y-%m-%d %H:%M:%S")
t1 = t1.utc_strftime("%Y-%m-%d %H:%M:%S")
# formatted TLE must have name%0Aline1%0Aline2
formatted_tle = f"{tle['name']}\n{tle['line1']}\n{tle['line2']}"
# Define the query parameters
params = {
"format": "json",
"COMMAND": "TLE",
"OBJ_DATA": "YES",
"MAKE_EPHEM": "YES",
"EPHEM_TYPE": "OBSERVER",
"CENTER": "X02",
"START_TIME": f"'{t0}'",
"STOP_TIME": f"'{t1}'",
"STEP_SIZE": "1",
"TIME_DIGITS": "SECONDS",
"TIME_TYPE": "UT",
"QUANTITIES": "'1,4,20,47'",
"OUT_UNITS": "KM-S",
"RANGE_UNITS": "KM",
"SKIP_DAYLT": "NO",
"CSV_FORMAT": "YES",
"TLE": formatted_tle,
}
# Assemble the query URL
base_url = "https://ssd.jpl.nasa.gov/api/horizons.api"
query_string = "&".join(f"{k}={quote(str(v))}" for k, v in params.items())
url = f"{base_url}?{query_string}"
# Send the request
response = requests.get(url, timeout=10)
try:
data = json.loads(response.text)
except json.JSONDecodeError:
print(f"Failed to decode JSON response. Status code: {response.status_code}")
print(response.text)
sys.exit(0)
if response.status_code != 200:
print(f"Failed to get response. Status code: {response.status_code}")
print(response.text)
print("Quitting.")
sys.exit(0)
jpl_response = data["result"]
# find rows starting with $$SOE and $$EOE
ephem_lines = []
in_ephem = False
for i, row in enumerate(jpl_response.split("\n")):
if row.strip().startswith("$$EOE"):
in_ephem = False
if in_ephem:
ephem_lines.append(row)
if row.strip().startswith("$$SOE"):
in_ephem = True
keys = [
"time",
"daylight",
"M",
"ra",
"dec",
"az",
"alt",
"dist",
"rad_vel",
"sky_mot",
"sky_mot_PA",
"los_ang",
"empty_column",
]
ephem_lines = [ephem_line.split(",") for ephem_line in ephem_lines]
df = pd.DataFrame(ephem_lines, columns=keys).drop("empty_column", axis=1)
df[["az", "alt", "dist", "rad_vel", "sky_mot", "sky_mot_PA", "los_ang"]] = (
df[["az", "alt", "dist", "rad_vel", "sky_mot", "sky_mot_PA", "los_ang"]]
.astype(float)
.round(2)
)
df["time"] = pd.to_datetime(df["time"].str.strip())
df["ra"] = df["ra"].str.strip()
df["dec"] = df["dec"].str.strip()
return df
print("-----------------------------------------------")
ts = load.timescale()
tle = {
"name": "LANDSAT 9",
"line1": "1 49260U 21088A 25216.90508576 .00000326 00000+0 82394-4 0 9996",
"line2": "2 49260 98.2249 286.7798 0001208 89.3134 270.8202 14.57094879204947",
}
satellite = EarthSatellite(tle["line1"], tle["line2"], tle["name"], ts)
print(satellite)
print("-----------------------------------------------")
el_sauce = wgs84.latlon(
longitude_degrees=289.23517, latitude_degrees=-30.302776, elevation_m=1490.0
)
eph = load("de421.bsp")
t0 = ts.utc(2025, 8, 5, 23, 44, 40)
t1 = ts.utc(2025, 8, 6, 10, 3, 3)
t, events = satellite.find_events(el_sauce, t0, t1, altitude_degrees=30.0)
event_names = "rise above 30°", "culminate", "set below 30°"
sunlit = satellite.at(t).is_sunlit(eph)
for ti, event, sunlit_flag in zip(t, events, sunlit):
name = event_names[event]
state = ("in shadow", "in sunlight")[int(sunlit_flag)]
print(f"{ti.utc_strftime('%Y %b %d %H:%M:%S'):22} {name:15} {state}")
print("-----------------------------------------------")
t0 = ts.utc(2025, 8, 6, 3, 14, 34)
t1 = ts.utc(2025, 8, 6, 3, 14, 35)
times = ts.linspace(t0, t1, 2) # times of interest
days = t0 - satellite.epoch
print(f"{days:.3f} days away from epoch")
print("-----------------------------------------------")
difference = satellite - el_sauce
topocentric = difference.at(times)
alt, az, distance = topocentric.altaz()
print(f"Skyfield J2000 Coordinates")
ra, dec, distance = topocentric.radec() # ICRF ("J2000")
df_skyfield_j2000 = pd.DataFrame(
{
"time": times.utc_strftime("%Y-%m-%d %H:%M:%S"),
"ra": ra.hstr(format="{0}{1} {2:02} {3:02}.{4}"),
"dec": dec.dstr(format="{0}{1} {2:02} {3:02}.{4}"),
"az": az.degrees,
"alt": alt.degrees,
"dist": distance.km,
}
).round(2)
print(df_skyfield_j2000.to_string(index=False))
print("-----------------------------------------------")
print(f"Skyfield 'Epoch' Coordinates")
ra, dec, distance = topocentric.radec(epoch="date") # Epoch
df_skyfield_epoch = pd.DataFrame(
{
"time": times.utc_strftime("%Y-%m-%d %H:%M:%S"),
"ra": ra.hstr(format="{0}{1} {2:02} {3:02}.{4}"),
"dec": dec.dstr(format="{0}{1} {2:02} {3:02}.{4}"),
"az": az.degrees,
"alt": alt.degrees,
"dist": distance.km,
}
).round(2)
print(df_skyfield_epoch.to_string(index=False))
print("-----------------------------------------------")
jpl_result = query_jpl_for_tle_and_time(t0, t1, tle)
print(f"JPL Coordinates")
print(jpl_result[["time", "ra", "dec", "az", "alt", "dist"]].to_string(index=False)) |
Beta Was this translation helpful? Give feedback.
Replies: 1 comment 1 reply
-
|
Good afternoon, @siarw! Thanks very much for the kind words about Skyfield. Whenever the output from HORIZONS does not seem to match up, I encourage that you switch to its raw In this case, a big difference seems to be observer location. In your Python code you say this: But your call to HORIZONS, instead of repeating these exact coordinates, uses the special code X02. The detailed HORIZONS output suggests that X02 is at a different latitude than the one you are using with Skyfield: Try adjusting your Skyfield code so that it uses that same latitude: 30.4706218 south. I think you will find that the results give much better agreement. |
Beta Was this translation helpful? Give feedback.
-
|
Hi @brandon-rhodes , thanks for the super-fast reply. Yes! The location was wrong - silly I missed that. It makes a huge difference for such a low orbit when projecting to RA/Dec. Indeed, the JPL output has lots of info which most astronomers tend to ignore in most situations. But when precision is needed that info is key. Once again, many thanks for Skyfield! |
Beta Was this translation helpful? Give feedback.
Good afternoon, @siarw! Thanks very much for the kind words about Skyfield.
Whenever the output from HORIZONS does not seem to match up, I encourage that you switch to its raw
"CSV_FORMAT": "NO"mode, save out the entire resulting document to a file, and then go through it very carefully word by word. It contains all kinds of detailed information about time scales, coordinate systems, and coordinates, that can be compared to your Skyfield code to find the difference.In this case, a big difference seems to be observer location. In your Python code you say this:
But your call to HORIZONS, instead of repeating the…