diff --git a/docs/index.rst b/docs/index.rst
index 2be718a71..3d3b7d869 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -30,6 +30,7 @@ Using specutils
    arithmetic
    smoothing
    fitting
+   spectral_regions
 
 .. toctree::
     :maxdepth: 1
diff --git a/docs/spectral_regions.rst b/docs/spectral_regions.rst
new file mode 100644
index 000000000..2d4f037e2
--- /dev/null
+++ b/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:
diff --git a/specutils/analysis/snr.py b/specutils/analysis/snr.py
index 5079ae603..940a0290b 100644
--- a/specutils/analysis/snr.py
+++ b/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
 
diff --git a/specutils/manipulation/__init__.py b/specutils/manipulation/__init__.py
index 7bc4d58b8..1f46be724 100644
--- a/specutils/manipulation/__init__.py
+++ b/specutils/manipulation/__init__.py
@@ -2,3 +2,4 @@
 
 from .smoothing import *  # noqa
 from .estimate_uncertainty import *  # noqa
+from .extract_spectral_region import *  # noqa
diff --git a/specutils/manipulation/extract_spectral_region.py b/specutils/manipulation/extract_spectral_region.py
new file mode 100644
index 000000000..b6118ff99
--- /dev/null
+++ b/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
diff --git a/specutils/spectra/spectral_region.py b/specutils/spectra/spectral_region.py
index 66c12d416..e83c76855 100644
--- a/specutils/spectra/spectral_region.py
+++ b/specutils/spectra/spectral_region.py
@@ -1,3 +1,6 @@
+import itertools
+import sys
+
 import numpy as np
 import astropy.units as u
 
@@ -10,6 +13,20 @@ class SpectralRegion:
     In the future, there might be more functionality added in here and there
     is some discussion that this might/could move to
     `Astropy Regions <http://astropy-regions.readthedocs.io/en/latest/>`_.
+
+    Parameters
+    ----------
+
+    lower: Scalar `~astropy.units.Quantity` with pixel or any valid ``spectral_axis`` unit
+       The lower bound of the region.
+
+    upper: Scalar `~astropy.units.Quantity` with pixel or any valid ``spectral_axis`` unit
+       The upper bound of the region.
+
+    Notes
+    -----
+    The subregions will be ordered based on the lower bound of each subregion.
+
     """
 
     @staticmethod
@@ -30,123 +47,254 @@ def from_center(center=None, width=None):
 
         return SpectralRegion(center - width, center + width)
 
-    def __init__(self, lower, upper):
+    def __init__(self, *args):
         """
         Lower and upper values for the interval.
+        """
 
-        Parameters
-        ----------
+        # Create instance variables
+        self._subregions = None
 
-        lower: Scalar `~astropy.units.Quantity` with pixel or any valid ``spectral_axis`` unit
-           The lower bound of the region.
+        #
+        # Set the values (using the setters for doing the proper checking)
+        #
+
+        if self._is_2_element(args):
+            self._subregions = [tuple(args)]
+        elif isinstance(args, (list, tuple)) and all([self._is_2_element(x) for x in args[0]]):
+            self._subregions = [tuple(x) for x in args[0]]
+        else:
+            raise ValueError('SpectralRegion input must be a 2-tuple or a list of 2-tuples.')
+
+        #
+        #  Check validity of the input sub regions.
+        #
+        if not self._valid():
+            raise ValueError('SpectralRegion 2-tuple lower extent must be less than upper extent.')
+
+        #
+        # The sub-regions are to always be ordered based on the lower bound.
+        #
+        self._reorder()
+
+    def _info(self):
+        """
+        Pretty print the sub-regions.
+        """
 
-        upper: Scalar `~astropy.units.Quantity` with pixel or any valid ``spectral_axis`` unit
-           The upper bound of the region.
+        toreturn = 'Spectral Region, {} sub-regions:\n'.format(len(self._subregions))
 
-        """
+        # Setup the subregion text.
+        subregion_text = []
+        for ii, subregion in enumerate(self._subregions):
+            subregion_text.append('  ({}, {})'.format(subregion[0], subregion[1]))
 
-        # Create instance variables
-        self._lower = None
-        self._upper = None
+        # Determine the length of the text boxes.
+        max_len = max(len(srt) for srt in subregion_text) + 1
+        ncols = 70 // max_len
 
