Different answers between .radec() and .altaz()

[ngc1535git]

I have recently started using Skyfield and it appears at first glance that the .radec() results do not agree with the .altaz() answer if I take the .altaz() answer and covert back to RA and Dec (or vice versa).

PyEphem says the satellite is at these coordinates at this time:
UTC Local AZ Alt RA Dec Airmass Velocity (""/sec) Elongation
2023/8/24 05:12:00 2023/8/23 19:12:00 229.44 56.25 229.745176 -3.326477 1.19 15.02 77.82

Skyfield outputs:
RA (degrees): 229.42112616876986, Dec (degrees): -3.2349986513594478, Az (degrees): 229.4459062645183, Alt (degrees): 56.24291331265545

A calculator with input of the Skyfield RA/Dec yeilds:
Julian date : 2460180.716666667
Local siderial time : 16.984115239341946
Altitude : 56.11780771891635
Azimuth : 229.88028053569792

My issue is that Pyephem and Skyfield agree on the Az/Alt. They do not agree on the RA/Dec. Furthermore, taking the Skyfield RA/Dec and converting to Alt/Az does not seem to yield the same answer. Is there some geometric term I need to take into account for the .radec() method?

The relevant part of my code is below.
Any help would be appreciated!

def find_nearby_satellites_by_az_alt(az_alt, all_satellites, ts, search_radius):
	# Suppress arctan warnings 
	warnings.filterwarnings("ignore", category=RuntimeWarning, module="skyfield")
	ref_az, ref_alt = az_alt #degrees
	satellites_within_threshold = []
	

	for _, (norad_id, sat_name, line1, line2) in all_satellites.items():
		satellite = EarthSatellite(line1, line2)
		observed_satellite = (satellite -ukirt).at(ts)
		alt,az, _ = observed_satellite.altaz()
		ra,dec,_ = observed_satellite.radec()

		angular_dist = angular_distance(ref_az, ref_alt, az._degrees, alt._degrees)

		if angular_dist <= search_radius:
			satellites_within_threshold.append((norad_id, sat_name, line1, line2,angular_dist))
			logging.debug("RA (degrees): %s, Dec (degrees): %s, Az (degrees): %s, Alt (degrees): %s",ra.hours * 15, dec.degrees, az.degrees, alt.degrees)
			
	# Sort the list in ascending order based on angular_dist
	satellites_within_threshold = sorted(satellites_within_threshold, key=lambda x: x[4])
	

	
	return satellites_within_threshold

[brandon-rhodes]

My issue is that Pyephem and Skyfield agree on the Az/Alt. They do not agree on the RA/Dec.

I don't think we can answer this without seeing your code that produced both numbers. Both PyEphem and Skyfield are capable of producing both J2000 coordinates, but also equinox-of-date coordinates, and there is additionally the question of apparent versus astrometric coordinates. If you made different choices in the two libraries, you would get slightly different angles.

The relevant part of my code is below.

To help Skyfield contributors dig deeper, try constructing a small stand-alone script that folks can run, that uses PyEphem and Skyfield to print out the two sets of coordinates that concern you, and maybe that then prints out how big the difference in coordinates is. Then folks will be able to see how big the difference is, and also will be able to reproduce the numbers on their own computers.

[ngc1535git]

I think my basic question is... why is the Pyephem and SkyField output for Az/Alt so close to agreement?
This confused me. This 5 minute video explains/demonstrates:

(first video on the page)
https://www.adamblockstudios.com/articles/galaxy-cellular-automata

The working code for the program is here:
https://gist.github.com/ngc1535git/6cee5084ecf4fddde4c91161626a6c27
You just need to put in your credentials (or disable the fetch and point to the satcat.txt file)

-adam

P.S. The code that prints out two sets of coordinates is here:

# -*- coding: utf-8 -*-
"""
Created on Thu Aug 24 00:43:31 2023

@author: ablock
"""

import ephem
import requests
from skyfield.api import Topos, load, EarthSatellite
from datetime import datetime, timedelta
import math

# Fetch TLE data
#url = "https://www.celestrak.com/NORAD/elements/active.txt"
#response = requests.get(url)
tle_all_data= []
filename = 'satcat.txt'
with open(filename, 'r') as file:
	lines = file.readlines()
