mdboom · Sep 22, 2014 · Aug 7, 2014 · Aug 7, 2014 · Aug 7, 2014 · Aug 7, 2014
diff --git a/CHANGES.rst b/CHANGES.rst
@@ -97,6 +97,11 @@ New Features
   - ``astropy.wcs.WCS.compare`` now supports a ``tolerance`` keyword argument
     to allow for approximate comparison of floating-point values. [#2503]
 
+  - added ``pixel_scale_matrix``, ``celestial``, ``is_celestial``, and
+    ``has_celestial`` convenience attributes.  Added ``celestial_pixel_scale``
+    and ``non_celestial_pixel_scales`` utility functions for retrieving WCS
+    pixel scale information [#2832]
+
 API Changes
 ^^^^^^^^^^^
 

diff --git a/astropy/wcs/tests/test_utils.py b/astropy/wcs/tests/test_utils.py
@@ -2,8 +2,13 @@
 from __future__ import absolute_import, division, print_function, unicode_literals
 from ...wcs import WCS
 from .. import utils
-from ...tests.helper import pytest
+from ..utils import celestial_pixel_scale, non_celestial_pixel_scales
+from ...tests.helper import pytest, catch_warnings
+from ...utils.exceptions import AstropyUserWarning
+from ... import units as u
+
 import numpy as np
+from numpy.testing import assert_almost_equal
 
 def test_wcs_dropping():
     wcs = WCS(naxis=4)
@@ -132,14 +137,20 @@ def test_invalid_slice():
 
 def test_axis_names():
     mywcs = WCS(naxis=4)
-    mywcs.wcs.ctype = ['RA---TAN','DEC---TAN','VOPT-LSR','STOKES']
+    mywcs.wcs.ctype = ['RA---TAN','DEC--TAN','VOPT-LSR','STOKES']
 
     assert mywcs.axis_type_names == ['RA','DEC','VOPT','STOKES']
 
-    mywcs.wcs.cname = ['RA','DEC','VOPT-LSR','STOKES']
+    mywcs.wcs.cname = ['RA','DEC','VOPT','STOKES']
 
-    assert mywcs.axis_type_names == ['RA','DEC','VOPT-LSR','STOKES']
+    assert mywcs.axis_type_names == ['RA','DEC','VOPT','STOKES']
 
+def test_celestial():
+    mywcs = WCS(naxis=4)
+    mywcs.wcs.ctype = ['RA---TAN','DEC--TAN','VOPT','STOKES']
+    cel = mywcs.celestial
+    assert list(cel.wcs.ctype) == ['RA---TAN','DEC--TAN']
+    assert cel.axis_type_names == ['RA','DEC']
 
 def test_wcs_to_celestial_frame():
 
@@ -219,3 +230,126 @@ def identify_offset(wcs):
     # Check that things are back to normal after the context manager
     with pytest.raises(ValueError):
         utils.wcs_to_celestial_frame(mywcs)
+
+def test_pixscale_nodrop():
+    mywcs = WCS(naxis=2)
+    mywcs.wcs.cdelt = [0.1,0.1]
+    mywcs.wcs.ctype = ['RA---TAN','DEC--TAN']
+    assert_almost_equal(celestial_pixel_scale(mywcs), 0.1*u.deg)
+
+    mywcs.wcs.cdelt = [-0.1,0.1]
+    assert_almost_equal(celestial_pixel_scale(mywcs), 0.1*u.deg)
+
+def test_pixscale_withdrop():
+    mywcs = WCS(naxis=3)
+    mywcs.wcs.cdelt = [0.1,0.1,1]
+    mywcs.wcs.ctype = ['RA---TAN','DEC--TAN','VOPT']
+    assert_almost_equal(celestial_pixel_scale(mywcs), 0.1*u.deg)
+
+    mywcs.wcs.cdelt = [-0.1,0.1,1]
+    assert_almost_equal(celestial_pixel_scale(mywcs), 0.1*u.deg)
+
+
+def test_pixscale_cd():
+    mywcs = WCS(naxis=2)
+    mywcs.wcs.cd = [[-0.1,0],[0,0.1]]
+    mywcs.wcs.ctype = ['RA---TAN','DEC--TAN']
+    assert_almost_equal(celestial_pixel_scale(mywcs), 0.1*u.deg)
+
+def test_pixscale_warning(recwarn):
+    mywcs = WCS(naxis=2)
+    mywcs.wcs.cd = [[-0.1,0],[0,0.1]]
+    mywcs.wcs.ctype = ['RA---TAN','DEC--TAN']
+
+    with catch_warnings(AstropyUserWarning) as warning_lines:
+
+        celestial_pixel_scale(mywcs)
+        assert ("Pixel sizes may very over the image for "
+                "projection class TAN"
+                in str(warning_lines[0].message))
+
+def test_pixscale_asymmetric():
+    mywcs = WCS(naxis=2)
+    mywcs.wcs.cd = [[-0.2,0],[0,0.1]]
+    mywcs.wcs.ctype = ['RA---TAN','DEC--TAN']
+
+    with pytest.raises(ValueError) as exc:
+        celestial_pixel_scale(mywcs)
+    assert exc.value.args[0] == "Pixels are not square: 'pixel scale' is ambiguous"
+
+@pytest.mark.parametrize('angle',
+                         (30,45,60,75))
+def test_pixscale_cd_rotated(angle):
+    mywcs = WCS(naxis=2)
+    rho = angle/180.*np.pi
+    scale = 0.1
+    mywcs.wcs.cd = [[scale*np.cos(rho), -scale*np.sin(rho)],
+                    [scale*np.sin(rho), scale*np.cos(rho)]]
+    mywcs.wcs.ctype = ['RA---TAN','DEC--TAN']
+    assert_almost_equal(celestial_pixel_scale(mywcs), 0.1*u.deg)
+
+@pytest.mark.parametrize('angle',
+                         (30,45,60,75))
+def test_pixscale_pc_rotated(angle):
+    mywcs = WCS(naxis=2)
+    rho = angle/180.*np.pi
+    scale = 0.1
+    mywcs.wcs.cdelt = [-scale, scale]
+    mywcs.wcs.pc = [[np.cos(rho), -np.sin(rho)],
+                    [np.sin(rho), np.cos(rho)]]
+    mywcs.wcs.ctype = ['RA---TAN','DEC--TAN']
+    assert_almost_equal(celestial_pixel_scale(mywcs), 0.1*u.deg)
+
+@pytest.mark.parametrize(('cdelt','pc','pccd'),
+                         (([0.1,0.2], np.eye(2), np.diag([0.1,0.2])),
+                          ([0.1,0.2,0.3], np.eye(3), np.diag([0.1,0.2,0.3])),
+                          ([1,1,1], np.diag([0.1,0.2,0.3]), np.diag([0.1,0.2,0.3])),
+                         ))
+def test_pixel_scale_matrix(cdelt, pc, pccd):
+
+    mywcs = WCS(naxis=(len(cdelt)))
+    mywcs.wcs.cdelt = cdelt
+    mywcs.wcs.pc = pc
+
+    assert_almost_equal(mywcs.pixel_scale_matrix, pccd)
+
+@pytest.mark.parametrize(('ctype', 'cel'),
+                         ((['RA---TAN','DEC--TAN'], True),
+                          (['RA---TAN','DEC--TAN','FREQ'], False),
+                          (['RA---TAN','FREQ'], False),))
+def test_is_celestial(ctype, cel):
+    mywcs = WCS(naxis=len(ctype))
+    mywcs.wcs.ctype = ctype
+
+    assert mywcs.is_celestial == cel
+
+@pytest.mark.parametrize(('ctype', 'cel'),
+                         ((['RA---TAN','DEC--TAN'], True),
+                          (['RA---TAN','DEC--TAN','FREQ'], True),
+                          (['RA---TAN','FREQ'], False),))
+def test_has_celestial(ctype, cel):
+    mywcs = WCS(naxis=len(ctype))
+    mywcs.wcs.ctype = ctype
+
+    assert mywcs.has_celestial == cel
+
+@pytest.mark.parametrize(('cdelt','pc','cd'),
+                         ((np.array([0.1,0.2]), np.eye(2), np.eye(2)),
+                          (np.array([1,1]), np.diag([0.1,0.2]), np.eye(2)),
+                          (np.array([0.1,0.2]), np.eye(2), None),
+                          (np.array([0.1,0.2]), None, np.eye(2)),
+                          )
+                        )
+def test_noncelestial_scale(cdelt, pc, cd):
+    mywcs = WCS(naxis=2)
+    if cd is not None:
+        mywcs.wcs.cd = cd
+    if pc is not None:
+        mywcs.wcs.pc = pc
+    mywcs.wcs.cdelt = cdelt
+
+    mywcs.wcs.ctype = ['RA---TAN','FREQ']
+
+    ps = non_celestial_pixel_scales(mywcs)
+
+    assert_almost_equal(ps, [0.1,0.2]*u.deg)
diff --git a/astropy/wcs/utils.py b/astropy/wcs/utils.py
@@ -2,9 +2,14 @@
 from __future__ import absolute_import, division, print_function, unicode_literals
 
 import numpy as np
+import warnings
+from .. import units as u
+from ..utils.exceptions import AstropyUserWarning
 
 __doctest_skip__ = ['wcs_to_celestial_frame']
 
+__all__ = ['add_stokes_axis_to_wcs', 'wcs_to_celestial_frame',
+           'celestial_pixel_scale', 'non_celestial_pixel_scales']
 
 def add_stokes_axis_to_wcs(wcs, add_before_ind):
     """
@@ -107,15 +112,15 @@ def wcs_to_celestial_frame(wcs):
     For a given WCS, return the coordinate frame that matches the celestial
     component of the WCS.
 
-    Paramters
-    ---------
+    Parameters
+    ----------
     wcs : :class:`~astropy.wcs.WCS` instance
         The WCS to find the frame for
 
     Returns
     -------
-    frame : :class:`~astropy.coordinates.base_frame.BaseFrame` subclass instance
-        An instance of a :class:`~astropy.coordinates.base_frame.BaseFrame`
+    frame : :class:`~astropy.coordinates.baseframe.BaseCoordinateFrame` subclass instance
+        An instance of a :class:`~astropy.coordinates.baseframe.BaseCoordinateFrame`
         subclass instance that best matches the specified WCS.
 
     Notes
@@ -126,9 +131,10 @@ def wcs_to_celestial_frame(wcs):
     instance and should return either an instance of a frame, or `None` if no
     matching frame was found. You can register this function temporarily with::
 
-    >>> from astropy.wcs.utils import wcs_to_celestial_frame, custom_frame_mappings
-    >>> with custom_frame_mappings(my_function):
-    ...     wcs_to_celestial_frame(...)
+        >>> from astropy.wcs.utils import wcs_to_celestial_frame, custom_frame_mappings
+        >>> with custom_frame_mappings(my_function):
+        ...     wcs_to_celestial_frame(...)
+
     """
     for mapping_set in WCS_FRAME_MAPPINGS:
         for func in mapping_set:
@@ -137,3 +143,67 @@ def wcs_to_celestial_frame(wcs):
                 return frame
     raise ValueError("Could not determine celestial frame corresponding to "
                      "the specified WCS object")
+
+
+def celestial_pixel_scale(inwcs, allow_nonsquare=False):
+    """
+    For a WCS, if the pixels are square, return the pixel scale in the spatial
+    dimensions
+
+    Parameters
+    ----------
+    inwcs: `astropy.wcs.WCS`
+        The world coordinate system object
+    allow_nonsquare : bool
+        Return the average of the X and Y scales if True
+
+    Returns
+    -------
+    scale : float
+        The square pixel scale
+
+    Raises
+    ------
+    ValueError if the pixels are nonsquare and ``allow_nonsquare==False``
+    """
+    cwcs = inwcs.celestial
+    if cwcs.wcs.ctype[0][-3:] != 'CAR':
+        warnings.warn("Pixel sizes may very over the image for "
+                      "projection class {0}".format(cwcs.wcs.ctype[0][-3:]),
+                      AstropyUserWarning)
+    scale = (cwcs.pixel_scale_matrix**2).sum(axis=0)**0.5
+    if not np.allclose(scale[0],scale[1]):
+        if allow_nonsquare:
+            warnings.warn("Pixels are not square, using an average pixel scale")
+            return np.mean(scale)*u.deg
+        else:
+            raise ValueError("Pixels are not square: 'pixel scale' is ambiguous")
+    # Return a quantity: WCS always stores in degrees
+    return scale[0]*u.deg
+
+
+def non_celestial_pixel_scales(inwcs):
+    """
+    For a non-celestial WCS, e.g. one with mixed spectral and spatial axes, it
+    is still sometimes possible to define a pixel scale.
+
+    Parameters
+    ----------
+    inwcs: `astropy.wcs.WCS`
+        The world coordinate system object
+
+    Returns
+    -------
+    scale : `numpy.ndarray`
+        The pixel scale along each axis
+    """
+
+    if inwcs.is_celestial:
+        raise ValueError("WCS is celestial, use celestial_pixel_scale instead")
+
+    pccd = inwcs.pixel_scale_matrix
+
+    if np.allclose(np.extract(1-np.eye(*pccd.shape), pccd), 0):
+        return np.abs(np.diagonal(pccd))*u.deg
+    else:
+        raise ValueError("WCS is rotated, cannot determine consistent pixel scales")
diff --git a/astropy/wcs/wcs.py b/astropy/wcs/wcs.py
@@ -2719,6 +2719,46 @@ def axis_type_names(self):
             names[i] = types[i].split('-')[0]
         return names
 
+    @property
+    def celestial(self):
+        """
+        A copy of the current WCS with only the celestial axes included
+        """
+        return self.sub([WCSSUB_CELESTIAL])
+
+    @property
+    def is_celestial(self):
+        return self.has_celestial and self.naxis==2
+
+    @property
+    def has_celestial(self):
+        try:
+            return self.celestial.naxis == 2
+        except InconsistentAxisTypesError:
+            return False
+
+    @property
+    def pixel_scale_matrix(self):
+
+        try:
+            cdelt = np.matrix(np.diag(self.wcs.get_cdelt()))
+            pc = np.matrix(self.wcs.get_pc())
+        except InconsistentAxisTypesError:
+            try:
+                # for non-celestial axes, get_cdelt doesnt work
+                cdelt = np.matrix(self.wcs.cd) * np.matrix(np.diag(self.wcs.cdelt))
+            except AttributeError:
+                cdelt = np.matrix(np.diag(self.wcs.cdelt))
+
+            try:
+                pc = np.matrix(self.wcs.pc)
+            except AttributeError:
+                pc = 1
+
+        pccd = np.array(cdelt * pc)
+
+        return pccd
+
 
 def __WCS_unpickle__(cls, dct, fits_data):
     """

diff --git a/docs/wcs/index.rst b/docs/wcs/index.rst
@@ -169,6 +169,19 @@ through `Wcsprm.cunit <astropy.wcs.Wcsprm.cunit>`), for example,
 - ``HPX``: HEALPix
 - ``XPH``: HEALPix polar, aka "butterfly"
 
+Subsetting and Pixel Scales
+===========================
+
+WCS objects can be broken apart into their constituent axes using the
+`~astropy.wcs.WCS.sub` function.  There is also a `~astropy.wcs.WCS.celestial`
+convenience function that will return a WCS object with only the celestial axes
+included.
+
+The pixel scale of a celestial image or the pixel dimensions of a non-celestial
+image can be extracted with the utility functions
+`~astropy.wcs.utils.celestial_pixel_scale` and
+`~astropy.wcs.utils.non_celestial_pixel_scales`.
+
 Other information
 =================
 
@@ -190,6 +203,7 @@ Reference/API
 
 .. automodapi:: astropy.wcs
 
+.. automodapi:: astropy.wcs.utils
 
 Acknowledgments and Licenses
 ============================