Add differential classes for angular systems to represent velocities in physical units - #20475

Thank you for your contribution to Astropy! 🌌 This checklist is meant to remind the package maintainers who will review this pull request of some common things to look for.

• Do the proposed changes actually accomplish desired goals?
• Do the proposed changes follow the Astropy coding guidelines?
• Are tests added/updated as required? If so, do they follow the Astropy testing guidelines?
• Are docs added/updated as required? If so, do they follow the Astropy documentation guidelines?
• Is rebase and/or squash necessary? If so, please provide the author with appropriate instructions. Also see instructions for rebase and squash.
• Did the CI pass? If no, are the failures related? If you need to run daily and weekly cron jobs as part of the PR, please apply the "Extra CI" label. Codestyle issues can be fixed by the bot.
• Is a change log needed? If yes, did the change log check pass? If no, add the "no-changelog-entry-needed" label. If this is a manual backport, use the "skip-changelog-checks" label unless special changelog handling is necessary.
• Is this a big PR that makes a "What's new?" entry worthwhile and if so, is (1) a "what's new" entry included in this PR and (2) the "whatsnew-needed" label applied?
• At the time of adding the milestone, if the milestone set requires a backport to release branch(es), apply the appropriate "backport-X.Y.x" label(s) before merge.

Thanks! When you think you are ready for the next step, might wanna add "Extra CI" label to make sure the tests won't fail in a more exotic setting, just in case.

Adding a RadialDifferential crashe

3. got cut off - can you try posting that point again?

EDIT: Resolved in comment above.

Oooh, this is nice to have. From just messing around I notice the output units can get slightly unusual when some components of the differentials are missing, e.g this looks great:

In [38]: f = ICRS(ra=1*u.deg, dec=2*u.deg, pm_ra_cosdec=1*u.mas/u.yr, pm_dec=2*u.mas/u.yr, distance=1*u.pc, radial_velocity=
       ⋮ 0*u.km/u.s)

In [39]: f.represent_as("spherical", s="sphericalphysical", in_frame_units=True).differentials["s"]
Out[39]: 
<SphericalPhysicalDifferential (d_lon, d_lat, d_distance) in km / s
    (0.00474047, 0.00948094, -2.05095566e-20)>

whereas this reports the same physical velocities along all the unit vectors, but with unusual units for some, but not all...

In [40]: f = ICRS(ra=1*u.deg, dec=2*u.deg, pm_ra_cosdec=1*u.mas/u.yr, pm_dec=2*u.mas/u.yr, distance=1*u.pc)

In [41]: f.represent_as("spherical", s="sphericalphysical", in_frame_units=True).differentials["s"]
Out[41]: 
<SphericalPhysicalDifferential (d_lon, d_lat, d_distance) in (mas pc / (rad yr), km / s, mas pc / (rad yr))
    (1., 0.00948094, -4.32648126e-18)>

the output units can get slightly unusual when some components of the differentials are missing

Thanks @StuartLittlefair - yes that's related to @nstarman 's point 2 above. So then the question is whether to error or warn in that situation. Thoughts?

Related but subtly different. In the missing distance case I’d expect it to error, but in the case of missing RVs the d_lon and d_lat attributes don’t depend on the RV and the values are right, so throwing an error seems wrong - I just expect different units.

I guess we could check if the missing attributes affect the scales and error if so?

Addressing the comments now. While I was fixing @nstarman's comment 2. above, I found another bug:

import astropy.units as u
import astropy.coordinates as coord
c = coord.SkyCoord(
    ra=1 * u.deg, dec=2 * u.deg, pm_ra_cosdec=1 * u.mas / u.yr, pm_dec=2 * u.mas / u.yr
)
c.set_representation_cls(s="spherical")
print(repr(c))

returns

<SkyCoord (ICRS): (ra, dec) in deg
    (1., 2.)
 (pm_ra_coslat, pm_dec) in mas / yr
    (1., 2.)>

Note the "pm_ra_coslat" instead of "pm_ra_cosdec"! I think this is a separate issue on main, but also affects this. I'll put a commit in here to fix it but fine if you think I should move it to a separate PR.

Adding a RadialDifferential crashes.

@nstarman - that's not unique to these classes. For example, this fails with that error on main:

import astropy.units as u
from astropy.coordinates import CylindricalDifferential, RadialDifferential
CylindricalDifferential(1 * u.km / u.s, 1 * u.mas / u.yr, 1 * u.km / u.s) + RadialDifferential(1 * u.km / u.s)

I guess we'd rather raise a more informative error, or maybe NotImplementedError? Can you make a separate issue to track that?

Cartesian round-trip for the common case. SphericalCosLat to SphericalPhysical is just "multiply by distance", but it goes through Cartesian and back, which is where the ~1e-17 residuals come from. A shortcut in from_representation would be cheaper and exact.

@nstarman I think the small 1e-17 is unfortunate but consistent with other re-representations in coordinates. Adding a shortcut, I think, would also mean special-casing the physical differentials in BaseSphericalDifferential and BaseSphericalCosLatDifferential, plus adding the custom from_representation. Do you think it's worth doing that here?

Those 1e-17 residuals cause issues elsewhere as well, particularly for differentials attached to representations near the poles, e.g the non-physical motion in ra_cosdec here:

b = SphericalRepresentation(0*u.deg, -90*u.deg, 1*u.kpc)
bd = b.with_differentials(SphericalDifferential(1*u.deg/u.s, 0*u.deg/u.s, 0*u.km/u.s))
ic = ICRS(bd)

<ICRS Coordinate: (ra, dec, distance) in (deg, deg, kpc)
    (0., -90., 1.)
 (pm_ra_cosdec, pm_dec, radial_velocity) in (mas / yr, mas / yr, km / s)
    (0.00695644, 0., 0.)>

which (I think) arises as:

In [19]: bd.differentials['s']._get_base_vectors(b)
Out[19]: 
({'lon': <CartesianRepresentation (x, y, z) [dimensionless]
      (-0., 1., 0.)>,
  'lat': <CartesianRepresentation (x, y, z) [dimensionless]
      (1., 0., 6.123234e-17)>,
  'distance': <CartesianRepresentation (x, y, z) [dimensionless]
      (6.123234e-17, 0., -1.)>},
 {'lon': <Quantity 6.123234e-17 kpc / rad>,
  'lat': <Quantity 1. kpc / rad>,
  'distance': <Quantity 1.>})

So I'd prefer raising this as a separate issue where we potentially special case a number of transforms

Thanks, yes I agree. I think other fixes have been implemented btw.