task: metadata: install: pip install -e .[test]
patch: diff --git a/CHANGES.rst b/CHANGES.rst
index 3ac05dd7..009685e7 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -9,6 +9,13 @@ Bug Fixes
 
 - Fixed cube transpose bug in moment analysis for non-square spatial arrays. [#912]
 
+- Fixed ``spectral_slab`` crashing when ``spectral_axis`` has unit of pixels and
+  the bounds are also defined in the unit of pixels. [#926]
+
+- Fixed resulting ``spectral_axis`` containing NaN when a cube is passed into
+  ``Spectrum1D`` without WCS nor spectral axis and the spatial-spatial dimension
+  is smaller than spectral dimension. [#926]
+
 Other Changes and Additions
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
diff --git a/specutils/manipulation/extract_spectral_region.py b/specutils/manipulation/extract_spectral_region.py
index 164de832..294737e9 100644
--- a/specutils/manipulation/extract_spectral_region.py
+++ b/specutils/manipulation/extract_spectral_region.py
@@ -70,8 +70,11 @@ def _subregion_to_edge_pixels(subregion, spectrum):
         right_func = min
 
     # Left/lower side of sub region
-    if subregion[0].unit.is_equivalent(u.pix) and not spectral_axis.unit.is_equivalent(u.pix):
-        left_index = floor(subregion[0].value)
+    if subregion[0].unit.is_equivalent(u.pix):
+        if not spectral_axis.unit.is_equivalent(u.pix):
+            left_index = floor(subregion[0].value)
+        else:
+            left_index = int(np.ceil(left_func(subregion).value))
     else:
         # Convert lower value to spectrum spectral_axis units
         left_reg_in_spec_unit = left_func(subregion).to(spectral_axis.unit,
@@ -79,8 +82,11 @@ def _subregion_to_edge_pixels(subregion, spectrum):
         left_index = _edge_value_to_pixel(left_reg_in_spec_unit, spectrum, order, "left")
 
     # Right/upper side of sub region
-    if subregion[1].unit.is_equivalent(u.pix) and not spectral_axis.unit.is_equivalent(u.pix):
-        right_index = ceil(subregion[1].value)
+    if subregion[1].unit.is_equivalent(u.pix):
+        if not spectral_axis.unit.is_equivalent(u.pix):
+            right_index = ceil(subregion[1].value)
+        else:
+            right_index = int(np.floor(right_func(subregion).value)) + 1
     else:
         # Convert upper value to spectrum spectral_axis units
         right_reg_in_spec_unit = right_func(subregion).to(spectral_axis.unit,
diff --git a/specutils/spectra/spectrum1d.py b/specutils/spectra/spectrum1d.py
index d0b22a0b..907b7aca 100644
--- a/specutils/spectra/spectrum1d.py
+++ b/specutils/spectra/spectrum1d.py
@@ -265,7 +265,7 @@ class Spectrum1D(OneDSpectrumMixin, NDCube, NDIOMixin, NDArithmeticMixin):
         elif wcs is None:
             # If no spectral axis or wcs information is provided, initialize
             # with an empty gwcs based on the flux.
-            size = len(flux) if not flux.isscalar else 1
+            size = flux.shape[-1] if not flux.isscalar else 1
             wcs = gwcs_from_array(np.arange(size) * u.Unit(""))
 
         super().__init__(

problem_statement: spectral_slab crashes when spectral axis is in the unit of pixel
I will have to try to get a minimally reproducible example later. Just posting the traceback now so I don't lose it. This is using spacetelescope/jdaviz#1040 (Collapse plugin) and glue-viz/glue-astronomy#61 on a data cube without WCS (uncalibrated data).

```
.../specutils/manipulation/extract_spectral_region.py in _edge_value_to_pixel(edge_value, spectrum, order, side)
     30         else:
---> 31             index = spectrum.wcs.world_to_pixel(edge_value)
     32 

.../gwcs/api.py in world_to_pixel(self, *world_objects)
    313         """
--> 314         result = self.invert(*world_objects, with_units=True)
    315         if not utils.isnumerical(result[0]):

.../gwcs/wcs.py in invert(self, *args, **kwargs)
    507             akwargs = {k: v for k, v in kwargs.items() if k not in _ITER_INV_KWARGS}
--> 508             result = self.backward_transform(*args, **akwargs)
    509         except (NotImplementedError, KeyError):

.../astropy/modeling/core.py in __call__(self, model_set_axis, with_bounding_box, fill_value, equivalencies, inputs_map, *inputs, **new_inputs)
    397 
--> 398             new_call = make_function_with_signature(
    399                 __call__, args, kwargs, varargs='inputs', varkwargs='new_inputs')

.../astropy/modeling/core.py in __call__(self, *inputs, **kwargs)
    376                 """Evaluate this model on the supplied inputs."""
--> 377                 return super(cls, self).__call__(*inputs, **kwargs)
    378 

.../astropy/modeling/core.py in __call__(self, *args, **kwargs)
   1068         # prepare for model evaluation (overridden in CompoundModel)
-> 1069         evaluate, inputs, broadcasted_shapes, kwargs = self._pre_evaluate(*args, **kwargs)
   1070 

.../astropy/modeling/core.py in _pre_evaluate(self, *args, **kwargs)
    924         # Broadcast inputs into common size
--> 925         inputs, broadcasted_shapes = self.prepare_inputs(*args, **kwargs)
    926 

.../astropy/modeling/core.py in prepare_inputs(self, model_set_axis, equivalencies, *inputs, **kwargs)
   2013 
-> 2014         inputs = self._validate_input_units(inputs, equivalencies, inputs_map)
   2015 

.../astropy/modeling/core.py in _validate_input_units(self, inputs, equivalencies, inputs_map)
   2081                         if input_unit is dimensionless_unscaled:
-> 2082                             raise UnitsError("{0}: Units of input '{1}', {2} ({3}),"
   2083                                              "could not be converted to "

UnitsError: SpectralTabular1D: Units of input 'x', pix (unknown),could not be converted to required dimensionless input

During handling of the above exception, another exception occurred:

...

.../jdaviz/configs/default/plugins/collapse/collapse.py in vue_collapse(self, *args, **kwargs)
    108         with warnings.catch_warnings():
    109             warnings.filterwarnings('ignore', message='No observer defined on WCS')
--> 110             spec = spectral_slab(self._selected_cube, spec_min, spec_max)
    111             # Spatial-spatial image only.
    112             collapsed_spec = spec.collapse(self.selected_func.lower(), axis=-1).T  # Quantity

.../specutils/manipulation/extract_spectral_region.py in spectral_slab(spectrum, lower, upper)
    224     region = SpectralRegion(lower, upper)
    225 
--> 226     return extract_region(spectrum, region)
    227 
    228 

.../specutils/manipulation/extract_spectral_region.py in extract_region(spectrum, region, return_single_spectrum)
    140     extracted_spectrum = []
    141     for subregion in region._subregions:
--> 142         left_index, right_index = _subregion_to_edge_pixels(subregion, spectrum)
    143 
    144         # If both indices are out of bounds then return an empty spectrum

.../specutils/manipulation/extract_spectral_region.py in _subregion_to_edge_pixels(subregion, spectrum)
     77         left_reg_in_spec_unit = left_func(subregion).to(spectral_axis.unit,
     78                                                         u.spectral())
---> 79         left_index = _edge_value_to_pixel(left_reg_in_spec_unit, spectrum, order, "left")
     80 
     81     # Right/upper side of sub region

.../specutils/manipulation/extract_spectral_region.py in _edge_value_to_pixel(edge_value, spectrum, order, side)
     39 
     40     except Exception as e:
---> 41         raise ValueError(f"Bound {edge_value}, could not be converted to pixel index"
     42                          f" using spectrum's WCS. Exception: {e}")
     43 

ValueError: Bound 0.0 pix, could not be converted to pixel index using spectrum's WCS.
Exception: SpectralTabular1D: Units of input 'x', pix (unknown),
could not be converted to required dimensionless input
```
repo: astropy/specutils
test_cmd: pytest --no-header -rA --tb=line --color=no -p no:cacheprovider -W ignore::DeprecationWarning
test_patch: diff --git a/specutils/tests/test_region_extract.py b/specutils/tests/test_region_extract.py
index 16389a2f..cd43bcca 100644
--- a/specutils/tests/test_region_extract.py
+++ b/specutils/tests/test_region_extract.py
@@ -3,7 +3,7 @@ import pytest
 
 import astropy.units as u
 from astropy.nddata import StdDevUncertainty
-from astropy.tests.helper import quantity_allclose
+from astropy.tests.helper import assert_quantity_allclose
 
 from ..spectra import Spectrum1D, SpectralRegion
 from ..manipulation import extract_region, extract_bounding_spectral_region, spectral_slab
@@ -34,7 +34,7 @@ def test_region_simple(simulated_spectra):
 
     sub_spectrum_flux_expected = np.array(FLUX_ARRAY)
 
-    assert quantity_allclose(sub_spectrum.flux.value, sub_spectrum_flux_expected)
+    assert_quantity_allclose(sub_spectrum.flux.value, sub_spectrum_flux_expected)
 
 
 def test_slab_simple(simulated_spectra):
@@ -48,7 +48,17 @@ def test_slab_simple(simulated_spectra):
 
     sub_spectrum_flux_expected = np.array(FLUX_ARRAY)
 
-    assert quantity_allclose(sub_spectrum.flux.value, sub_spectrum_flux_expected)
+    assert_quantity_allclose(sub_spectrum.flux.value, sub_spectrum_flux_expected)
+
+
+def test_slab_pixels():
+    range_iter = range(5)
+    spectrum = Spectrum1D(flux=np.stack(
+        [np.zeros((2, 2)) + i for i in range_iter], axis=-1) * u.nJy)
+    for i in range_iter:
+        sub_spectrum = spectral_slab(spectrum, i * u.pix, (i + 0.5) * u.pix)
+        assert_quantity_allclose(sub_spectrum.flux, i * u.nJy)
+        assert_quantity_allclose(sub_spectrum.spectral_axis, i * u.pix)
 
 
 def test_region_ghz(simulated_spectra):
@@ -61,7 +71,7 @@ def test_region_ghz(simulated_spectra):
 
     sub_spectrum_flux_expected = FLUX_ARRAY * u.mJy
 
-    assert quantity_allclose(sub_spectrum.flux, sub_spectrum_flux_expected)
+    assert_quantity_allclose(sub_spectrum.flux, sub_spectrum_flux_expected)
 
 
 def test_region_simple_check_ends(simulated_spectra):
@@ -140,7 +150,7 @@ def test_region_descending(simulated_spectra):
 
     sub_spectrum_flux_expected = np.array(FLUX_ARRAY)
 
-    assert quantity_allclose(sub_spectrum.flux.value, sub_spectrum_flux_expected)
+    assert_quantity_allclose(sub_spectrum.flux.value, sub_spectrum_flux_expected)
 
 
 def test_descending_spectral_axis(simulated_spectra):
@@ -151,12 +161,12 @@ def test_descending_spectral_axis(simulated_spectra):
     region = SpectralRegion(0.8*u.um, 0.6*u.um)
     sub_spectrum = extract_region(spectrum, region)
 
-    assert quantity_allclose(sub_spectrum.flux.value, sub_spectrum_flux_expected)
+    assert_quantity_allclose(sub_spectrum.flux.value, sub_spectrum_flux_expected)
 
     region = SpectralRegion(0.6*u.um, 0.8*u.um)
     sub_spectrum = extract_region(spectrum, region)
 
-    assert quantity_allclose(sub_spectrum.flux.value, sub_spectrum_flux_expected)
+    assert_quantity_allclose(sub_spectrum.flux.value, sub_spectrum_flux_expected)
 
 
 def test_region_two_sub(simulated_spectra):
@@ -171,10 +181,10 @@ def test_region_two_sub(simulated_spectra):
     sub_spectra = extract_region(spectrum, region)
 
     # Confirm the end points of the subspectra are correct
-    assert quantity_allclose(sub_spectra[0].spectral_axis[[0, -1]],
+    assert_quantity_allclose(sub_spectra[0].spectral_axis[[0, -1]],
                        [0.6035353535353536, 0.793939393939394]*u.um)
 
-    assert quantity_allclose(sub_spectra[1].spectral_axis[[0, -1]],
+    assert_quantity_allclose(sub_spectra[1].spectral_axis[[0, -1]],
                        [0.8661616161616162, 0.8858585858585859]*u.um)
 
     sub_spectrum_0_flux_expected = FLUX_ARRAY * u.mJy
@@ -182,15 +192,15 @@ def test_region_two_sub(simulated_spectra):
     sub_spectrum_1_flux_expected = [1337.65312465, 1263.48914109,
                                     1589.81797876, 1548.46068415]*u.mJy
 
-    assert quantity_allclose(sub_spectra[0].flux, sub_spectrum_0_flux_expected)
-    assert quantity_allclose(sub_spectra[1].flux, sub_spectrum_1_flux_expected)
+    assert_quantity_allclose(sub_spectra[0].flux, sub_spectrum_0_flux_expected)
+    assert_quantity_allclose(sub_spectra[1].flux, sub_spectrum_1_flux_expected)
 
     # also ensure this works if the multi-region is expressed as a single
     # Quantity
     region2 = SpectralRegion([(0.6, 0.8), (0.86, 0.89)]*u.um)
     sub_spectra2 = extract_region(spectrum, region2)
-    assert quantity_allclose(sub_spectra[0].flux, sub_spectra2[0].flux)
-    assert quantity_allclose(sub_spectra[1].flux, sub_spectra2[1].flux)
+    assert_quantity_allclose(sub_spectra[0].flux, sub_spectra2[0].flux)
+    assert_quantity_allclose(sub_spectra[1].flux, sub_spectra2[1].flux)
 
     # Check that the return_single_spectrum argument works properly
     concatenated_spectrum = extract_region(spectrum, region2, return_single_spectrum=True)
@@ -210,24 +220,24 @@ def test_bounding_region(simulated_spectra):
     extracted_spectrum = extract_bounding_spectral_region(spectrum, region)
 
     # Confirm the end points are correct
-    assert quantity_allclose(extracted_spectrum.spectral_axis[[0, -1]],
+    assert_quantity_allclose(extracted_spectrum.spectral_axis[[0, -1]],
                        [0.6035353535353536, 0.8858585858585859]*u.um)
 
-    flux_expected = [FLUX_ARRAY + [948.81864554, 1197.84859443, 1069.75268943,
+    flux_expected = (FLUX_ARRAY + [948.81864554, 1197.84859443, 1069.75268943,
                                    1118.27269184, 1301.7695563, 1206.62880648,
                                    1518.16549319, 1256.84259015, 1638.76791267,
                                    1562.05642302, 1337.65312465, 1263.48914109,
-                                   1589.81797876, 1548.46068415]]*u.mJy
+                                   1589.81797876, 1548.46068415])*u.mJy
 
-    assert quantity_allclose(extracted_spectrum.flux, flux_expected)
+    assert_quantity_allclose(extracted_spectrum.flux, flux_expected)
 
     # also ensure this works if the multi-region is expressed as a single
     # Quantity
     region2 = SpectralRegion([(0.6, 0.8), (0.86, 0.89)]*u.um)
     extracted_spectrum2 = extract_bounding_spectral_region(spectrum, region2)
-    assert quantity_allclose(extracted_spectrum2.spectral_axis[[0, -1]],
+    assert_quantity_allclose(extracted_spectrum2.spectral_axis[[0, -1]],
                        [0.6035353535353536, 0.8858585858585859]*u.um)
-    assert quantity_allclose(extracted_spectrum2.flux, flux_expected)
+    assert_quantity_allclose(extracted_spectrum2.flux, flux_expected)
 
 
 def test_extract_region_pixels():
@@ -237,7 +247,7 @@ def test_extract_region_pixels():
 
     extracted = extract_region(spectrum, region)
 
-    assert quantity_allclose(extracted.flux, [10, 11]*u.Jy)
+    assert_quantity_allclose(extracted.flux, [10, 11]*u.Jy)
 
 
 def test_extract_region_mismatched_units():
@@ -248,7 +258,7 @@ def test_extract_region_mismatched_units():
 
     extracted = extract_region(spectrum, region)
 
-    assert quantity_allclose(extracted.flux, [10, 11]*u.Jy)
+    assert_quantity_allclose(extracted.flux, [10, 11]*u.Jy)
 
 
 def test_linear_excise_invert_from_spectrum():
@@ -261,9 +271,9 @@ def test_linear_excise_invert_from_spectrum():
 
     excised_spec = linear_exciser(spec, exc_regs)
 
-    assert quantity_allclose(np.diff(excised_spec[50:60].flux),
+    assert_quantity_allclose(np.diff(excised_spec[50:60].flux),
                              np.diff(excised_spec[51:61].flux))
-    assert quantity_allclose(np.diff(excised_spec[80:90].flux),
+    assert_quantity_allclose(np.diff(excised_spec[80:90].flux),
                              np.diff(excised_spec[81:91].flux))