Spectral Region Enhancement

docs/index.rst

@@ -30,6 +30,7 @@ Using specutils
    arithmetic
    smoothing
    fitting
+   spectral_regions
 
 .. toctree::
     :maxdepth: 1

docs/spectral_regions.rst

@@ -0,0 +1,196 @@
+================
+Spectral Regions
+================
+
+A spectral region may be defined and may encompass one, or more,
+sub-regions. They are independent of a `~specutils.Spectrum1D`
+object.  There, one may define a spectral region and it may be used
+on `~specutils.Spectrum1D` objects of different dispersion
+sampling.
+
+Spectral regions can be defined either as a single region by passing 
+two `~astropy.units.Quantity`'s or by passing a list of 2-tuples.  By
+definition the `~specutils.SpectralRegion` will be ordered
+based on the lower bound of each 2-tuple.
+
+.. code-block:: python
+
+    >>> from astropy import units as u
+    >>> from specutils.spectra import SpectralRegion
+
+    >>> sr = SpectralRegion(0.45*u.um, 0.6*u.um)
+    >>> sr_two = SpectralRegion([(0.45*u.um, 0.6*u.um), (0.8*u.um, 0.9*u.um)])
+
+`~specutils.SpectralRegion` can be combined by using the '+' operator:
+
+.. code-block:: python
+
+    >>> from astropy import units as u
+    >>> from specutils.spectra import SpectralRegion
+    >>> sr = SpectralRegion(0.45*u.um, 0.6*u.um) + SpectralRegion(0.8*u.um, 0.9*u.um)
+
+Regions can also be added in place:
+
+.. code-block:: python
+
+    >>> from astropy import units as u
+    >>> from specutils.spectra import SpectralRegion
+
+    >>> sr1 = SpectralRegion(0.45*u.um, 0.6*u.um)
+    >>> sr2 = SpectralRegion(0.8*u.um, 0.9*u.um)
+    >>> sr1 += sr2
+
+Regions can be sliced by indexing by an integer or by a range:
+
+.. code-block:: python
+
+    >>> from astropy import units as u
+    >>> from specutils.spectra import SpectralRegion
+
+    >>> sr = SpectralRegion(0.15*u.um, 0.2*u.um) + SpectralRegion(0.3*u.um, 0.4*u.um) +\
+    ...      SpectralRegion(0.45*u.um, 0.6*u.um) + SpectralRegion(0.8*u.um, 0.9*u.um) +\
+    ...      SpectralRegion(1.0*u.um, 1.2*u.um) + SpectralRegion(1.3*u.um, 1.5*u.um)
+
+    >>> # Get on spectral region (returns a SpectralRegion instance)
+    >>> sone = sr1[0]
+
+    >>> # Slice spectral region.
+    >>> subsr = sr[3:5]
+    >>> # SpectralRegion: 0.8 um - 0.9 um, 1.0 um - 1.2 um
+
+The lower and upper bounds on a region are accessed by calling lower
+or upper. The lower bound of a `~specutils.SpectralRegion` is the
+minimum of the lower bounds of each sub-region and the upper bound is the
+maximum of the upper bounds:
+
+.. code-block:: python
+
+    >>> from astropy import units as u
+    >>> from specutils.spectra import SpectralRegion
+
+    >>> sr = SpectralRegion(0.15*u.um, 0.2*u.um) + SpectralRegion(0.3*u.um, 0.4*u.um) +\
+    ...      SpectralRegion(0.45*u.um, 0.6*u.um) + SpectralRegion(0.8*u.um, 0.9*u.um) +\
+    ...      SpectralRegion(1.0*u.um, 1.2*u.um) + SpectralRegion(1.3*u.um, 1.5*u.um)
+
+    >>> # Bounds on the spectral region (most minimum and maximum bound)
+    >>> print(sr.bounds) #doctest:+SKIP
+    (<Quantity 0.15 um>, <Quantity 1.5 um>)
+
+    >>> # Lower bound on the spectral region (most minimum)
+    >>> sr.lower #doctest:+SKIP
+    <Quantity 0.15 um>
+
+    >>> sr.upper #doctest:+SKIP
+    <Quantity 1.5 um>
+
+    >>> # Lower bound on one element of the spectral region.
+    >>> sr[3].lower #doctest:+SKIP
+    <Quantity 0.8 um>
+
+One can also delete a sub-region:
+.. code-block:: python
+
+    >>> from astropy import units as u
+    >>> from specutils.spectra import SpectralRegion
+
+    >>> sr = SpectralRegion(0.15*u.um, 0.2*u.um) + SpectralRegion(0.3*u.um, 0.4*u.um) +\
+    ...      SpectralRegion(0.45*u.um, 0.6*u.um) + SpectralRegion(0.8*u.um, 0.9*u.um) +\
+    ...      SpectralRegion(1.0*u.um, 1.2*u.um) + SpectralRegion(1.3*u.um, 1.5*u.um)
+
+    >>> del sr[1]
+    >>> sr #doctest:+SKIP
+    [(<Quantity 0.15 um>, <Quantity 0.2 um>),
+     (<Quantity 0.45 um>, <Quantity 0.6 um>),
+     (<Quantity 0.8 um>, <Quantity 0.9 um>),
+     (<Quantity 1. um>, <Quantity 1.2 um>),
+     (<Quantity 1.3 um>, <Quantity 1.5 um>)]
+
+There is also the ability to iterate:
+
+.. code-block:: python
+
+    >>> from astropy import units as u
+    >>> from specutils.spectra import SpectralRegion
+
+    >>> sr = SpectralRegion(0.15*u.um, 0.2*u.um) + SpectralRegion(0.3*u.um, 0.4*u.um) +\
+    ...      SpectralRegion(0.45*u.um, 0.6*u.um) + SpectralRegion(0.8*u.um, 0.9*u.um) +\
+    ...      SpectralRegion(1.0*u.um, 1.2*u.um) + SpectralRegion(1.3*u.um, 1.5*u.um)
+
+    >>> for s in sr:
+    ...     print(s.lower) #doctest:+SKIP
+    SpectralRegion: 0.15 um - 0.2 um
+    SpectralRegion: 0.3 um - 0.4 um
+    SpectralRegion: 0.45 um - 0.6 um
+    SpectralRegion: 0.8 um - 0.9 um
+    SpectralRegion: 1.0 um - 1.2 um
+    SpectralRegion: 1.3 um - 1.5 um
+
+
+And, lastly, there is the ability to invert a `~specutils.SpectralRegion` given a 
+lower and upper bound. For example, if a set of ranges are defined each defining a range
+around lines, then calling invert will return a `~specutils.SpectralRegion` that
+defines the baseline/noise regions:
+
+.. code-block:: python
+
+    >>> from astropy import units as u
+    >>> from specutils.spectra import SpectralRegion
+
+    >>> sr = SpectralRegion(0.15*u.um, 0.2*u.um) + SpectralRegion(0.3*u.um, 0.4*u.um) +\
+    ...      SpectralRegion(0.45*u.um, 0.6*u.um) + SpectralRegion(0.8*u.um, 0.9*u.um) +\
+    ...      SpectralRegion(1.0*u.um, 1.2*u.um) + SpectralRegion(1.3*u.um, 1.5*u.um)
+
+    >>> sr_inverted = sr.invert(0.05*u.um, 3*u.um)
+    >>> sr_inverted #doctest:+SKIP
+    SpectralRegion: 0.05 um - 0.15 um, 0.2 um - 0.3 um, 0.4 um - 0.45 um, 0.6 um - 0.8 um, 0.9 um - 1.0 um, 1.2 um - 1 .3 um, 1.5 um - 3.0 um
+
+
+Region Extraction
+-----------------
+
+Given a `~specutils.SpectralRegion`, one can extract a sub-spectrum
+from a `~specutils.Spectrum1D` object. If the `~specutils.SpectralRegion`
+has multiple sub-regions then a list of `~specutils.Spectrum1D` objects will
+be returned.
+
+An example of a single sub-region `~specutils.SpectralRegion`:
+
+.. code-block:: python
+
+    >>> from astropy import units as u
+    >>> import numpy as np
+    >>> from specutils import Spectrum1D, SpectralRegion
+    >>> from specutils.manipulation import extract_region
+
+    >>> region = SpectralRegion(8*u.nm, 22*u.nm)
+    >>> spectrum = Spectrum1D(spectral_axis=np.arange(1, 50) * u.nm, flux=np.random.sample(49)*u.Jy)
+    >>> sub_spectrum = extract_region(spectrum, region)
+    >>> sub_spectrum.spectral_axis
+    <Quantity [ 8.,  9., 10., 11., 12., 13., 14., 15., 16., 17., 18., 19., 20.,
+               21., 22.] nm>
+
+
+An example of a multiple sub-region `~specutils.SpectralRegion`:
+
+.. code-block:: python
+
+    >>> from astropy import units as u
+    >>> import numpy as np
+    >>> from specutils import Spectrum1D, SpectralRegion
+    >>> from specutils.manipulation import extract_region
+
+    >>> region = SpectralRegion([(8*u.nm, 22*u.nm), (34*u.nm, 40*u.nm)])
+    >>> spectrum = Spectrum1D(spectral_axis=np.arange(1, 50) * u.nm, flux=np.random.sample(49)*u.Jy)
+    >>> sub_spectra = extract_region(spectrum, region)
+    >>> sub_spectra[0].spectral_axis
+    <Quantity [ 8.,  9., 10., 11., 12., 13., 14., 15., 16., 17., 18., 19., 20.,
+               21., 22.] nm>
+    >>> sub_spectra[1].spectral_axis
+    <Quantity [34., 35., 36., 37., 38., 39., 40.] nm>
+
+
+Reference/API
+-------------
+
+.. automodapi:: specutils.spectra.spectral_region
+    :no-heading:

