diff --git a/spectral_cube/analysis_utilities.py b/spectral_cube/analysis_utilities.py
index d92ec2cd3..8496141f6 100644
--- a/spectral_cube/analysis_utilities.py
+++ b/spectral_cube/analysis_utilities.py
@@ -232,15 +232,22 @@ def stack_spectra(cube, velocity_surface, v0=None,
         # Find max +/- pixel shifts, rounding up to the nearest integer
         max_pos_shift = np.ceil(np.nanmax(pix_shifts)).astype(int)
         max_neg_shift = np.ceil(np.nanmin(pix_shifts)).astype(int)
+        if max_neg_shift > 0:
+            # if there are no negative shifts, we can ignore them and just
+            # use the positive shift
+            max_neg_shift = 0
+        if max_pos_shift < 0:
+            # same for positive
+            max_pos_shift = 0
 
         # The total pixel size of the new spectral axis
         num_vel_pix = cube.spectral_axis.size + max_pos_shift - max_neg_shift
         new_header['NAXIS1'] = num_vel_pix
 
         # Adjust CRPIX in header
-        new_header['CRPIX1'] += max_pos_shift
+        new_header['CRPIX1'] += -max_neg_shift
 
-        pad_size = (max_pos_shift, -max_neg_shift)
+        pad_size = (-max_neg_shift, max_pos_shift)
 
     else:
         pad_size = None
@@ -269,7 +276,10 @@ def stack_spectra(cube, velocity_surface, v0=None,
 
             all_shifted_spectra.extend([out for out in shifted_spectra])
 
-    stacked = stack_function(all_shifted_spectra, axis=0)
+    shifted_spectra_array = np.array(all_shifted_spectra)
+    assert shifted_spectra_array.ndim == 2
+
+    stacked = stack_function(shifted_spectra_array, axis=0)
 
     stack_spec = \
         OneDSpectrum(stacked, unit=cube.unit, wcs=WCS(new_header),
diff --git a/spectral_cube/tests/test_analysis_functions.py b/spectral_cube/tests/test_analysis_functions.py
index 4d46b2e77..d113aa213 100644
--- a/spectral_cube/tests/test_analysis_functions.py
+++ b/spectral_cube/tests/test_analysis_functions.py
@@ -164,6 +164,41 @@ def test_stacking_wpadding():
         or (stacked.size == stack_shape + 1)
 
 
+def test_padding_direction():
+
+    amp = 1.
+    sigma = 8.
+    v0 = 0. * u.km / u.s
+    noise = None
+    shape = (100, 2, 2)
+
+    vel_surface = np.array([[0, 5], [5, 10]])
+
+    test_cube, test_vels = \
+        generate_gaussian_cube(shape=shape, amp=amp, sigma=sigma, noise=noise,
+                               vel_surface=vel_surface)
+
+    # 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)
+
+    # now check that the stacked spectral axis is right
+    # (all shifts are negative, so vmin < -50 km/s, should be -60?)
+    assert stacked.spectral_axis.min() == -60*u.km/u.s
+    assert stacked.spectral_axis.max() == 49*u.km/u.s
+
+    # Calculate residuals
+    resid = np.abs(stacked.value - true_spectrum)
+    assert np.std(resid) <= 1e-3
+
+
 def test_stacking_woffset():
     '''
     Use a set of identical Gaussian profiles randomly offset to ensure the
diff --git a/spectral_cube/tests/utilities.py b/spectral_cube/tests/utilities.py
index 370bebbf8..7c6aec328 100644
--- a/spectral_cube/tests/utilities.py
+++ b/spectral_cube/tests/utilities.py
@@ -89,9 +89,12 @@ def generate_gaussian_cube(shape=(100, 25, 25), sigma=8., amp=1.,
         spat_inds = np.indices(shape[1:])
         for y, x in zip(spat_inds[0].flatten(), spat_inds[1].flatten()):
             # Lock the mean to within 25% from the centre
-            mean_pos = \
-                np.random.uniform(low=spec_inds[spec_quarter],
-                                  high=spec_inds[spec_quarter + spec_middle])
+            if vel_surface is not None:
+                mean_pos = vel_surface[y,x]
+            else:
+                mean_pos = \
+                    np.random.uniform(low=spec_inds[spec_quarter],
+                                      high=spec_inds[spec_quarter + spec_middle])
             test_cube[:, y, x] = gaussian(spec_inds, amp, mean_pos, sigma)
             mean_positions[y, x] = mean_pos + v0
             if noise is not None: