Various modifications, mostly as requested during review.

astropy/coordinates/builtin_frames/__init__.py

@@ -48,6 +48,7 @@
 from . import icrs_cirs_transforms
 from . import cirs_observed_transforms
 from . import icrs_observed_transforms
+from . import itrs_observed_transforms
 from . import intermediate_rotation_transforms
 from . import ecliptic_transforms
 

astropy/coordinates/builtin_frames/itrs.py

@@ -45,12 +45,13 @@ class ITRS(BaseCoordinateFrame):
     position to a topocentric ITRS frame and add the observing site's EarthLocation geocentric
     ITRS coordinates to yield the object's geocentric ITRS coordinates.
 
-    On the other hand, direct transforms between geocentric ITRS coordinates and observed
-    `~astropy.coordinates.AltAz` or `~astropy.coordinates.HADec` coordinates include the
-    difference between stellar aberration from the point of view of an observer at the
-    geocenter and stellar aberration from the point of view of an observer on the surface
-    of the Earth. If the geocentric ITRS coordinates of the object include stellar aberration
-    at the geocenter, then this is the way to go.
+    On the other hand, using ``transform_to`` to transform geocentric ITRS coordinates to
+    topocentric ITRS, observed `~astropy.coordinates.AltAz`, or observed
+    `~astropy.coordinates.HADec` coordinates includes the difference between stellar aberration
+    from the point of view of an observer at the geocenter and stellar aberration from the
+    point of view of an observer on the surface of the Earth. If the geocentric ITRS
+    coordinates of the object include stellar aberration at the geocenter (e.g. certain ILRS
+    ephemerides), then this is the way to go.
 
     Note to ILRS ephemeris users: Astropy does not currently consider relativistic
     effects of the Earth's gravatational field. Nor do the `~astropy.coordinates.AltAz`

astropy/coordinates/builtin_frames/itrs_to_observed_transforms.py → astropy/coordinates/builtin_frames/itrs_observed_transforms.py

@@ -8,45 +8,42 @@
 from .altaz import AltAz
 from .hadec import HADec
 from .itrs import ITRS
-from .utils import PIOVER2
 
 # Minimum cos(alt) and sin(alt) for refraction purposes
 CELMIN = 1e-6
 SELMIN = 0.05
+# Latitude of the north pole.
+NORTH_POLE = 90.0*u.deg
 
 
 def itrs_to_altaz_mat(lon, lat):
     # form ITRS to AltAz matrix
-    elat = lat.to_value(u.radian)
-    elong = lon.to_value(u.radian)
     # AltAz frame is left handed
     minus_x = np.eye(3)
     minus_x[0][0] = -1.0
     mat = (minus_x
-           @ rotation_matrix(PIOVER2 - elat, 'y', unit=u.radian)
-           @ rotation_matrix(elong, 'z', unit=u.radian))
+           @ rotation_matrix(NORTH_POLE - lat, 'y')
+           @ rotation_matrix(lon, 'z'))
     return mat
 
 
 def itrs_to_hadec_mat(lon):
     # form ITRS to HADec matrix
-    elong = lon.to_value(u.radian)
     # HADec frame is left handed
     minus_y = np.eye(3)
     minus_y[1][1] = -1.0
     mat = (minus_y
-           @ rotation_matrix(elong, 'z', unit=u.radian))
+           @ rotation_matrix(lon, 'z'))
     return mat
 
 
 def altaz_to_hadec_mat(lat):
     # form AltAz to HADec matrix
-    elat = lat.to_value(u.radian)
     z180 = np.eye(3)
     z180[0][0] = -1.0
     z180[1][1] = -1.0
     mat = (z180
-           @ rotation_matrix(PIOVER2 - elat, 'y', unit=u.radian))
+           @ rotation_matrix(NORTH_POLE - lat, 'y'))
     return mat
 
 
@@ -58,7 +55,7 @@ def add_refraction(aa_crepr, observed_frame):
         observed_frame.relative_humidity.value,
         observed_frame.obswl.to_value(u.micron)
     )
-    #reference: erfa.atioq()
+    # reference: erfa.atioq()
     norm, uv = erfa.pn(aa_crepr.get_xyz(xyz_axis=-1).to_value())
     # Cosine and sine of altitude, with precautions.
     sel = np.maximum(uv[..., 2], SELMIN)
