diff --git a/docs/conf.py b/docs/conf.py
index 0e2f1bc59..261036949 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -70,6 +70,7 @@
 """
 
 intersphinx_mapping['astroquery'] = ('http://astroquery.readthedocs.org/en/latest/', None)
+intersphinx_mapping['regions'] = ('http://astropy-regions.readthedocs.org/en/latest/', None)
 intersphinx_mapping['radio_beam'] = ('https://radio-beam.readthedocs.io/en/latest/', None)
 
 # -- Project information ------------------------------------------------------
diff --git a/docs/installing.rst b/docs/installing.rst
index ce81ae0da..e657d3066 100644
--- a/docs/installing.rst
+++ b/docs/installing.rst
@@ -13,6 +13,9 @@ This package has the following dependencies:
   reading in spectral cubes that use the BMAJ/BMIN convention for specifying the beam size.
 * `Bottleneck <http://berkeleyanalytics.com/bottleneck/>`_, optional (speeds
   up median and percentile operations on cubes with missing data)
+* `Regions <https://astropy-regions.readthedocs.io/en/latest>`_ >=0.3dev, optional
+  (Serialises/Deserialises DS9/CRTF region files and handles them. Used when
+  extracting a subcube from region)
 
 Installation
 ------------
diff --git a/docs/manipulating.rst b/docs/manipulating.rst
index c230de10b..412ef5a54 100644
--- a/docs/manipulating.rst
+++ b/docs/manipulating.rst
@@ -79,31 +79,69 @@ Numpy slicing notation::
 This returns a new :class:`~spectral_cube.SpectralCube` object
 with updated WCS information.
 
-Extracting a subcube from a ds9 region
---------------------------------------
+.. _reg:
 
-Starting with spectral_cube v0.2, you can use ds9 regions to extract subcubes.
-The minimal enclosing subcube will be extracted with a two-dimensional mask
-corresponding to the ds9 region.  `pyregion
-<http://leejjoon.github.io/pyregion/>`_ is required for region parsing::
+Extracting a subcube from a DS9/CRTF region
+-------------------------------------------
 
-    >>> region_list = pyregion.open('file.reg')  # doctest: +SKIP
-    >>> sub_cube = cube.subcube_from_ds9region(region_list)  # doctest: +SKIP
+You can use `DS9
+<http://ds9.si.edu/doc/ref/region.html>`_/`CRTF
+<https://casaguides.nrao.edu/index.php/CASA_Region_Format>`_ regions to extract
+subcubes. The minimal enclosing subcube will be extracted with a two-dimensional
+mask corresponding to the DS9/CRTF region.  `Regions
+<https://astropy-regions.readthedocs.io/en/latest/>`_ is required for region
+parsing.  CRTF regions may also contain spectral cutout information.
 
-If you want to loop over individual regions with a single region file, you need to convert the individual
-region to a shape list due to limitations in pyregion::
+This example shows extraction of a subcube from a ds9 region file ``file.reg``.
+`~regions.read_ds9` parses the ds9 file and converts it to a list of
+`~regions.Region` objects::
 
-  >>> region_list = pyregion.open('file.reg')  #doctest: +SKIP
-  >>> for region in region_list: #doctest: +SKIP
-  >>>     sub_cube = cube.subcube_from_ds9region(pyregion.ShapeList([region])) #doctest: +SKIP
-    
-You can also create a region on the fly using ds9 region syntax.  This extracts
-a 0.1 degree circle around the Galactic Center::
+    >>> import regions # doctest: +SKIP
+    >>> region_list = regions.read_ds9('file.reg')  # doctest: +SKIP
+    >>> sub_cube = cube.subcube_from_regions(region_list)  # doctest: +SKIP
+
+This one shows extraction of a subcube from a CRTF region file ``file.crtf``,
+parsed using `~regions.read_crtf`::
+
+    >>> import regions # doctest: +SKIP
+    >>> region_list = regions.read_crtf('file.reg')  # doctest: +SKIP
+    >>> sub_cube = cube.subcube_from_regions(region_list)  # doctest: +SKIP
 
-    >>> region_list = pyregion.parse("galactic; circle(0,0,0.1)")  # doctest: +SKIP
-    >>> sub_cube = cube.subcube_from_ds9region(region_list)  # doctest: +SKIP
+If you want to loop over individual regions with a single region file, you need
+to convert the individual regions to lists of that region::
 
+    >>> region_list = regions.read_ds9('file.reg')  #doctest: +SKIP
+    >>> for region in region_list: #doctest: +SKIP
+    >>>     sub_cube = cube.subcube_from_regions([region]) #doctest: +SKIP
     
+You can also directly use a ds9 region string.  This example extracts a 0.1
+degree circle around the Galactic Center::
+
+    >>> region_str = "galactic; circle(0, 0, 0.1)"  # doctest: +SKIP
+    >>> sub_cube = cube.subcube_from_ds9region(region_str)  # doctest: +SKIP
+
+Similarly, you can also use a CRTF region string::
+
+    >>> region_str = "circle[[0deg, 0deg], 0.1deg], coord=galactic, range=[150km/s, 300km/s]"  # doctest: +SKIP
+    >>> sub_cube = cube.subcube_from_crtfregion(region_str)  # doctest: +SKIP
+
+CRTF regions that specify a subset in the spectral dimension can be used to
+produce full 3D cutouts.  The ``meta`` attribute of a `regions.Region` object
+contains the spectral information for that region in the three special keywords
+``range``, ``restfreq``, and ``veltype``::
+
+    >>> import regions # doctest: +SKIP
+    >>> from astropy import units as u
+
+    >>> regpix = regions.RectanglePixelRegion(regions.PixCoord(0.5, 1), width=4, height=2)  # doctest: +SKIP
+    >>> regpix.meta['range'] = [150 * u.km/u.s, 300 * u.km/u.s] # spectral range # doctest: +SKIP
+    >>> regpix.meta['restfreq'] = [100 * u.GHz] # rest frequency # doctest: +SKIP
+    >>> regpix.meta['veltype'] = 'OPTICAL' # velocity convention # doctest: +SKIP
+    >>> subcube = cube.subcube_from_regions([regpix])  # doctest: +SKIP
+
+If ``range`` is specified, but the other two keywords are not, the code will
+likely crash.
+
 Extract the minimal valid subcube
 ---------------------------------
 
diff --git a/docs/spectral_extraction.rst b/docs/spectral_extraction.rst
index a7e9d3d55..d9560e75a 100644
--- a/docs/spectral_extraction.rst
+++ b/docs/spectral_extraction.rst
@@ -40,22 +40,34 @@ mask from scratch and apply it to the data.::
     >>> maskedsubcube = subcube.with_mask(mask)  # doctest: +SKIP
     >>> spectrum = maskedsubcube.mean(axis=(1,2))  # doctest: +SKIP
 
-Aperture extraction using regions
----------------------------------
+Aperture and spectral extraction using regions
+----------------------------------------------
 
-Spectral-cube supports ds9 regions, so you can use the ds9 region to create a
-mask.  The ds9 region support relies on `pyregion
-<https://pyregion.readthedocs.io/en/latest/>`_, which supports most shapes in
-ds9, so you are not limited to circular apertures.
+Spectral-cube supports ds9 and crtf regions, so you can use them to create a
+mask.  The ds9/crtf region support relies on `regions
+<https://astropy-regions.readthedocs.io/en/latest/>`, which supports most
+shapes in ds9 and crtf, so you are not limited to circular apertures.
 
-In this example, we'll create a region "from scratch", but you can also use a
-predefined region file using `pyregion.open
-<http://pyregion.readthedocs.io/en/latest/api/pyregion.open.html>`_.::
+In this example, we'll extract a subcube from ds9 region string using
+`~spectral_cube.spectral_cube.SpectralCube.subcube_from_ds9region`::
 
-    >>> shapelist = pyregion.parse("fk5; circle(19:23:43.907,+14:30:34.66, 3\")")  # doctest: +SKIP
-    >>> subcube = cube.subcube_from_ds9region(shapelist)  # doctest: +SKIP
-    >>> spectrum = subcube.mean(axis=(1,2))  # doctest: +SKIP
+    >>> ds9_str = 'fk5; circle(19:23:43.907, +14:30:34.66, 3")'  # doctest: +SKIP
+    >>> subcube = cube.subcube_from_ds9region(ds9_str)  # doctest: +SKIP
+    >>> spectrum = subcube.mean(axis=(1, 2))  # doctest: +SKIP
+
+Similarly, one can extract a subcube from a crtf region string using
+`~spectral_cube.spectral_cube.SpectralCube.subcube_from_crtfregion`::
+
+    >>> crtf_str = 'circle[[19:23:43.907, +14:30:34.66], 3"], coord=fk5, range=[150km/s, 300km/s]'  # doctest: +SKIP
+    >>> subcube = cube.subcube_from_crtfregion(crtf_str)  # doctest: +SKIP
+    >>> spectrum = subcube.mean(axis=(1, 2))  # doctest: +SKIP
+
+You can also use a _list_ of `~regions.Region` objects to extract a subcube using
+`~spectral_cube.spectral_cube.SpectralCube.subcube_from_regions`::
+
+    >>> import regions # doctest: +SKIP
+    >>> regpix = regions.RectanglePixelRegion(regions.PixCoord(0.5, 1), width=4, height=2)  # doctest: +SKIP
+    >>> subcube = cube.subcube_from_regions([regpix])  # doctest: +SKIP
+    >>> spectrum = subcube.mean(axis=(1, 2))  # doctest: +SKIP
 
-Eventually, we hope to change the region support from pyregion to `astropy
-regions <http://astropy-regions.readthedocs.io/en/latest/>`_, so the
-above example may become obsolete.
+To learn more, go to :ref:`reg`.
diff --git a/spectral_cube/spectral_cube.py b/spectral_cube/spectral_cube.py
index 8aced3c5f..55bea3a05 100644
--- a/spectral_cube/spectral_cube.py
+++ b/spectral_cube/spectral_cube.py
@@ -1842,7 +1842,7 @@ def subcube_from_ds9region(self, ds9_region, allow_empty=False):
 
         return self.subcube_from_regions(region_list, allow_empty)
 
-    def subcube_from_crtf(self, crtf_region, allow_empty=False):
+    def subcube_from_crtfregion(self, crtf_region, allow_empty=False):
         """
         Extract a masked subcube from a CRTF region.
 
