Fix pixel region index retrieval

specutils/manipulation/extract_spectral_region.py

@@ -75,7 +75,19 @@ def _subregion_to_edge_pixels(subregion, spectrum):
         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))
+            if order != "ascending":
+                raise ValueError("u.pix handling should always be ascending")
+            # If the lower bound is larger than the largest value, immediately return nothing
+            # Assuming ascending, both bounds are "out of bounds"
+            if subregion[0] > spectral_axis[-1]:
+                return None, None
+            else:
+                # Get index of closest value
+                # See https://stackoverflow.com/a/26026189 if performance becomes an issue
+                left_index = np.nanargmin((np.abs(spectral_axis - subregion[0])))
+                # Ensure index is inclusive of region bounds
+                if (spectral_axis[left_index] > subregion[0]) and (left_index >= 1):
+                    left_index -= 1
     else:
         # Convert lower value to spectrum spectral_axis units
         left_reg_in_spec_unit = left_func(subregion).to(spectral_axis.unit,
@@ -87,7 +99,19 @@ def _subregion_to_edge_pixels(subregion, spectrum):
         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
+            if order != "ascending":
+                raise ValueError("u.pix handling should always be ascending")
+            # If upper bound is smaller than smallest value, immediately return nothing
+            # Assuming ascending, both bounds are "out of bounds"
+            if subregion[1] < spectral_axis[0]:
+                return None, None
+            else:
+                # Get index of closest value
+                # See https://stackoverflow.com/a/26026189 if performance becomes an issue
+                right_index = np.nanargmin((np.abs(spectral_axis - subregion[1])))
+                # Ensure index is inclusive of region bounds
+                if (spectral_axis[right_index] < subregion[1]) and (right_index < len(spectral_axis)):
+                    right_index += 1
     else:
         # Convert upper value to spectrum spectral_axis units
         right_reg_in_spec_unit = right_func(subregion).to(spectral_axis.unit,

specutils/tests/test_region_extract.py

@@ -33,6 +33,84 @@ def test_region_simple(simulated_spectra):
     assert_quantity_allclose(sub_spectrum.flux.value, sub_spectrum_flux_expected)
 
 
+def test_pixel_spectralaxis_extraction():
+    """
+    Tests a region extraction on a Spectrum1D with a u.pix defined spectral axis
+    """
+    flux_unit = u.dimensionless_unscaled
+    spec_unit = u.pix
+
+    spec1d = Spectrum1D(spectral_axis=np.arange(5100, 5300)*spec_unit,
+                        flux=np.random.randn(200)*flux_unit)
+
+    # Case 1: Region is safely within the bounds of a continuously-defined Spec1D
+    # Region is intentionally defined "between values" to test for index rounding
+    region = SpectralRegion.from_center(center=5200.5*spec_unit,
+                                        width=100*spec_unit)
+    extracted_spec1d = extract_region(spec1d, region)
+
+    # Extraction is of the right shape and value
+    assert extracted_spec1d.shape == (101,)
+    assert_quantity_allclose(extracted_spec1d.spectral_axis,
+                             spec1d.spectral_axis[50:151])
+    assert_quantity_allclose(extracted_spec1d.flux,
+                             spec1d.flux[50:151])
+
+    # Follow Python slicing conventions:
+    # Lower slice is inclusive
+    assert extracted_spec1d.spectral_axis[0].quantity <= region.lower
+    # Upper slice is exclusive
+    assert extracted_spec1d.spectral_axis[-1].quantity < region.upper
+
+    # Case 2: Region is outside the lower bounds of the Spec1D
+    region2 = SpectralRegion.from_center(center=spec1d.spectral_axis[0].quantity,
+                                         width=100*spec_unit)
+    extracted_spec1d_2 = extract_region(spec1d, region2)
+
+    assert extracted_spec1d_2.shape == (50,)  # Included upper bound is exclusive
+    assert_quantity_allclose(extracted_spec1d_2.spectral_axis,
+                             spec1d.spectral_axis[0:50])
+    assert_quantity_allclose(extracted_spec1d_2.flux,
+                             spec1d.flux[0:50])
+
+    # Case 3: Region is outside the upper bounds of the Spec1D
+    region3 = SpectralRegion.from_center(center=spec1d.spectral_axis[-1].quantity,
+                                         width=100*spec_unit)
+    extracted_spec1d_3 = extract_region(spec1d, region3)
+
+    assert extracted_spec1d_3.shape == (51,)  # Included lower bound is inclusive (+1)
+    assert_quantity_allclose(extracted_spec1d_3.spectral_axis,
+                             spec1d.spectral_axis[149:])
+    assert_quantity_allclose(extracted_spec1d_3.flux,
+                             spec1d.flux[149:])
+
+    # Case 4: Compound region with the two definitions above
+    extracted_spec1d_4 = extract_region(spec1d, (region2+region3))
+
+    # Each extracted part of compound region should be identical to the case 2 and 3 extractions
+    # If they're identical, then the previous value checks for case 2/3 should cover this as well
+    assert len(extracted_spec1d_4) == 2
+
+    for original_case_spectra, compound_spectra in [(extracted_spec1d_2, extracted_spec1d_4[0]),
+                                                    (extracted_spec1d_3, extracted_spec1d_4[1])]:
+        assert original_case_spectra.shape == compound_spectra.shape
+        assert_quantity_allclose(original_case_spectra.spectral_axis,
+                                 compound_spectra.spectral_axis)
+        assert_quantity_allclose(original_case_spectra.flux,
+                                 compound_spectra.flux)
+
+    # Case 5: Region is entirely outside bounds of the spectra should return nothing
+    upper_bound = spec1d.spectral_axis[-1].quantity
+    region = SpectralRegion((upper_bound + 10*spec_unit), (upper_bound + 100*spec_unit))
+    extracted_spec1d = extract_region(spec1d, region)
+    assert extracted_spec1d.shape == (0,)
+
+    lower_bound = spec1d.spectral_axis[0].quantity
+    region = SpectralRegion((lower_bound - 100*spec_unit), (lower_bound - 10*spec_unit))
+    extracted_spec1d = extract_region(spec1d, region)
+    assert extracted_spec1d.shape == (0,)
+
+
 def test_slab_simple(simulated_spectra):
     np.random.seed(42)