specutils/analysis/snr.py

@@ -1,5 +1,6 @@
 import numpy as np
 from ..spectra import SpectralRegion
+from ..manipulation import extract_region
 
 __all__ = ['snr']
 
@@ -72,7 +73,7 @@ def _snr_single_region(spectrum, region=None):
     """
 
     if region is not None:
-        calc_spectrum = region.extract(spectrum)
+        calc_spectrum = extract_region(spectrum, region)
     else:
         calc_spectrum = spectrum
 

specutils/manipulation/__init__.py

@@ -2,3 +2,4 @@
 
 from .smoothing import *  # noqa
 from .estimate_uncertainty import *  # noqa
+from .extract_spectral_region import *  # noqa

specutils/manipulation/extract_spectral_region.py

@@ -0,0 +1,99 @@
+import numpy as np
+
+__all__ = ['extract_region']
+
+
+def to_pixel(subregion, spectrum):
+    """
+    Calculate and return the left and right indices defined
+    by the lower and upper bounds and based on the input
+    `~specutils.spectra.spectrum1d.Spectrum1D`. The left and right indices will
+    be returned.
+
+    Parameters
+    ----------
+    spectrum: `~specutils.spectra.spectrum1d.Spectrum1D`
+        The spectrum object from which the region will be extracted.
+
+    Returns
+    -------
+    left_index, right_index: int, int
+        Left and right indices defined by the lower and upper bounds.
+
+    """
+
+    try:
+        left_index = int(np.ceil(spectrum.wcs.world_to_pixel([subregion[0]])))
+    except Exception as e:
+        left_index = None
+
+    try:
+        right_index = int(np.floor(spectrum.wcs.world_to_pixel([subregion[1]]))) + 1
+    except Exception as e:
+        right_index = None
+
+    return left_index, right_index
+
+
+def extract_region(spectrum, region):
+    """
+    Extract a region from the input `~specutils.spectra.spectrum1d.Spectrum1D`
+    defined by the lower and upper bounds defined by this SpectralRegion
+    instance.  The extracted region will be returned.
+
+    Parameters
+    ----------
+    spectrum: `~specutils.spectra.spectrum1d.Spectrum1D`
+        The spectrum object from which the region will be extracted.
+
+    Returns
+    -------
+    spectrum: `~specutils.spectra.spectrum1d.Spectrum1D`
+        Excised spectrum.
+
+    Notes
+    -----
+    The region extracted is a discrete subset of the input spectrum. No interpolation is done
+    on the left and right side of the spectrum.
+
+    The region is assumed to be a closed interval (as opposed to Python which is open
+    on the upper end).  For example:
+
+        Given:
+           A ``spectrum`` with spectral_axis of ``[0.1, 0.2, 0.3, 0.4, 0.5, 0.6]*u.um``.
+
+           A ``region`` defined as ``SpectralRegion(0.2*u.um, 0.5*u.um)``
+
+        And we calculate ``sub_spectrum = extract_region(spectrum, region)``, then the ``sub_spectrum``
+        spectral axis will be ``[0.2, 0.3, 0.4, 0.5] * u.um``.
+
+    """
+
+    extracted_spectrum = []
+    for subregion in region._subregions:
+        left_index, right_index = to_pixel(subregion, spectrum)
+
+        # If both indices are out of bounds then return None
+        if left_index is None and right_index is None:
+            extracted_spectrum.append(None)
+        else:
+
+            # If only one index is out of bounds then set it to
+            # the lower or upper extent
+            if left_index is None:
+                left_index = 0
+
+            if right_index is None:
+                right_index = len(spectrum.spectral_axis)
+
+            if left_index > right_index:
+                left_index, right_index = right_index, left_index
+
+            extracted_spectrum.append(spectrum[left_index:right_index])
+
+    # If there is only one subregion in the region then we will
+    # just return a spectrum.
+    if len(region) == 1:
+        extracted_spectrum = extracted_spectrum[0]
+
+    return extracted_spectrum