Fix mixed unit spectral region extraction (#1187)

* Fix extraction when only one of region/spectrum are in a reversed spectral axis unit

* Changelog

* Add test

* Get rid of more pointless random calls

* Update specutils/tests/test_region_extract.py

Co-authored-by: P. L. Lim <[email protected]>

* Better test improvement

---------

Co-authored-by: P. L. Lim <[email protected]>

CHANGES.rst

@@ -13,6 +13,9 @@ Bug Fixes
 - Fixed ``Spectrum1D.with_flux_unit()`` not converting uncertainty along
   with flux unit. [#1181]
 
+- Fixed extracting a spectral region when one of spectrum/region is in wavelength
+  and the other is in frequency units. [#1187]
+
 Other Changes and Additions
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^
 

specutils/manipulation/extract_spectral_region.py

@@ -103,7 +103,11 @@ def _subregion_to_edge_pixels(subregion, spectrum):
 
         right_index = _edge_value_to_pixel(right_reg_in_spec_unit, spectrum, order, "right")
 
-    return left_index, right_index
+    # If the spectrum is in wavelength and region is in Hz (for example), these still might be reversed
+    if left_index < right_index:
+        return left_index, right_index
+    else:
+        return right_index, left_index
 
 
 def extract_region(spectrum, region, return_single_spectrum=False):

specutils/tests/test_region_extract.py

@@ -18,10 +18,9 @@
 
 
 def test_region_simple(simulated_spectra):
-    np.random.seed(42)
 
     spectrum = simulated_spectra.s1_um_mJy_e1
-    uncertainty = StdDevUncertainty(0.1*np.random.random(len(spectrum.flux))*u.mJy)
+    uncertainty = StdDevUncertainty(np.full(spectrum.flux.shape, 0.1) * u.mJy)
     spectrum.uncertainty = uncertainty
 
     region = SpectralRegion(0.6*u.um, 0.8*u.um)
@@ -112,10 +111,8 @@ def test_pixel_spectralaxis_extraction():
 
 
 def test_slab_simple(simulated_spectra):
-    np.random.seed(42)
-
     spectrum = simulated_spectra.s1_um_mJy_e1
-    uncertainty = StdDevUncertainty(0.1*np.random.random(len(spectrum.flux))*u.mJy)
+    uncertainty = StdDevUncertainty(np.full(spectrum.flux.shape, 0.1) * u.mJy)
     spectrum.uncertainty = uncertainty
 
     sub_spectrum = spectral_slab(spectrum, 0.6*u.um, 0.8*u.um)
@@ -148,6 +145,16 @@ def test_region_ghz(simulated_spectra):
     assert_quantity_allclose(sub_spectrum.flux, sub_spectrum_flux_expected)
 
 
+def test_region_ghz_spectrum_wave(simulated_spectra):
+    spectrum = simulated_spectra.s1_um_mJy_e1
+    region = SpectralRegion(499654.09666667*u.GHz, 374740.5725*u.GHz)
+
+    sub_spectrum = extract_region(spectrum, region)
+
+    sub_spectrum_flux_expected = FLUX_ARRAY * u.mJy
+    assert_quantity_allclose(sub_spectrum.flux, sub_spectrum_flux_expected)
+
+
 def test_region_simple_check_ends(simulated_spectra):
     np.random.seed(42)
 
@@ -168,13 +175,11 @@ def test_region_simple_check_ends(simulated_spectra):
     assert sub_spectrum.spectral_axis.value[-1] == 25
 
 
-def test_region_empty(simulated_spectra):
-    np.random.seed(42)
-
+def test_region_empty():
     empty_spectrum = Spectrum1D(spectral_axis=[]*u.um, flux=[]*u.Jy)
 
     # Region past upper range of spectrum
-    spectrum = Spectrum1D(spectral_axis=np.linspace(1, 25, 25)*u.um, flux=np.random.random(25)*u.Jy)
+    spectrum = Spectrum1D(spectral_axis=np.linspace(1, 25, 25)*u.um, flux=np.ones(25)*u.Jy)
     region = SpectralRegion(28*u.um, 30*u.um)
     sub_spectrum = extract_region(spectrum, region)
 
@@ -185,7 +190,7 @@ def test_region_empty(simulated_spectra):
     assert sub_spectrum.flux.unit == empty_spectrum.flux.unit
 
     # Region below lower range of spectrum
-    spectrum = Spectrum1D(spectral_axis=np.linspace(1, 25, 25)*u.um, flux=np.random.random(25)*u.Jy)
+    spectrum = Spectrum1D(spectral_axis=np.linspace(1, 25, 25)*u.um, flux=np.ones(25)*u.Jy)
     region = SpectralRegion(0.1*u.um, 0.3*u.um)
     sub_spectrum = extract_region(spectrum, region)
 
@@ -212,10 +217,8 @@ def test_region_empty(simulated_spectra):
 
 
 def test_region_descending(simulated_spectra):
-    np.random.seed(42)
-
     spectrum = simulated_spectra.s1_um_mJy_e1
-    uncertainty = StdDevUncertainty(0.1*np.random.random(len(spectrum.flux))*u.mJy)
+    uncertainty = StdDevUncertainty(np.full(spectrum.flux.shape, 0.1) * u.mJy)
     spectrum.uncertainty = uncertainty
 
     region = SpectralRegion(0.8*u.um, 0.6*u.um)
@@ -244,10 +247,8 @@ def test_descending_spectral_axis(simulated_spectra):
 
 
 def test_region_two_sub(simulated_spectra):
-    np.random.seed(42)
-
     spectrum = simulated_spectra.s1_um_mJy_e1
-    uncertainty = StdDevUncertainty(0.1*np.random.random(len(spectrum.flux))*u.mJy)
+    uncertainty = StdDevUncertainty(np.full(spectrum.flux.shape, 0.1) * u.mJy)
     spectrum.uncertainty = uncertainty
 
     region = SpectralRegion([(0.6*u.um, 0.8*u.um), (0.86*u.um, 0.89*u.um)])
@@ -284,10 +285,8 @@ def test_region_two_sub(simulated_spectra):
 
 
 def test_bounding_region(simulated_spectra):
-    np.random.seed(42)
-
     spectrum = simulated_spectra.s1_um_mJy_e1
-    uncertainty = StdDevUncertainty(0.1*np.random.random(len(spectrum.flux))*u.mJy)
+    uncertainty = StdDevUncertainty(np.full(spectrum.flux.shape, 0.1) * u.mJy)
     spectrum.uncertainty = uncertainty
 
     region = SpectralRegion([(0.6*u.um, 0.8*u.um), (0.86*u.um, 0.89*u.um)])