AltAz radial_velocity transformation fails

[tvwenger]

Description

Transforming coordinates that include a radial_velocity to or from AltAz fails with an obscure error. This is related to #12593.

This functionality is useful for observation planning and Doppler corrections. For example, if I know that an object has a given radial velocity in the GalacticLSR frame, I might want to know what the position is in the AltAz frame, so I know where to point the telescope, and what the radial velocity is in the topocentric frame, so that I can set my receiver system to observe the correct frequency.

Steps to Reproduce

import astropy.units as u
from astropy.time import Time
from astropy.coordinates import EarthLocation, ICRS, AltAz
import astropy
print(astropy.__version__)
# 5.0 (also tried dev version: 5.1.dev251+gcd25e6329)

location = EarthLocation.of_site("ALMA")
obstime = Time('2021-12-13T21:10:00.624')

coord = ICRS(356.38*u.deg, -71.154*u.deg, radial_velocity=0.0*u.km/u.s)
coord = coord.transform_to(AltAz(location=location, obstime=obstime))
# UnitConversionError: Can only apply 'add' function to quantities with compatible dimensions

I suspect that the issue is that AltAz and ICRS have different origins (the former at location, and the later at the solar system barycenter). Sure enough, there is no error if you supply additional information when creating the ICRS coordinate:

coord = ICRS(356.38*u.deg, -71.154*u.deg, radial_velocity=0.0*u.km/u.s,
             pm_ra_cosdec=0.0*u.mas/u.yr, pm_dec=0.0*u.mas/u.yr,
             distance=1.0*u.Mpc)
coord = coord.transform_to(AltAz(location=location, obstime=obstime))
print(coord)
# <AltAz Coordinate (obstime=2021-12-13T21:10:00.624, location=(2225015.30883296, -5440016.41799762, -2481631.27428014) m, pressure=0.0 hPa, temperature=0.0 deg_C, relative_humidity=0.0, obswl=1.0 micron): (az, alt, distance) in (deg, deg, Mpc)
#    (169.99955542, 40.00005221, 1.)
# (pm_az_cosalt, pm_alt, radial_velocity) in (mas / yr, mas / yr, km / s)
#    (1.34282337e+11, 7.58579968e+10, 1151.58325584)>

It doesn't give an error, but it gives unrealistic values for the radial velocity and proper motions.

Expected behavior

We should be able to convert between ICRS and AltAz coordinates for coordinates that include a radial_velocity.

Actual behavior

ICRS to AltAz fails with an obscure error message when only a radial_velocity is supplied. If radial_velocity and proper motions are supplied, the conversion does not have an error, but the resulting kinematics are wildly incorrect.

System Details

Linux-5.10.60.1-microsoft-standard-WSL2-x86_64-with-glibc2.10
Python 3.8.5 (default, Sep 4 2020, 07:30:14)
[GCC 7.3.0]
Numpy 1.21.4
pyerfa 2.0.0.1
astropy 5.1.dev251+gcd25e6329

[tvwenger]

Basically, I would expect that a transformation from AltAz to ICRS (or another barycentric frame) would result in the coordinate radial_velocity to be modified by the coordinate's radial_velocity_correction. For example:

import astropy
from astropy.coordinates import EarthLocation, SpectralCoord, SkyCoord, AltAz, ICRS
from astropy.time import Time
import astropy.units as u

location = EarthLocation.of_site("ALMA")
obstime = Time('2021-12-13T21:10:00.624')

target = AltAz(az=170.0*u.deg, alt=40.0*u.deg, location=location,
               obstime=obstime, radial_velocity=0.0*u.km/u.s)
target_icrs = target.transform_to(ICRS())
# Fails with noted error

# I would expect an ICRS coordinate with a radial velocity given
# by `radial_velocity_correction`. For example:
target = SkyCoord(az=170.0*u.deg, alt=40.0*u.deg, location=location,
                  obstime=obstime, frame="altaz")
target_radial_velocity = 0.0*u.km/u.s
barycorr = target.radial_velocity_correction()
expected_target_icrs = SkyCoord(
    target.icrs, radial_velocity=target_radial_velocity + barycorr)

[adrn]

For the first error you note, your interpretation is right: ICRS is defined with an origin at the SS barycenter whereas AltAz is defined from a location on the WGS84 ellipsoid (~Earth surface). In order to transform a barycentric, spherical velocity between these frames, astropy.coordinates must know the full 3D velocity and position of the source. I think the error message is confusing and we should probably improve that!

The huge proper motions you get even after specifying the full 3d position and velocity I think are correct, and they're telling you what Earth's rotation is (2π rad/day ~= 1e11 mas/yr). The radial velocity looks wrong though... that could be a bug!

Pinging @eteq who probably knows more about the AltAz implementation off the top of his head.