-        # Set the values (using the setters for doing the proper checking)
-        self.lower = lower
-        self.upper = upper
+        # Add sub region info to the output text.
+        fmt = '{' + ':<{}'.format(max_len) + '}'
+        for ii, srt in enumerate(subregion_text):
+            toreturn += fmt.format(srt)
+            if ii % ncols == (ncols-1):
+                toreturn += '\n'
 
-    @property
-    def lower(self):
-        return self._lower
+        return toreturn
 
-    @lower.setter
-    def lower(self, value):
-        """
-        Lower bound for the interval.
+    def __str__(self):
+        return self._info()
 
-        Parameters
-        ----------
+    def __repr__(self):
+        return self._info()
 
-        value: Scalar `~astropy.units.Quantity` with pixel or any valid ``spectral_axis`` unit
-           The lower bound of the region.
+    def __add__(self, other):
+        """
+        Ability to add two SpectralRegion classes together.
         """
+        return SpectralRegion(self._subregions + other._subregions)
 
-        if not (value.unit.is_equivalent(u.pixel) or
-           value.unit.is_equivalent(u.angstrom, equivalencies=u.spectral())):
-            raise u.UnitError('Lower bound of region is not a spectral region unit')
+    def __iadd__(self, other):
+        """
+        Ability to add one SpectralRegion to another using +=.
+        """
+        self._subregions += other._subregions
+        self._reorder()
+        return self
 
-        self._lower = value
+    def __len__(self):
+        """
+        Number of spectral regions.
+        """
+        return len(self._subregions)
 
-    @property
-    def upper(self):
-        return self._upper
+    def __getslice__(self, item):
+        """
+        Enable slicing of the SpectralRegion list.
+        """
+        return SpectralRegion(self._subregions[item])
+
+    def __getitem__(self, item):
+        """
+        Enable slicing or extracting the SpectralRegion.
+        """
+        if isinstance(item, slice):
+            return self.__getslice__(item)
+        else:
+            return SpectralRegion([self._subregions[item]])
 
-    @upper.setter
-    def upper(self, value):
+    def __delitem__(self, item):
+        """
+        Delete a specific item from the list.
         """
-        Upper bound for the interval.
+        del self._subregions[item]
 
-        Parameters
-        ----------
+    def _valid(self):
 
-        value: Scalar `~astropy.units.Quantity` with pixel or any valid ``spectral_axis`` unit
-           The upper bound of the region.
+        # Lower bound < Upper bound for all sub regions
+        if any(x[0] >= x[1] for x in self._subregions):
+            raise ValueError('Lower bound must be strictly less than the upper bound')
+
+        return True
+
+    def _is_2_element(self, value):
+        """
+        Helper function to check a variable to see if it
+        is a 2-tuple.
         """
+        return len(value) == 2 and \
+               isinstance(value[0], u.Quantity) and \
+               isinstance(value[1], u.Quantity)
 
-        if not (value.unit.is_equivalent(u.pixel) or
-           value.unit.is_equivalent(u.angstrom, equivalencies=u.spectral())):
-            raise u.UnitError('Upper bound of region is not a spectral region unit')
+    def _reorder(self):
+        """
+        Re-order the  list based on lower bounds.
+        """
+        self._subregions.sort(key=lambda k: k[0])
 
-        self._upper = value
+    @property
+    def subregions(self):
+        return self._subregions
 
-    def to_pixel(self, spectrum):
+    @property
+    def bounds(self):
+        """
+        Compute the lower and upper extent of the SpectralRegion.
         """
-        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.
+        return (self.lower, self.upper)
 
-        Parameters
-        ----------
-        spectrum: `~specutils.spectra.spectrum1d.Spectrum1D`
-            The spectrum object from which the region will be extracted.
+    @property
+    def lower(self):
+        """
+        The most minimum value of the sub-regions.
 
-        Return
-        ------
-        left_index, right_index: int, int
-            Left and right indices defined by the lower and upper bounds.
+        The sub-regions are ordered based on the lower bound, so the
+        lower bound for this instance is the lower bound of the first
+        sub-region.
+        """
+        return self._subregions[0][0]
 
+    @property
+    def upper(self):
         """
+        The most maximum value of the sub-regions.
 