for line in lines:
	tle_all_data.append(line)

# Find the TLE corresponding to NORAD ID 38978
for i in range(0, len(tle_all_data), 3):
	if "YAMAL 300K" in tle_all_data[i]:
		print(tle_all_data[i])
		tle_line1 = tle_all_data[i + 1]
		tle_line2 = tle_all_data[i + 2]
		print(tle_line1,tle_line2)
		break

# Define Observer location
obs = ephem.Observer()
obs.lat, obs.lon = '19.8224', '-155.4703'  #UKIRT

# Use pyephem
sat = ephem.readtle(tle_all_data[i], tle_line1, tle_line2)

# Use Skyfield
ts = load.timescale()
satellite = EarthSatellite(tle_line1, tle_line2, tle_all_data[i], ts)
topos = Topos(latitude_degrees=19.8224, longitude_degrees=-155.4703,elevation_m = 4800)

# Initialize time
current_time = datetime.utcnow()

# Time increment
time_increment = timedelta(minutes=1)

# Header for the output
print("Date-Time,                  RA (PE) RA(SF) Dec(PE), Dec(SF)  Az(PE) Az(SF)   Alt(PE)  Alt(SF)")

# Loop to increment time and calculate positions
for i in range(10):  # 10 increments; you can adjust the number as needed
	# pyephem calculations
	obs.date = ephem.Date(current_time)
	sat.compute(obs)
	ra_pyephem_decimal = math.degrees(float(sat.ra))
	dec_pyephem_decimal = math.degrees(float(sat.dec))
	az_pyephem_decimal = math.degrees(float(sat.az))
	alt_pyephem_decimal = math.degrees(float(sat.alt))
	# Skyfield calculations
	t = ts.utc(current_time.year, current_time.month, current_time.day, current_time.hour, current_time.minute, current_time.second)
	astrometric = (satellite - topos).at(t)
	ra , dec, _ = astrometric.radec() 
	alt, az, d = astrometric.altaz()
	ra = ra.hours *15

	# Print output
	print(f"{current_time}, {ra_pyephem_decimal:.3f}, {ra:.3f}, {dec_pyephem_decimal:.3f}, {dec.degrees:.3f}, {az_pyephem_decimal:.3f}, {az.degrees:.3f}, {alt_pyephem_decimal:.3f}, {alt.degrees:.3f}")

	# Increment time
	current_time += time_increment

[brandon-rhodes]

Thank you for the example script, which made it easier to test some possibilities! (And also for the video; it was almost kind of startling, for the very first time ever, to hear one of my Skyfield users as an actual voice instead of simply text on a screen.)

