Fix stacking when the spectral axis is decreasing

spectral_cube/analysis_utilities.py

@@ -205,7 +205,11 @@ def stack_spectra(cube, velocity_surface, v0=None,
                              "axis of the cube.")
 
     # Calculate the pixel shifts that will be applied.
-    vdiff = np.abs(np.diff(cube.spectral_axis[:2])[0])
+    spec_size = np.diff(cube.spectral_axis[:2])[0]
+    # Assign the correct +/- for pixel shifts based on whether the spectral
+    # axis is increasing (-1) or decreasing (+1)
+    vdiff_sign = -1. if spec_size.value > 0. else 1.
+    vdiff = np.abs(spec_size)
     vel_unit = vdiff.unit
 
     # Check to make sure vdiff doesn't change more than the allowed tolerance
@@ -214,8 +218,8 @@ def stack_spectra(cube, velocity_surface, v0=None,
     if not np.isclose(vdiff2.value, vdiff.value, rtol=vdiff_tol):
         raise ValueError("Cannot shift spectra on a non-linear axes")
 
-    pix_shifts = -((velocity_surface.to(vel_unit) -
-                    v0.to(vel_unit)) / vdiff).value[xy_posns]
+    pix_shifts = vdiff_sign * ((velocity_surface.to(vel_unit) -
+                                v0.to(vel_unit)) / vdiff).value[xy_posns]
 
     # May a header copy so we can start altering
     new_header = cube[:, 0, 0].header.copy()

spectral_cube/tests/test_analysis_functions.py

@@ -85,6 +85,42 @@ def test_stacking():
                                test_cube.spectral_axis.value)
 
 
+def test_stacking_reversed_specaxis():
+    '''
+    Use a set of identical Gaussian profiles randomly offset to ensure the
+    shifted spectrum has the correct properties.
+    '''
+
+    amp = 1.
+    v0 = 0. * u.km / u.s
+    sigma = 8.
+    noise = None
+    shape = (100, 25, 25)
+
+    test_cube, test_vels = \
+        generate_gaussian_cube(amp=amp, sigma=sigma, noise=noise,
+                               shape=shape, spec_scale=-1. * u.km / u.s)
+
+    true_spectrum = gaussian(test_cube.spectral_axis.value,
+                             amp, v0.value, sigma)
+
+    # Stack the spectra in the cube
+    stacked = \
+        stack_spectra(test_cube, test_vels, v0=v0,
+                      stack_function=np.nanmean,
+                      xy_posns=None, num_cores=1,
+                      chunk_size=-1,
+                      progressbar=False, pad_edges=False)
+
+    # Calculate residuals
+    resid = np.abs(stacked.value - true_spectrum)
+    assert np.std(resid) <= 1e-3
+
+    # The stacked spectrum should have the same spectral axis
+    np.testing.assert_allclose(stacked.spectral_axis.value,
+                               test_cube.spectral_axis.value)
+
+
 def test_stacking_wpadding():
     '''
     Use a set of identical Gaussian profiles randomly offset to ensure the
@@ -128,6 +164,47 @@ def test_stacking_wpadding():
         or (stacked.size == stack_shape + 1)
 
 
+def test_stacking_woffset():
+    '''
+    Use a set of identical Gaussian profiles randomly offset to ensure the
+    shifted spectrum has the correct properties.
+
+    Make sure the operations aren't affected by absolute velocity offsets
+
+    '''
+
+    amp = 1.
+    sigma = 8.
+    v0 = 100. * u.km / u.s
+    noise = None
+    shape = (100, 25, 25)
+
+    test_cube, test_vels = \
+        generate_gaussian_cube(shape=shape, amp=amp, sigma=sigma, noise=noise,
+                               v0=v0.value)
+
+    # Stack the spectra in the cube
+    stacked = \
+        stack_spectra(test_cube, test_vels, v0=v0,
+                      stack_function=np.nanmean,
+                      xy_posns=None, num_cores=1,
+                      chunk_size=-1,
+                      progressbar=False, pad_edges=True)
+
+    true_spectrum = gaussian(stacked.spectral_axis.value,
+                             amp, v0.value, sigma)
+
+    # Calculate residuals
+    resid = np.abs(stacked.value - true_spectrum)
+    assert np.std(resid) <= 1e-3
+
+    # The spectral axis should be padded by ~25% on each side
+    stack_shape = int(test_cube.shape[0] * 1.5)
+    # This is rounded, so the shape could be +/- 1
+    assert (stacked.size == stack_shape) or (stacked.size == stack_shape - 1) \
+        or (stacked.size == stack_shape + 1)
+
+
 def test_stacking_noisy():
 
     # Test stack w/ S/N of 0.2

spectral_cube/tests/utilities.py

@@ -55,6 +55,7 @@ def generate_gaussian_cube(shape=(100, 25, 25), sigma=8., amp=1.,
                            noise=None, spec_scale=1 * u.km / u.s,
                            pixel_scale=1 * u.arcsec,
                            beamfwhm=3 * u.arcsec,
+                           v0=None,
                            vel_surface=None,
                            seed=247825498):
     '''
@@ -79,6 +80,9 @@ def generate_gaussian_cube(shape=(100, 25, 25), sigma=8., amp=1.,
     spec_middle = int(shape[0] / 2)
     spec_quarter = int(shape[0] / 4)
 
+    if v0 is None:
+        v0 = 0
+
     with NumpyRNGContext(seed):
 
         spec_inds = np.mgrid[-spec_middle:spec_middle] * spec_scale.value
@@ -88,13 +92,13 @@ def generate_gaussian_cube(shape=(100, 25, 25), sigma=8., amp=1.,
             mean_pos = \
                 np.random.uniform(low=spec_inds[spec_quarter],
                                   high=spec_inds[spec_quarter + spec_middle])
-            mean_positions[y, x] = mean_pos
             test_cube[:, y, x] = gaussian(spec_inds, amp, mean_pos, sigma)
+            mean_positions[y, x] = mean_pos + v0
             if noise is not None:
                 test_cube[:, y, x] += np.random.normal(0, noise, shape[0])
 
     test_hdu = generate_hdu(test_cube, pixel_scale, spec_scale, beamfwhm,
-                            spec_inds[0])
+                            spec_inds[0] + v0)
 
     spec_cube = SpectralCube.read(test_hdu)