From f3a72dc279cbe4e61cd2f8ab09e6d63870a9849b Mon Sep 17 00:00:00 2001
From: "Adam Ginsburg (keflavich)" <keflavich@gmail.com>
Date: Sat, 3 Jun 2017 16:44:09 -0600
Subject: [PATCH 01/12] slice the third axis, don't drop it

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

diff --git a/spectral_cube/spectral_cube.py b/spectral_cube/spectral_cube.py
index 7c5420a42..c23700fe9 100644
--- a/spectral_cube/spectral_cube.py
+++ b/spectral_cube/spectral_cube.py
@@ -1373,7 +1373,10 @@ def moment(self, order=0, axis=0, how='auto'):
         if order == 1 and axis == 0:
             out += self.world[0, :, :][0]
 
-        new_wcs = wcs_utils.drop_axis(self._wcs, np2wcs[axis])
+        view = [slice(None) if ax!=axis else slice(self.shape[axis]//2,
+                                                   self.shape[axis]//2+1)
+                for ax in range(self.ndim)]
+        new_wcs = wcs_utils.slice_wcs(self._wcs, view)
 
         meta = {'moment_order': order,
                 'moment_axis': axis,

From 14f9186fa3de3626f3d45cf824b95fb3018c1733 Mon Sep 17 00:00:00 2001
From: "Adam Ginsburg (keflavich)" <keflavich@gmail.com>
Date: Sat, 3 Jun 2017 17:13:10 -0600
Subject: [PATCH 02/12] add a new tool to "drop" an axis by slicing it to 1
 pixel wide

---
 spectral_cube/spectral_cube.py | 29 +++++++++++++++---
 spectral_cube/wcs_utils.py     | 55 ++++++++++++++++++++++++++++++++++
 2 files changed, 80 insertions(+), 4 deletions(-)

diff --git a/spectral_cube/spectral_cube.py b/spectral_cube/spectral_cube.py
index c23700fe9..0ced72dad 100644
--- a/spectral_cube/spectral_cube.py
+++ b/spectral_cube/spectral_cube.py
@@ -268,6 +268,8 @@ def apply_numpy_function(self, function, fill=np.nan,
                              unit=None,
                              check_endian=False,
                              progressbar=False,
+                             dropped_axis_slice_position='middle',
+                             dropped_axis_cdelt='same',
                              **kwargs):
         """
         Apply a numpy function to the cube
@@ -306,6 +308,15 @@ def apply_numpy_function(self, function, fill=np.nan,
         progressbar : bool
             Show a progressbar while iterating over the slices through the
             cube?
+        dropped_axis_slice_position : 'middle', 'start', 'end'
+            If an axis is being dropped, where should the WCS say the
+            projection is?  It can be at the start, middle, or end of the
+            axis.
+        dropped_axis_cdelt : 'same', 'full_range', or value
+            If an axis is being dropped, what should the new CDELT be?  For an
+            integral, for example, one might want the value to be the full
+            range.  For a slice, it should stay the same.  For something like
+            min or max, it might be zero.
         kwargs : dict
             Passed to the numpy function.
 
@@ -375,7 +386,14 @@ def apply_numpy_function(self, function, fill=np.nan,
                     return out
 
             else:
-                new_wcs = wcs_utils.drop_axis(self._wcs, np2wcs[axis])
+
+                new_wcs = wcs_utils.drop_axis_by_slicing(self._wcs,
+                                                         self.shape,
+                                                         axis,
+                                                         dropped_axis_cdelt=dropped_axis_cdelt,
+                                                         dropped_axis_slice_position=dropped_axis_slice_position,
+                                                        )
+
                 header = self._nowcs_header
 
                 return Projection(out, copy=False, wcs=new_wcs, meta=meta,
@@ -1373,8 +1391,8 @@ def moment(self, order=0, axis=0, how='auto'):
         if order == 1 and axis == 0:
             out += self.world[0, :, :][0]
 
-        view = [slice(None) if ax!=axis else slice(self.shape[axis]//2,
-                                                   self.shape[axis]//2+1)
+        view = [slice(None) if ax!=axis else slice(self.shape[ax]//2,
+                                                   self.shape[ax]//2+1)
                 for ax in range(self.ndim)]
         new_wcs = wcs_utils.slice_wcs(self._wcs, view)
 
@@ -2667,7 +2685,10 @@ def __getitem__(self, view):
                                     meta=meta)
 
             # only one element, so drop an axis
-            newwcs = wcs_utils.drop_axis(self._wcs, intslices[0])
+            view = [slice(None) if ax!=intslices[0]
+                    else slice(self.shape[ax]//2, self.shape[ax]//2+1)
+                    for ax in range(self.ndim)]
+            newwcs = wcs_utils.slice_wcs(self._wcs, view)
             header = self._nowcs_header
 
             # Slice objects know how to parse Beam objects stored in the
diff --git a/spectral_cube/wcs_utils.py b/spectral_cube/wcs_utils.py
index 2b780119d..7d5666592 100644
--- a/spectral_cube/wcs_utils.py
+++ b/spectral_cube/wcs_utils.py
@@ -420,3 +420,58 @@ def diagonal_wcs_to_cdelt(mywcs):
         del mywcs.wcs.cd
         mywcs.wcs.cdelt = cdelt
     return mywcs
+
+def drop_axis_by_slicing(mywcs, shape, dropped_axis,
+                         dropped_axis_slice_position='middle',
+                         dropped_axis_cdelt='same',
+                         convert_misaligned_to_offset=True,
+                        ):
+    """
+    Parameters
+    ----------
+        dropped_axis_slice_position : 'middle', 'start', 'end'
+            If an axis is being dropped, where should the WCS say the
+            projection is?  It can be at the start, middle, or end of the
+            axis.
+        dropped_axis_cdelt : 'same', 'full_range', or value
+            If an axis is being dropped, what should the new CDELT be?  For an
+            integral, for example, one might want the value to be the full
+            range.  For a slice, it should stay the same.  For something like
+            min or max, it might be zero.
+        convert_misaligned_to_offset : bool
+            If the axes are misaligned, it is not possible to "drop" an axis.
+            In this case, a generic "offset axis" will be returned.
+    """
+    ndim = len(shape)
+
+    if dropped_axis_slice_position == 'middle':
+        dropped_axis_slice_position = shape[dropped_axis]//2
+    elif dropped_axis_slice_position == 'start':
+        dropped_axis_slice_position = 0
+    elif dropped_axis_slice_position == 'end':
+        dropped_axis_slice_position = shape[dropped_axis] - 1
+
+    dropax_slice = slice(dropped_axis_slice_position,
+                         dropped_axis_slice_position+1)
+
+    view = [slice(None) if ax!=dropped_axis else dropax_slice
+            for ax in range(ndim)]
+
+    crpix_new = [0 if ax!=dropped_axis else dropped_axis_slice_position
+                 for ax in range(ndim)]
+    new_crval = mywcs.wcs_pix2world(crpix_new, 0)[dropped_axis]
+
+    result = slice_wcs(mywcs, view, shape=shape)
+
+    result.wcs.crval[dropped_axis] = new_crval
+    result.wcs.crpix[dropped_axis] = 1
+
+    if dropped_axis_cdelt == 'same':
+        dropped_axis_cdelt = mywcs.wcs.cdelt[dropped_axis]
+    elif dropped_axis_cdelt == 'full_range':
+        ref_pixels = [(0,0) if ax!=dropped_axis else (0, shape[dropped_axis]-1)
+                      for ax in range(ndim)]
+        dropped_axis_cdelt = mywcs.wcs_pix2world(ref_pixels, 0)[dropped_axis]
+    result.wcs.cdelt[dropped_axis] = dropped_axis_cdelt
+
+    return result

From ab2be8c1150fc8664766bbe1a5ba0dbff3c5049a Mon Sep 17 00:00:00 2001
From: "Adam Ginsburg (keflavich)" <keflavich@gmail.com>
Date: Sat, 3 Jun 2017 17:19:38 -0600
Subject: [PATCH 03/12] bugfix because I indexed the WCS incorrectly, and first
 test...

---
 spectral_cube/tests/test_wcs_utils.py | 12 ++++++++++++
 spectral_cube/wcs_utils.py            |  2 +-
 2 files changed, 13 insertions(+), 1 deletion(-)

diff --git a/spectral_cube/tests/test_wcs_utils.py b/spectral_cube/tests/test_wcs_utils.py
index 241a560e2..d4cb49226 100644
--- a/spectral_cube/tests/test_wcs_utils.py
+++ b/spectral_cube/tests/test_wcs_utils.py
@@ -128,3 +128,15 @@ def test_strip_wcs():
     header2_stripped = strip_wcs_from_header(header2)
 
     assert header1_stripped == header2_stripped
+
+def test_drop_by_slice():
+
+    wcs = WCS(naxis=3)
+    wcs.wcs.crpix = [1., 1., 1.]
+    wcs.wcs.crval = [0., 0., 0.]
+    wcs.wcs.cdelt = [1., 1., 1.]
+    wcs.wcs.ctype = ['RA---TAN', 'DEC--TAN', 'FREQ']
+
+    newwcs = drop_axis_by_slicing(wcs, shape=[10,12,14], dropped_axis=0)
+
+    assert newwcs.wcs.crval[0] == 5
diff --git a/spectral_cube/wcs_utils.py b/spectral_cube/wcs_utils.py
index 7d5666592..844fb4b21 100644
--- a/spectral_cube/wcs_utils.py
+++ b/spectral_cube/wcs_utils.py
@@ -459,7 +459,7 @@ def drop_axis_by_slicing(mywcs, shape, dropped_axis,
 
     crpix_new = [0 if ax!=dropped_axis else dropped_axis_slice_position
                  for ax in range(ndim)]
-    new_crval = mywcs.wcs_pix2world(crpix_new, 0)[dropped_axis]
+    new_crval = mywcs.wcs_pix2world([crpix_new], 0)[0, dropped_axis]
 
     result = slice_wcs(mywcs, view, shape=shape)
 

From bec73849eb53aaefce15c918943c6caeb6f9a902 Mon Sep 17 00:00:00 2001
From: "Adam Ginsburg (keflavich)" <keflavich@gmail.com>
Date: Sat, 3 Jun 2017 17:34:49 -0600
Subject: [PATCH 04/12] used tests to debug new tool

---
 spectral_cube/tests/test_wcs_utils.py | 21 ++++++++++++++++++---
 spectral_cube/wcs_utils.py            |  8 +++++---
 2 files changed, 23 insertions(+), 6 deletions(-)

diff --git a/spectral_cube/tests/test_wcs_utils.py b/spectral_cube/tests/test_wcs_utils.py
index d4cb49226..da3479677 100644
--- a/spectral_cube/tests/test_wcs_utils.py
+++ b/spectral_cube/tests/test_wcs_utils.py
@@ -129,14 +129,29 @@ def test_strip_wcs():
 
     assert header1_stripped == header2_stripped
 
-def test_drop_by_slice():
+def test_drop_by_slice_middle():
 
     wcs = WCS(naxis=3)
     wcs.wcs.crpix = [1., 1., 1.]
     wcs.wcs.crval = [0., 0., 0.]
-    wcs.wcs.cdelt = [1., 1., 1.]
+    wcs.wcs.cdelt = [1e-5, 1e-5, 1e-5]
     wcs.wcs.ctype = ['RA---TAN', 'DEC--TAN', 'FREQ']
 
     newwcs = drop_axis_by_slicing(wcs, shape=[10,12,14], dropped_axis=0)
 
-    assert newwcs.wcs.crval[0] == 5
+    np.testing.assert_almost_equal(newwcs.wcs.crval[0], 5e-5)
+    assert all(newwcs.wcs.cdelt == [1e-5,1e-5,1e-5])
+
+def test_drop_by_slice_middle_fullrange():
+
+    wcs = WCS(naxis=3)
+    wcs.wcs.crpix = [1., 1., 1.]
+    wcs.wcs.crval = [0., 0., 0.]
+    wcs.wcs.cdelt = [1e-5, 1e-5, 1e-5]
+    wcs.wcs.ctype = ['RA---TAN', 'DEC--TAN', 'FREQ']
+
+    newwcs = drop_axis_by_slicing(wcs, shape=[10,12,14], dropped_axis=0,
+                                  dropped_axis_cdelt='full_range')
+
+    np.testing.assert_almost_equal(newwcs.wcs.crval[0], 5e-5)
+    np.testing.assert_almost_equal(newwcs.wcs.cdelt[0], 10e-5)
diff --git a/spectral_cube/wcs_utils.py b/spectral_cube/wcs_utils.py
index 844fb4b21..d31a5db5d 100644
--- a/spectral_cube/wcs_utils.py
+++ b/spectral_cube/wcs_utils.py
@@ -469,9 +469,11 @@ def drop_axis_by_slicing(mywcs, shape, dropped_axis,
     if dropped_axis_cdelt == 'same':
         dropped_axis_cdelt = mywcs.wcs.cdelt[dropped_axis]
     elif dropped_axis_cdelt == 'full_range':
-        ref_pixels = [(0,0) if ax!=dropped_axis else (0, shape[dropped_axis]-1)
-                      for ax in range(ndim)]
-        dropped_axis_cdelt = mywcs.wcs_pix2world(ref_pixels, 0)[dropped_axis]
+        ref_pixels = np.array(
+            [(0,0) if ax!=dropped_axis else (0, shape[dropped_axis])
+             for ax in range(ndim)])
+        refvals = mywcs.wcs_pix2world(ref_pixels.T, 0)[:, dropped_axis]
+        dropped_axis_cdelt = refvals[1]-refvals[0]
     result.wcs.cdelt[dropped_axis] = dropped_axis_cdelt
 
     return result

From f901a40ec0cecea266277fa57e2a62ecdd488b5a Mon Sep 17 00:00:00 2001
From: "Adam Ginsburg (keflavich)" <keflavich@gmail.com>
Date: Sat, 3 Jun 2017 17:48:28 -0600
Subject: [PATCH 05/12] improve test coverage and fix another edge case
 (literally - the end case)

---
 spectral_cube/tests/test_wcs_utils.py | 13 ++++++++++---
 spectral_cube/wcs_utils.py            | 27 ++++++++++++++-------------
 2 files changed, 24 insertions(+), 16 deletions(-)

diff --git a/spectral_cube/tests/test_wcs_utils.py b/spectral_cube/tests/test_wcs_utils.py
index da3479677..ddeb0c58f 100644
--- a/spectral_cube/tests/test_wcs_utils.py
+++ b/spectral_cube/tests/test_wcs_utils.py
@@ -2,6 +2,8 @@
 
 from astropy.io import fits
 
+import pytest
+
 from ..wcs_utils import *
 from . import path
 
@@ -129,7 +131,11 @@ def test_strip_wcs():
 
     assert header1_stripped == header2_stripped
 
-def test_drop_by_slice_middle():
+@pytest.mark.parametrize(('position', 'result'),
+                    (('start', 0.),
+                     ('middle', 5e-5),
+                     ('end', 10e-5)))
+def test_drop_by_slice(position, result):
 
     wcs = WCS(naxis=3)
     wcs.wcs.crpix = [1., 1., 1.]
@@ -137,9 +143,10 @@ def test_drop_by_slice_middle():
     wcs.wcs.cdelt = [1e-5, 1e-5, 1e-5]
     wcs.wcs.ctype = ['RA---TAN', 'DEC--TAN', 'FREQ']
 
-    newwcs = drop_axis_by_slicing(wcs, shape=[10,12,14], dropped_axis=0)
+    newwcs = drop_axis_by_slicing(wcs, shape=[10,12,14], dropped_axis=0,
+                                  dropped_axis_slice_position=position)
 
-    np.testing.assert_almost_equal(newwcs.wcs.crval[0], 5e-5)
+    np.testing.assert_almost_equal(newwcs.wcs.crval[0], result)
     assert all(newwcs.wcs.cdelt == [1e-5,1e-5,1e-5])
 
 def test_drop_by_slice_middle_fullrange():
diff --git a/spectral_cube/wcs_utils.py b/spectral_cube/wcs_utils.py
index d31a5db5d..a28fa56f0 100644
--- a/spectral_cube/wcs_utils.py
+++ b/spectral_cube/wcs_utils.py
@@ -429,18 +429,18 @@ def drop_axis_by_slicing(mywcs, shape, dropped_axis,
     """
     Parameters
     ----------
-        dropped_axis_slice_position : 'middle', 'start', 'end'
-            If an axis is being dropped, where should the WCS say the
-            projection is?  It can be at the start, middle, or end of the
-            axis.
-        dropped_axis_cdelt : 'same', 'full_range', or value
-            If an axis is being dropped, what should the new CDELT be?  For an
-            integral, for example, one might want the value to be the full
-            range.  For a slice, it should stay the same.  For something like
-            min or max, it might be zero.
-        convert_misaligned_to_offset : bool
-            If the axes are misaligned, it is not possible to "drop" an axis.
-            In this case, a generic "offset axis" will be returned.
+    dropped_axis_slice_position : 'middle', 'start', 'end'
+        If an axis is being dropped, where should the WCS say the
+        projection is?  It can be at the start, middle, or end of the
+        axis.
+    dropped_axis_cdelt : 'same', 'full_range', or value
+        If an axis is being dropped, what should the new CDELT be?  For an
+        integral, for example, one might want the value to be the full
+        range.  For a slice, it should stay the same.  For something like
+        min or max, it might be zero.
+    convert_misaligned_to_offset : bool
+        If the axes are misaligned, it is not possible to "drop" an axis.
+        In this case, a generic "offset axis" will be returned.
     """
     ndim = len(shape)
 
@@ -449,7 +449,7 @@ def drop_axis_by_slicing(mywcs, shape, dropped_axis,
     elif dropped_axis_slice_position == 'start':
         dropped_axis_slice_position = 0
     elif dropped_axis_slice_position == 'end':
-        dropped_axis_slice_position = shape[dropped_axis] - 1
+        dropped_axis_slice_position = shape[dropped_axis]
 
     dropax_slice = slice(dropped_axis_slice_position,
                          dropped_axis_slice_position+1)
@@ -474,6 +474,7 @@ def drop_axis_by_slicing(mywcs, shape, dropped_axis,
              for ax in range(ndim)])
         refvals = mywcs.wcs_pix2world(ref_pixels.T, 0)[:, dropped_axis]
         dropped_axis_cdelt = refvals[1]-refvals[0]
+
     result.wcs.cdelt[dropped_axis] = dropped_axis_cdelt
 
     return result

From 95f5eba69830606832481ea70cd622c24cf6fe06 Mon Sep 17 00:00:00 2001
From: "Adam Ginsburg (keflavich)" <keflavich@gmail.com>
Date: Sat, 3 Jun 2017 18:07:53 -0600
Subject: [PATCH 06/12] start dealing with some more error cases and move the
 dropped axis to the end

---
 spectral_cube/lower_dimensional_structures.py |  6 ------
 spectral_cube/tests/test_wcs_utils.py         | 17 +++++++++--------
 spectral_cube/wcs_utils.py                    |  4 ++++
 3 files changed, 13 insertions(+), 14 deletions(-)

diff --git a/spectral_cube/lower_dimensional_structures.py b/spectral_cube/lower_dimensional_structures.py
index e7a077919..b1f4bf74d 100644
--- a/spectral_cube/lower_dimensional_structures.py
+++ b/spectral_cube/lower_dimensional_structures.py
@@ -210,9 +210,6 @@ def __new__(cls, value, unit=None, dtype=None, copy=True, wcs=None,
         if np.asarray(value).ndim != 2:
             raise ValueError("value should be a 2-d array")
 
-        if wcs is not None and wcs.wcs.naxis != 2:
-            raise ValueError("wcs should have two dimension")
-
         self = u.Quantity.__new__(cls, value, unit=unit, dtype=dtype,
                                   copy=copy).view(cls)
         self._wcs = wcs
@@ -466,9 +463,6 @@ def __new__(cls, value, unit=None, dtype=None, copy=True, wcs=None,
         if np.asarray(value).ndim != 1:
             raise ValueError("value should be a 1-d array")
 
-        if wcs is not None and wcs.wcs.naxis != 1:
-            raise ValueError("wcs should have two dimension")
-
         self = u.Quantity.__new__(cls, value, unit=unit, dtype=dtype,
                                   copy=copy).view(cls)
         self._wcs = wcs
diff --git a/spectral_cube/tests/test_wcs_utils.py b/spectral_cube/tests/test_wcs_utils.py
index ddeb0c58f..43279a46a 100644
--- a/spectral_cube/tests/test_wcs_utils.py
+++ b/spectral_cube/tests/test_wcs_utils.py
@@ -132,22 +132,23 @@ def test_strip_wcs():
     assert header1_stripped == header2_stripped
 
 @pytest.mark.parametrize(('position', 'result'),
-                    (('start', 0.),
-                     ('middle', 5e-5),
-                     ('end', 10e-5)))
+                         (('start', 0.),
+                          ('middle', 5e-5),
+                          ('end', 10e-5)))
 def test_drop_by_slice(position, result):
 
     wcs = WCS(naxis=3)
     wcs.wcs.crpix = [1., 1., 1.]
     wcs.wcs.crval = [0., 0., 0.]
-    wcs.wcs.cdelt = [1e-5, 1e-5, 1e-5]
+    wcs.wcs.cdelt = [1e-5, 2e-5, 3e-5]
     wcs.wcs.ctype = ['RA---TAN', 'DEC--TAN', 'FREQ']
 
     newwcs = drop_axis_by_slicing(wcs, shape=[10,12,14], dropped_axis=0,
                                   dropped_axis_slice_position=position)
 
-    np.testing.assert_almost_equal(newwcs.wcs.crval[0], result)
-    assert all(newwcs.wcs.cdelt == [1e-5,1e-5,1e-5])
+    # drop-by-slicing moves axis to be last
+    np.testing.assert_almost_equal(newwcs.wcs.crval[2], result)
+    assert all(newwcs.wcs.cdelt == [2e-5,3e-5,1e-5])
 
 def test_drop_by_slice_middle_fullrange():
 
@@ -160,5 +161,5 @@ def test_drop_by_slice_middle_fullrange():
     newwcs = drop_axis_by_slicing(wcs, shape=[10,12,14], dropped_axis=0,
                                   dropped_axis_cdelt='full_range')
 
-    np.testing.assert_almost_equal(newwcs.wcs.crval[0], 5e-5)
-    np.testing.assert_almost_equal(newwcs.wcs.cdelt[0], 10e-5)
+    np.testing.assert_almost_equal(newwcs.wcs.crval[2], 5e-5)
+    np.testing.assert_almost_equal(newwcs.wcs.cdelt[2], 10e-5)
diff --git a/spectral_cube/wcs_utils.py b/spectral_cube/wcs_utils.py
index a28fa56f0..99bfda121 100644
--- a/spectral_cube/wcs_utils.py
+++ b/spectral_cube/wcs_utils.py
@@ -477,4 +477,8 @@ def drop_axis_by_slicing(mywcs, shape, dropped_axis,
 
     result.wcs.cdelt[dropped_axis] = dropped_axis_cdelt
 
+    new_inds = np.array([ii for ii in range(ndim) if ii != dropped_axis] +
+                        [dropped_axis])
+    result = reindex_wcs(result, new_inds)
+
     return result

From 4176dd880b924c25f8703978d22932ce3912e2ac Mon Sep 17 00:00:00 2001
From: "Adam Ginsburg (keflavich)" <keflavich@gmail.com>
Date: Sat, 3 Jun 2017 18:19:19 -0600
Subject: [PATCH 07/12] working through some more test cases

---
 spectral_cube/spectral_cube.py |  8 ++++----
 spectral_cube/wcs_utils.py     | 35 +++++++++++++++++++++++++++++++++-
 2 files changed, 38 insertions(+), 5 deletions(-)

diff --git a/spectral_cube/spectral_cube.py b/spectral_cube/spectral_cube.py
index 0ced72dad..cc423edc1 100644
--- a/spectral_cube/spectral_cube.py
+++ b/spectral_cube/spectral_cube.py
@@ -1391,10 +1391,10 @@ def moment(self, order=0, axis=0, how='auto'):
         if order == 1 and axis == 0:
             out += self.world[0, :, :][0]
 
-        view = [slice(None) if ax!=axis else slice(self.shape[ax]//2,
-                                                   self.shape[ax]//2+1)
-                for ax in range(self.ndim)]
-        new_wcs = wcs_utils.slice_wcs(self._wcs, view)
+        new_wcs = wcs_utils.drop_axis_by_slicing(self.wcs, self.shape,
+                                                 dropped_axis=axis,
+                                                 dropped_axis_slice_position='middle',
+                                                 dropped_axis_cdelt='full_range')
 
         meta = {'moment_order': order,
                 'moment_axis': axis,
diff --git a/spectral_cube/wcs_utils.py b/spectral_cube/wcs_utils.py
index 99bfda121..4a5cfb9aa 100644
--- a/spectral_cube/wcs_utils.py
+++ b/spectral_cube/wcs_utils.py
@@ -1,7 +1,7 @@
 from __future__ import print_function, absolute_import, division
 
 import numpy as np
-from astropy.wcs import WCS
+from astropy.wcs import WCS, InconsistentAxisTypesError
 import warnings
 from astropy import units as u
 from astropy import log
@@ -18,6 +18,28 @@
                          'WAVELENG':'WAVE',
                          }
 
+class WCSWrapper(WCS):
+    """
+    Wrapper of WCS to deal with some of the special cases we face within
+    spectral_cube
+    """
+    @property
+    def has_celestial(self):
+        if hasattr(self, '_has_celestial'):
+            return self._has_celestial
+        try:
+            return self.celestial.naxis == 2
+        except InconsistentAxisTypesError:
+            return False
+
+    @has_celestial.setter
+    def has_celestial(self, value):
+        if value is not False:
+            warnings.warn("_has_celestial is being set to {0}, "
+                          "which may not be what you want."
+                          .format(value))
+        self._has_celestial = value
+
 def drop_axis(wcs, dropax):
     """
     Drop the ax on axis dropax
@@ -444,6 +466,11 @@ def drop_axis_by_slicing(mywcs, shape, dropped_axis,
     """
     ndim = len(shape)
 
+    if mywcs.get_axis_types()[dropped_axis]['coordinate_type'] == 'celestial':
+        dropping_celestial = True
+    else:
+        dropping_celestial = False
+
     if dropped_axis_slice_position == 'middle':
         dropped_axis_slice_position = shape[dropped_axis]//2
     elif dropped_axis_slice_position == 'start':
@@ -481,4 +508,10 @@ def drop_axis_by_slicing(mywcs, shape, dropped_axis,
                         [dropped_axis])
     result = reindex_wcs(result, new_inds)
 
+    if dropping_celestial:
+        new_result = WCSWrapper()
+        new_result.wcs = result.wcs
+        new_result.has_celestial = False
+        result = new_result
+
     return result

From d67854d48ab0533672ab4135d559505660894d3e Mon Sep 17 00:00:00 2001
From: "Adam Ginsburg (keflavich)" <keflavich@gmail.com>
Date: Sat, 3 Jun 2017 18:35:45 -0600
Subject: [PATCH 08/12] make the wrapper do something sensible if a smaller
 array is passed than wcs_pix2world allows (so if you've 'dropped' a
 dimension, allow indexing with ndim-1 coordinates)

---
 spectral_cube/spectral_cube.py |  8 ++++----
 spectral_cube/wcs_utils.py     | 12 ++++++++++++
 2 files changed, 16 insertions(+), 4 deletions(-)

diff --git a/spectral_cube/spectral_cube.py b/spectral_cube/spectral_cube.py
index cc423edc1..8a2b87081 100644
--- a/spectral_cube/spectral_cube.py
+++ b/spectral_cube/spectral_cube.py
@@ -2685,10 +2685,10 @@ def __getitem__(self, view):
                                     meta=meta)
 
             # only one element, so drop an axis
-            view = [slice(None) if ax!=intslices[0]
-                    else slice(self.shape[ax]//2, self.shape[ax]//2+1)
-                    for ax in range(self.ndim)]
-            newwcs = wcs_utils.slice_wcs(self._wcs, view)
+            newwcs = wcs_utils.drop_axis_by_slicing(self.wcs, self.shape,
+                                                    dropped_axis=intslices[0],
+                                                    dropped_axis_slice_position='middle',
+                                                    dropped_axis_cdelt='full_range')
             header = self._nowcs_header
 
             # Slice objects know how to parse Beam objects stored in the
diff --git a/spectral_cube/wcs_utils.py b/spectral_cube/wcs_utils.py
index 4a5cfb9aa..284eee188 100644
--- a/spectral_cube/wcs_utils.py
+++ b/spectral_cube/wcs_utils.py
@@ -40,6 +40,17 @@ def has_celestial(self, value):
                           .format(value))
         self._has_celestial = value
 
+    def wcs_pix2world(self, pixels, reference):
+        if ((pixels.shape[1] < self.naxis and
+            hasattr(self, 'active_dimensions') and
+             len(self.active_dimensions) < self.naxis)):
+            pixels = np.asarray(pixels)
+            pixels = np.c_[pixels, np.zeros(pixels.shape[0])]
+            result = super(WCSWrapper, self).wcs_pix2world(pixels, reference)
+            return result[:, :len(self.active_dimensions) - self.naxis]
+        else:
+            return super(WCSWrapper, self).wcs_pix2world(pixels, reference)
+
 def drop_axis(wcs, dropax):
     """
     Drop the ax on axis dropax
@@ -512,6 +523,7 @@ def drop_axis_by_slicing(mywcs, shape, dropped_axis,
         new_result = WCSWrapper()
         new_result.wcs = result.wcs
         new_result.has_celestial = False
+        new_result.active_dimensions = list(range(ndim-1))
         result = new_result
 
     return result

From acffbdb7f12d81385d7d10dc1771235243ded7af Mon Sep 17 00:00:00 2001
From: "Adam Ginsburg (keflavich)" <keflavich@gmail.com>
Date: Sun, 4 Jun 2017 17:22:47 -0600
Subject: [PATCH 09/12] add some more tests that are failing.  For some reason,
 the new slicing operation makes a non-self-consistent WCS object, i.e., one
 where `wcs.naxis != wcs.wcs.naxis`, which is supposed to be impossible, and
 unfortunately is beyond my immediate capability to debug because that part is
 not written in python

---
 spectral_cube/lower_dimensional_structures.py | 24 ++++++-------
 spectral_cube/spectral_cube.py                | 34 ++++++++++++++++---
 spectral_cube/tests/test_projection.py        |  7 ++++
 spectral_cube/tests/test_visualization.py     |  2 +-
 spectral_cube/wcs_utils.py                    |  8 +++++
 5 files changed, 57 insertions(+), 18 deletions(-)

diff --git a/spectral_cube/lower_dimensional_structures.py b/spectral_cube/lower_dimensional_structures.py
index b1f4bf74d..d8667b7d9 100644
--- a/spectral_cube/lower_dimensional_structures.py
+++ b/spectral_cube/lower_dimensional_structures.py
@@ -8,7 +8,7 @@
 from astropy import convolution
 from astropy import units as u
 from astropy import wcs
-#from astropy import log
+from astropy import log
 from astropy.io.fits import Header, Card, HDUList, PrimaryHDU
 
 from .io.core import determine_format
@@ -41,10 +41,18 @@ def header(self):
         for keyword in header:
             if 'NAXIS' in keyword:
                 del header[keyword]
-        header.insert(2, Card(keyword='NAXIS', value=self.ndim))
+
+        header.insert(self.wcs.naxis, Card(keyword='NAXIS', value=self.ndim))
         for ind,sh in enumerate(self.shape[::-1]):
             header.insert(3+ind, Card(keyword='NAXIS{0:1d}'.format(ind+1),
                                       value=sh))
+        if self.wcs.naxis > self.ndim:
+            assert ind != 0
+            for ii in range(self.wcs.naxis - self.ndim):
+                log.debug('Adding NAXIS{0}'.format(ind+ii+2))
+                header.insert(3+ind+ii+1,
+                              Card(keyword='NAXIS{0:1d}'.format(ind+ii+2),
+                                   value=1))
 
         if 'beam' in self.meta:
             header.update(self.meta['beam'].to_header_keywords())
@@ -518,16 +526,8 @@ def __repr__(self):
 
     @property
     def header(self):
-        header = self._header
-        # This inplace update is OK; it's not bad to overwrite WCS in this
-        # header
-        if self.wcs is not None:
-            header.update(self.wcs.to_header())
-        header['BUNIT'] = self.unit.to_string(format='fits')
-        header.insert(2, Card(keyword='NAXIS', value=self.ndim))
-        for ind,sh in enumerate(self.shape[::-1]):
-            header.insert(3+ind, Card(keyword='NAXIS{0:1d}'.format(ind+1),
-                                      value=sh))
+        # use the generic LowerDimensionalObject header, but more
+        header = super(OneDSpectrum, self).header
 
         # Preserve the spectrum's spectral units
         if 'CUNIT1' in header and self._spectral_unit != u.Unit(header['CUNIT1']):
diff --git a/spectral_cube/spectral_cube.py b/spectral_cube/spectral_cube.py
index 8a2b87081..030dfac29 100644
--- a/spectral_cube/spectral_cube.py
+++ b/spectral_cube/spectral_cube.py
@@ -490,7 +490,10 @@ def mean(self, axis=None, how='cube'):
                                             projection=False)
             out = ttl / counts
             if projection:
-                new_wcs = wcs_utils.drop_axis(self._wcs, np2wcs[axis])
+                new_wcs = wcs_utils.drop_axis_by_slicing(self.wcs, self.shape,
+                                                         dropped_axis=axis,
+                                                         dropped_axis_slice_position='middle',
+                                                         dropped_axis_cdelt='full_range')
                 meta = {'collapse_axis': axis}
                 meta.update(self._meta)
                 return Projection(out, copy=False, wcs=new_wcs,
@@ -558,7 +561,10 @@ def std(self, axis=None, how='cube', ddof=0):
             out = (result/(counts-ddof))**0.5
 
             if projection:
-                new_wcs = wcs_utils.drop_axis(self._wcs, np2wcs[axis])
+                new_wcs = wcs_utils.drop_axis_by_slicing(self.wcs, self.shape,
+                                                         dropped_axis=axis,
+                                                         dropped_axis_slice_position='middle',
+                                                         dropped_axis_cdelt='full_range')
                 meta = {'collapse_axis': axis}
                 meta.update(self._meta)
                 return Projection(out, copy=False, wcs=new_wcs,
@@ -730,7 +736,10 @@ def _cube_on_cube_operation(self, function, cube, equivalencies=[], **kwargs):
         return self._new_cube_with(data=data, unit=unit)
 
     def apply_function(self, function, axis=None, weights=None, unit=None,
-                       projection=False, progressbar=False, **kwargs):
+                       projection=False, progressbar=False,
+                       dropped_axis_slice_position='middle',
+                       dropped_axis_cdelt='same',
+                       **kwargs):
         """
         Apply a function to valid data along the specified axis or to the whole
         cube, optionally using a weight array that is the same shape (or at
@@ -754,6 +763,15 @@ def apply_function(self, function, axis=None, weights=None, unit=None,
         progressbar : bool
             Show a progressbar while iterating over the slices/rays through the
             cube?
+        dropped_axis_slice_position : 'middle', 'start', 'end'
+            If an axis is being dropped, where should the WCS say the
+            projection is?  It can be at the start, middle, or end of the
+            axis.
+        dropped_axis_cdelt : 'same', 'full_range', or value
+            If an axis is being dropped, what should the new CDELT be?  For an
+            integral, for example, one might want the value to be the full
+            range.  For a slice, it should stay the same.  For something like
+            min or max, it might be zero.
 
         Returns
         -------
@@ -796,7 +814,10 @@ def apply_function(self, function, axis=None, weights=None, unit=None,
             pbu()
 
         if projection and axis in (0,1,2):
-            new_wcs = wcs_utils.drop_axis(self._wcs, np2wcs[axis])
+            new_wcs = wcs_utils.drop_axis_by_slicing(self.wcs, self.shape,
+                                                     dropped_axis=axis,
+                                                     dropped_axis_slice_position=dropped_axis_slice_position,
+                                                     dropped_axis_cdelt=dropped_axis_cdelt)
 
             meta = {'collapse_axis': axis}
             meta.update(self._meta)
@@ -1041,7 +1062,10 @@ def __getitem__(self, view):
                                    )
 
             # only one element, so drop an axis
-            newwcs = wcs_utils.drop_axis(self._wcs, intslices[0])
+            newwcs = wcs_utils.drop_axis_by_slicing(self.wcs, self.shape,
+                                                    dropped_axis=intslices[0],
+                                                    dropped_axis_slice_position=view[2-intslices[0]],
+                                                    dropped_axis_cdelt='same')
             header = self._nowcs_header
 
             if intslices[0] == 0:
diff --git a/spectral_cube/tests/test_projection.py b/spectral_cube/tests/test_projection.py
index 648359f43..bae8cb73f 100644
--- a/spectral_cube/tests/test_projection.py
+++ b/spectral_cube/tests/test_projection.py
@@ -511,3 +511,10 @@ def test_1d_slice_round():
 
     assert 'OneDSpectrum' in sp.round().__repr__()
     assert 'OneDSpectrum' in sp[1:-1].round().__repr__()
+
+def test_slice_wcs():
+
+    cube, data = cube_and_raw('255_delta.fits')
+
+    mom0 = cube.moment0(axis=0)
+    assert mom0.header['NAXIS3'] == 1
diff --git a/spectral_cube/tests/test_visualization.py b/spectral_cube/tests/test_visualization.py
index b13790af4..833bd69e9 100644
--- a/spectral_cube/tests/test_visualization.py
+++ b/spectral_cube/tests/test_visualization.py
@@ -38,7 +38,7 @@ def test_to_pvextractor():
 
 
 @pytest.mark.skipif("not MATPLOTLIB_INSTALLED")
-def test_projvis():
+def test_projvis_noaplpy():
 
     cube, data = cube_and_raw('vda_Jybeam_lower.fits')
 
diff --git a/spectral_cube/wcs_utils.py b/spectral_cube/wcs_utils.py
index 284eee188..4d01eb5cc 100644
--- a/spectral_cube/wcs_utils.py
+++ b/spectral_cube/wcs_utils.py
@@ -183,6 +183,8 @@ def reindex_wcs(wcs, inds):
                 ps_cards.append((i, m, v))
     outwcs.wcs.set_ps(ps_cards)
 
+    assert outwcs.naxis == len(inds)
+
     return outwcs
 
 
@@ -475,6 +477,9 @@ def drop_axis_by_slicing(mywcs, shape, dropped_axis,
         If the axes are misaligned, it is not possible to "drop" an axis.
         In this case, a generic "offset axis" will be returned.
     """
+    log.debug("Dropping axis by slicing with args: {0}, {1}, {2}, {3}, {4}"
+              .format(mywcs, shape, dropped_axis, dropped_axis_slice_position,
+                      dropped_axis_cdelt))
     ndim = len(shape)
 
     if mywcs.get_axis_types()[dropped_axis]['coordinate_type'] == 'celestial':
@@ -526,4 +531,7 @@ def drop_axis_by_slicing(mywcs, shape, dropped_axis,
         new_result.active_dimensions = list(range(ndim-1))
         result = new_result
 
+    assert result.naxis == result.wcs.naxis
+    assert result.naxis == mywcs.naxis
+
     return result

From 2322ecb603981dd622e1c3d308a2a26d4bf0b045 Mon Sep 17 00:00:00 2001
From: "Adam Ginsburg (keflavich)" <keflavich@gmail.com>
Date: Sun, 4 Jun 2017 17:37:14 -0600
Subject: [PATCH 10/12] fixed naxis disagreement issue

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

diff --git a/spectral_cube/wcs_utils.py b/spectral_cube/wcs_utils.py
index 4d01eb5cc..fd3f6fa2b 100644
--- a/spectral_cube/wcs_utils.py
+++ b/spectral_cube/wcs_utils.py
@@ -23,6 +23,20 @@ class WCSWrapper(WCS):
     Wrapper of WCS to deal with some of the special cases we face within
     spectral_cube
     """
+    @staticmethod
+    def from_wcs(otherwcs):
+        """
+        Create a WCSWrapper class from another WCS object
+        """
+        new_wcs = WCSWrapper()
+
+        new_wcs.wcs = otherwcs.wcs
+        new_wcs.naxis = otherwcs.naxis
+
+        assert new_wcs.wcs.naxis == otherwcs.wcs.naxis == otherwcs.naxis == new_wcs.naxis
+
+        return new_wcs
+
     @property
     def has_celestial(self):
         if hasattr(self, '_has_celestial'):
@@ -184,6 +198,7 @@ def reindex_wcs(wcs, inds):
     outwcs.wcs.set_ps(ps_cards)
 
     assert outwcs.naxis == len(inds)
+    assert outwcs.wcs.naxis == len(inds)
 
     return outwcs
 
@@ -525,10 +540,10 @@ def drop_axis_by_slicing(mywcs, shape, dropped_axis,
     result = reindex_wcs(result, new_inds)
 
     if dropping_celestial:
-        new_result = WCSWrapper()
-        new_result.wcs = result.wcs
+        new_result = WCSWrapper.from_wcs(result)
         new_result.has_celestial = False
         new_result.active_dimensions = list(range(ndim-1))
+        assert new_result.naxis == new_result.wcs.naxis
         result = new_result
 
     assert result.naxis == result.wcs.naxis

From 9b1b5e87df4dc25d24831aafbb46aa422c266d46 Mon Sep 17 00:00:00 2001
From: "Adam Ginsburg (keflavich)" <keflavich@gmail.com>
Date: Sun, 4 Jun 2017 17:44:01 -0600
Subject: [PATCH 11/12] add some debug statements (to be removed eventually)
 and fix the NAXIS inserting

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

diff --git a/spectral_cube/lower_dimensional_structures.py b/spectral_cube/lower_dimensional_structures.py
index d8667b7d9..3b9f1dcd2 100644
--- a/spectral_cube/lower_dimensional_structures.py
+++ b/spectral_cube/lower_dimensional_structures.py
@@ -42,15 +42,18 @@ def header(self):
             if 'NAXIS' in keyword:
                 del header[keyword]
 
-        header.insert(self.wcs.naxis, Card(keyword='NAXIS', value=self.ndim))
+        header.insert(3, Card(keyword='NAXIS', value=self.ndim))
         for ind,sh in enumerate(self.shape[::-1]):
-            header.insert(3+ind, Card(keyword='NAXIS{0:1d}'.format(ind+1),
+            log.debug('Adding NAXIS{0} at position {1}'.format(ind+1,
+                                                               3+ind+1))
+            header.insert(3+ind+1, Card(keyword='NAXIS{0:1d}'.format(ind+1),
                                       value=sh))
         if self.wcs.naxis > self.ndim:
             assert ind != 0
             for ii in range(self.wcs.naxis - self.ndim):
-                log.debug('Adding NAXIS{0}'.format(ind+ii+2))
-                header.insert(3+ind+ii+1,
+                log.debug('Adding NAXIS{0} at position {1}'.format(ind+ii+2,
+                                                                   3+ind+ii+2))
+                header.insert(3+ind+ii+2,
                               Card(keyword='NAXIS{0:1d}'.format(ind+ii+2),
                                    value=1))
 

From 8254a7e8da5bfb5d790b56043d8161c134d58a53 Mon Sep 17 00:00:00 2001
From: "Adam Ginsburg (keflavich)" <keflavich@gmail.com>
Date: Sun, 4 Jun 2017 17:56:21 -0600
Subject: [PATCH 12/12] minor whitespace

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

diff --git a/spectral_cube/lower_dimensional_structures.py b/spectral_cube/lower_dimensional_structures.py
index 3b9f1dcd2..92d2906bd 100644
--- a/spectral_cube/lower_dimensional_structures.py
+++ b/spectral_cube/lower_dimensional_structures.py
@@ -47,7 +47,7 @@ def header(self):
             log.debug('Adding NAXIS{0} at position {1}'.format(ind+1,
                                                                3+ind+1))
             header.insert(3+ind+1, Card(keyword='NAXIS{0:1d}'.format(ind+1),
-                                      value=sh))
+                                        value=sh))
         if self.wcs.naxis > self.ndim:
             assert ind != 0
             for ii in range(self.wcs.naxis - self.ndim):