@@ -108,6 +105,7 @@ def remove_refraction(aa_crepr, observed_frame):
 def itrs_to_observed(itrs_coo, observed_frame):
     if (np.any(itrs_coo.location != observed_frame.location) or
             np.any(itrs_coo.obstime != observed_frame.obstime)):
+        # This transform will go through the CIRS and alter stellar aberration.
         itrs_coo = itrs_coo.transform_to(ITRS(obstime=observed_frame.obstime,
                                               location=observed_frame.location))
 
@@ -142,5 +140,6 @@ def observed_to_itrs(observed_coo, itrs_frame):
 
     itrs_at_obs_time = ITRS(crepr, obstime=observed_coo.obstime,
                             location=observed_coo.location)
-    # this final transform may be a no-op if the obstimes and locations are the same
+    # This final transform may be a no-op if the obstimes and locations are the same.
+    # Otherwise, this transform will go through the CIRS and alter stellar aberration.
     return itrs_at_obs_time.transform_to(itrs_frame)

astropy/coordinates/tests/test_intermediate_transformations.py

@@ -203,37 +203,44 @@ def test_itrs_topo_to_altaz_with_refraction():
     altaz_frame1 = AltAz(obstime = 'J2000', location=loc)
     altaz_frame2 = AltAz(obstime = 'J2000', location=loc, pressure=1000.0 * u.hPa,
                          relative_humidity=0.5)
+    cirs_frame = CIRS(obstime = 'J2000', location=loc)
     itrs_frame = ITRS(location=loc)
 
-    #Normal route
-    #No Refraction
+    # Normal route
+    # No Refraction
     altaz1 = icrs.transform_to(altaz_frame1)
 
-    #Refraction added
+    # Refraction added
     altaz2 = icrs.transform_to(altaz_frame2)
 
-    #Refraction removed
-    cirs_frame = CIRS(obstime = 'J2000', location=loc)
+    # Refraction removed
     cirs = altaz2.transform_to(cirs_frame)
     altaz3 = cirs.transform_to(altaz_frame1)
 
-    #Through ITRS
-    #No Refraction
+    # Through ITRS
+    # No Refraction
     itrs = icrs.transform_to(itrs_frame)
     altaz11 = itrs.transform_to(altaz_frame1)
 
     assert_allclose(altaz11.az - altaz1.az, 0*u.mas, atol=0.1*u.mas)
     assert_allclose(altaz11.alt - altaz1.alt, 0*u.mas, atol=0.1*u.mas)
     assert_allclose(altaz11.distance - altaz1.distance, 0*u.cm, atol=10.0*u.cm)
 
-    #Refraction added
+    # Round trip
+    itrs11 = altaz11.transform_to(itrs_frame)
+
+    assert_allclose(itrs11.x, itrs.x)
+    assert_allclose(itrs11.y, itrs.y)
+    assert_allclose(itrs11.z, itrs.z)
+
+    # Refraction added
     altaz22 = itrs.transform_to(altaz_frame2)
 
     assert_allclose(altaz22.az - altaz2.az, 0*u.mas, atol=0.1*u.mas)
     assert_allclose(altaz22.alt - altaz2.alt, 0*u.mas, atol=0.1*u.mas)
     assert_allclose(altaz22.distance - altaz2.distance, 0*u.cm, atol=10.0*u.cm)
 
-    #Refraction removed
+    # Refraction removed
     itrs = altaz22.transform_to(itrs_frame)
     altaz33 = itrs.transform_to(altaz_frame1)
 
@@ -251,44 +258,52 @@ def test_itrs_topo_to_hadec_with_refraction():
     hadec_frame1 = HADec(obstime = 'J2000', location=loc)
     hadec_frame2 = HADec(obstime = 'J2000', location=loc, pressure=1000.0 * u.hPa,
                          relative_humidity=0.5)
+    cirs_frame = CIRS(obstime = 'J2000', location=loc)
     itrs_frame = ITRS(location=loc)
 
-    #Normal route
-    #No Refraction
+    # Normal route
+    # No Refraction
     hadec1 = icrs.transform_to(hadec_frame1)
 
-    #Refraction added
+    # Refraction added
     hadec2 = icrs.transform_to(hadec_frame2)
 