[mkbrewer]

One thing that I would note is that "proper_motion" in AltAz is just the rate of change of the coordinates. More proper units would be arcsec per second, the conversion factor being 1 / 3.15576e+10. This yields the right order of magnitude, but the numbers still seem a bit off as I would expect the square root of the sum of the squares to be 15 arcsec per second.

[mkbrewer]

Regarding the radial velocity, if the target's radial velocity wrt the SSB is 0, the most it can be is ~30 km/sec wrt an observer on Earth.

[mkbrewer]

AltAz seemed like a rather strange frame to be doing this in, so I tried TETE and got:

>>> location = EarthLocation.of_site("ALMA")
>>> obstime = Time('2021-12-13T21:10:00.624')
>>> coord = ICRS(356.38*u.deg, -71.154*u.deg, radial_velocity=0.0*u.km/u.s,
...              pm_ra_cosdec=0.0*u.mas/u.yr, pm_dec=0.0*u.mas/u.yr,
...              distance=1.0*u.Mpc)
>>> coord = coord.transform_to(TETE(location=location, obstime=obstime))
>>> coord
<TETE Coordinate (obstime=2021-12-13T21:10:00.624, location=(2225015.30883296, -5440016.41799762, -2481631.27428014) m): (ra, dec, distance) in (deg, deg, Mpc)
    (356.68330707, -71.03910289, 1.)
 (pm_ra_cosdec, pm_dec, radial_velocity) in (mas / yr, mas / yr, km / s)
    (159345.92360131, -596898.48095559, -5577.83235682)>

So there is definitely a bug here.

[mkbrewer]

And CIRS yields something completely different, which it shouldn't.

>>> coord = coord.transform_to(CIRS(location=location, obstime=obstime))
>>> coord
<CIRS Coordinate (obstime=2021-12-13T21:10:00.624, location=(2225015.30883296, -5440016.41799762, -2481631.27428014) m): (ra, dec, distance) in (deg, deg, Mpc)
    (356.40595924, -71.03910289, 1.)
 (pm_ra_cosdec, pm_dec, radial_velocity) in (mas / yr, mas / yr, km / s)
    (147989.52564203, -596897.14656101, 217.77449531)>

[mkbrewer]

And CIRS at the geocenter:

>>> coord = coord.transform_to(CIRS(obstime=obstime))
>>> coord
<CIRS Coordinate (obstime=2021-12-13T21:10:00.624, location=(0., 0., 0.) km): (ra, dec, distance) in (deg, deg, Mpc)
    (356.40572924, -71.03913461, 1.)
 (pm_ra_cosdec, pm_dec, radial_velocity) in (mas / yr, mas / yr, km / s)
    (-129835.56232232, -11380.7407748, 202.5673867)>

No way is ALMA moving 15.2 km/s wrt the geocenter. Not to mention the totally bogus difference in proper motion.

[mkbrewer]

I have never looked at velocities in Astropy as they are done with differentials, which seems over complicated, but I guess I should.

[mkbrewer]

Looking into this, I see that the issue is documented here: https://docs.astropy.org/en/stable/coordinates/velocities.html#finite-difference-transformations, so it's a known problem.

[mkbrewer]

It seems to me that the problem here is the way that velocities are computed and intimately coupled to space motion. Using finite differences is rather compute intensive and clearly wildly inaccurate for anything beyond the solar system. Maybe I'm missing something, but it seems that finite differences are only needed for calculating the velocities of solar system objects that don't have an ephemeris such as comets, solar flares, etc. In fact, this is what stymied my ill-fated attempt at getting an astrometric frame implemented.

Many stars have a known distance along with components of proper motion and radial velocity. For these, position and velocity vectors can be constructed and full 3D space motion can be derived. Parallax can be computed by subtracting the observer's position vector wrt the SSB from that of the target body and redshift can be computed by adding in the component of the observer's velocity wrt the SSB along the line of site to the radial velocity. Please note here that components of the observer's velocity wrt the SSB should NOT be added in to the components of proper motion. Proper motion is defined as the secular rates of change in right ascension and declination (Explanatory Supplement, 3rd Edition, pg. 28).

On the other hand, many objects have no given distance, but still have a radial velocity and perhaps proper motion associated with them. For these, no true position and velocity vectors can be constructed. Parallax cannot be computed, but 2D space motion can be derived from the proper motion (if any) applied to the position unit vectors and redshift can be derived from the radial velocity.

Solar system bodies such as planets, asteroids and comets can be included in the first case. In the case of planets, position and velocity vectors can be read from the ephemerides. For asteroids and comets, orbital elements can be used to construct a model which yields short term position and velocity predictions. Radial velocities can be computed in order to observe lines in the atmospheres of planets and the head and tail of comets. All of this without a single finite difference.

