Map ITRS frames to terrestrial WCS coordinates

This will make it possible to use WCSAxes to make figures that
combine both celestial and terrestrial features. An example is
plotting the coordinates of an astronomical transient over an all-
sky satellite image to illustrate the position relative to the
Earth at the time of the event.

The ITRS frame is identified with WCSs that use the `TLON-` and
`TLAT-` coordinate types. There are several examples of WCSs where
this syntax is used to describe terrestrial coordinate systems:
Section 7.4.1 of the WCS in FITS "Paper II", and the WCSTools docs:
http://tdc-www.harvard.edu/software/wcstools/wcstools.multiwcs.html

This is the technique that I am using to make images such as this one:
http://ligo.org/detections/GW170817/images-GW170817/skymap+earth.jpg

CHANGES.rst

@@ -71,6 +71,17 @@ astropy.wcs
 ^^^^^^^^^^^
 - ``wcslib`` was updated to v 5.18. [#7066]
 
+- Map ITRS frames to terrestrial WCS coordinates. This will make it possible to
+  use WCSAxes to make figures that combine both celestial and terrestrial
+  features. An example is plotting the coordinates of an astronomical transient
+  over an all- sky satellite image to illustrate the position relative to the
+  Earth at the time of the event. The ITRS frame is identified with WCSs that
+  use the ``TLON-`` and ``TLAT-`` coordinate types. There are several examples
+  of WCSs where this syntax is used to describe terrestrial coordinate systems:
+  Section 7.4.1 of `WCS in FITS "Paper II" <http://adsabs.harvard.edu/abs/2002A%26A...395.1077C>`_
+  and the `WCSTools documentation <http://tdc-www.harvard.edu/software/wcstools/wcstools.multiwcs.html>`_.
+  [#6990]
+
 API Changes
 -----------
 

astropy/wcs/tests/test_utils.py

@@ -187,7 +187,7 @@ def test_celestial():
 def test_wcs_to_celestial_frame():
 
     # Import astropy.coordinates here to avoid circular imports
-    from ...coordinates.builtin_frames import ICRS, FK5, FK4, Galactic
+    from ...coordinates.builtin_frames import ICRS, ITRS, FK5, FK4, Galactic
 
     mywcs = WCS(naxis=2)
     with pytest.raises(ValueError) as exc:
@@ -217,6 +217,12 @@ def test_wcs_to_celestial_frame():
     frame = wcs_to_celestial_frame(mywcs)
     assert isinstance(frame, Galactic)
 
+    mywcs.wcs.ctype = ['TLON-CAR', 'TLAT-CAR']
+    mywcs.wcs.dateobs = '2017-08-17T12:41:04.430'
+    frame = wcs_to_celestial_frame(mywcs)
+    assert isinstance(frame, ITRS)
+    assert frame.obstime == Time('2017-08-17T12:41:04.430')
+
     mywcs.wcs.ctype = ['RA---TAN', 'DEC--TAN']
     mywcs.wcs.radesys = 'ICRS'
 
@@ -264,7 +270,7 @@ def identify_offset(wcs):
 def test_celestial_frame_to_wcs():
 
     # Import astropy.coordinates here to avoid circular imports
-    from ...coordinates import ICRS, FK5, FK4, FK4NoETerms, Galactic, BaseCoordinateFrame
+    from ...coordinates import ICRS, ITRS, FK5, FK4, FK4NoETerms, Galactic, BaseCoordinateFrame
 
     class FakeFrame(BaseCoordinateFrame):
         pass
@@ -320,6 +326,18 @@ class FakeFrame(BaseCoordinateFrame):
     mywcs.wcs.set()
     assert_allclose((mywcs.wcs.lonpole, mywcs.wcs.latpole), (180, 60))
 
+    frame = ITRS(obstime=Time('2017-08-17T12:41:04.43'))
+    mywcs = celestial_frame_to_wcs(frame, projection='CAR')
+    assert tuple(mywcs.wcs.ctype) == ('TLON-CAR', 'TLAT-CAR')
+    assert mywcs.wcs.radesys == 'ITRS'
+    assert mywcs.wcs.dateobs == '2017-08-17T12:41:04.430'
+
+    frame = ITRS()
+    mywcs = celestial_frame_to_wcs(frame, projection='CAR')
+    assert tuple(mywcs.wcs.ctype) == ('TLON-CAR', 'TLAT-CAR')
+    assert mywcs.wcs.radesys == 'ITRS'
+    assert mywcs.wcs.dateobs == Time('J2000').utc.isot
+
 
 def test_celestial_frame_to_wcs_extend():
 

astropy/wcs/utils.py

@@ -45,7 +45,7 @@ def add_stokes_axis_to_wcs(wcs, add_before_ind):
 def _wcs_to_celestial_frame_builtin(wcs):
 
     # Import astropy.coordinates here to avoid circular imports
-    from ..coordinates import FK4, FK4NoETerms, FK5, ICRS, Galactic
+    from ..coordinates import FK4, FK4NoETerms, FK5, ICRS, ITRS, Galactic
 
     # Import astropy.time here otherwise setup.py fails before extensions are compiled
     from ..time import Time
@@ -92,6 +92,8 @@ def _wcs_to_celestial_frame_builtin(wcs):
     else:
         if xcoord == 'GLON' and ycoord == 'GLAT':
             frame = Galactic()
+        elif xcoord == 'TLON' and ycoord == 'TLAT':
+            frame = ITRS(obstime=wcs.wcs.dateobs or None)
         else:
             frame = None
 
@@ -101,7 +103,7 @@ def _wcs_to_celestial_frame_builtin(wcs):
 def _celestial_frame_to_wcs_builtin(frame, projection='TAN'):
 
     # Import astropy.coordinates here to avoid circular imports
-    from ..coordinates import BaseRADecFrame, FK4, FK4NoETerms, FK5, ICRS, Galactic
+    from ..coordinates import BaseRADecFrame, FK4, FK4NoETerms, FK5, ICRS, ITRS, Galactic
 
     # Create a 2-dimensional WCS
     wcs = WCS(naxis=2)
@@ -126,6 +128,11 @@ def _celestial_frame_to_wcs_builtin(frame, projection='TAN'):
     elif isinstance(frame, Galactic):
         xcoord = 'GLON'
         ycoord = 'GLAT'
+    elif isinstance(frame, ITRS):
+        xcoord = 'TLON'
+        ycoord = 'TLAT'
+        wcs.wcs.radesys = 'ITRS'
+        wcs.wcs.dateobs = frame.obstime.utc.isot
     else:
         return None