-        left_index = int(np.ceil(spectrum.wcs.world_to_pixel(np.array(self._lower))))
-        right_index = int(np.floor(spectrum.wcs.world_to_pixel(np.array(self._upper))))
+        The sub-regions are ordered based on the lower bound, but the
+        upper bound might not be the upper bound of the last sub-region
+        so we have to look for it.
+        """
+        return max(x[1] for x in self._subregions)
+
+    def invert_from_spectrum(self, spectrum):
+        """
+        Invert a SpectralRegion based on the extent of the
+        input spectrum.
+
+        See notes in SpectralRegion.invert() method.
+        """
+        return self.invert(spectrum.spectral_axis[0], spectrum.spectral_axis[-1])
 
-        return left_index, right_index
+    def _in_range(self, value, lower, upper):
+        return (value >= lower) and (value <= upper)
 
-    def extract(self, spectrum):
+    def invert(self, lower_bound, upper_bound):
         """
-        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.
+        Given a set of sub-regions this SpectralRegion defines,
+        create a new SpectralRegion such that the sub-regions
+        are defined in the new one as regions *not* in this
+        SpectralRegion.
 
         Parameters
         ----------
-        spectrum: `~specutils.spectra.spectrum1d.Spectrum1D`
-            The spectrum object from which the region will be extracted.
 
-        Return
-        ------
-        spectrum: `~specutils.spectra.spectrum1d.Spectrum1D`
-            Excised spectrum.
+        lower: Scalar `~astropy.units.Quantity` with pixel or any valid ``spectral_axis`` unit
+           The lower bound of the region.
+
+        upper: Scalar `~astropy.units.Quantity` with pixel or any valid ``spectral_axis`` unit
+           The upper bound of the region.
+
+        Returns
+        -------
+        spectral_region: `~specutils.SpectralRegion`
+           Spectral region of the non-selected regions
 
         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.
 
-        """
+        For example, if this SpectralRegion is defined as:
+           sr = SpectralRegion([(0.45*u.um, 0.6*u.um), (0.8*u.um, 0.9*u.um)])
 
-        left_index, right_index = self.to_pixel(spectrum)
+        and we call ``sr_invert = sr.invert(0.3*u.um, 1.0*u.um)`` then
+        ``sr_invert`` will be SpectralRegion([(0.3*u.um, 0.45*u.um), (0.6*u.um, 0.8*u.um), (0.9*u.um, 1*u.um)])
 
-        if left_index >= right_index:
-            raise ValueError('Lower region, {}, appears to be greater than the upper region, {}.'.format(self._lower, self._upper))
+        This could be useful if a SpectralRegion has sub-regions defined
+        for peaks in a spectrum and then one wants to create a SpectralRegion
+        defined as all the *non*-peaks, then one could use this function.
+
+        """
 
-        return spectrum[left_index:right_index]
+        #
+        # Create 'rs' region list with left and right extra ranges.
+        #
+        min_num = -sys.maxsize-1
+        max_num = sys.maxsize
+        rs = self._subregions + [(min_num*u.um, lower_bound), (upper_bound, max_num*u.um)]
+
+        #
+        # Sort the region list based on lower bound.
+        #
+
+        sorted_regions = sorted(rs, key=lambda k: k[0])
+
+        #
+        # Create new region list that has overlapping regions merged
+        #
+
+        merged = []
+        for higher in sorted_regions:
+            if not merged:
+                merged.append(higher)
+            else:
+                lower = merged[-1]
+                # test for intersection between lower and higher:
+                # we know via sorting that lower[0] <= higher[0]
+                if higher[0] <= lower[1]:
+                    upper_bound = max(lower[1], higher[1])
+                    merged[-1] = (lower[0], upper_bound)  # replace by merged interval
+                else:
+                    merged.append(higher)
+
+        #
+        # Create new list and drop first and last (the maxsize ones).
+        # We go from -inf, upper1, lower2, upper2....
+        # and remap to     lower1, upper1, lower2, ...
+        #
+
+        newlist = list(itertools.chain.from_iterable(merged))
+        newlist = newlist[1:-1]
+
+        #
+        # Now create new Spectrum region
+        #
+
+        return SpectralRegion([(x, y) for x, y in zip(newlist[0::2], newlist[1::2])])
diff --git a/specutils/tests/test_analysis.py b/specutils/tests/test_analysis.py
index fd638d80c..d01f09b9c 100644
--- a/specutils/tests/test_analysis.py
+++ b/specutils/tests/test_analysis.py
@@ -1,14 +1,11 @@
 import astropy.units as u
