Repository: astropy/astropy

Issue:
SkyCoord should support velocities
The velocity machinery added in v2.0 should get SkyCoord support.

cc @adrn

Discussion:
astropy-bot[bot]: Hi there @eteq :wave: - thanks for the pull request! I'm just a friendly :robot: that checks for issues related to the changelog and making sure that this pull request is milestoned and labeled correctly. This is mainly intended for the maintainers, so if you are not a maintainer you can ignore this, and a maintainer will let you know if any action is required on your part :smiley:.

Everything looks good from my point of view! :+1:

*If there are any issues with this message, please report them [here](https://github.com/astropy/astropy-bot/issues).*

Patch:
diff --git a/CHANGES.rst b/CHANGES.rst
index e54f01820e..f3cc6b3eab 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -26,6 +26,9 @@ astropy.coordinates
   a ``CartesianDifferential``, the 2D proper motion as a ``Quantity``, and the
   radial or line-of-sight velocity as a ``Quantity``. [#6869]
 
+- SkyCoord objects now support storing and tranforming differentials - i.e.,
+  both radial velocities and proper motions. [#6944]
+
 astropy.cosmology
 ^^^^^^^^^^^^^^^^^
 
@@ -302,6 +305,10 @@ astropy.convolution
 astropy.coordinates
 ^^^^^^^^^^^^^^^^^^^
 
+- Frame objects now use the default differential even if the representation is
+  explicitly provided as long as the representation provided is the same type as
+  the default representation. [#6944]
+
 astropy.cosmology
 ^^^^^^^^^^^^^^^^^
 
diff --git a/astropy/coordinates/baseframe.py b/astropy/coordinates/baseframe.py
index ebbf05a52e..2999e3bae1 100644
--- a/astropy/coordinates/baseframe.py
+++ b/astropy/coordinates/baseframe.py
@@ -61,6 +61,24 @@ def _get_repr_cls(value):
     return value
 
 
+def _get_diff_cls(value):
+    """
+    Return a valid differential class from ``value`` or raise exception.
+
+    As originally created, this is only used in the SkyCoord initializer, so if
+    that is refactored, this function my no longer be necessary.
+    """
+
+    if value in r.DIFFERENTIAL_CLASSES:
+        value = r.DIFFERENTIAL_CLASSES[value]
+    elif (not isinstance(value, type) or
+          not issubclass(value, r.BaseDifferential)):
+        raise ValueError(
+            'Differential is {0!r} but must be a BaseDifferential class '
+            'or one of the string aliases {1}'.format(
+                value, list(r.DIFFERENTIAL_CLASSES)))
+    return value
+
 def _get_repr_classes(base, **differentials):
     """Get valid representation and differential classes.
 
@@ -284,7 +302,10 @@ class BaseCoordinateFrame(ShapedLikeNDArray, metaclass=FrameMeta):
                 differential_cls = {'s': differential_cls}
 
             elif differential_cls is None:
-                differential_cls = {'s': 'base'} # see set_representation_cls()
+                if representation == self.default_representation:
+                    differential_cls = {'s': self.default_differential}
+                else:
+                    differential_cls = {'s': 'base'}  # see set_representation_cls()
 
             self.set_representation_cls(representation, **differential_cls)
 
diff --git a/astropy/coordinates/sky_coordinate.py b/astropy/coordinates/sky_coordinate.py
index 3e1e1910fc..4d05fbc192 100644
--- a/astropy/coordinates/sky_coordinate.py
+++ b/astropy/coordinates/sky_coordinate.py
@@ -17,7 +17,7 @@ from ..utils import ShapedLikeNDArray
 
 from .distances import Distance
 from .angles import Angle
-from .baseframe import BaseCoordinateFrame, frame_transform_graph, GenericFrame, _get_repr_cls
+from .baseframe import BaseCoordinateFrame, frame_transform_graph, GenericFrame, _get_repr_cls, _get_diff_cls
 from .builtin_frames import ICRS, SkyOffsetFrame
 from .representation import (BaseRepresentation, SphericalRepresentation,
                              UnitSphericalRepresentation)
@@ -146,6 +146,13 @@ class SkyCoord(ShapedLikeNDArray):
 
       >>> c = SkyCoord([ICRS(ra=1*u.deg, dec=2*u.deg), ICRS(ra=3*u.deg, dec=4*u.deg)])
 
+    Velocity components (proper motions or radial velocities) can also be
+    provided in a similar manner::
+
+      >>> c = SkyCoord(ra=1*u.deg, dec=2*u.deg, radial_velocity=10*u.km/u.s)
+
+      >>> c = SkyCoord(ra=1*u.deg, dec=2*u.deg, pm_ra_cosdec=2*u.mas/u.yr, pm_dec=1*u.mas/u.yr)
+
     As shown, the frame can be a `~astropy.coordinates.BaseCoordinateFrame`
     class or the corresponding string alias.  The frame classes that are built in
     to astropy are `ICRS`, `FK5`, `FK4`, `FK4NoETerms`, and `Galactic`.
@@ -181,14 +188,22 @@ class SkyCoord(ShapedLikeNDArray):
             RA and Dec for frames where ``ra`` and ``dec`` are keys in the
             frame's ``representation_component_names``, including `ICRS`,
             `FK5`, `FK4`, and `FK4NoETerms`.
+        pm_ra_cosdec, pm_dec  : `~astropy.units.Quantity`, optional
+            Proper motion components, in angle per time units.
         l, b : valid `~astropy.coordinates.Angle` initializer, optional
             Galactic ``l`` and ``b`` for for frames where ``l`` and ``b`` are
             keys in the frame's ``representation_component_names``, including
             the `Galactic` frame.
+        pm_l_cosb, pm_b : `~astropy.units.Quantity`, optional
+            Proper motion components in the `Galactic` frame, in angle per time
+            units.
         x, y, z : float or `~astropy.units.Quantity`, optional
             Cartesian coordinates values
         u, v, w : float or `~astropy.units.Quantity`, optional
             Cartesian coordinates values for the Galactic frame.
+        radial_velocity : `~astropy.units.Quantity`, optional
+            The component of the velocity along the line-of-sight (i.e., the
+            radial direction), in velocity units.
     """
 
     # Declare that SkyCoord can be used as a Table column by defining the
@@ -217,11 +232,15 @@ class SkyCoord(ShapedLikeNDArray):
                 setattr(self, attr, kwargs[attr])
 
         coord_kwargs = {}
+        component_names = frame.representation_component_names
+        component_names.update(frame.get_representation_component_names('s'))
         if 'representation' in kwargs:
             coord_kwargs['representation'] = _get_repr_cls(kwargs['representation'])
+        if 'differential_cls' in kwargs:
+            coord_kwargs['differential_cls'] = _get_diff_cls(kwargs['differential_cls'])
         for attr, value in kwargs.items():
-            if value is not None and (attr in frame.representation_component_names
-                                      or attr in frame.get_frame_attr_names()):
+            if value is not None and (attr in component_names
+                                      or attr in frame_attr_names):
                 coord_kwargs[attr] = value
 
         # Finally make the internal coordinate object.
@@ -321,15 +340,20 @@ class SkyCoord(ShapedLikeNDArray):
         # by keyword args or else get a None default.  Pop them off of kwargs
         # in the process.
         frame = valid_kwargs['frame'] = _get_frame(args, kwargs)
+        # TODO: possibly remove the below.  The representation/differential
+        # information should *already* be stored in the frame object, as it is
+        # extracted in _get_frame.  So it may be redundent to include it below.
         if 'representation' in kwargs:
             valid_kwargs['representation'] = _get_repr_cls(kwargs.pop('representation'))
+        if 'differential_cls' in kwargs:
+            valid_kwargs['differential_cls'] = _get_diff_cls(kwargs.pop('differential_cls'))
 
         for attr in frame_transform_graph.frame_attributes:
             if attr in kwargs:
                 valid_kwargs[attr] = kwargs.pop(attr)
 
         # Get units
-        units = _get_units(args, kwargs)
+        units = _get_representation_component_units(args, kwargs)
 
         # Grab any frame-specific attr names like `ra` or `l` or `distance` from kwargs
         # and migrate to valid_kwargs.
@@ -337,22 +361,14 @@ class SkyCoord(ShapedLikeNDArray):
 
         # Error if anything is still left in kwargs
         if kwargs:
-
-            # TODO: remove this when velocities are supported in SkyCoord
-            vel_url = 'http://docs.astropy.org/en/stable/coordinates/velocities.html'
-            for k in kwargs:
-                if k.startswith('pm_') or k == 'radial_velocity':
-                    raise ValueError('Velocity data is currently only supported'
-                                     ' in the coordinate frame objects, not in '
-                                     'SkyCoord. See the velocities '
-                                     'documentation page for more information: '
-                                     '{0}'.format(vel_url))
-
             raise ValueError('Unrecognized keyword argument(s) {0}'
                              .format(', '.join("'{0}'".format(key) for key in kwargs)))
 
-        # Finally deal with the unnamed args.  This figures out what the arg[0] is
-        # and returns a dict with appropriate key/values for initializing frame class.
+        # Finally deal with the unnamed args.  This figures out what the arg[0]
+        # is and returns a dict with appropriate key/values for initializing
+        # frame class. Note that differentials are *never* valid args, only
+        # kwargs.  So they are not accounted for here (unless they're in a frame
+        # or SkyCoord object)
         if args:
             if len(args) == 1:
                 # One arg which must be a coordinate.  In this case
@@ -1192,7 +1208,17 @@ class SkyCoord(ShapedLikeNDArray):
         """
         from .funcs import get_constellation
 
-        return get_constellation(self, short_name, constellation_list)
+        # because of issue #7028, the conversion to a PrecessedGeocentric
+        # system fails in some cases.  Work around is to  drop the velocities.
+        # they are not needed here since only position infromation is used
+        extra_frameattrs = {nm: getattr(self, nm)
+                            for nm in self._extra_frameattr_names}
+        novel = SkyCoord(self.realize_frame(self.data.without_differentials()),
+                         **extra_frameattrs)
+        return get_constellation(novel, short_name, constellation_list)
+
+        # the simpler version below can be used when gh-issue #7028 is resolved
+        #return get_constellation(self, short_name, constellation_list)
 
     # WCS pixel to/from sky conversions
     def to_pixel(self, wcs, origin=0, mode='all'):
@@ -1612,11 +1638,15 @@ def _get_frame(args, kwargs):
              len(args) == 1 and len(arg) > 0):
             arg = arg[0]
 
-        coord_frame_cls = None
+        coord_frame_obj = coord_frame_cls = None
         if isinstance(arg, BaseCoordinateFrame):
-            coord_frame_cls = arg.__class__
+            coord_frame_obj = arg
         elif isinstance(arg, SkyCoord):
-            coord_frame_cls = arg.frame.__class__
+            coord_frame_obj = arg.frame
+        if coord_frame_obj is not None:
+            coord_frame_cls = coord_frame_obj.__class__
+            kwargs.setdefault('differential_cls',
+                              coord_frame_obj.get_representation_cls('s'))
 
         if coord_frame_cls is not None:
             if not frame_specified_explicitly:
@@ -1626,18 +1656,19 @@ def _get_frame(args, kwargs):
                                  "new frame='{1}'.  Instead transform the coordinate."
                                  .format(coord_frame_cls.__name__, frame_cls.__name__))
 
+    frame_cls_kwargs = {}
     if 'representation' in kwargs:
-        frame = frame_cls(representation=_get_repr_cls(kwargs['representation']))
-    else:
-        frame = frame_cls()
+        frame_cls_kwargs['representation'] = _get_repr_cls(kwargs['representation'])
+    if 'differential_cls' in kwargs:
+        frame_cls_kwargs['differential_cls'] = _get_diff_cls(kwargs['differential_cls'])
 
-    return frame
+    return frame_cls(**frame_cls_kwargs)
 
 
-def _get_units(args, kwargs):
+def _get_representation_component_units(args, kwargs):
     """
-    Get the longitude unit and latitude unit from kwargs.  Possible enhancement
-    is to allow input from args as well.
+    Get the unit from kwargs for the *representation* components (not the
+    differentials).
     """
     if 'unit' not in kwargs:
         units = [None, None, None]
@@ -1682,9 +1713,9 @@ def _parse_coordinate_arg(coords, frame, units, init_kwargs):
     is_scalar = False  # Differentiate between scalar and list input
     valid_kwargs = {}  # Returned dict of lon, lat, and distance (optional)
 
-    frame_attr_names = frame.representation_component_names.keys()
-    repr_attr_names = frame.representation_component_names.values()
-    repr_attr_classes = frame.representation.attr_classes.values()
+    frame_attr_names = list(frame.representation_component_names.keys())
+    repr_attr_names = list(frame.representation_component_names.values())
+    repr_attr_classes = list(frame.representation.attr_classes.values())
     n_attr_names = len(repr_attr_names)
 
     # Turn a single string into a list of strings for convenience
@@ -1702,15 +1733,38 @@ def _parse_coordinate_arg(coords, frame, units, init_kwargs):
         data = coords.data.represent_as(frame.representation)
 
         values = []  # List of values corresponding to representation attrs
+        repr_attr_name_to_drop = []
         for repr_attr_name in repr_attr_names:
             # If coords did not have an explicit distance then don't include in initializers.
             if (isinstance(coords.data, UnitSphericalRepresentation) and
                     repr_attr_name == 'distance'):
+                repr_attr_name_to_drop.append(repr_attr_name)
                 continue
 
             # Get the value from `data` in the eventual representation
             values.append(getattr(data, repr_attr_name))
 
+        # drop the ones that were skipped because they were distances
+        for nametodrop in repr_attr_name_to_drop:
+            nameidx = repr_attr_names.index(nametodrop)
+            del repr_attr_names[nameidx]
+            del units[nameidx]
+            del frame_attr_names[nameidx]
+            del repr_attr_classes[nameidx]
+
+        if coords.data.differentials and 's' in coords.data.differentials:
+            orig_vel = coords.data.differentials['s']
+            vel = coords.data.represent_as(frame.representation, frame.get_representation_cls('s')).differentials['s']
+            for frname, reprname in frame.get_representation_component_names('s').items():
+                if (reprname == 'd_distance' and not hasattr(orig_vel, reprname) and
+                    'unit' in orig_vel.get_name()):
+                    continue
+                values.append(getattr(vel, reprname))
+                units.append(None)
+                frame_attr_names.append(frname)
+                repr_attr_names.append(reprname)
+                repr_attr_classes.append(vel.attr_classes[reprname])
+
         for attr in frame_transform_graph.frame_attributes:
             value = getattr(coords, attr, None)
             use_value = (isinstance(coords, SkyCoord)
@@ -1719,8 +1773,20 @@ def _parse_coordinate_arg(coords, frame, units, init_kwargs):
                 valid_kwargs[attr] = value
 
     elif isinstance(coords, BaseRepresentation):
-        data = coords.represent_as(frame.representation)
-        values = [getattr(data, repr_attr_name) for repr_attr_name in repr_attr_names]
+        if coords.differentials and 's' in coords.differentials:
+            diffs = frame.get_representation_cls('s')
+            data = coords.represent_as(frame.representation, diffs)
+            values = [getattr(data, repr_attr_name) for repr_attr_name in repr_attr_names]
+            for frname, reprname in frame.get_representation_component_names('s').items():
+                values.append(getattr(data.differentials['s'], reprname))
+                units.append(None)
+                frame_attr_names.append(frname)
+                repr_attr_names.append(reprname)
+                repr_attr_classes.append(data.differentials['s'].attr_classes[reprname])
+
+        else:
+            data = coords.represent_as(frame.representation)
+            values = [getattr(data, repr_attr_name) for repr_attr_name in repr_attr_names]
 
     elif (isinstance(coords, np.ndarray) and coords.dtype.kind in 'if'
           and coords.ndim == 2 and coords.shape[1] <= 3):
@@ -1839,6 +1905,9 @@ def _get_representation_attrs(frame, units, kwargs):
     for the underlying data values in the representation, e.g. "ra" for "lon"
     for many equatorial spherical representations, or "w" for "x" in the
     cartesian representation of Galactic.
+
+    This also gets any *differential* kwargs, because they go into the same
+    frame initializer later on.
     """
     frame_attr_names = frame.representation_component_names.keys()
     repr_attr_classes = frame.representation.attr_classes.values()
@@ -1849,6 +1918,17 @@ def _get_representation_attrs(frame, units, kwargs):
         if value is not None:
             valid_kwargs[frame_attr_name] = repr_attr_class(value, unit=unit)
 
+    # also check the differentials.  They aren't included in the units keyword,
+    # so we only look for the names.
+
+    # TODO: change _representation['s'] to `differential_type` when that exists
+    differential_type = frame._representation['s'].attr_classes
+    for frame_name, repr_name in frame.get_representation_component_names('s').items():
+        diff_attr_class = differential_type[repr_name]
+        value = kwargs.pop(frame_name, None)
+        if value is not None:
+            valid_kwargs[frame_name] = diff_attr_class(value)
+
     return valid_kwargs
 
 
diff --git a/docs/coordinates/index.rst b/docs/coordinates/index.rst
index 83b86b53dd..f1ce767859 100644
--- a/docs/coordinates/index.rst
+++ b/docs/coordinates/index.rst
@@ -291,10 +291,21 @@ high-level |skycoord| method - see :ref:`astropy-coordinates-rv-corrs`)::
 Velocities (Proper Motions and Radial Velocities)
 -------------------------------------------------
 
-New in Astropy v2.0, the :doc:`coordinate frame classes <frames>` can now store
-and transform velocities along with positional coordinate information. This
-is not available from the |skycoord| class as this new functionality is
-experimental, but is accessible  for more information see the :doc:`velocities` page.
+In addition to positional coordinates, `~astropy.coordinates` supports storing
+and transforming velocities.  These are available both via the lower-level
+:doc:`coordinate frame classes <frames>`, and (new in v3.0) via  |skycoord|
+objects::
+
+    >>> sc = SkyCoord(1*u.deg, 2*u.deg, radial_velocity=20*u.km/u.s)
+    >>> sc  # doctest: +SKIP
+    <SkyCoord (ICRS): (ra, dec) in deg
+        ( 1.,  2.)
+     (radial_velocity) in km / s
+        ( 20.,)>
+
+.. the SKIP above in the ``sc`` line is because numpy has a subtly different output in versions < 12 - the trailing comma is missing.  If a NPY_LT_1_12 comes in to being this can switch to that.  But don't forget to *also* change this in the velocities.rst file
+
+For more details on velocity support (and limitations), see the :doc:`velocities` page.
 
 .. _astropy-coordinates-overview:
 
diff --git a/docs/coordinates/velocities.rst b/docs/coordinates/velocities.rst
index d773828e8a..331e6ea31c 100644
--- a/docs/coordinates/velocities.rst
+++ b/docs/coordinates/velocities.rst
@@ -5,19 +5,52 @@
 Working with velocities in Astropy coordinates
 **********************************************
 
-.. warning::
-    Velocities support, new in Astropy v2.0, is an experimental feature and is
-    subject to change based on user feedback.  While we do not expect major API
-    changes, the possibility exists based on the precedent of earlier changes
-    in the ``coordinates`` subpackage based on user feedback from previous
-    versions of Astropy.
+Using velocities with ``SkyCoord``
+==================================
+
+The easiest way of getting a coordinate object with velocities is to use the
+|skycoord| interface.  For example, a |skycoord| to represent a star with a
+measured radial velocity but unknown proper motion and distance could be
+created as::
+
+    >>> from astropy.coordinates import SkyCoord
+    >>> import astropy.units as u
+    >>> sc = SkyCoord(1*u.deg, 2*u.deg, radial_velocity=20*u.km/u.s)
+    >>> sc  # doctest: +SKIP
+    <SkyCoord (ICRS): (ra, dec) in deg
+        ( 1.,  2.)
+     (radial_velocity) in km / s
+        ( 20.,)>
+    >>> sc.radial_velocity  # doctest: +FLOAT_CMP
+    <Quantity 20.0 km / s>
+
+.. the SKIP above in the ``sc`` line is because numpy has a subtly different output in versions < 12 - the trailing comma is missing.  If a NPY_LT_1_12 comes in to being this can switch to that.  But don't forget to *also* change this in the coordinates/index.rst file
+
+|skycoord| objects created in this manner follow all the same transformation
+rules, and will correctly update their velocities when transformed to other
+frames.  For example, to determine proper motions in Galactic coordinates for
+a star with proper motions measured in ICRS::
+
+    >>> sc = SkyCoord(1*u.deg, 2*u.deg, pm_ra_cosdec=.2*u.mas/u.yr, pm_dec=.1*u.mas/u.yr)
+    >>> sc.galactic  # doctest: +FLOAT_CMP
+    <SkyCoord (Galactic): (l, b) in deg
+      ( 99.63785528, -58.70969293)
+    (pm_l_cosb, pm_b) in mas / yr
+      ( 0.22240398,  0.02316181)>
+
+For more details on valid operations and limitations of velocity support in
+`astropy.coordinates` (particularly the :ref:`current accuracy limitations
+<astropy-coordinate-finite-difference-velocities>` ), see the more detailed
+discussions below of velocity support in the lower-level frame objects.  All
+these same rules apply for |skycoord| objects, as they are built directly on top
+of the frame classes' velocity functionality detailed here.
 
 .. _astropy-coordinate-custom-frame-with-velocities:
 
 Creating frame objects with velocity data
 =========================================
 
-The coordinate frame classes now support storing and transforming velocity data
+The coordinate frame classes support storing and transforming velocity data
 (along side the positional coordinate data). Similar to the positional data ---
 that use the ``Representation`` classes to abstract away the particular
 representation and allow re-representing from, e.g., Cartesian to Spherical
@@ -36,7 +69,6 @@ term. For example, the proper motion components for the ``ICRS`` frame are
 (``pm_ra_cosdec``, ``pm_dec``)::
 
     >>> from astropy.coordinates import ICRS
-    >>> import astropy.units as u
     >>> ICRS(ra=8.67*u.degree, dec=53.09*u.degree,
     ...      pm_ra_cosdec=4.8*u.mas/u.yr, pm_dec=-15.16*u.mas/u.yr)  # doctest: +FLOAT_CMP
     <ICRS Coordinate: (ra, dec) in deg
@@ -198,6 +230,8 @@ for example, `~astropy.coordinates.ICRS` to `~astropy.coordinates.LSR`::
      (pm_ra_cosdec, pm_dec, radial_velocity) in (mas / yr, mas / yr, km / s)
         (-24.51315607, -2.67935501, 27.07339176)>
 
+.. _astropy-coordinate-finite-difference-velocities:
+
 Finite Difference Transformations
 ---------------------------------
 
@@ -300,15 +334,6 @@ a particular direction changes dramatically over the course of one year).
 Future versions of Astropy will improve on this algorithm to make the results
 more numerically stable and practical for use in these (not unusual) use cases.
 
-
-``SkyCoord`` support for Velocities
-===================================
-
-|skycoord| currently does *not* support velocities as of Astropy v2.0.  This is
-an intentional choice, allowing the "power-user" community to provide feedback
-on the API and functionality in the frame-level classes before it is adopted in
-|skycoord| (currently planned for the next Astropy version, v3.0).
-
 .. _astropy-coordinates-rv-corrs:
 
 Radial Velocity Corrections


Test patch:
diff --git a/astropy/coordinates/tests/test_sky_coord_velocities.py b/astropy/coordinates/tests/test_sky_coord_velocities.py
new file mode 100644
index 0000000000..601d7ee35c
--- /dev/null
+++ b/astropy/coordinates/tests/test_sky_coord_velocities.py
@@ -0,0 +1,205 @@
+# -*- coding: utf-8 -*-
+# Licensed under a 3-clause BSD style license - see LICENSE.rst
+
+"""
+Tests for putting velocity differentials into SkyCoord objects.
+"""
+
+import pytest
+
+import numpy as np
+
+from ... import units as u
+from ...tests.helper import assert_quantity_allclose
+from .. import (SkyCoord, ICRS, SphericalRepresentation, SphericalDifferential,
+                SphericalCosLatDifferential, CartesianRepresentation,
+                CartesianDifferential, Galactic, PrecessedGeocentric)
+
+try:
+    import scipy
+    HAS_SCIPY = True
+except ImportError:
+    HAS_SCIPY = False
+
+def test_creation_frameobjs():
+    i = ICRS(1*u.deg, 2*u.deg, pm_ra_cosdec=.2*u.mas/u.yr, pm_dec=.1*u.mas/u.yr)
+    sc = SkyCoord(i)
+
+    for attrnm in ['ra', 'dec', 'pm_ra_cosdec', 'pm_dec']:
+        assert_quantity_allclose(getattr(i, attrnm), getattr(sc, attrnm))
+
+    sc_nod = SkyCoord(ICRS(1*u.deg, 2*u.deg))
+
+    for attrnm in ['ra', 'dec']:
+        assert_quantity_allclose(getattr(sc, attrnm), getattr(sc_nod, attrnm))
+
+
+def test_creation_attrs():
+    sc1 = SkyCoord(1*u.deg, 2*u.deg,
+                   pm_ra_cosdec=.2*u.mas/u.yr, pm_dec=.1*u.mas/u.yr,
+                   frame='fk5')
+    assert_quantity_allclose(sc1.ra, 1*u.deg)
+    assert_quantity_allclose(sc1.dec, 2*u.deg)
+    assert_quantity_allclose(sc1.pm_ra_cosdec, .2*u.arcsec/u.kyr)
+    assert_quantity_allclose(sc1.pm_dec, .1*u.arcsec/u.kyr)
+
+    sc2 = SkyCoord(1*u.deg, 2*u.deg,
+                   pm_ra=.2*u.mas/u.yr, pm_dec=.1*u.mas/u.yr,
+                   differential_cls=SphericalDifferential)
+    assert_quantity_allclose(sc2.ra, 1*u.deg)
+    assert_quantity_allclose(sc2.dec, 2*u.deg)
+    assert_quantity_allclose(sc2.pm_ra, .2*u.arcsec/u.kyr)
+    assert_quantity_allclose(sc2.pm_dec, .1*u.arcsec/u.kyr)
+
+    sc3 = SkyCoord('1:2:3 4:5:6',
+                   pm_ra_cosdec=.2*u.mas/u.yr, pm_dec=.1*u.mas/u.yr,
+                   unit=(u.hour, u.deg))
+
+    assert_quantity_allclose(sc3.ra, 1*u.hourangle + 2*u.arcmin*15 + 3*u.arcsec*15)
+    assert_quantity_allclose(sc3.dec, 4*u.deg + 5*u.arcmin + 6*u.arcsec)
+    # might as well check with sillier units?
+    assert_quantity_allclose(sc3.pm_ra_cosdec, 1.2776637006616473e-07 * u.arcmin / u.fortnight)
+    assert_quantity_allclose(sc3.pm_dec, 6.388318503308237e-08 * u.arcmin / u.fortnight)
+
+
+def test_creation_copy_basic():
+    i = ICRS(1*u.deg, 2*u.deg, pm_ra_cosdec=.2*u.mas/u.yr, pm_dec=.1*u.mas/u.yr)
+    sc = SkyCoord(i)
+    sc_cpy = SkyCoord(sc)
+
+    for attrnm in ['ra', 'dec', 'pm_ra_cosdec', 'pm_dec']:
+        assert_quantity_allclose(getattr(sc, attrnm), getattr(sc_cpy, attrnm))
+
+
+def test_creation_copy_rediff():
+    sc = SkyCoord(1*u.deg, 2*u.deg,
+                  pm_ra=.2*u.mas/u.yr, pm_dec=.1*u.mas/u.yr,
+                  differential_cls=SphericalDifferential)
+
+    sc_cpy = SkyCoord(sc)
+    for attrnm in ['ra', 'dec', 'pm_ra', 'pm_dec']:
+        assert_quantity_allclose(getattr(sc, attrnm), getattr(sc_cpy, attrnm))
+
+    sc_newdiff = SkyCoord(sc, differential_cls=SphericalCosLatDifferential)
+    reprepr = sc.represent_as(SphericalRepresentation, SphericalCosLatDifferential)
+    assert_quantity_allclose(sc_newdiff.pm_ra_cosdec,
+                             reprepr.differentials['s'].d_lon_coslat)
+
+
+def test_creation_cartesian():
+    rep = CartesianRepresentation([10, 0., 0.]*u.pc)
+    dif = CartesianDifferential([0, 100, 0.]*u.pc/u.Myr)
+    rep = rep.with_differentials(dif)
+    c = SkyCoord(rep)
+
+    sdif = dif.represent_as(SphericalCosLatDifferential, rep)
+    assert_quantity_allclose(c.pm_ra_cosdec, sdif.d_lon_coslat)
+
+
+def test_useful_error_missing():
+    sc_nod = SkyCoord(ICRS(1*u.deg, 2*u.deg))
+    try:
+        sc_nod.l
+    except AttributeError as e:
+        # this is double-checking the *normal* behavior
+        msg_l = e.args[0]
+
+    try:
+        sc_nod.pm_dec
+    except Exception as e:
+        msg_pm_dec = e.args[0]
+
+    assert "has no attribute" in msg_l
+    assert "has no associated differentials" in msg_pm_dec
+
+
+# ----------------------Operations on SkyCoords w/ velocities-------------------
+
+# define some fixtires to get baseline coordinates to try operations with
+@pytest.fixture(scope="module", params=[(False, False),
+                                        (True, False),
+                                        (False, True),
+                                        (True, True)])
+def sc(request):
+    incldist, inclrv = request.param
+
+    args = [1*u.deg, 2*u.deg]
+    kwargs = dict(pm_dec=1*u.mas/u.yr, pm_ra_cosdec=2*u.mas/u.yr)
+    if incldist:
+        kwargs['distance'] = 213.4*u.pc
+    if inclrv:
+        kwargs['radial_velocity'] = 61*u.km/u.s
+
+    return SkyCoord(*args, **kwargs)
+
+@pytest.fixture(scope="module")
+def scmany():
+    return SkyCoord(ICRS(ra=[1]*100*u.deg, dec=[2]*100*u.deg,
+                     pm_ra_cosdec=np.random.randn(100)*u.mas/u.yr,
+                     pm_dec=np.random.randn(100)*u.mas/u.yr,))
+
+@pytest.fixture(scope="module")
+def sc_for_sep():
+    return SkyCoord(1*u.deg, 2*u.deg,
+                    pm_dec=1*u.mas/u.yr, pm_ra_cosdec=2*u.mas/u.yr)
+
+
+def test_separation(sc, sc_for_sep):
+    sc.separation(sc_for_sep)
+
+
+def test_accessors(sc, scmany):
+    sc.data.differentials['s']
+    sph = sc.spherical
+    gal = sc.galactic
+
+    if (sc.data.get_name().startswith('unit') and not
+        sc.data.differentials['s'].get_name().startswith('unit')):
+        # this xfail can be eliminated when issue #7028 is resolved
+        pytest.xfail('.velocity fails if there is an RV but not distance')
+    sc.velocity
+
+    assert isinstance(sph, SphericalRepresentation)
+    assert gal.data.differentials is not None
+
+    scmany[0]
+    sph = scmany.spherical
+    gal = scmany.galactic
+
+    assert isinstance(sph, SphericalRepresentation)
+    assert gal.data.differentials is not None
+
+def test_transforms(sc):
+    trans = sc.transform_to('galactic')
+    assert isinstance(trans.frame, Galactic)
+
+
+@pytest.mark.xfail
+def test_transforms_diff(sc):
+    # note that arguably this *should* fail for the no-distance cases: 3D
+    # information is necessary to truly solve this, hence the xfail
+    trans = sc.transform_to(PrecessedGeocentric(equinox='B1975'))
+    assert isinstance(trans.frame, PrecessedGeocentric)
+
+
+@pytest.mark.skipif(str('not HAS_SCIPY'))
+def test_matching(sc, scmany):
+    # just check that it works and yields something
+    idx, d2d, d3d = sc.match_to_catalog_sky(scmany)
+
+
+def test_position_angle(sc, sc_for_sep):
+    sc.position_angle(sc_for_sep)
+
+
+def test_constellations(sc):
+    const = sc.get_constellation()
+    assert const == 'Pisces'
+
+
+# at the time of this writing, sky offset frames don't support velocities
+# TODO: make them actually work
+@pytest.mark.xfail
+def test_skyoffset_frame(sc):
+
+    sc.skyoffset_frame()