-    #Refraction removed
-    cirs_frame = CIRS(obstime = 'J2000', location=loc)
+    # Refraction removed
     cirs = hadec2.transform_to(cirs_frame)
     hadec3 = cirs.transform_to(hadec_frame1)
 
-    #Through ITRS
-    #No Refraction
+    # Through ITRS
+    # No Refraction
     itrs = icrs.transform_to(itrs_frame)
     hadec11 = itrs.transform_to(hadec_frame1)
 
     assert_allclose(hadec11.ha - hadec1.ha, 0*u.mas, atol=0.1*u.mas)
     assert_allclose(hadec11.dec - hadec1.dec, 0*u.mas, atol=0.1*u.mas)
     assert_allclose(hadec11.distance - hadec1.distance, 0*u.cm, atol=10.0*u.cm)
 
-    #Refraction added
+    # Round trip
+    itrs11 = hadec11.transform_to(itrs_frame)
+
+    assert_allclose(itrs11.x, itrs.x)
+    assert_allclose(itrs11.y, itrs.y)
+    assert_allclose(itrs11.z, itrs.z)
+
+    # Refraction added
     hadec22 = itrs.transform_to(hadec_frame2)
 
     assert_allclose(hadec22.ha - hadec2.ha, 0*u.mas, atol=0.1*u.mas)
     assert_allclose(hadec22.dec - hadec2.dec, 0*u.mas, atol=0.1*u.mas)
     assert_allclose(hadec22.distance - hadec2.distance, 0*u.cm, atol=10.0*u.cm)
 
-    #Refraction removed
+    # Refraction removed
     itrs = hadec22.transform_to(itrs_frame)
     hadec33 = itrs.transform_to(hadec_frame1)
 
     assert_allclose(hadec33.ha - hadec3.ha, 0*u.mas, atol=0.1*u.mas)
     assert_allclose(hadec33.dec - hadec3.dec, 0*u.mas, atol=0.1*u.mas)
     assert_allclose(hadec33.distance - hadec3.distance, 0*u.cm, atol=10.0*u.cm)
 
