Merge branch 'main' into mwmvisit-boss-fix

CHANGES.rst

@@ -1,4 +1,4 @@
-1.18.0 (unreleased)
+1.19.0 (unreleased)
 -------------------
 
 New Features
@@ -7,13 +7,33 @@ New Features
 Bug Fixes
 ^^^^^^^^^
 
+Other Changes and Additions
+^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+1.18.0 (2024-10-16)
+-------------------
+
+New Features
+^^^^^^^^^^^^
+
+- New ``Spectrum1D.with_spectral_axis_and_flux_units`` method to convert both
+  spectral axis and flux units at the same time. [#1184]
+
+Bug Fixes
+^^^^^^^^^
+
 - Fixed ``Spectrum1D.with_flux_unit()`` not converting uncertainty along
   with flux unit. [#1181]
 
-- Fixed ``mwmVisit`` SDSS-V ``Spectrum1D`` and ``SpectrumList`` default loader being unable to load files containing only BOSS instrument spectra. [#1185]
+- Fixed ``mwmVisit`` SDSS-V ``Spectrum1D`` and ``SpectrumList`` default loader 
+  being unable to load files containing only BOSS instrument spectra. [#1185]
 
 - Fixed automatic format detection for SDSS-V ``SpectrumList`` default loaders. [#1185]
 
+- 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/manipulation/resample.py

@@ -461,8 +461,8 @@ def resample1d(self, orig_spectrum, fin_spec_axis):
         if self.extrapolation_treatment == 'zero_fill':
             fill_val = 0
 
-        origedges = orig_spectrum.spectral_axis.bin_edges
-        off_edges = (fin_spec_axis < origedges[0]) | (origedges[-1] < fin_spec_axis)
+        orig_edges = orig_spectrum.spectral_axis.bin_edges
+        off_edges = (fin_spec_axis < np.min(orig_edges)) | (np.max(orig_edges) < fin_spec_axis)
         out_flux_val[off_edges] = fill_val
         if new_unc is not None:
             new_unc.array[off_edges] = fill_val

specutils/spectra/spectrum_mixin.py

@@ -86,6 +86,26 @@ def new_flux_unit(self, unit, equivalencies=None, suppress_conversion=False):
         return self.with_flux_unit(unit, equivalencies=equivalencies,
                                   suppress_conversion=suppress_conversion)
 
+    def _convert_flux(self, unit, equivalencies=None, suppress_conversion=False):
+        """This is always done in-place.
+        Also see :meth:`with_flux_unit`."""
+
+        if not suppress_conversion:
+            if equivalencies is None:
+                equivalencies = eq.spectral_density(self.spectral_axis)
+
+            new_data = self.flux.to(unit, equivalencies=equivalencies)
+
+            self._data = new_data.value
+            self._unit = new_data.unit
+        else:
+            self._unit = u.Unit(unit)
+
+        if self.uncertainty is not None:
+            self.uncertainty = StdDevUncertainty(
+                self.uncertainty.represent_as(StdDevUncertainty).quantity.to(
+                    unit, equivalencies=equivalencies))
+
     def with_flux_unit(self, unit, equivalencies=None, suppress_conversion=False):
         """Returns a new spectrum with a different flux unit.
         If uncertainty is defined, it will be converted to
@@ -107,30 +127,14 @@ def with_flux_unit(self, unit, equivalencies=None, suppress_conversion=False):
 
         Returns
         -------
-        `~specutils.Spectrum1D`
+        new_spec : `~specutils.Spectrum1D`
             A new spectrum with the converted flux array
             (and uncertainty, if applicable).
 
         """
         new_spec = deepcopy(self)
-
-        if not suppress_conversion:
-            if equivalencies is None:
-                equivalencies = eq.spectral_density(self.spectral_axis)
-
-            new_data = self.flux.to(
-                unit, equivalencies=equivalencies)
-
-            new_spec._data = new_data.value
-            new_spec._unit = new_data.unit
-        else:
-            new_spec._unit = u.Unit(unit)
-
-        if self.uncertainty is not None:
-            new_spec.uncertainty = StdDevUncertainty(
-                self.uncertainty.represent_as(StdDevUncertainty).quantity.to(
-                    unit, equivalencies=equivalencies))
-
+        new_spec._convert_flux(
+            unit, equivalencies=equivalencies, suppress_conversion=suppress_conversion)
         return new_spec
 
     @property
@@ -175,7 +179,7 @@ def velocity(self):
 
         Returns
         -------
-        ~`astropy.units.Quantity`
+        new_data : `~astropy.units.Quantity`
             The converted dispersion array in the new dispersion space.
         """
         if self.rest_value is None:
@@ -202,8 +206,7 @@ def with_spectral_unit(self, unit, velocity_convention=None,
         self.with_spectral_axis_unit(unit, velocity_convention=velocity_convention,
                                      rest_value=rest_value)
 
-    def with_spectral_axis_unit(self, unit, velocity_convention=None,
-                           rest_value=None):
+    def with_spectral_axis_unit(self, unit, velocity_convention=None, rest_value=None):
         """
         Returns a new spectrum with a different spectral axis unit. Note that this creates a new
         object using the converted spectral axis and thus drops the original WCS, if it existed,
@@ -230,11 +233,9 @@ def with_spectral_axis_unit(self, unit, velocity_convention=None,
                      even if your spectrum has air wavelength units
 
         """
-
-        velocity_convention = velocity_convention if velocity_convention is not None else self.velocity_convention # noqa
+        velocity_convention = velocity_convention if velocity_convention is not None else self.velocity_convention  # noqa
         rest_value = rest_value if rest_value is not None else self.rest_value
-        unit = self._new_wcs_argument_validation(unit, velocity_convention,
-                                                 rest_value)
+        unit = self._new_wcs_argument_validation(unit, velocity_convention, rest_value)
 
         # Store the original unit information and WCS for posterity
         meta = deepcopy(self._meta)
@@ -252,6 +253,24 @@ def with_spectral_axis_unit(self, unit, velocity_convention=None,
         return self.__class__(flux=self.flux, spectral_axis=new_spectral_axis, meta=meta,
                               uncertainty=self.uncertainty, mask=self.mask)
 
+    def with_spectral_axis_and_flux_units(self, spectral_axis_unit, flux_unit,
+                                          velocity_convention=None, rest_value=None,
+                                          flux_equivalencies=None, suppress_flux_conversion=False):
+        """Perform :meth:`with_spectral_axis_unit` and :meth:`with_flux_unit` together.
+        See the respective methods for input and output definitions.
+
+        Returns
+        -------
+        new_spec : `~specutils.Spectrum1D`
+            Spectrum in requested units.
+
+        """
+        new_spec = self.with_spectral_axis_unit(
+            spectral_axis_unit, velocity_convention=velocity_convention, rest_value=rest_value)
+        new_spec._convert_flux(
+            flux_unit, equivalencies=flux_equivalencies, suppress_conversion=suppress_flux_conversion)
+        return new_spec
+
     def _new_wcs_argument_validation(self, unit, velocity_convention,
                                      rest_value):
         # Allow string specification of units, for example

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)])