diff --git a/spectral_cube/tests/test_subcubes.py b/spectral_cube/tests/test_subcubes.py
index c58eae82d..a9f36e2b0 100644
--- a/spectral_cube/tests/test_subcubes.py
+++ b/spectral_cube/tests/test_subcubes.py
@@ -1,6 +1,7 @@
 from __future__ import print_function, absolute_import, division
 
 import pytest
+from distutils.version import LooseVersion
 
 from astropy import units as u
 from astropy import wcs
@@ -14,8 +15,9 @@
 try:
     import regions
     regionsOK = True
+    REGIONS_GT_03 = LooseVersion(regions.__version__) >= LooseVersion('0.3')
 except ImportError:
-    regionsOK = False
+    regionsOK = REGIONS_GT_03 = False
 
 try:
     import scipy
@@ -55,6 +57,7 @@ def test_subcube():
 
 @pytest.mark.skipif('not scipyOK', reason='Could not import scipy')
 @pytest.mark.skipif('not regionsOK', reason='Could not import regions')
+@pytest.mark.skipif('not REGIONS_GT_03', reason='regions version should be >= 0.3')
 @pytest.mark.parametrize('regfile',
                          ('255-fk5.reg', '255-pixel.reg'),
                         )
@@ -71,6 +74,7 @@ def test_ds9region_255(regfile):
 
 @pytest.mark.skipif('not scipyOK', reason='Could not import scipy')
 @pytest.mark.skipif('not regionsOK', reason='Could not import regions')
+@pytest.mark.skipif('not REGIONS_GT_03', reason='regions version should be >= 0.3')
 @pytest.mark.parametrize(('regfile', 'result'),
                              (('fk5.reg', [slice(None), 1, 1]),
                               ('image.reg', [slice(None), 1, slice(None)]),
@@ -108,6 +112,7 @@ def test_ds9region_new(regfile, result):
 
 @pytest.mark.skipif('not scipyOK', reason='Could not import scipy')
 @pytest.mark.skipif('not regionsOK', reason='Could not import regions')
+@pytest.mark.skipif('not REGIONS_GT_03', reason='regions version should be >= 0.3')
 def test_regions_spectral():
     cube, data = cube_and_raw('adv.fits')
     rf_cube = get_rest_value_from_wcs(cube.wcs).to("GHz",