diff --git a/spectral_cube/analysis_utilities.py b/spectral_cube/analysis_utilities.py
index 28aebcb67..25f732bb6 100644
--- a/spectral_cube/analysis_utilities.py
+++ b/spectral_cube/analysis_utilities.py
@@ -245,7 +245,7 @@ def stack_spectra(cube, velocity_surface, v0=None,
     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
+    # Make a header copy so we can start altering
     new_header = cube[:, 0, 0].header.copy()
 
     if pad_edges:
diff --git a/spectral_cube/spectral_cube.py b/spectral_cube/spectral_cube.py
index 2d45cdb57..d9b88f010 100644
--- a/spectral_cube/spectral_cube.py
+++ b/spectral_cube/spectral_cube.py
@@ -2050,7 +2050,8 @@ def subcube_from_crtfregion(self, crtf_region, allow_empty=False):
 
         return self.subcube_from_regions(region_list, allow_empty)
 
-    def subcube_from_regions(self, region_list, allow_empty=False):
+    def subcube_from_regions(self, region_list, allow_empty=False,
+                             minimize=False):
         """
         Extract a masked subcube from a list of ``regions.Region`` object
         (only functions on celestial dimensions)
@@ -2062,6 +2063,10 @@ def subcube_from_regions(self, region_list, allow_empty=False):
         allow_empty: bool, optional
             If this is False, an exception will be raised if the region
             contains no overlap with the cube. Default is False.
+        minimize : bool
+            Run `minimal_subcube`?  This should be redundant, since the
+            bounding box of the region is already used, but it might slice off
+            an extra few pixels depending on the details of the region shape.
         """
         import regions
 
@@ -2123,10 +2128,14 @@ def subcube_from_regions(self, region_list, allow_empty=False):
         maskarray = np.zeros(subcube.shape[1:], dtype='bool')
         maskarray[:] = mask.data[slices_small]
 
-        masked_subcube = subcube.with_mask(BooleanArrayMask(maskarray, subcube.wcs, shape=subcube.shape))
+        BAM = BooleanArrayMask(maskarray, subcube.wcs.celestial, shape=subcube.shape[1:])
+        masked_subcube = subcube.with_mask(BAM)
         # by using ceil / floor above, we potentially introduced a NaN buffer
         # that we can now crop out
-        return masked_subcube.minimal_subcube(spatial_only=True)
+        if minimize:
+            return masked_subcube.minimal_subcube(spatial_only=True)
+        else:
+            return masked_subcube
 
     def _velocity_freq_conversion_regions(self, ranges, veltypes, restfreqs):
         """