From 143190df480712c2866dfd0cc93853afef6cfae1 Mon Sep 17 00:00:00 2001
From: e-koch <koch.eric.w@gmail.com>
Date: Thu, 22 Oct 2020 14:43:04 -0400
Subject: [PATCH 01/23] Return 2D planes in WCS coordsinates given a 2D plane
 of pixel positions

---
 spectral_cube/cube_utils.py | 40 +++++++++++++++++++++++++++++++++++++
 1 file changed, 40 insertions(+)

diff --git a/spectral_cube/cube_utils.py b/spectral_cube/cube_utils.py
index 76b70ff30..6863d039c 100644
--- a/spectral_cube/cube_utils.py
+++ b/spectral_cube/cube_utils.py
@@ -486,3 +486,43 @@ def convert_bunit(bunit):
             unit = None
 
     return unit
+
+
+def world_take_along_axis(cube, posn_plane, axis):
+    '''
+    Apply a 2D plane of pixel positions along collapsed along a third axis
+    to the equivalent WCS coordinates. For example, this will convert `argmax`
+    along the spectral axis to the equivalent spectral value (e.g., velocity at
+    peak intensity).
+    Parameters
+    ----------
+    cube : SpectralCube
+        A spectral cube.
+    posn_plane : 2D numpy.ndarray
+        2D array of pixel dimensions along `axis`.
+    axis : int
+        The axis that `posn_plane` is collapsed along.
+    Returns
+    -------
+    out : astropy.units.Quantity
+        2D array of WCS coordinates.
+    '''
+    # Get 1D slice along that axis.
+    world_slice = [0, 0]
+    world_slice.insert(axis, slice(None))
+
+    world_coords = cube.world[tuple(world_slice)][axis]
+
+    world_newaxis = [np.newaxis] * 2
+    world_newaxis.insert(axis, slice(None))
+    world_newaxis = tuple(world_newaxis)
+
+    plane_newaxis = [slice(None), slice(None)]
+    plane_newaxis.insert(axis, np.newaxis)
+    plane_newaxis = tuple(plane_newaxis)
+
+    out = np.take_along_axis(world_coords[world_newaxis],
+                             posn_plane[plane_newaxis], axis=axis)
+    out = out.squeeze()
+
+    return out

From 8723522aa05c35e689c883d4fd468ec7ba14688d Mon Sep 17 00:00:00 2001
From: e-koch <koch.eric.w@gmail.com>
Date: Thu, 22 Oct 2020 14:44:35 -0400
Subject: [PATCH 02/23] Add functions to return argmax/min along one dimension
 in WCS coordinates

---
 spectral_cube/spectral_cube.py | 60 ++++++++++++++++++++++++++++++++++
 1 file changed, 60 insertions(+)

diff --git a/spectral_cube/spectral_cube.py b/spectral_cube/spectral_cube.py
index 797a12117..e508f0fae 100644
--- a/spectral_cube/spectral_cube.py
+++ b/spectral_cube/spectral_cube.py
@@ -812,6 +812,66 @@ def argmin(self, axis=None, how='auto', **kwargs):
                                          reduce=False, projection=False,
                                          how=how, axis=axis, **kwargs)
 
+    @warn_slow
+    def argmax_world(self, axis, **kwargs):
+        '''
+        Return the spatial or spectral index of the maximum value
+        along a line of sight.
+        Parameters
+        ----------
+        axis : int
+            The axis to return the peak location along. e.g., `axis=0`
+            will return the value of the spectral axis at the peak value.
+        '''
+
+        argmax_plane = self.argmax(axis=axis, **kwargs)
+
+        # Convert to WCS coordinates.
+        out = cube_utils.world_take_along_axis(self, argmax_plane, axis)
+
+        # Compute whether the mask has any valid data along `axis`
+        collapsed_mask = self.mask.include().any(axis=axis)
+        out[~collapsed_mask] = np.NaN
+
+        # Return a Projection.
+        new_wcs = wcs_utils.drop_axis(self._wcs, np2wcs[axis])
+
+        meta = {'collapse_axis': axis}
+        meta.update(self._meta)
+
+        return Projection(out, copy=False, wcs=new_wcs, meta=meta,
+                          unit=out.unit, header=self._nowcs_header)
+
+    @warn_slow
+    def argmin_world(self, axis, **kwargs):
+        '''
+        Return the spatial or spectral index of the minimum value
+        along a line of sight.
+        Parameters
+        ----------
+        axis : int
+            The axis to return the peak location along. e.g., `axis=0`
+            will return the value of the spectral axis at the peak value.
+        '''
+
+        argmin_plane = self.argmin(axis=axis, **kwargs)
+
+        # Convert to WCS coordinates.
+        out = cube_utils.world_take_along_axis(self, argmin_plane, axis)
+
+        # Compute whether the mask has any valid data along `axis`
+        collapsed_mask = self.mask.include().any(axis=axis)
+        out[~collapsed_mask] = np.NaN
+
+        # Return a Projection.
+        new_wcs = wcs_utils.drop_axis(self._wcs, np2wcs[axis])
+
+        meta = {'collapse_axis': axis}
+        meta.update(self._meta)
+
+        return Projection(out, copy=False, wcs=new_wcs, meta=meta,
+                          unit=out.unit, header=self._nowcs_header)
+
     def chunked(self, chunksize=1000):
         """
         Not Implemented.

From 3779d4fa7d76f7c08208d234c4c6981b4241d90e Mon Sep 17 00:00:00 2001
From: e-koch <koch.eric.w@gmail.com>
Date: Thu, 22 Oct 2020 14:46:13 -0400
Subject: [PATCH 03/23] Add world argmin/max to tests

---
 spectral_cube/tests/test_spectral_cube.py | 13 +++++++++++++
 1 file changed, 13 insertions(+)

diff --git a/spectral_cube/tests/test_spectral_cube.py b/spectral_cube/tests/test_spectral_cube.py
index 7cc5f68b3..822c3173a 100644
--- a/spectral_cube/tests/test_spectral_cube.py
+++ b/spectral_cube/tests/test_spectral_cube.py
@@ -719,6 +719,19 @@ def test_transpose(self, method, data_adv, data_vad, use_dask):
                                 getattr(c2, method)(axis=axis, progressbar=True))
         self.c = self.d = None
 
+    @pytest.mark.parametrize('method', ('argmax_world', 'argmin_world'))
+    def test_arg_world(self, method, data_adv, use_dask):
+        c1, d1 = cube_and_raw(data_adv, use_dask=use_dask)
+
+        # Pixel operation is same name with "_world" removed.
+        arg0_pixel = getattr(c1, method.split("_")[0])(axis=0)
+
+        arg0_world = np.take_along_axis(c1.spectral_axis[:, np.newaxis, np.newaxis],
+                                        arg0_pixel[np.newaxis, :, :], axis=0).squeeze()
+
+        assert_allclose(getattr(c1, method)(axis=0), arg0_world)
+
+        self.c = self.d = None
 
 class TestSlab(BaseTest):
 

From abb0fed87d78b906ae4f574e172a65bb97ade772 Mon Sep 17 00:00:00 2001
From: e-koch <koch.eric.w@gmail.com>
Date: Thu, 22 Oct 2020 14:47:35 -0400
Subject: [PATCH 04/23] Also test WCS argmin/max in existing tests. But need to
 account for axis=None cases that aren't defined for these operations

---
 spectral_cube/tests/test_spectral_cube.py | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/spectral_cube/tests/test_spectral_cube.py b/spectral_cube/tests/test_spectral_cube.py
index 822c3173a..7d1a2aaad 100644
--- a/spectral_cube/tests/test_spectral_cube.py
+++ b/spectral_cube/tests/test_spectral_cube.py
@@ -706,11 +706,14 @@ def test_percentile(self, pct, iterate_rays, use_dask):
         self.c = self.d = None
 
     @pytest.mark.parametrize('method', ('sum', 'min', 'max', 'std', 'mad_std',
-                                        'median', 'argmin', 'argmax'))
+                                        'median', 'argmin', 'argmax',
+                                        'argmax_world', 'argmin_world'))
     def test_transpose(self, method, data_adv, data_vad, use_dask):
         c1, d1 = cube_and_raw(data_adv, use_dask=use_dask)
         c2, d2 = cube_and_raw(data_vad, use_dask=use_dask)
         for axis in [None, 0, 1, 2]:
+            if 'world' in method and axis == None:
+                continue
             assert_allclose(getattr(c1, method)(axis=axis),
                             getattr(c2, method)(axis=axis))
             if not use_dask:

From c6c4ca84dc405c5bf2eab1e0b4e54281a6768e16 Mon Sep 17 00:00:00 2001
From: e-koch <koch.eric.w@gmail.com>
Date: Thu, 22 Oct 2020 16:50:03 -0400
Subject: [PATCH 05/23] Add to the moments doc page with description of other
 common operations to make 2D maps (max/min, argmin/max)

---
 docs/moments.rst | 53 +++++++++++++++++++++++++++++++++++++++---------
 1 file changed, 43 insertions(+), 10 deletions(-)

diff --git a/docs/moments.rst b/docs/moments.rst
index 1b5464bc2..acbb472da 100644
--- a/docs/moments.rst
+++ b/docs/moments.rst
@@ -22,25 +22,25 @@ axis::
           moment is computed as the variance. For the actual formulas used for
           the moments, please see :class:`~spectral_cube.SpectralCube.moment`.
           For linewidth maps, see the `Linewidth maps`_ section.
-          
+
 You may also want to convert the unit of the datacube into a velocity one before
-you can obtain a genuine velocity map via a 1st moment map. So first it will be necessary to 
+you can obtain a genuine velocity map via a 1st moment map. So first it will be necessary to
 apply the :class:`~spectral_cube.SpectralCube.with_spectral_unit` method from this package with the proper attribute settings::
 
     >>> nii_cube = cube.with_spectral_unit(u.km/u.s,
                                            velocity_convention='optical',
                                            rest_value=6584*u.AA)  # doctest: +SKIP
 
-Note that the ``rest_value`` in the above code refers to the wavelength of the targeted line 
-in the 1D spectrum corresponding to the 3rd dimension. Also, since not all velocity values are relevant, 
-next we will use the :class:`~spectral_cube.SpectralCube.spectral_slab` method to slice out the chunk of 
+Note that the ``rest_value`` in the above code refers to the wavelength of the targeted line
+in the 1D spectrum corresponding to the 3rd dimension. Also, since not all velocity values are relevant,
+next we will use the :class:`~spectral_cube.SpectralCube.spectral_slab` method to slice out the chunk of
 the cube that actually contains the line::
 
     >>> nii_cube = cube.with_spectral_unit(u.km/u.s,
                                            velocity_convention='optical',
                                            rest_value=6584*u.AA)  # doctest: +SKIP
     >>> nii_subcube = nii_cube.spectral_slab(-60*u.km/u.s,-20*u.km/u.s)  # doctest: +SKIP
-    
+
 Finally, we can now generate the 1st moment map containing the expected velocity structure::
 
     >>> moment_1 = nii_subcube.moment(order=1)  # doctest: +SKIP
@@ -73,7 +73,7 @@ will create the quicklook grayscale image and save it to a png all in one go.
 
 Moment map equations
 ^^^^^^^^^^^^^^^^^^^^
- 
+
 The moments are defined below, using :math:`v` for the spectral (velocity,
 frequency, wavelength, or energy) axis and :math:`I_v` as the intensity,
 or otherwise measured flux, value in a given spectral channel.
@@ -85,7 +85,7 @@ The equation for the 0th moment is:
 The equation for the 1st moment is:
 
 .. math:: M_1 = \frac{\int v I_v  dv}{\int I_v dv} = \frac{\int v I_v dv}{M_0}
-   
+
 Higher-order moments (:math:`N\geq2`) are defined as:
 
 .. math:: M_N = \frac{\int I_v (v - M_1)^N dv}{M_0}
@@ -112,14 +112,47 @@ with either of these two commands::
     >>> fwhm_map = cube.linewidth_fwhm()  # doctest: +SKIP
 
 ``~spectral_cube.SpectralCube.linewidth_sigma`` computes a sigma linewidth map
-along the spectral axis, where sigma is the width of a Gaussian, while 
+along the spectral axis, where sigma is the width of a Gaussian, while
 ``~spectral_cube.SpectralCube.linewidth_fwhm`` computes a FWHM
 linewidth map along the same spectral axis.
 
 The linewidth maps are related to the second moment by
 
 .. math:: \sigma = \sqrt{M_2} \\
-          FWHM = \sigma \sqrt{8 ln{2}} 
+          FWHM = \sigma \sqrt{8 ln{2}}
 
 These functions return :class:`~spectral_cube.lower_dimensional_structures.Projection` instances as for the
 `Moment maps`_.
+
+Additional 2D maps
+------------------
+
+Other common 2D views of a spectral cube include the maximum/minimum along a dimension and
+the location of the maximum or the minimum along that dimension.
+
+To produce a 2D map of the maximum along a cube dimension, use `~spectral_cube.SpectralCube.max`
+
+    >>> max_map = cube.max(axis=0)  # doctest: +SKIP
+
+Along the spectral axis, this will return the peak intensity map.
+
+Similarly we can use `~spectral_cube.SpectralCube.min` to make a 2D minimum map:
+
+    >>> min_map = cube.min(axis=0)  # doctest: +SKIP
+
+The `~spectral_cube.SpectralCube.argmax` and `~spectral_cube.SpectralCube.argmin` operations will
+return the pixel positions of the max/min along that axis:
+
+    >>> argmax_map = cube.argmax(axis=0)  # doctest: +SKIP
+    >>> argmin_map = cube.argmin(axis=0)  # doctest: +SKIP
+
+These maps are useful for identifying where signal is located within the spectral cube, however,
+it is more useful to return the WCS values of those pixels for comparisons with other data sets
+or for modeling. The `~spectral_cube.SpectralCube.argmax_world` and
+`~spectral_cube.SpectralCube.argmin_world` should be used in this case:
+
+    >>> world_argmax_map = cube.argmax_world(axis=0)  # doctest: +SKIP
+    >>> world_argmin_map = cube.argmin_world(axis=0)  # doctest: +SKIP
+
+Along the spectral axis, `~spectral_cube.SpectralCube.argmax_world` creates the often used
+"velocity at peak intensity," which may also be called the "peak velocity."

From 056b88ac9ca02b817e1a078fd909ed556d9acd1a Mon Sep 17 00:00:00 2001
From: e-koch <koch.eric.w@gmail.com>
Date: Thu, 22 Oct 2020 16:51:48 -0400
Subject: [PATCH 06/23] Finish off docstrings

---
 spectral_cube/spectral_cube.py | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/spectral_cube/spectral_cube.py b/spectral_cube/spectral_cube.py
index e508f0fae..7312eb5d8 100644
--- a/spectral_cube/spectral_cube.py
+++ b/spectral_cube/spectral_cube.py
@@ -817,11 +817,14 @@ def argmax_world(self, axis, **kwargs):
         '''
         Return the spatial or spectral index of the maximum value
         along a line of sight.
+
         Parameters
         ----------
         axis : int
             The axis to return the peak location along. e.g., `axis=0`
             will return the value of the spectral axis at the peak value.
+        kwargs : dict
+             Passed to `~SpectralCube.argmax`.
         '''
 
         argmax_plane = self.argmax(axis=axis, **kwargs)
@@ -852,6 +855,8 @@ def argmin_world(self, axis, **kwargs):
         axis : int
             The axis to return the peak location along. e.g., `axis=0`
             will return the value of the spectral axis at the peak value.
+        kwargs : dict
+             Passed to `~SpectralCube.argmax`.
         '''
 
         argmin_plane = self.argmin(axis=axis, **kwargs)

From b32cea839e77ca46ae9cfed3b42cba9dac07eef8 Mon Sep 17 00:00:00 2001
From: e-koch <koch.eric.w@gmail.com>
Date: Sat, 24 Oct 2020 17:05:14 -0400
Subject: [PATCH 07/23] Fix docstring issue for sphinx build

---
 spectral_cube/spectral_cube.py | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/spectral_cube/spectral_cube.py b/spectral_cube/spectral_cube.py
index 7312eb5d8..8bfb4fa04 100644
--- a/spectral_cube/spectral_cube.py
+++ b/spectral_cube/spectral_cube.py
@@ -824,7 +824,7 @@ def argmax_world(self, axis, **kwargs):
             The axis to return the peak location along. e.g., `axis=0`
             will return the value of the spectral axis at the peak value.
         kwargs : dict
-             Passed to `~SpectralCube.argmax`.
+            Passed to `~SpectralCube.argmax`.
         '''
 
         argmax_plane = self.argmax(axis=axis, **kwargs)
@@ -850,13 +850,14 @@ def argmin_world(self, axis, **kwargs):
         '''
         Return the spatial or spectral index of the minimum value
         along a line of sight.
+
         Parameters
         ----------
         axis : int
             The axis to return the peak location along. e.g., `axis=0`
             will return the value of the spectral axis at the peak value.
         kwargs : dict
-             Passed to `~SpectralCube.argmax`.
+            Passed to `~SpectralCube.argmax`.
         '''
 
         argmin_plane = self.argmin(axis=axis, **kwargs)

From bdebe9df838d431f9d76a17331071123fd932eba Mon Sep 17 00:00:00 2001
From: e-koch <koch.eric.w@gmail.com>
Date: Thu, 12 Nov 2020 10:43:22 -0500
Subject: [PATCH 08/23] Fix docstring issues from @keflavich

---
 spectral_cube/cube_utils.py | 13 ++++++++-----
 1 file changed, 8 insertions(+), 5 deletions(-)

diff --git a/spectral_cube/cube_utils.py b/spectral_cube/cube_utils.py
index 6863d039c..de3f0bc90 100644
--- a/spectral_cube/cube_utils.py
+++ b/spectral_cube/cube_utils.py
@@ -488,20 +488,23 @@ def convert_bunit(bunit):
     return unit
 
 
-def world_take_along_axis(cube, posn_plane, axis):
+def world_take_along_axis(cube, position_plane, axis):
     '''
     Apply a 2D plane of pixel positions along collapsed along a third axis
     to the equivalent WCS coordinates. For example, this will convert `argmax`
     along the spectral axis to the equivalent spectral value (e.g., velocity at
     peak intensity).
+
     Parameters
     ----------
     cube : SpectralCube
         A spectral cube.
-    posn_plane : 2D numpy.ndarray
-        2D array of pixel dimensions along `axis`.
+    position_plane : 2D numpy.ndarray
+        2D array of pixel positions along `axis`. For example, `position_plane` can
+        be the output of `argmax` or `argmin` along an axis.
     axis : int
-        The axis that `posn_plane` is collapsed along.
+        The axis that `position_plane` is collapsed along.
+
     Returns
     -------
     out : astropy.units.Quantity
@@ -522,7 +525,7 @@ def world_take_along_axis(cube, posn_plane, axis):
     plane_newaxis = tuple(plane_newaxis)
 
     out = np.take_along_axis(world_coords[world_newaxis],
-                             posn_plane[plane_newaxis], axis=axis)
+                             position_plane[plane_newaxis], axis=axis)
     out = out.squeeze()
 
     return out

From f6427de21055246ebc3fbcd88531cc6916b32b53 Mon Sep 17 00:00:00 2001
From: e-koch <koch.eric.w@gmail.com>
Date: Thu, 12 Nov 2020 10:54:12 -0500
Subject: [PATCH 09/23] Require the celestial axes to be aligned when using
 world_take_along_axis

---
 spectral_cube/cube_utils.py | 8 +++++++-
 1 file changed, 7 insertions(+), 1 deletion(-)

diff --git a/spectral_cube/cube_utils.py b/spectral_cube/cube_utils.py
index de3f0bc90..8764d5edb 100644
--- a/spectral_cube/cube_utils.py
+++ b/spectral_cube/cube_utils.py
@@ -12,9 +12,10 @@
 import numpy as np
 from astropy.wcs import (WCSSUB_SPECTRAL, WCSSUB_LONGITUDE, WCSSUB_LATITUDE)
 from . import wcs_utils
-from .utils import FITSWarning, AstropyUserWarning
+from .utils import FITSWarning, AstropyUserWarning, WCSCelestialError
 from astropy import log
 from astropy.io import fits
+from astropy.wcs.utils import is_proj_plane_distorted
 from astropy.io.fits import BinTableHDU, Column
 from astropy import units as u
 import itertools
@@ -510,6 +511,11 @@ def world_take_along_axis(cube, position_plane, axis):
     out : astropy.units.Quantity
         2D array of WCS coordinates.
     '''
+
+    if is_proj_plane_distorted(cube.wcs):
+        raise WCSCelestialError("world_take_along_axis requires the celestial axes"
+                                " to be aligned along image axes.")
+
     # Get 1D slice along that axis.
     world_slice = [0, 0]
     world_slice.insert(axis, slice(None))

From 94c999dbe8ea67735e54c5e3789eed8465c5d577 Mon Sep 17 00:00:00 2001
From: e-koch <koch.eric.w@gmail.com>
Date: Thu, 12 Nov 2020 10:55:40 -0500
Subject: [PATCH 10/23] Test is_proj_plane_distorted before computing the
 argmax/argmin for the cube

---
 spectral_cube/spectral_cube.py | 8 ++++++++
 1 file changed, 8 insertions(+)

diff --git a/spectral_cube/spectral_cube.py b/spectral_cube/spectral_cube.py
index 8bfb4fa04..53c683d9b 100644
--- a/spectral_cube/spectral_cube.py
+++ b/spectral_cube/spectral_cube.py
@@ -827,6 +827,10 @@ def argmax_world(self, axis, **kwargs):
             Passed to `~SpectralCube.argmax`.
         '''
 
+        if is_proj_plane_distorted(self.wcs):
+            raise WCSCelestialError("argmax_world requires the celestial axes"
+                                    " to be aligned along image axes.")
+
         argmax_plane = self.argmax(axis=axis, **kwargs)
 
         # Convert to WCS coordinates.
@@ -860,6 +864,10 @@ def argmin_world(self, axis, **kwargs):
             Passed to `~SpectralCube.argmax`.
         '''
 
+        if is_proj_plane_distorted(self.wcs):
+            raise WCSCelestialError("argmax_world requires the celestial axes"
+                                    " to be aligned along image axes.")
+
         argmin_plane = self.argmin(axis=axis, **kwargs)
 
         # Convert to WCS coordinates.

From fb9bd05ce045e92ed895653f299e14c29411a278 Mon Sep 17 00:00:00 2001
From: e-koch <koch.eric.w@gmail.com>
Date: Thu, 12 Nov 2020 10:56:10 -0500
Subject: [PATCH 11/23] max -> min

---
 spectral_cube/spectral_cube.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/spectral_cube/spectral_cube.py b/spectral_cube/spectral_cube.py
index 53c683d9b..0a68d9306 100644
--- a/spectral_cube/spectral_cube.py
+++ b/spectral_cube/spectral_cube.py
@@ -861,11 +861,11 @@ def argmin_world(self, axis, **kwargs):
             The axis to return the peak location along. e.g., `axis=0`
             will return the value of the spectral axis at the peak value.
         kwargs : dict
-            Passed to `~SpectralCube.argmax`.
+            Passed to `~SpectralCube.argmin`.
         '''
 
         if is_proj_plane_distorted(self.wcs):
-            raise WCSCelestialError("argmax_world requires the celestial axes"
+            raise WCSCelestialError("argmin_world requires the celestial axes"
                                     " to be aligned along image axes.")
 
         argmin_plane = self.argmin(axis=axis, **kwargs)

From 6e1e6d50de3f2950579ff7d656f7a4708bc4e389 Mon Sep 17 00:00:00 2001
From: e-koch <koch.eric.w@gmail.com>
Date: Thu, 12 Nov 2020 12:49:48 -0500
Subject: [PATCH 12/23] Handle finding where image -> WCS axes are correlated
 and argmax/min_world cannot handle

---
 spectral_cube/cube_utils.py    |  2 +-
 spectral_cube/spectral_cube.py |  4 ++--
 spectral_cube/wcs_utils.py     | 21 +++++++++++++++++++++
 3 files changed, 24 insertions(+), 3 deletions(-)

diff --git a/spectral_cube/cube_utils.py b/spectral_cube/cube_utils.py
index 8764d5edb..e38feb984 100644
--- a/spectral_cube/cube_utils.py
+++ b/spectral_cube/cube_utils.py
@@ -512,7 +512,7 @@ def world_take_along_axis(cube, position_plane, axis):
         2D array of WCS coordinates.
     '''
 
-    if is_proj_plane_distorted(cube.wcs):
+    if wcs_utils.is_pixel_axis_to_wcs_correlated(cube.wcs):
         raise WCSCelestialError("world_take_along_axis requires the celestial axes"
                                 " to be aligned along image axes.")
 
diff --git a/spectral_cube/spectral_cube.py b/spectral_cube/spectral_cube.py
index 0a68d9306..11a1a6bfa 100644
--- a/spectral_cube/spectral_cube.py
+++ b/spectral_cube/spectral_cube.py
@@ -827,7 +827,7 @@ def argmax_world(self, axis, **kwargs):
             Passed to `~SpectralCube.argmax`.
         '''
 
-        if is_proj_plane_distorted(self.wcs):
+        if wcs_utils.is_pixel_axis_to_wcs_correlated(self.wcs):
             raise WCSCelestialError("argmax_world requires the celestial axes"
                                     " to be aligned along image axes.")
 
@@ -864,7 +864,7 @@ def argmin_world(self, axis, **kwargs):
             Passed to `~SpectralCube.argmin`.
         '''
 
-        if is_proj_plane_distorted(self.wcs):
+        if wcs_utils.is_pixel_axis_to_wcs_correlated(self.wcs):
             raise WCSCelestialError("argmin_world requires the celestial axes"
                                     " to be aligned along image axes.")
 
diff --git a/spectral_cube/wcs_utils.py b/spectral_cube/wcs_utils.py
index 28594712e..2d451ef59 100644
--- a/spectral_cube/wcs_utils.py
+++ b/spectral_cube/wcs_utils.py
@@ -483,6 +483,27 @@ def diagonal_wcs_to_cdelt(mywcs):
     return mywcs
 
 
+def is_pixel_axis_to_wcs_correlated(cube, axis):
+    """
+    Check if the chosen pixel axis correlates to other WCS axes. This tests
+    whether the pixel axis is correlated only to 1 WCS axis and can be
+    considered independent of the others.
+    """
+
+    axis_corr_matrix = cube.wcs.axis_correlation_matrix
+
+    wcs_axis = cube.ndim - (axis + 1)
+
+    wcs_axis_correlations = axis_corr_matrix[:, wcs_axis]
+
+    # The trace is always one. Correlations with other axes will give
+    # a sum > 1
+    if wcs_axis_correlations.sum() > 1:
+        return True
+
+    return False
+
+
 def find_spatial_pixel_index(cube, xlo, xhi, ylo, yhi):
     '''
     Given low and high cuts, return the pixel coordinates for a rectangular

From 5e8428e49e63483bc4ea3d88e8753e206bccf50a Mon Sep 17 00:00:00 2001
From: e-koch <koch.eric.w@gmail.com>
Date: Thu, 12 Nov 2020 14:29:47 -0500
Subject: [PATCH 13/23] Fix calls to check for correlated WCS axes from image

---
 spectral_cube/cube_utils.py    | 2 +-
 spectral_cube/spectral_cube.py | 7 ++++---
 spectral_cube/wcs_utils.py     | 6 +++---
 3 files changed, 8 insertions(+), 7 deletions(-)

diff --git a/spectral_cube/cube_utils.py b/spectral_cube/cube_utils.py
index e38feb984..b4cd6b1ea 100644
--- a/spectral_cube/cube_utils.py
+++ b/spectral_cube/cube_utils.py
@@ -512,7 +512,7 @@ def world_take_along_axis(cube, position_plane, axis):
         2D array of WCS coordinates.
     '''
 
-    if wcs_utils.is_pixel_axis_to_wcs_correlated(cube.wcs):
+    if wcs_utils.is_pixel_axis_to_wcs_correlated(cube.wcs, axis):
         raise WCSCelestialError("world_take_along_axis requires the celestial axes"
                                 " to be aligned along image axes.")
 
diff --git a/spectral_cube/spectral_cube.py b/spectral_cube/spectral_cube.py
index 11a1a6bfa..2b129d742 100644
--- a/spectral_cube/spectral_cube.py
+++ b/spectral_cube/spectral_cube.py
@@ -51,7 +51,8 @@
                     UnsupportedIterationStrategyWarning, WCSMismatchWarning,
                     NotImplementedWarning, SliceWarning, SmoothingWarning,
                     StokesWarning, ExperimentalImplementationWarning,
-                    BeamAverageWarning, NonFiniteBeamsWarning, BeamWarning)
+                    BeamAverageWarning, NonFiniteBeamsWarning, BeamWarning,
+                    WCSCelestialError)
 from .spectral_axis import (determine_vconv_from_ctype, get_rest_value_from_wcs,
                             doppler_beta, doppler_gamma, doppler_z)
 from .io.core import SpectralCubeRead, SpectralCubeWrite
@@ -827,7 +828,7 @@ def argmax_world(self, axis, **kwargs):
             Passed to `~SpectralCube.argmax`.
         '''
 
-        if wcs_utils.is_pixel_axis_to_wcs_correlated(self.wcs):
+        if wcs_utils.is_pixel_axis_to_wcs_correlated(self.wcs, axis):
             raise WCSCelestialError("argmax_world requires the celestial axes"
                                     " to be aligned along image axes.")
 
@@ -864,7 +865,7 @@ def argmin_world(self, axis, **kwargs):
             Passed to `~SpectralCube.argmin`.
         '''
 
-        if wcs_utils.is_pixel_axis_to_wcs_correlated(self.wcs):
+        if wcs_utils.is_pixel_axis_to_wcs_correlated(self.wcs, axis):
             raise WCSCelestialError("argmin_world requires the celestial axes"
                                     " to be aligned along image axes.")
 
diff --git a/spectral_cube/wcs_utils.py b/spectral_cube/wcs_utils.py
index 2d451ef59..0b6c4458c 100644
--- a/spectral_cube/wcs_utils.py
+++ b/spectral_cube/wcs_utils.py
@@ -483,16 +483,16 @@ def diagonal_wcs_to_cdelt(mywcs):
     return mywcs
 
 
-def is_pixel_axis_to_wcs_correlated(cube, axis):
+def is_pixel_axis_to_wcs_correlated(mywcs, axis):
     """
     Check if the chosen pixel axis correlates to other WCS axes. This tests
     whether the pixel axis is correlated only to 1 WCS axis and can be
     considered independent of the others.
     """
 
-    axis_corr_matrix = cube.wcs.axis_correlation_matrix
+    axis_corr_matrix = mywcs.axis_correlation_matrix
 
-    wcs_axis = cube.ndim - (axis + 1)
+    wcs_axis = mywcs.world_n_dim - (axis + 1)
 
     wcs_axis_correlations = axis_corr_matrix[:, wcs_axis]
 

From 5667fe314265750c28ad8b7767001a9cf51ac56e Mon Sep 17 00:00:00 2001
From: e-koch <koch.eric.w@gmail.com>
Date: Thu, 12 Nov 2020 14:37:31 -0500
Subject: [PATCH 14/23] Add more description to the correlated WCS axes check

---
 spectral_cube/wcs_utils.py | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/spectral_cube/wcs_utils.py b/spectral_cube/wcs_utils.py
index 0b6c4458c..f8e090627 100644
--- a/spectral_cube/wcs_utils.py
+++ b/spectral_cube/wcs_utils.py
@@ -492,11 +492,14 @@ def is_pixel_axis_to_wcs_correlated(mywcs, axis):
 
     axis_corr_matrix = mywcs.axis_correlation_matrix
 
+    # Map from numpy axis to WCS axis
     wcs_axis = mywcs.world_n_dim - (axis + 1)
 
+    # Grab the row along the given spatial axis. This slice is along the WCS axes
     wcs_axis_correlations = axis_corr_matrix[:, wcs_axis]
 
-    # The trace is always one. Correlations with other axes will give
+    # The image axis should always be correlated to at least 1 WCS axis.
+    # i.e., the diagonal term is one in the matrix. Correlations with other axes will give
     # a sum > 1
     if wcs_axis_correlations.sum() > 1:
         return True

From f37b2b3646f5134050dc3c0122f67f43fda92ae5 Mon Sep 17 00:00:00 2001
From: e-koch <koch.eric.w@gmail.com>
Date: Thu, 19 Nov 2020 10:13:44 -0500
Subject: [PATCH 15/23] Split arg world tests into own function; check for
 exception when taking along celestial axis

---
 spectral_cube/tests/test_spectral_cube.py | 29 +++++++++++++++++++----
 1 file changed, 25 insertions(+), 4 deletions(-)

diff --git a/spectral_cube/tests/test_spectral_cube.py b/spectral_cube/tests/test_spectral_cube.py
index 7d1a2aaad..384a51d85 100644
--- a/spectral_cube/tests/test_spectral_cube.py
+++ b/spectral_cube/tests/test_spectral_cube.py
@@ -706,14 +706,11 @@ def test_percentile(self, pct, iterate_rays, use_dask):
         self.c = self.d = None
 
     @pytest.mark.parametrize('method', ('sum', 'min', 'max', 'std', 'mad_std',
-                                        'median', 'argmin', 'argmax',
-                                        'argmax_world', 'argmin_world'))
+                                        'median', 'argmin', 'argmax'))
     def test_transpose(self, method, data_adv, data_vad, use_dask):
         c1, d1 = cube_and_raw(data_adv, use_dask=use_dask)
         c2, d2 = cube_and_raw(data_vad, use_dask=use_dask)
         for axis in [None, 0, 1, 2]:
-            if 'world' in method and axis == None:
-                continue
             assert_allclose(getattr(c1, method)(axis=axis),
                             getattr(c2, method)(axis=axis))
             if not use_dask:
@@ -722,6 +719,30 @@ def test_transpose(self, method, data_adv, data_vad, use_dask):
                                 getattr(c2, method)(axis=axis, progressbar=True))
         self.c = self.d = None
 
+    @pytest.mark.parametrize('method', ('argmax_world', 'argmin_world'))
+    def test_transpose_arg_world(self, method, data_adv, data_vad, use_dask):
+        c1, d1 = cube_and_raw(data_adv, use_dask=use_dask)
+        c2, d2 = cube_and_raw(data_vad, use_dask=use_dask)
+
+        # The spectral axis should work in all of these test cases.
+        axis = 0
+        assert_allclose(getattr(c1, method)(axis=axis),
+                        getattr(c2, method)(axis=axis))
+        if not use_dask:
+            # check that all these accept progressbar kwargs
+            assert_allclose(getattr(c1, method)(axis=axis, progressbar=True),
+                            getattr(c2, method)(axis=axis, progressbar=True))
+
+        # But the spatial axes should fail since the pixel axes are correlated to
+        # the WCS celestial axes. Currently this will happen for ALL celestial axes.
+        for axis in [1, 2]:
+
+            with pytest.raises(utils.WCSCelestialError,
+                               match=re.escape("world_take_along_axis requires the celestial axes")):
+
+                assert_allclose(getattr(c1, method)(axis=axis),
+                                getattr(c2, method)(axis=axis))
+
     @pytest.mark.parametrize('method', ('argmax_world', 'argmin_world'))
     def test_arg_world(self, method, data_adv, use_dask):
         c1, d1 = cube_and_raw(data_adv, use_dask=use_dask)

From eb5651296d0c722a811fa13381b5df10e5c5e943 Mon Sep 17 00:00:00 2001
From: e-koch <koch.eric.w@gmail.com>
Date: Thu, 19 Nov 2020 10:37:57 -0500
Subject: [PATCH 16/23] Error message changes with the method

---
 spectral_cube/tests/test_spectral_cube.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/spectral_cube/tests/test_spectral_cube.py b/spectral_cube/tests/test_spectral_cube.py
index 384a51d85..a9d0af989 100644
--- a/spectral_cube/tests/test_spectral_cube.py
+++ b/spectral_cube/tests/test_spectral_cube.py
@@ -738,7 +738,7 @@ def test_transpose_arg_world(self, method, data_adv, data_vad, use_dask):
         for axis in [1, 2]:
 
             with pytest.raises(utils.WCSCelestialError,
-                               match=re.escape("world_take_along_axis requires the celestial axes")):
+                               match=re.escape(f"{method} requires the celestial axes")):
 
                 assert_allclose(getattr(c1, method)(axis=axis),
                                 getattr(c2, method)(axis=axis))

From 4c686317151ad377c653aff0f5235df1e6e4cc29 Mon Sep 17 00:00:00 2001
From: e-koch <koch.eric.w@gmail.com>
Date: Fri, 27 Nov 2020 11:52:20 -0500
Subject: [PATCH 17/23] Missing setting attributes back to None; causes
 InvocationError

---
 spectral_cube/tests/test_spectral_cube.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/spectral_cube/tests/test_spectral_cube.py b/spectral_cube/tests/test_spectral_cube.py
index a9d0af989..858ec463b 100644
--- a/spectral_cube/tests/test_spectral_cube.py
+++ b/spectral_cube/tests/test_spectral_cube.py
@@ -743,6 +743,8 @@ def test_transpose_arg_world(self, method, data_adv, data_vad, use_dask):
                 assert_allclose(getattr(c1, method)(axis=axis),
                                 getattr(c2, method)(axis=axis))
 
+        self.c = self.d = None
+
     @pytest.mark.parametrize('method', ('argmax_world', 'argmin_world'))
     def test_arg_world(self, method, data_adv, use_dask):
         c1, d1 = cube_and_raw(data_adv, use_dask=use_dask)

From 04f8256616dd81234dd268bcaebe25514f509b84 Mon Sep 17 00:00:00 2001
From: e-koch <koch.eric.w@gmail.com>
Date: Fri, 27 Nov 2020 12:23:38 -0500
Subject: [PATCH 18/23] Add explanation for why spatial axes cannot be used for
 argmax/min_world

---
 docs/moments.rst | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/docs/moments.rst b/docs/moments.rst
index acbb472da..05a2c8ca5 100644
--- a/docs/moments.rst
+++ b/docs/moments.rst
@@ -156,3 +156,7 @@ or for modeling. The `~spectral_cube.SpectralCube.argmax_world` and
 
 Along the spectral axis, `~spectral_cube.SpectralCube.argmax_world` creates the often used
 "velocity at peak intensity," which may also be called the "peak velocity."
+
+.. note:: `cube.argmax_world` and `cube.argmin_world` are currently only defined along the spectral axis,
+as the example above shows. This is because `argmax_world` and `argmin_world` operate along the pixel axes,
+but they are not independent in WCS coordinates due to the curvature of the sky.

From 184c51a1d0fc10898a40327251a5e296551c453f Mon Sep 17 00:00:00 2001
From: e-koch <koch.eric.w@gmail.com>
Date: Fri, 27 Nov 2020 12:25:42 -0500
Subject: [PATCH 19/23] Update CHANGE

---
 CHANGES.rst | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/CHANGES.rst b/CHANGES.rst
index 3d4e29f4f..726fe13e0 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -1,5 +1,8 @@
 0.6 (unreleased)
 ----------------
+- Add ``argmax_world`` and ``argmin_world`` to return the argmin/max position
+  in WCS coordinates. This is ONLY defined for independent WCS axes
+  (e.g., spectral) #680
 - Bugfix: subcube producing spatial offsets for large images #666
 - Switch to using standalone casa-formats-io package for reading in CASA
   .image files (this was split out from spectral-cube). #684

From 3b7a75d17f95f87f6613375eb2ebfe04115bbde6 Mon Sep 17 00:00:00 2001
From: e-koch <koch.eric.w@gmail.com>
Date: Fri, 27 Nov 2020 12:33:03 -0500
Subject: [PATCH 20/23] Fix docs formatting

---
 docs/moments.rst | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/docs/moments.rst b/docs/moments.rst
index 05a2c8ca5..0e92b1582 100644
--- a/docs/moments.rst
+++ b/docs/moments.rst
@@ -158,5 +158,5 @@ Along the spectral axis, `~spectral_cube.SpectralCube.argmax_world` creates the
 "velocity at peak intensity," which may also be called the "peak velocity."
 
 .. note:: `cube.argmax_world` and `cube.argmin_world` are currently only defined along the spectral axis,
-as the example above shows. This is because `argmax_world` and `argmin_world` operate along the pixel axes,
-but they are not independent in WCS coordinates due to the curvature of the sky.
+          as the example above shows. This is because `argmax_world` and `argmin_world` operate along the
+          pixel axes, but they are not independent in WCS coordinates due to the curvature of the sky.

From 3fe4b83a6675a02b66521f7527593d2b1057cb10 Mon Sep 17 00:00:00 2001
From: Eric Koch <koch.eric.w@gmail.com>
Date: Thu, 3 Dec 2020 16:15:12 -0500
Subject: [PATCH 21/23] Update spectral_cube/cube_utils.py

Co-authored-by: Adam Ginsburg <keflavich@gmail.com>
---
 spectral_cube/cube_utils.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/spectral_cube/cube_utils.py b/spectral_cube/cube_utils.py
index b4cd6b1ea..9e1656ae4 100644
--- a/spectral_cube/cube_utils.py
+++ b/spectral_cube/cube_utils.py
@@ -491,8 +491,8 @@ def convert_bunit(bunit):
 
 def world_take_along_axis(cube, position_plane, axis):
     '''
-    Apply a 2D plane of pixel positions along collapsed along a third axis
-    to the equivalent WCS coordinates. For example, this will convert `argmax`
+    Convert a 2D plane of pixel positions to the equivalent WCS coordinates.
+    For example, this will convert `argmax`
     along the spectral axis to the equivalent spectral value (e.g., velocity at
     peak intensity).
 

From 2119822a2632825c4992f4e99695b642ce031105 Mon Sep 17 00:00:00 2001
From: Eric Koch <koch.eric.w@gmail.com>
Date: Thu, 3 Dec 2020 16:15:21 -0500
Subject: [PATCH 22/23] Update docs/moments.rst

Co-authored-by: Adam Ginsburg <keflavich@gmail.com>
---
 docs/moments.rst | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/moments.rst b/docs/moments.rst
index 0e92b1582..9af2fea6f 100644
--- a/docs/moments.rst
+++ b/docs/moments.rst
@@ -147,7 +147,7 @@ return the pixel positions of the max/min along that axis:
     >>> argmin_map = cube.argmin(axis=0)  # doctest: +SKIP
 
 These maps are useful for identifying where signal is located within the spectral cube, however,
-it is more useful to return the WCS values of those pixels for comparisons with other data sets
+it is more useful to return the WCS values of those pixels for comparison with other data sets
 or for modeling. The `~spectral_cube.SpectralCube.argmax_world` and
 `~spectral_cube.SpectralCube.argmin_world` should be used in this case:
 

From 1affb34201d7ec5035af1cb8abda07b2cfaa6257 Mon Sep 17 00:00:00 2001
From: e-koch <koch.eric.w@gmail.com>
Date: Thu, 3 Dec 2020 18:17:48 -0500
Subject: [PATCH 23/23] Remove duplicated code per @keflavich's comment

---
 spectral_cube/spectral_cube.py | 72 ++++++++++++++++------------------
 1 file changed, 34 insertions(+), 38 deletions(-)

diff --git a/spectral_cube/spectral_cube.py b/spectral_cube/spectral_cube.py
index 2b129d742..ee0abc120 100644
--- a/spectral_cube/spectral_cube.py
+++ b/spectral_cube/spectral_cube.py
@@ -813,29 +813,28 @@ def argmin(self, axis=None, how='auto', **kwargs):
                                          reduce=False, projection=False,
                                          how=how, axis=axis, **kwargs)
 
-    @warn_slow
-    def argmax_world(self, axis, **kwargs):
+    def _argmaxmin_world(self, axis, method, **kwargs):
         '''
-        Return the spatial or spectral index of the maximum value
-        along a line of sight.
-
-        Parameters
-        ----------
-        axis : int
-            The axis to return the peak location along. e.g., `axis=0`
-            will return the value of the spectral axis at the peak value.
-        kwargs : dict
-            Passed to `~SpectralCube.argmax`.
+        Return the spatial or spectral index of the maximum or minimum value.
+        Use `argmax_world` and `argmin_world` directly.
         '''
 
-        if wcs_utils.is_pixel_axis_to_wcs_correlated(self.wcs, axis):
-            raise WCSCelestialError("argmax_world requires the celestial axes"
-                                    " to be aligned along image axes.")
+        operation_name = '{}_world'.format(method)
 
-        argmax_plane = self.argmax(axis=axis, **kwargs)
+        if wcs_utils.is_pixel_axis_to_wcs_correlated(self.wcs, axis):
+            raise WCSCelestialError("{} requires the celestial axes"
+                                    " to be aligned along image axes."
+                                    .format(operation_name))
+
+        if method == 'argmin':
+            arg_pixel_plane = self.argmin(axis=axis, **kwargs)
+        elif method == 'argmax':
+            arg_pixel_plane = self.argmax(axis=axis, **kwargs)
+        else:
+            raise ValueError("`method` must be 'argmin' or 'argmax'")
 
         # Convert to WCS coordinates.
-        out = cube_utils.world_take_along_axis(self, argmax_plane, axis)
+        out = cube_utils.world_take_along_axis(self, arg_pixel_plane, axis)
 
         # Compute whether the mask has any valid data along `axis`
         collapsed_mask = self.mask.include().any(axis=axis)
@@ -851,9 +850,9 @@ def argmax_world(self, axis, **kwargs):
                           unit=out.unit, header=self._nowcs_header)
 
     @warn_slow
-    def argmin_world(self, axis, **kwargs):
+    def argmax_world(self, axis, **kwargs):
         '''
-        Return the spatial or spectral index of the minimum value
+        Return the spatial or spectral index of the maximum value
         along a line of sight.
 
         Parameters
@@ -862,30 +861,27 @@ def argmin_world(self, axis, **kwargs):
             The axis to return the peak location along. e.g., `axis=0`
             will return the value of the spectral axis at the peak value.
         kwargs : dict
-            Passed to `~SpectralCube.argmin`.
+            Passed to `~SpectralCube.argmax`.
         '''
 
-        if wcs_utils.is_pixel_axis_to_wcs_correlated(self.wcs, axis):
-            raise WCSCelestialError("argmin_world requires the celestial axes"
-                                    " to be aligned along image axes.")
-
-        argmin_plane = self.argmin(axis=axis, **kwargs)
-
-        # Convert to WCS coordinates.
-        out = cube_utils.world_take_along_axis(self, argmin_plane, axis)
-
-        # Compute whether the mask has any valid data along `axis`
-        collapsed_mask = self.mask.include().any(axis=axis)
-        out[~collapsed_mask] = np.NaN
+        return self._argmaxmin_world(axis, 'argmax', **kwargs)
 
-        # Return a Projection.
-        new_wcs = wcs_utils.drop_axis(self._wcs, np2wcs[axis])
+    @warn_slow
+    def argmin_world(self, axis, **kwargs):
+        '''
+        Return the spatial or spectral index of the minimum value
+        along a line of sight.
 
-        meta = {'collapse_axis': axis}
-        meta.update(self._meta)
+        Parameters
+        ----------
+        axis : int
+            The axis to return the peak location along. e.g., `axis=0`
+            will return the value of the spectral axis at the peak value.
+        kwargs : dict
+            Passed to `~SpectralCube.argmin`.
+        '''
 
-        return Projection(out, copy=False, wcs=new_wcs, meta=meta,
-                          unit=out.unit, header=self._nowcs_header)
+        return self._argmaxmin_world(axis, 'argmin', **kwargs)
 
     def chunked(self, chunksize=1000):
         """