+
 def test_gcrs_itrs():
     """
     Check basic GCRS<->ITRS transforms for round-tripping.
@@ -939,8 +954,7 @@ def test_itrs_straight_overhead():
     itrs_repr = itrs_geo - obsrepr
 
     # create a ITRS object that appears straight overhead for a TOPOCENTRIC OBSERVER
-    topocentric_itrs_frame = ITRS(obstime=t, location=home)
-    itrs_topo = topocentric_itrs_frame.realize_frame(itrs_repr)
+    itrs_topo = ITRS(itrs_repr, obstime=t, location=home)
 
     # Check AltAz (though Azimuth can be anything so is not tested).
     aa = itrs_topo.transform_to(AltAz(obstime=t, location=home))
@@ -951,6 +965,7 @@ def test_itrs_straight_overhead():
     assert_allclose(hd.ha, 0*u.hourangle, atol=1*u.uas, rtol=0)
     assert_allclose(hd.dec, 52*u.deg, atol=1*u.uas, rtol=0)
 
+
 def jplephem_ge(minversion):
     """Check if jplephem is installed and has version >= minversion."""
     # This is a separate routine since somehow with pyinstaller the stanza

docs/changes/coordinates/13398.feature.rst

@@ -1 +1,6 @@
-Adds new topocentric ITRS frame and direct transforms to and from the observed frames ``AltAz`` and ``HADec`` with the ability to add or remove refraction corrections as required. Since these frames are all within the ITRS, there are no corrections applied other than refraction in the transforms. This makes the topocentric ITRS frame and these transforms convenient for observers of near Earth objects where stellar aberration should be omitted.
+Adds new topocentric ITRS frame and direct transforms to and from the observed
+frames ``AltAz`` and ``HADec`` with the ability to add or remove refraction
+corrections as required. Since these frames are all within the ITRS, there are
+no corrections applied other than refraction in the transforms. This makes the
+topocentric ITRS frame and these transforms convenient for observers of near
+Earth objects where stellar aberration should be omitted.

docs/coordinates/common_errors.rst

@@ -44,8 +44,8 @@ One might expect that the following code snippet would produce an altitude of ex
     >>> t = Time('J2010')
     >>> obj = EarthLocation(-1*u.deg, 52*u.deg, height=10.*u.km)
     >>> home = EarthLocation(-1*u.deg, 52*u.deg, height=0.*u.km)
-    >>> altaz_frame = AltAz(obstime=t, location=home)
-    >>> obj.get_itrs(t).transform_to(altaz_frame).alt # doctest: +FLOAT_CMP
+    >>> aa = obj.get_itrs(t).transform_to(AltAz(obstime=t, location=home))
+    >>> aa.alt # doctest: +FLOAT_CMP
     <Latitude 86.32878441 deg>
 
 Why is the result over three degrees away from the zenith? It is only possible to understand by taking a very careful
@@ -59,28 +59,24 @@ around 600 metres away from where it appears to be. This 600 metre shift, for an
 is an angular difference of just over three degrees - which is why this object does not appear overhead for a topocentric
 observer - one on the surface of the Earth.
 
-The correct way to construct a |SkyCoord| object for a source that is directly overhead a topocentric observer is
+The correct way to construct a |SkyCoord| object for a source that is directly overhead for a topocentric observer is
 as follows::
 
-    >>> from astropy.coordinates import EarthLocation, AltAz, ITRS, CIRS
+    >>> from astropy.coordinates import EarthLocation, AltAz, ITRS
     >>> from astropy.time import Time
     >>> from astropy import units as u
 
     >>> t = Time('J2010')
     >>> obj = EarthLocation(-1*u.deg, 52*u.deg, height=10.*u.km)
     >>> home = EarthLocation(-1*u.deg, 52*u.deg, height=0.*u.km)
 
-    >>> # Now we make a ITRS vector of a straight overhead object
+    >>> # First we make an ITRS vector of a straight overhead object
     >>> itrs_vec = obj.get_itrs(t).cartesian - home.get_itrs(t).cartesian
 
-    >>> # Now we create a topocentric coordinate with this data
-    >>> # Any topocentric frame will work, we use CIRS
-    >>> # Start by transforming the ITRS vector to CIRS
-    >>> cirs_vec = ITRS(itrs_vec, obstime=t).transform_to(CIRS(obstime=t)).cartesian
-    >>> # Finally, make CIRS frame object with the correct data
-    >>> cirs_topo = CIRS(cirs_vec, obstime=t, location=home)
+    >>> # Now we create a topocentric ITRS frame with this data
+    >>> itrs_topo = ITRS(itrs_vec, obstime=t, location=home)
 
     >>> # convert to AltAz
-    >>> aa = cirs_topo.transform_to(AltAz(obstime=t, location=home))
+    >>> aa = itrs_topo.transform_to(AltAz(obstime=t, location=home))
     >>> aa.alt # doctest: +FLOAT_CMP
     <Latitude 90. deg>

docs/coordinates/satellites.rst

@@ -5,8 +5,8 @@ Working with Earth Satellites Using Astropy Coordinates
 This document discusses Two-Line Element ephemerides and the True Equator, Mean Equinox frame.
 For satellite ephemerides given directly in geocentric ITRS coordinates
 (e.g. `ILRS ephemeris format <https://ilrs.gsfc.nasa.gov/data_and_products/formats/cpf.html>`_)
-please see the documentation of the ITRS frame and the example transform to `~astropy.coordinates.AltAz`
-below starting with the geocentric ITRS coordinate frame.
+please see the example transform to `~astropy.coordinates.AltAz` below starting with
+the geocentric ITRS coordinate frame.
 
 Satellite data is normally provided in the Two-Line Element (TLE) format
 (see `here <https://www.celestrak.com/NORAD/documentation/tle-fmt.php>`_
@@ -89,7 +89,7 @@ For example, to find the overhead latitude, longitude, and height of the satelli
 .. doctest-requires:: sgp4
 
     >>> from astropy.coordinates import ITRS
-    >>> itrs_geo = teme.transform_to(ITRS(obstime=t))  # doctest: +IGNORE_WARNINGS
+    >>> itrs_geo = teme.transform_to(ITRS(obstime=t))
     >>> location = itrs_geo.earth_location
     >>> location.geodetic  # doctest: +FLOAT_CMP
     GeodeticLocation(lon=<Longitude 160.34199789 deg>, lat=<Latitude -24.6609379 deg>, height=<Quantity 420.17927591 km>)
@@ -101,16 +101,33 @@ For example, to find the overhead latitude, longitude, and height of the satelli
 
 Or, if you want to find the altitude and azimuth of the satellite from a particular location:
 
+.. note ::
+    In this example, the intermediate step of manually setting up a topocentric `~astropy.coordinates.ITRS`
+    frame is necessary in order to avoid the change in stellar aberration that would occur if a direct
+    transform from geocentric to topocentric coordinates using ``transform_to`` was used. Please see
+    the documentation of the `~astropy.coordinates.ITRS` frame for more details.
+
 .. doctest-requires:: sgp4
 
     >>> from astropy.coordinates import EarthLocation, AltAz
     >>> siding_spring = EarthLocation.of_site('aao')  # doctest: +SKIP
-    >>> topo_itrs_repr = itrs_geo.cartesian - siding_spring.get_itrs(t).cartesian  # doctest: +IGNORE_WARNINGS
-    >>> itrs_topo = ITRS(topo_itrs_repr, obstime = t, location=siding_spring)  # doctest: +IGNORE_WARNINGS
-    >>> aa = itrs_topo.transform_to(AltAz(obstime=t, location=siding_spring))  # doctest: +IGNORE_WARNINGS
+    >>> topo_itrs_repr = itrs_geo.cartesian.without_differentials() - siding_spring.get_itrs(t).cartesian
+    >>> itrs_topo = ITRS(topo_itrs_repr, obstime = t, location=siding_spring)
+    >>> aa = itrs_topo.transform_to(AltAz(obstime=t, location=siding_spring))
     >>> aa.alt  # doctest: +FLOAT_CMP
     <Latitude 10.94799670 deg>
     >>> aa.az  # doctest: +FLOAT_CMP
     <Longitude 59.28803392 deg>
 
+For a stationary observer, velocity in the `~astropy.coordinates.ITRS` is independent of location,
+so if you want to carry the velocity to the topocentric frame, you can do so as follows:
+
+.. doctest-requires:: sgp4
+
+    >>> itrs_geo_p = itrs_geo.cartesian.without_differentials()
+    >>> itrs_geo_v = itrs_geo.cartesian.differentials['s']
+    >>> topo_itrs_p = itrs_geo_p - siding_spring.get_itrs(t).cartesian
+    >>> topo_itrs_repr = topo_itrs_p.with_differentials(itrs_geo_v)
+    >>> itrs_topo = ITRS(topo_itrs_repr, obstime = t, location=siding_spring)
+
 .. EXAMPLE END

docs/whatsnew/5.2.rst

@@ -14,6 +14,7 @@ In particular, this release includes:
 
 * :ref:`whatsnew-5.2-quantity-dtype`
 * :ref:`whatsnew-5.2-cosmology`
+* :ref:`whatsnew-5.2-coordinates`
 
 
 .. _whatsnew-5.2-quantity-dtype:
@@ -46,6 +47,38 @@ cosmologies with their non-flat equivalents.
     True
 
 
+.. _whatsnew-5.2-coordinates:
+
+Topocentric ITRS Frame
+======================
+
+A topocentric ITRS frame has been added that makes dealing with near-Earth objects
+easier and more intuitive.::
+
+    >>> from astropy.coordinates import EarthLocation, AltAz, ITRS
+    >>> from astropy.time import Time
+    >>> from astropy import units as u
+
+    >>> t = Time('J2010')
+    >>> obj = EarthLocation(-1*u.deg, 52*u.deg, height=10.*u.km)
+    >>> home = EarthLocation(-1*u.deg, 52*u.deg, height=0.*u.km)
+
+    >>> # Direction of object from GEOCENTER
+    >>> itrs_geo = obj.get_itrs(t).cartesian
+
+    >>> # now get the Geocentric ITRS position of observatory
+    >>> obsrepr = home.get_itrs(t).cartesian
+
+    >>> # topocentric ITRS position of a straight overhead object
+    >>> itrs_repr = itrs_geo - obsrepr
+
+    >>> # create an ITRS object that appears straight overhead for a TOPOCENTRIC OBSERVER
+    >>> itrs_topo = ITRS(itrs_repr, obstime=t, location=home)
+
+    >>> # convert to AltAz
+    >>> aa = itrs_topo.transform_to(AltAz(obstime=t, location=home))
+
+
 Full change log
 ===============