diff --git a/docs/creating_reading.rst b/docs/creating_reading.rst
index f1ef02a49..dbfecebc3 100644
--- a/docs/creating_reading.rst
+++ b/docs/creating_reading.rst
@@ -35,6 +35,60 @@ accessing other Stokes parameters, see :doc:`stokes`.
           files. See the :doc:`big_data` page for more details about dealing
           with large data sets.
 
+
+Reading images from file
+------------------------
+
+While spectral-cube is designed for cube analysis, in the course of normal
+analysis procedures you are likely to need to load up one- and two-dimensional
+subsets or views of the data.
+
+You can load :class:`~spectral_cube.lower_dimensional_objects.Projection`
+objects from 2D FITS HDU objects with
+:meth:`~spectral_cube.lower_dimensional_objects.Projection.from_hdu`.  Only
+FITS reading is currently supported::
+
+    >>> from astropy.io import fits  # doctest: +SKIP
+    >>> hdul = fits.open('file.fits')
+    >>> projection = Projection.from_hdu(hdul)
+
+Note that if you pass in a :class:`astropy.io.fits.HDUList` object, by default the data will be loaded
+from the first HDU. To load a different HDU in the list, the index can be passed to the ``ext`` keyword (e.g., `ext=1` to load the second HDU in the list).
+
+The resulting :class:`~spectral_cube.lower_dimensional_objects.Projection`
+object will have ``.unit``, ``.wcs``, and (if available) ``.beam`` attributes.
+
+If you are working with two dimensional data that have "dummy" third dimensions,
+you may load them using the normal :meth:`~spectral_cube.SpectralCube.read` method.
+This case is common as such files are normally output from CASA using
+``exportfits`` with no additional keywords.  To get a 2D slice, you simply index the
+result::
+
+    >>> flat_cube = SpectralCube.read('casa_exported_file.fits')  # doctest: +SKIP
+    >>> image = flat_cube[0]
+
+
+Reading spectra from file
+-------------------------
+
+Similar to 2D objects (images), you may want to load 1D slices - spectra - from disk.
+
+
+You can load :class:`~spectral_cube.lower_dimensional_objects.OneDSpectrum`
+objects from FITS HDU objects with
+:meth:`~spectral_cube.lower_dimensional_objects.OneDSpectrum.from_hdu`.  As with
+projections, only
+FITS reading is supported::
+
+    >>> from astropy.io import fits  # doctest: +SKIP
+    >>> hdul = fits.open('file.fits')
+    >>> projection = OneDSpectrum.from_hdu(hdul)
+
+The spectrum loader only works for 1D spectra with valid FITS WCS in their
+headers.  For other types of spectra, you may want to use `specutils
+<https://specutils.readthedocs.io/en/stable/>`_ instead.
+
+
 Direct Initialization
 ---------------------
 
diff --git a/spectral_cube/lower_dimensional_structures.py b/spectral_cube/lower_dimensional_structures.py
index 6c603fc17..f531cea1e 100644
--- a/spectral_cube/lower_dimensional_structures.py
+++ b/spectral_cube/lower_dimensional_structures.py
@@ -367,13 +367,20 @@ def _new_projection_with(self, data=None, wcs=None, mask=None, meta=None,
         return newproj
 
     @staticmethod
-    def from_hdu(hdu):
+    def from_hdu(hdu, ext=0):
         '''
         Return a projection from a FITS HDU.
+
+        Parameters
+        -----------
+        ext : int
+            The integer index to load when given an :class:`astropy.io.fits.HDUList`.
+            Default is 0 (the first HDU in the list.
+
         '''
         if isinstance(hdu, HDUList):
             hdul = hdu
-            hdu = hdul[0]
+            hdu = hdul[ext]
 
         if not len(hdu.data.shape) == 2:
             raise ValueError("HDU must contain two-dimensional data.")
@@ -625,14 +632,21 @@ def __repr__(self):
         return '{0}{1}{2:s}>'.format(prefixstr, arrstr, self._unitstr)
 
     @staticmethod
-    def from_hdu(hdu):
+    def from_hdu(hdu, ext=0):
         '''
         Return a OneDSpectrum from a FITS HDU or HDU list.
+
+        Parameters
+        -----------
+        ext : int
+            The integer index to load when given an :class:`astropy.io.fits.HDUList`.
+            Default is 0 (the first HDU in the list.
+
         '''
 
         if isinstance(hdu, HDUList):
             hdul = hdu
-            hdu = hdul[0]
+            hdu = hdul[ext]
         else:
             hdul = HDUList([hdu])