SpectralCoord gives unexpected results with AltAz target

[tvwenger]

Description

Using SpectralCoord with a target specified in the AltAz frame yields unexpected results. There are some invalid value warnings, NaNs, and incorrect transformation results. This does not happen when the same target is passed in the ICRS frame.

Steps to Reproduce

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

freq_topo = np.arange(1.0, 2.0, 0.1)*u.GHz

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

observer = location.get_itrs(obstime=obstime)
target = SkyCoord(170.0*u.deg, 40.0*u.deg, location=location, obstime=obstime, frame="altaz")

spec_altaz = SpectralCoord(freq_topo, observer=observer, target=target)
spec_altaz_lsrk = spec_altaz.with_observer_stationary_relative_to('lsrk')

spec_icrs = SpectralCoord(freq_topo, observer=observer, target=target.icrs)
spec_icrs_lsrk = spec_icrs.with_observer_stationary_relative_to('lsrk')

Expected behavior

SpectralCoord should yield the same results in both the AltAz and ICRS frames.

Actual behavior

print(spec_altaz.redshift)
# .../astropy/units/quantity.py:614: RuntimeWarning: invalid value encountered in sqrt
#  result = super().__array_ufunc__(function, method, *arrays, **kwargs)
# nan
print(spec_icrs.redshift)
# 3.746574655116852e-05
print(spec_altaz_lsrk)
# [1.0000701  1.10007711 1.20008412 1.30009113 1.40009814 1.50010515
#  1.60011216 1.70011918 1.80012619 1.9001332 ] GHz
print(spec_icrs_lsrk)
# [1.00006755 1.1000743  1.20008106 1.30008781 1.40009457 1.50010132
#  1.60010808 1.70011483 1.80012159 1.90012834] GHz

System Details

Linux-3.10.0-1160.49.1.el7.x86_64-x86_64-with-glibc2.10
Python 3.8.3 (default, Jul 2 2020, 16:21:59)
[GCC 7.3.0]
Numpy 1.21.4
pyerfa 2.0.0.1
astropy 5.0
Scipy 1.6.2
Matplotlib 3.3.4

[github-actions]

Welcome to Astropy 👋 and thank you for your first issue!

A project member will respond to you as soon as possible; in the meantime, please double-check the guidelines for submitting issues and make sure you've provided the requested details.

GitHub issues in the Astropy repository are used to track bug reports and feature requests; If your issue poses a question about how to use Astropy, please instead raise your question in the Astropy Discourse user forum and close this issue.

If you feel that this issue has not been responded to in a timely manner, please leave a comment mentioning our software support engineer @embray, or send a message directly to the development mailing list. If the issue is urgent or sensitive in nature (e.g., a security vulnerability) please send an e-mail directly to the private e-mail feedback@astropy.org.

[tvwenger]

I've pinned down at least one issue to the following function:

astropy/astropy/coordinates/spectral_coordinate.py

Line 482 in cd25e63

def _calculate_radial_velocity(observer, target, as_scalar=False):

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

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

observer = location.get_itrs(obstime=obstime)
target = SkyCoord(170.0*u.deg, 40.0*u.deg, location=location,
                  obstime=obstime, frame="altaz")

observer = SpectralCoord._validate_coordinate(observer)
target = SpectralCoord._validate_coordinate(target)
print(SpectralCoord._calculate_radial_velocity(observer, target, as_scalar=True))
# 2032008.296875 km / s

observer = location.get_itrs(obstime=obstime)
target = SkyCoord(170.0*u.deg, 40.0*u.deg, location=location,
                  obstime=obstime, frame="altaz")
target = target.icrs

observer = SpectralCoord._validate_coordinate(observer)
target = SpectralCoord._validate_coordinate(target)
print(SpectralCoord._calculate_radial_velocity(observer, target, as_scalar=True))
# 11.231737842710979 km / s

Since the first step in _calculate_radial_velocity is to convert the coordinate to ICRS, I suspect that the issue is in _validate_coordinate.

Indeed, if I try to set a radial_velocity on the target at all, I get a whole new error:

target = SkyCoord(170.0*u.deg, 40.0*u.deg, location=location,
                  obstime=obstime, frame="altaz", radial_velocity=0.0*u.km/u.s)
target = SpectralCoord._validate_coordinate(target)
# ValueError: alt has unit "deg" with physical type "angle", but pm_alt has incompatible unit "km rad / s" with physical type "unknown" instead of the expected "angular frequency/angular speed/angular velocity".

I will look into this further and, if I find the problem, I'll submit a PR!

[tvwenger]

I think I've pinned down the problem to simply how radial_velocity is handled in AltAz. Since _validate_coordinate assigns a zero velocity to the coordinate in the original frame, before it is transformed to ICRS in _calculate_radial_velocity, this bug manifests itself there.

For example, this fails and I don't think it should:

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

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))

I will submit another issue related to this bug.