Merge pull request #478 from keflavich/wcs_reversal

MAJOR: WCS reversal is non-functional

CHANGES.rst

@@ -5,6 +5,9 @@
    ``apply_function_parallel_spectral`` to make this general.  Added
    ``joblib`` as a dependency.
    (https://github.com/radio-astro-tools/spectral-cube/pull/474)
+ - Bugfix: Reversing a cube's spectral axis should now do something reasonable
+   instead of unreasonable
+   (https://github.com/radio-astro-tools/spectral-cube/pull/478)
 
 0.4.2 (2018-02-21)
 ------------------

spectral_cube/base_class.py

@@ -142,8 +142,9 @@ def world_extrema(self):
         _lat_min = lat.min()
         _lat_max = lat.max()
 
-        return ((_lon_min, _lon_max),
-                (_lat_min, _lat_max))
+        return u.Quantity(((_lon_min.to(u.deg).value, _lon_max.to(u.deg).value),
+                           (_lat_min.to(u.deg).value, _lat_max.to(u.deg).value)),
+                          u.deg)
 
     @property
     @cached

spectral_cube/tests/test_spectral_cube.py

@@ -932,6 +932,7 @@ def test_slicing():
                           ((slice(None), slice(None), 1), 2),
                           ((slice(None), slice(None), slice(1)), 3),
                           ((slice(1), slice(1), slice(1)), 3),
+                          ((slice(None, None, -1), slice(None), slice(None)), 3),
                          ])
 def test_slice_wcs(view, naxis):
 
@@ -940,6 +941,28 @@ def test_slice_wcs(view, naxis):
     sl = cube[view]
     assert sl.wcs.naxis == naxis
 
+def test_slice_wcs_reversal():
+    cube, data = cube_and_raw('advs.fits')
+    view = (slice(None,None,-1), slice(None), slice(None))
+
+    rcube = cube[view]
+    rrcube = rcube[view]
+
+    np.testing.assert_array_equal(np.diff(cube.spectral_axis),
+                                  -np.diff(rcube.spectral_axis))
+
+    np.testing.assert_array_equal(rrcube.spectral_axis.value,
+                                  cube.spectral_axis.value)
+    np.testing.assert_array_equal(rcube.spectral_axis.value,
+                                  cube.spectral_axis.value[::-1])
+    np.testing.assert_array_equal(rrcube.world_extrema.value,
+                                  cube.world_extrema.value)
+    # check that the lon, lat arrays are *entirely* unchanged
+    np.testing.assert_array_equal(rrcube.spatial_coordinate_map[0].value,
+                                  cube.spatial_coordinate_map[0].value)
+    np.testing.assert_array_equal(rrcube.spatial_coordinate_map[1].value,
+                                  cube.spatial_coordinate_map[1].value)
+
 def test_spectral_slice_preserve_units():
     cube, data = cube_and_raw('advs.fits')
     cube = cube.with_spectral_unit(u.km/u.s)
@@ -1607,7 +1630,8 @@ def test_caching():
     world_extrema_function = base_class.SpatialCoordMixinClass.world_extrema.fget.wrapped_function
 
     assert cube.world_extrema is cube._cache[(world_extrema_function, ())]
-    assert worldextrema == cube.world_extrema
+    np.testing.assert_almost_equal(worldextrema.value,
+                                   cube.world_extrema.value)
 
 def test_spatial_smooth_g2d():
     cube, data = cube_and_raw('adv.fits')

spectral_cube/tests/test_wcs_utils.py

@@ -93,6 +93,37 @@ def test_wcs_slice_reversal():
 
     np.testing.assert_allclose(spaxis, new_spaxis[::-1])
 
+def test_reversal_roundtrip():
+    wcs = WCS(naxis=3)
+    wcs.wcs.crpix = [50., 45., 30.]
+    wcs.wcs.crval = [0., 0., 0.]
+    wcs.wcs.cdelt = [1., 1., 1.]
+    wcs_new = slice_wcs(wcs, (slice(None, None, -1), slice(None), slice(None)),
+                        shape=[100., 150., 200.])
+    spaxis = wcs.sub([0]).wcs_pix2world(np.arange(100), 0)
+
+    new_spaxis = wcs_new.sub([0]).wcs_pix2world(np.arange(100), 0)
+
+    np.testing.assert_allclose(spaxis, new_spaxis[::-1])
+
+    re_reverse = slice_wcs(wcs_new, (slice(None, None, -1), slice(None), slice(None)),
+                           shape=[100., 150., 200.])
+    new_spaxis = re_reverse.sub([0]).wcs_pix2world(np.arange(100), 0)
+
+    np.testing.assert_allclose(spaxis, new_spaxis[::-1])
+
+    #These are NOT equal, but they are equivalent: CRVAL and CRPIX are shifted
+    #by an acceptable amount
+    # assert check_equality(wcs, re_reverse)
+
+    re_re_reverse = slice_wcs(re_reverse, (slice(None, None, -1), slice(None),
+                                           slice(None)),
+                              shape=[100., 150., 200.])
+    re_re_re_reverse = slice_wcs(re_re_reverse, (slice(None, None, -1),
+                                                 slice(None), slice(None)),
+                                 shape=[100., 150., 200.])
+
+    assert check_equality(re_re_re_reverse, re_reverse)
 
 def test_wcs_comparison():
     wcs1 = WCS(naxis=3)

spectral_cube/wcs_utils.py

@@ -246,10 +246,13 @@ def slice_wcs(mywcs, view, shape=None, numpy_order=True,
             # this will raise an inconsistent axis type error if slicing over
             # celestial axes is attempted
             # wcs_index+1 is required because sub([0]) = sub([all])
-            crval = mywcs.sub([wcs_index+1]).wcs_pix2world([refpix], 0)[0]
+            crval = mywcs.sub([wcs_index+1]).wcs_pix2world([refpix-1], 0)[0]
             crpix = 1
+            cdelt = mywcs.wcs.cdelt[wcs_index]
+
             wcs_new.wcs.crpix[wcs_index] = crpix
             wcs_new.wcs.crval[wcs_index] = crval
+            wcs_new.wcs.cdelt[wcs_index] = -cdelt
 
         elif iview.start is not None: