Coordinates transformation in Astropy

Question

I am doing a transformation of coordinates in Astropy. But this is not working correctly. The situation is the following, I have a Coordinate of a star in J2000, and I want to transform it into actual epoch (or another specific date). I am verifying the code with the coordinates providing by Stellarium. The problem is that the results are different in 20 minutes approx and that is a great problem because I need a precision of some tens of milliseconds. The utilized code is the following:

from astropy.coordinates import SkyCoord
from astropy.coordinates import FK5

c = SkyCoord(20.398617733743833, 38.466348612533892,
         unit='deg', frame='icrs')

c_fk5 = c.transform_to(FK5(equinox='2018-10-19 00:19:41'))

I would like to know if my method is wrong and how I can get better results. Please, everyone be free to answer and post comments and suggestions. Thanks.

Answer 1

You are mixing up equinox (the classical zeropoint for right ascension) and epoch (when something is being observed). Also FK5 is the old celestial reference frame which has been replaced by ICRS and the two frames are very closely aligned so it doesn't make sense to transform between the two.

Astropy uses ERFA which is a re-implementation of the SOFA code which implements the IAU2000/2006 resolutions to calculate positions and is using the modern CIO-based transformations. Stellarium is using the older equinox-based approach. You can see the difference between the two approaches in Figure 2 of the SOFA Tools for Earth Attitude Cookbook where Astropy goes down the right-hand CIO branch and Stellarium goes down the left-hand equinox-based branch.

You can calculate what you are asking i.e. "where in the sky is my object which had these ICRS co-ordinates at J2000.0 at my time of 2018-10-19 00:19:41 UTC" which would be the Celestial Intermediate Reference System (CIRS) using astropy, but this is not comparable to geocentric apparent positions reported in Stellarium. This is because the models for where the Earth's pole is (the precession-nutation models), the origins of right ascension and the model of Earth rotation are different between the two approaches. This makes it harder to compare things until the two approaches have met up again at the local apparent [h, delta] step in Figure 2 (corresponding to HA/Dec (apparent) in Stellarium). It is not helped by the fact that Astropy doesn't really have a HA/Dec frame as it is usually an intermediate step to an Alt/Az one i.e. how far West/East and above the horizon an object will be for an Earth-based observer.

The following code should let you calculate a local apparent HA, Dec for comparison with Stellarium (providing you sent the date and place in Stellarium to be the same). This will work for a distant object for which proper motion and parallax are negligible; otherwise you will need to add these in when you declare the SkyCoord - see Using velocities with SkyCoord for more details.

import astropy.units as u
from astropy.coordinates import SkyCoord, ITRS, EarthLocation
from astropy.time import Time

c = SkyCoord(20.398617733743833, 38.466348612533892, unit='deg', frame='icrs')
t = Time('2018-10-19 00:19:41', scale='utc')
loc = EarthLocation(lon=30*u.deg, lat=30*u.deg, height=0*u.m)

c_ITRS = c.transform_to(ITRS(obstime=t))
# Calculate local apparent Hour Angle (HA), wrap at 0/24h
local_ha = loc.lon - c_ITRS.spherical.lon
local_ha.wrap_at(24*u.hourangle, inplace=True)
# Calculate local apparent Declination
local_dec = c_ITRS.spherical.lat
print("Local apparent HA, Dec={} {}".format(local_ha.to_string(unit=u.hourangle, sep=':'), local_dec.to_string(unit=u.deg, sep=':', alwayssign=True) ))

There will be some difference due to the different models used, accounting for the Earth orientation parameters (UT1-UTC, polar motion) etc but they should compare at the sub-second level.