The answer to your question is that you are asking for coordinates in two different RA/dec systems. The two are explained here (maybe you'll suggest a way to make this section even clearer for folks facing your situation):

https://rhodesmill.org/skyfield/coordinates.html#icrs-and-equatorial-right-ascension-and-declination

One RA/dec system is the J2000 (or, officially, "ICRS" system) which never moves and is permanent. It's not PyEphem's usual system, but you can print it out with:

	ra_pyephem_decimal = math.degrees(float(sat.a_ra))
	dec_pyephem_decimal = math.degrees(float(sat.a_dec))

You'll see that the results agree very well with Skyfield.

The other RA/dec system is based on the Earth's actual pole and equator as precession gradually skews the coordinate system across the sky over thousands of years. It's the one PyEphem is using in your script. To have Skyfield use the same system, try this:

	ra , dec, _ = astrometric.radec(epoch='date')

Agreement doesn't look to be quite as high when both PyEphem and Skyfield are using this "equinox-of-date" system, but it's still closer than the difference between one system to the other. The discrepancy might simply be that PyEphem and Skyfield use different models of the Earth's precession and nutation.

Anyway, try these two adjustments, and see if you can confirm that Skyfield and PyEphem think the satellite is in just about the same place in the sky, but that they are naming that place differently in the two slightly different RA/dec coordinate systems.

[ngc1535git]

Thank you Brandon. I hope I can impose on your patience a little more. I still do not quite understand. 1. If I do not specify an epoch, do I understand correctly we are computing (propagating) for epoch.now() (or whatever the method)? It seems like from the below you are suggesting that if I put both pyephem and Skyfield at the same epoch- I will get answers that ae closer. I understand there will be differences (see #2 below). What is not clear to me is that why would specifying "2000" do this but letting them both be epoch.now would not? You say PyEphem usual system is not J2000... what is it? 2. My main concern was that by switching between RA/Dec and Alt/AZ I get different answers. If you covert the AZ/Alt propagated values to RA/Dec you do not get the propagated RA/Dec values. How could you? Both PyEphem and Skyfield agree on he Az/Alt but not the RA/Dec (which again, I get why that is true). What is special about Az/Alt? I will modify my code and check it out shortly.

…

-adam

On 08/24/2023 12:34 PM MST Brandon Rhodes ***@***.***> wrote: Thank you for the example script, which made it easier to test some possibilities! (And also for the video; it was almost kind of startling, for the very first time ever, to hear one of my Skyfield users as an actual voice instead of simply text on a screen.) The answer to your question is that you are asking for coordinates in two different RA/dec systems. The two are explained here (maybe you'll suggest a way to make this section even clearer for folks facing your situation): https://rhodesmill.org/skyfield/coordinates.html#icrs-and-equatorial-right-ascension-and-declination One RA/dec system is the J2000 (or, officially, "ICRS" system) which never moves and is permanent. It's not PyEphem's usual system, but you can print it out with: ra_pyephem_decimal = math.degrees(float(sat.a_ra)) dec_pyephem_decimal = math.degrees(float(sat.a_dec)) You'll see that the results agree very well with Skyfield. The other RA/dec system is based on the Earth's actual pole and equator as precession gradually skews the coordinate system across the sky over thousands of years. It's the one PyEphem is using in your script. To have Skyfield use the same system, try this: ra , dec, _ = astrometric.radec(epoch='date') Agreement doesn't look to be quite as high when both PyEphem and Skyfield are using this "equinox-of-date" system, but it's still closer than the difference between one system to the other. The discrepancy might simply be that PyEphem and Skyfield use different models of the Earth's precession and nutation. Anyway, try these two adjustments, and see if you can confirm that Skyfield and PyEphem think the satellite is in just about the same place in the sky, but that they are naming that place differently in the two slightly different RA/dec coordinate systems. — Reply to this email directly, view it on GitHub #894 (reply in thread), or unsubscribe https://github.com/notifications/unsubscribe-auth/ALJOA7JBKJZSI5X5E3BRBC3XW6UFBANCNFSM6AAAAAA34FETMM. You are receiving this because you authored the thread.Message ID: ***@***.***>

[brandon-rhodes]

1. I think your idea is perfect to "modify my code and check it out shortly." Before answering in more detail, I'm going to give you time to try out the code changes, and especially to read the documentation that I linked to; I think that should hopefully clear up your confusion about what an epoch/equinox date is and how it affects coordinates (but not positions!).
2. Converting from alt/az to ra/dec is fraught with design decisions that different calculators make in different ways; for example, should they account for atmospheric refraction or not? So I'm not confident, unless I'd written the code for the calculator myself and knew which effects it was accounting for, whether its results could be compared directly to Skyfield or not. But the big question to ask is whether it's returning J2000 RA/dec, or whether it's returning RA/dec for the equinox of the particular date of the observation. Then smaller questions arise: whether it's using the same Earth rotation/position model, and so forth.

[ngc1535git]

OK. I got it.
The default PyEphem is precessed (representing "now") and SkyField is not ...it is the ICRS fixed coordinates which roughly corresponds to J2000. (this is the sentence I think I needed to hear to get me on track) I understand there are additional terms of nutation and other things that may or may not be included in each.

The two answers refer to the same position in the sky- it is just that with one the time element of precession means I am not comparing the same coordinates.

If I took the coordinates that SkyField is giving me... and precessed them forward in time... the angular distance "error" I was measuring will likely be of similar size.

Thank you for your patience. This was an important exercise.

-adam

[brandon-rhodes]

The two answers refer to the same position in the sky- it is just that with one the time element of precession means I am not comparing the same coordinates.

Yes, that's exactly correct! I'm glad that we tracked down the problem, and I hope your project goes forward more smoothly now.