Cannot use SC.apply_space_motion twice on the same SC

Hi everyone,

I’m noticing some strange behaviour when using apply_space_motion to a SkyCoord. First of all, it doesn’t only return a new SC, but the original SC gets altered, which is surprising. And when I try to run it again, it fails with an error message which I can’t quite parse. I’ve created a minimal viable example, posted below. This occurs both with Astropy 5.0.2, and 5.3.4.

from astropy.coordinates import SkyCoord
from astropy.time import Time
import astropy.units as u
import astropy

print('Version:', astropy.__version__)

T = Time([2020, 2021, 2021], format = 'decimalyear')

SC = SkyCoord(15 * u.deg, 15 * u.deg,
    pm_ra_cosdec = 750 * u.mas/u.year, pm_dec = -200 * u.mas/u.year,
    distance = 21 * u.pc, obstime  = T[0])

print('Original SC:')
print(SC)

SC1 = SC.apply_space_motion(new_obstime = T[1])

print('Original SC, now mysteriously has a radial velocity added')
print(SC)
print('New SC, looks fine')
print(SC1)

# This second call to apply_space_motion fails
print('Now calling SC.apply_space_motion again...')
SC2 = SC.apply_space_motion(new_obstime = T[2])

The output including error message stack is:

Astropy Version: 5.3.4
Original SC:
<SkyCoord (ICRS): (ra, dec, distance) in (deg, deg, pc)
    (15., 15., 21.)
 (pm_ra_cosdec, pm_dec) in mas / yr
    (750., -200.)>
Original SC, now mysteriously has a radial velocity added:
<SkyCoord (ICRS): (ra, dec, distance) in (deg, deg, pc)
    (15., 15., 21.)
 (pm_ra_cosdec, pm_dec, radial_velocity) in (mas / yr, mas / yr, km / s)
    (750., -200., -1.07785819e-15)>
New SC, looks fine:
<SkyCoord (ICRS): (ra, dec, distance) in (deg, deg, pc)
    (15.00021613, 14.99994433, 20.99999999)
 (pm_ra_cosdec, pm_dec, radial_velocity) in (mas / yr, mas / yr, km / s)
    (749.99980473, -200.00073221, 0.00029138)>
Now calling SC.apply_space_motion again...
Traceback (most recent call last):
  File "/home/boven/./SC-bug.py", line 24, in <module>
    SC2 = SC.apply_space_motion(new_obstime = T[2])
  File "/home/boven/venv/apy/lib/python3.10/site-packages/astropy/coordinates/sky_coordinate.py", line 811, in apply_space_motion
    icrsvel.d_lon.to_value(u.radian / u.yr),
AttributeError: 'SphericalCosLatDifferential' object has no attribute 'd_lon'. Did you mean: '_d_lon'?

The SC still clearly has a velocity (proper motion) in the longitudinal direction, despite the error message.

It’s a Heisenbug: if I remove the print statements for the SkyCoords, the code works!

To be more precise: it’s the print statement for the original SC, after calling ‘apply_space_motion’ that breaks things, the other print statements don’t seem to matter.

We suspect that the problem may be in line 1234 of baseframe.py in represent_as():

# Here we have to bypass using with_differentials() because
                    # it has a validation check. But because
                    # .representation_type and .differential_type don't point to
                    # the original classes, if the input differential is a
                    # RadialDifferential, it usually gets turned into a
                    # SphericalCosLatDifferential (or whatever the default is)
                    # with strange units for the d_lon and d_lat attributes.
                    # This then causes the dictionary key check to fail (i.e.
                    # comparison against `diff._get_deriv_key()`)
                    data._differentials.update({"s": diff})

The statement above gets called when you do a print() on a SkyCoord with differentials, and just goes and updates those.