-import astropy.wcs as fitswcs
-import gwcs
 import numpy as np
 import pytest
 
 from ..tests.spectral_examples import simulated_spectra
 from ..spectra.spectrum1d import Spectrum1D
-import astropy.units as u
 from ..analysis import equivalent_width, snr, centroid
-from ..manipulation import noise_region_uncertainty
+from ..manipulation import noise_region_uncertainty, extract_region
 from ..spectra import SpectralRegion
 from astropy.nddata import StdDevUncertainty
 
@@ -128,8 +125,9 @@ def test_snr_single_region(simulated_spectra):
     wavelengths = spectrum.spectral_axis
     flux = spectrum.flux
 
+    # +1 because we want to include it in the calculation
     l = np.nonzero(wavelengths>region.lower)[0][0]
-    r = np.nonzero(wavelengths<region.upper)[0][-1]
+    r = np.nonzero(wavelengths<region.upper)[0][-1]+1
 
     spec_snr_expected = np.mean(flux[l:r] / (uncertainty.array[l:r]*uncertainty.unit))
 
@@ -170,7 +168,7 @@ def test_snr_two_regions(simulated_spectra):
     for region in regions:
 
         l = np.nonzero(wavelengths >= region.lower)[0][0]
-        r = np.nonzero(wavelengths <= region.upper)[0][-1]
+        r = np.nonzero(wavelengths <= region.upper)[0][-1]+1
 
         spec_snr_expected.append(np.mean(flux[l:r] / (uncertainty.array[l:r]*uncertainty.unit)))
 
@@ -183,43 +181,6 @@ def test_snr_two_regions(simulated_spectra):
     assert np.allclose(spec_snr, spec_snr_expected)
 
 
