diff --git a/specutils/manipulation/extract_spectral_region.py b/specutils/manipulation/extract_spectral_region.py index a01c5b041..18290bb7d 100644 --- a/specutils/manipulation/extract_spectral_region.py +++ b/specutils/manipulation/extract_spectral_region.py @@ -72,10 +72,21 @@ def _subregion_to_edge_pixels(subregion, spectrum): # Left/lower side of sub region if subregion[0].unit.is_equivalent(u.pix): + # Pixel handling assumes always ascending 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 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, @@ -84,10 +95,21 @@ def _subregion_to_edge_pixels(subregion, spectrum): # Right/upper side of sub region if subregion[1].unit.is_equivalent(u.pix): + # Pixel handling assumes always ascending 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 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, diff --git a/specutils/spectra/spectral_axis.py b/specutils/spectra/spectral_axis.py index 5ab61956c..8dee94783 100644 --- a/specutils/spectra/spectral_axis.py +++ b/specutils/spectra/spectral_axis.py @@ -26,6 +26,13 @@ class SpectralAxis(SpectralCoord): def __new__(cls, value, *args, bin_specification="centers", **kwargs): + # Enforce pixel axes are ascending + if ((type(value) is u.quantity.Quantity) and + (value.size > 1) and + (value.unit is u.pix) and + (value[-1] <= value[0])): + raise ValueError("u.pix spectral axes should always be ascending") + # Convert to bin centers if bin edges were given, since SpectralCoord # only accepts centers if bin_specification == "edges": diff --git a/specutils/tests/test_region_extract.py b/specutils/tests/test_region_extract.py index 32a8426f0..f40eb5b9f 100644 --- a/specutils/tests/test_region_extract.py +++ b/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) diff --git a/specutils/tests/test_spectral_axis.py b/specutils/tests/test_spectral_axis.py index 65832aba6..d5013f0f2 100644 --- a/specutils/tests/test_spectral_axis.py +++ b/specutils/tests/test_spectral_axis.py @@ -217,3 +217,16 @@ def test_no_change_redshift(): assert type(spec.spectral_axis) == SpectralAxis assert_quantity_allclose(spec.wavelength, wave) + + +def test_pixel_descending_error(): + ''' + Spectral axes of pixel units must always be ascending. + This test checks that an error is thrown if one is provided descending + ''' + flux_unit = u.dimensionless_unscaled + spec_unit = u.pix + + with pytest.raises(ValueError, match="u.pix spectral axes should always be ascending"): + Spectrum1D(spectral_axis=(np.arange(5100, 5300)[::-1])*spec_unit, + flux=np.random.randn(200)*flux_unit)