In radio astronomy, we normally only deal with the second case. We are given a set of RA/Dec coordinates along with a VLSR (LSR usually being what Astropy calls LSRK) for each source. From this, we would like to compute the radial velocity of the source in the rest frame of the observer and the sky frequency for a given line in order to tune our receiver. It is a bit of a peeve of mine that the SpectralCoord class doesn't have a method for doing this. Now, I realize that most astronomers these days are mainly consumers of data and aren't involved too much with the hardware. Hence the focus of the SpectralCoord class on data analysis. Historically, however, tuning the receiver was very much the job of the observer. Back then we had to physically tune the Klystron tube by adjusting its micrometer head amongst other things before any observing session could get underway. Yep, I'm sounding like the old fogey that I am reminiscing about the bad old days.

Anyway, I experimented around with this and it can be done. First I defined a SkyCoord in the LSRK frame with a VLSR of 10 km / s:

>>> location = EarthLocation.of_site("ALMA")
>>> obstime = Time('2021-12-13T21:10:00.624')
>>> target = SkyCoord(271.0*u.deg, 30.0*u.deg, radial_velocity=10.0*u.km/u.s,
...                   pm_ra_cosdec=0.0*u.mas/u.yr, pm_dec=0.0*u.mas/u.yr, 
...                   distance=1.0*u.Mpc, location=location, obstime=obstime,
...                   frame='lsrk')

Then I defined a SpectralCoord using this target.

>>> observer = location.get_itrs(obstime=obstime)
>>> sp = SpectralCoord(100.0 * u.GHz, observer=observer, target=target)
WARNING: NoVelocityWarning: No velocity defined on frame, assuming (0., 0., 0.) km / s. [astropy.coordinates.spectral_coordinate]

This gave a warning, but it worked.

>>> sp
<SpectralCoord 
   (observer: <ITRS Coordinate (obstime=2021-12-13T21:10:00.624): (x, y, z) in m
                  (2225015.30883296, -5440016.41799762, -2481631.27428014)
               (v_x, v_y, v_z) in km / s
                  (0., 0., 0.)>
    target: <LSRK Coordinate: (ra, dec, distance) in (deg, deg, Mpc)
                (271., 30., 1.)
             (pm_ra_cosdec, pm_dec, radial_velocity) in (mas / yr, mas / yr, km / s)
                (0., 0., 10.)>
    observer to target (computed from above):
      radial_velocity=-6.767342719936956 km / s
      redshift=-2.2573170738726667e-05)
  100. GHz>

I chose RA 271, Dec 30 for the target as that is the direction that the SSB is moving in relative to the LSRK frame at J2000. This makes it easy to verify the radial velocity of the target wrt the SSB.

>>> sp.with_observer_stationary_relative_to('icrs')
<SpectralCoord 
   (observer: <ICRS Coordinate: (x, y, z) in m
                  (1.97765916e+10, 1.34195248e+11, 5.81997515e+10)
               (v_x, v_y, v_z) in km / s
                  (0., 0., 0.)>
    target: <LSRK Coordinate: (ra, dec, distance) in (deg, deg, Mpc)
                (271., 30., 1.)
             (pm_ra_cosdec, pm_dec, radial_velocity) in (mas / yr, mas / yr, km / s)
                (0., 0., 10.)>
    observer to target (computed from above):
      radial_velocity=-9.999996166444316 km / s
      redshift=-3.3355840426407823e-05)
  100.0010783 GHz>

Close enough to the -10 km / s that I expected. The documentation for the LSRK frame is incorrect, btw.

This frame is defined as having a velocity of 20 km/s towards RA=270 Dec=30 (B1900) relative to the solar system Barycenter. 

should read:

This frame is defined such that the solar system Barycenter has a velocity of 20 km/s towards RA=270 Dec=30 (B1900) relative to it. 

Then I did this:

>>> from astropy.coordinates.spectral_coordinate import _apply_relativistic_doppler_shift
>>> _apply_relativistic_doppler_shift(sp, sp.radial_velocity)
<SpectralQuantity 100.00225737 GHz>

This is the correct sky frequency for a line with a rest frequency of 100 GHz in a source with a velocity of 10 km / s wrt LSRK.

So, what I would propose is adding a from_rest() method to the SpectralCoord class. Whereas to_rest() is defined as transforming a SpectralQuantity representing a sky frequency in the rest frame of the observer to a rest frequency in the rest frame of the target, from_rest() would be defined as transforming a SpectralQuantity representing a rest frequency in the rest frame of the target to a sky frequency in the rest frame of the observer.

I will raise a separate issue regarding Astropy's computation of velocities and cross post much of what I have posted here. Sorry about the "wall of words", but I wanted to thoroughly address your issues.