-def test_snr_single_region_with_noise_region(simulated_spectra):
-    """
-    Test the simple version of the spectral SNR over a region of the spectrum.
-    """
-
-    np.random.seed(42)
-
-    region = SpectralRegion(0.52*u.um, 0.59*u.um)
-    noise_region = SpectralRegion(0.40*u.um, 0.45*u.um)
-
-    #
-    #  Set up the data
-    #
-
-    spectrum = simulated_spectra.s1_um_mJy_e1
-    spectrum_uncertainty = noise_region_uncertainty(spectrum, noise_region)
-
-    wavelengths = spectrum.spectral_axis
-    flux = spectrum.flux
-
-    l = np.nonzero(wavelengths>=region.lower)[0][0]
-    r = np.nonzero(wavelengths<region.upper)[0][-1]
-
-    noise_l = np.nonzero(wavelengths>=noise_region.lower)[0][0]
-    noise_r = np.nonzero(wavelengths<noise_region.upper)[0][-1]
-
-    spec_snr_expected = np.mean(flux[l:r] / np.std(flux[noise_l:noise_r]))
-
-    #
-    # SNR of the whole spectrum
-    #
-
-    spec_snr = snr(spectrum_uncertainty, region)
-
-    assert np.allclose(spec_snr.value, spec_snr_expected.value)
-
-
 def test_centroid(simulated_spectra):
     """
     Test the simple version of the spectral centroid.
diff --git a/specutils/tests/test_region_extract.py b/specutils/tests/test_region_extract.py
new file mode 100644
index 000000000..ec3a95795
--- /dev/null
+++ b/specutils/tests/test_region_extract.py
@@ -0,0 +1,158 @@
+import numpy as np
+import pytest
+
+import astropy.units as u
+from astropy.nddata import StdDevUncertainty
+
+from ..spectra import Spectrum1D, SpectralRegion
+from ..manipulation import extract_region
+from .spectral_examples import simulated_spectra
+
+
+def test_region_simple(simulated_spectra):
+    """
+    Test the simple version of the spectral SNR.
+    """
+
+    np.random.seed(42)
+
+    #
+    #  Set up the data and add the uncertainty and calculate the expected SNR
+    #
+
+    spectrum = simulated_spectra.s1_um_mJy_e1
+    uncertainty = StdDevUncertainty(0.1*np.random.random(len(spectrum.flux))*u.mJy)
+    spectrum.uncertainty = uncertainty
+
+    region = SpectralRegion(0.6*u.um, 0.8*u.um)
+
+    sub_spectrum = extract_region(spectrum, region)
+
+    sub_spectrum_flux_expected = np.array(
+            [1605.71612173, 1651.41650744, 2057.65798618, 2066.73502361, 1955.75832537,
+             1670.52711471, 1491.10034446, 1637.08084112, 1471.28982259, 1299.19484483,
+             1423.11195734, 1226.74494917, 1572.31888312, 1311.50503403, 1474.05051673,
+             1335.39944397, 1420.61880528, 1433.18623759, 1290.26966668, 1605.67341284,
+             1528.52281708, 1592.74392861, 1568.74162534, 1435.29407808, 1536.68040935,
+             1157.33825995, 1136.12679394,  999.92394692, 1038.61546167, 1011.60297294])
+
+    assert np.allclose(sub_spectrum.flux.value, sub_spectrum_flux_expected)
+
+def test_region_ghz(simulated_spectra):
+    """
+    Test the simple version of the spectral SNR.
+    """
+
+    np.random.seed(42)
+
+    #
+    #  Set up the data and add the uncertainty and calculate the expected SNR
+    #
+
+    spectrum = simulated_spectra.s1_um_mJy_e1
+    spectrum = Spectrum1D(flux=simulated_spectra.s1_um_mJy_e1,
+                          spectral_axis=simulated_spectra.s1_um_mJy_e1.frequency)
+
+    region = SpectralRegion(374740.5725*u.GHz, 499654.09666667*u.GHz)
+
+    sub_spectrum = extract_region(spectrum, region)
+
+    sub_spectrum_flux_expected = np.array(
+            [1605.71612173, 1651.41650744, 2057.65798618, 2066.73502361, 1955.75832537,
+             1670.52711471, 1491.10034446, 1637.08084112, 1471.28982259, 1299.19484483,
+             1423.11195734, 1226.74494917, 1572.31888312, 1311.50503403, 1474.05051673,
+             1335.39944397, 1420.61880528, 1433.18623759, 1290.26966668, 1605.67341284,
+             1528.52281708, 1592.74392861, 1568.74162534, 1435.29407808, 1536.68040935,
+             1157.33825995, 1136.12679394,  999.92394692, 1038.61546167, 1011.60297294])
+
+    assert np.allclose(sub_spectrum.flux.value, sub_spectrum_flux_expected)
+
+def test_region_simple_check_ends(simulated_spectra):
+    """
+    Test the simple version of the spectral SNR.
+    """
+
+    np.random.seed(42)
+
+    #
+    #  Set up the data and add the uncertainty and calculate the expected SNR
+    #
+
+    spectrum = Spectrum1D(spectral_axis=np.linspace(1, 25, 25)*u.um, flux=np.random.random(25)*u.Jy)
+    region = SpectralRegion(8*u.um, 15*u.um)
+
+    sub_spectrum = extract_region(spectrum, region)
+
+    assert sub_spectrum.spectral_axis.value[0] == 8
+    assert sub_spectrum.spectral_axis.value[-1] == 15
+
+
+def test_region_empty(simulated_spectra):
+    """
+    Test the simple version of the spectral SNR.
+    """
+
+    np.random.seed(42)
+
+    # Region past upper range of spectrum
+    spectrum = Spectrum1D(spectral_axis=np.linspace(1, 25, 25)*u.um, flux=np.random.random(25)*u.Jy)
+    region = SpectralRegion(28*u.um, 30*u.um)
+    sub_spectrum = extract_region(spectrum, region)
+    assert sub_spectrum is None
+
+    # Region below lower range of spectrum
+    spectrum = Spectrum1D(spectral_axis=np.linspace(1, 25, 25)*u.um, flux=np.random.random(25)*u.Jy)
+    region = SpectralRegion(0.1*u.um, 0.3*u.um)
+    sub_spectrum = extract_region(spectrum, region)
+    assert sub_spectrum is None
+
+    # Region below lower range of spectrum and upper range in the spectrum.
+    spectrum = Spectrum1D(spectral_axis=np.linspace(1, 25, 25)*u.um, flux=np.random.random(25)*u.Jy)
+    region = SpectralRegion(0.1*u.um, 3.3*u.um)
+    sub_spectrum = extract_region(spectrum, region)
+    assert sub_spectrum is not None
+
+    # Region has lower and upper bound the same
+    with pytest.raises(Exception) as e_info:
+        region = SpectralRegion(3*u.um, 3*u.um)
+
+
+def test_region_two_sub(simulated_spectra):
+    """
+    Test the simple version of the spectral SNR.
+    """
+
+    np.random.seed(42)
+
+    #
+    #  Set up the data and add the uncertainty and calculate the expected SNR
+    #
+
+    spectrum = simulated_spectra.s1_um_mJy_e1
+    uncertainty = StdDevUncertainty(0.1*np.random.random(len(spectrum.flux))*u.mJy)
+    spectrum.uncertainty = uncertainty
+
+    region = SpectralRegion([(0.6*u.um, 0.8*u.um), (0.86*u.um, 0.89*u.um)])
+
+    sub_spectra = extract_region(spectrum, region)
+
+    # Confirm the end points of the subspectra are correct
+    assert np.allclose(sub_spectra[0].spectral_axis.value[[0, -1]],
+                       np.array([0.6035353535353536, 0.793939393939394]))
+
+    assert np.allclose(sub_spectra[1].spectral_axis.value[[0, -1]],
+                       np.array([0.8661616161616162, 0.8858585858585859]))
+
+    sub_spectrum_0_flux_expected = np.array(
+            [1605.71612173, 1651.41650744, 2057.65798618, 2066.73502361, 1955.75832537,
+             1670.52711471, 1491.10034446, 1637.08084112, 1471.28982259, 1299.19484483,
+             1423.11195734, 1226.74494917, 1572.31888312, 1311.50503403, 1474.05051673,
+             1335.39944397, 1420.61880528, 1433.18623759, 1290.26966668, 1605.67341284,
+             1528.52281708, 1592.74392861, 1568.74162534, 1435.29407808, 1536.68040935,
+             1157.33825995, 1136.12679394,  999.92394692, 1038.61546167, 1011.60297294])
+
+    sub_spectrum_1_flux_expected = np.array(
+            [1337.65312465, 1263.48914109, 1589.81797876, 1548.46068415])
+
+    assert np.allclose(sub_spectra[0].flux.value, sub_spectrum_0_flux_expected)
+    assert np.allclose(sub_spectra[1].flux.value, sub_spectrum_1_flux_expected)
diff --git a/specutils/tests/test_regions.py b/specutils/tests/test_regions.py
new file mode 100644
index 000000000..cea211471
--- /dev/null
+++ b/specutils/tests/test_regions.py
@@ -0,0 +1,136 @@
+import numpy as np
+
+import astropy.units as u
+from astropy.nddata import StdDevUncertainty
+
+from ..spectra import Spectrum1D, SpectralRegion
+
+
+def test_lower_upper():
+    # Spectral region with just one range (lower and upper bound)
+    sr = SpectralRegion(0.45*u.um, 0.6*u.um)
+
+    assert sr.lower == 0.45*u.um
+    assert sr.upper == 0.6*u.um
+
+    # Spectral region with just two ranges
+    sr = SpectralRegion([(0.45*u.um, 0.6*u.um), (0.8*u.um, 0.9*u.um)])
+
+    assert sr.lower == 0.45*u.um
+    assert sr.upper == 0.9*u.um
+
+    # Spectral region with multiple ranges and not ordered
+    sr = SpectralRegion([(0.3*u.um, 1.0*u.um), (0.45*u.um, 0.6*u.um), (0.04*u.um, 0.05*u.um),
+                         (0.8*u.um, 0.9*u.um)])
+
+    assert sr.lower == 0.04*u.um
+    assert sr.upper == 1.0*u.um
+
+    # Get lower bound of a single sub-region:
+    assert sr[0].lower == 0.04*u.um
+    assert sr[0].upper == 0.05*u.um
+
+
+def test_adding_spectral_regions():
+
+    # Combine two Spectral regions into one:
+    sr = SpectralRegion(0.45*u.um, 0.6*u.um) + SpectralRegion(0.8*u.um, 0.9*u.um)
+
+    assert set(sr.subregions) == set([(0.45*u.um, 0.6*u.um), (0.8*u.um, 0.9*u.um)])
+
+    # In-place adding spectral regions:
+    sr1 = SpectralRegion(0.45*u.um, 0.6*u.um)
+    sr2 = SpectralRegion(0.8*u.um, 0.9*u.um)
+    sr1 += sr2
+
+    assert set(sr1.subregions) == set([(0.45*u.um, 0.6*u.um), (0.8*u.um, 0.9*u.um)])
+
+
+def test_getitem():
+
+    sr = SpectralRegion([(0.8*u.um, 0.9*u.um), (0.3*u.um, 1.0*u.um),
+                         (0.45*u.um, 0.6*u.um), (0.04*u.um, 0.05*u.um)])
+
+    assert sr[0].subregions == [(0.04*u.um, 0.05*u.um)]
+    assert sr[1].subregions == [(0.3*u.um, 1.0*u.um)]
+    assert sr[2].subregions == [(0.45*u.um, 0.6*u.um)]
+    assert sr[3].subregions == [(0.8*u.um, 0.9*u.um)]
+    assert sr[-1].subregions == [(0.8*u.um, 0.9*u.um)]
+
+
+def test_bounds():
+
+    # Single subregion
+    sr = SpectralRegion(0.45*u.um, 0.6*u.um)
+    assert sr.bounds == (0.45*u.um, 0.6*u.um)
+
+    # Multiple subregions
+    sr = SpectralRegion([(0.8*u.um, 0.9*u.um), (0.3*u.um, 1.0*u.um),
+                         (0.45*u.um, 0.6*u.um), (0.04*u.um, 0.05*u.um)])
+    assert sr.bounds == (0.04*u.um, 1.0*u.um)
+
+
+def test_delitem():
+
+    # Single subregion
+    sr = SpectralRegion(0.45*u.um, 0.6*u.um)
+
+    del sr[0]
+    assert sr.subregions == []
+
+    # Multiple sub-regions
+    sr = SpectralRegion([(0.8*u.um, 0.9*u.um), (0.3*u.um, 1.0*u.um),
+                         (0.45*u.um, 0.6*u.um), (0.04*u.um, 0.05*u.um)])
+    del sr[1]
+
+    assert sr[0].subregions == [(0.04*u.um, 0.05*u.um)]
+    assert sr[1].subregions == [(0.45*u.um, 0.6*u.um)]
+    assert sr[2].subregions == [(0.8*u.um, 0.9*u.um)]
+
+
+def test_iterate():
+
+    # Create the Spectral region
+    subregions = [(0.8*u.um, 0.9*u.um), (0.3*u.um, 1.0*u.um),
+                  (0.45*u.um, 0.6*u.um), (0.04*u.um, 0.05*u.um)]
+    sr = SpectralRegion(subregions)
+
+    # For testing, sort our subregion list.
+    subregions.sort(key=lambda k: k[0])
+
+    for ii, s in enumerate(sr):
+        assert s.subregions[0] == subregions[ii]
+
+
+def test_slicing():
+
+    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)
+
+    subsr = sr[3:5]
+
+    assert subsr[0].subregions == [(0.8*u.um, 0.9*u.um)]
+    assert subsr[1].subregions == [(1.0*u.um, 1.2*u.um)]
+
+
+def test_invert():
+    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_expected = [(0.05*u.um, 0.15*u.um), (0.2*u.um, 0.3*u.um), (0.4*u.um, 0.45*u.um),
+                            (0.6*u.um, 0.8*u.um), (0.9*u.um, 1.0*u.um), (1.2*u.um, 1.3*u.um),
+                            (1.5*u.um, 3.0*u.um)]
+
+    # Invert from range.
+    sr_inverted = sr.invert(0.05*u.um, 3*u.um)
+
+    for ii, expected in enumerate(sr_inverted_expected):
+        assert sr_inverted.subregions[ii] == sr_inverted_expected[ii]
+
+    # Invert from spectrum.
+    spectrum = Spectrum1D(spectral_axis=np.linspace(0.05, 3, 20)*u.um, flux=np.random.random(20)*u.Jy)
+    sr_inverted = sr.invert_from_spectrum(spectrum)
+    for ii, expected in enumerate(sr_inverted_expected):
+        assert sr_inverted.subregions[ii] == sr_inverted_expected[ii]