diff --git a/astrodendro/__init__.py b/astrodendro/__init__.py
index 37934fd..d88039b 100644
--- a/astrodendro/__init__.py
+++ b/astrodendro/__init__.py
@@ -3,4 +3,5 @@
 from .dendrogram import Dendrogram
 from .analysis import ppv_catalog, pp_catalog
 from .plot import DendrogramPlotter
-from .viewer import BasicDendrogramViewer
\ No newline at end of file
+from .viewer import BasicDendrogramViewer
+from .datasets import load_perseus
diff --git a/docs/PerA_Extn2MASS_F_Gal.fits b/astrodendro/datasets/PerA_Extn2MASS_F_Gal.fits
similarity index 100%
rename from docs/PerA_Extn2MASS_F_Gal.fits
rename to astrodendro/datasets/PerA_Extn2MASS_F_Gal.fits
diff --git a/astrodendro/datasets/__init__.py b/astrodendro/datasets/__init__.py
new file mode 100644
index 0000000..614b1f6
--- /dev/null
+++ b/astrodendro/datasets/__init__.py
@@ -0,0 +1,14 @@
+import os
+
+from astropy.io import fits
+
+path = os.path.dirname(__file__)
+
+def load_perseus():
+    """ Load the Primary HDU for the Perseus extinction map
+    used throughout the documentation.
+
+    For more details about the data, see http://hdl.handle.net/10904/10080
+    """
+    pth = os.path.join(path, 'PerA_Extn2MASS_F_Gal.fits')
+    return fits.open(pth)[0]
diff --git a/astrodendro/test/test_datasets.py b/astrodendro/test/test_datasets.py
new file mode 100644
index 0000000..2963b2f
--- /dev/null
+++ b/astrodendro/test/test_datasets.py
@@ -0,0 +1,5 @@
+from astrodendro import load_perseus
+
+def test_load_perseus():
+    p = load_perseus()
+    assert p.data.shape == (313, 217)
diff --git a/docs/catalog.rst b/docs/catalog.rst
index 8aa86fa..d45f03b 100644
--- a/docs/catalog.rst
+++ b/docs/catalog.rst
@@ -22,9 +22,8 @@ where ``structure`` is a :class:`~astrodendro.structure.Structure` instance
 from a dendrogram. The resulting object then has methods to compute various
 statistics. Using the example data from :doc:`using`::
 
-    >>> from astrodendro import Dendrogram
-    >>> from astropy.io import fits
-    >>> image = fits.getdata('PerA_Extn2MASS_F_Gal.fits')
+    >>> from astrodendro import Dendrogram, load_perseus
+    >>> image = load_perseus().data
     >>> d = Dendrogram.compute(image, min_value=2.0, min_delta=1., min_npix=10)
 
 we can get statistics for the first structure in the trunk, which is a leaf::
@@ -200,13 +199,13 @@ approximating the structures on top of the structures themselves:
 
     from astropy.io import fits
 
-    from astrodendro import Dendrogram
+    from astrodendro import Dendrogram, load_perseus
     from astrodendro.analysis import PPStatistic
 
     import matplotlib.pyplot as plt
     from matplotlib.patches import Ellipse
 
-    hdu = fits.open('PerA_Extn2MASS_F_Gal.fits')[0]
+    hdu = load_perseus()
 
     d = Dendrogram.compute(hdu.data, min_value=2.0, min_delta=1., min_npix=10)
     p = d.plotter()
diff --git a/docs/plotting.rst b/docs/plotting.rst
index 1928947..a71d5ee 100644
--- a/docs/plotting.rst
+++ b/docs/plotting.rst
@@ -63,10 +63,9 @@ main branches:
    :include-source:
 
     import matplotlib.pyplot as plt
-    from astropy.io import fits
-    from astrodendro import Dendrogram
+    from astrodendro import Dendrogram, load_perseus
 
-    image = fits.getdata('PerA_Extn2MASS_F_Gal.fits')
+    image = load_perseus().data
     d = Dendrogram.compute(image, min_value=2.0, min_delta=1., min_npix=10)
     p = d.plotter()
 
@@ -98,10 +97,9 @@ shown.
    :include-source:
 
     import matplotlib.pyplot as plt
-    from astropy.io import fits
-    from astrodendro import Dendrogram
+    from astrodendro import Dendrogram, load_perseus
 
-    image = fits.getdata('PerA_Extn2MASS_F_Gal.fits')
+    image = load_perseus().data
     d = Dendrogram.compute(image, min_value=2.0, min_delta=1., min_npix=10)
     p = d.plotter()
 
diff --git a/docs/using.rst b/docs/using.rst
index f452e2c..608e891 100644
--- a/docs/using.rst
+++ b/docs/using.rst
@@ -58,15 +58,15 @@ These options are illustrated graphically in :doc:`algorithm`.
 As an example, we can use a publicly available extinction map of the Perseus
 star-formation region from the The COordinated Molecular Probe Line Extinction
 Thermal Emission (COMPLETE) Survey of Star Forming Regions
-(:download:`PerA_Extn2MASS_F_Gal.fits`, originally obtained from
+(originally obtained from
 `<http://hdl.handle.net/10904/10080>`_). The units of the map are magnitudes of
 extinction, and we want to make a dendrogram of all structures above a minimum
 value of 2 magnitudes, and we only consider leaves with at least 10 pixels and
 which have a peak to base different larger than one magnitude of extinction::
 
-    >>> from astrodendro import Dendrogram
+    >>> from astrodendro import Dendrogram, load_perseus
     >>> from astropy.io import fits
-    >>> image = fits.getdata('PerA_Extn2MASS_F_Gal.fits')
+    >>> image = load_perseus().data
     >>> d = Dendrogram.compute(image, min_value=2.0, min_delta=1., min_npix=10)
 
 By default, the computation will be silent, but for large dendrograms, it can
diff --git a/setup.py b/setup.py
index 6023531..ee2bb9b 100755
--- a/setup.py
+++ b/setup.py
@@ -49,8 +49,10 @@ def run(self):
       description='Python package for computation of astronomical dendrograms',
       author='Thomas Robitaille, Chris Beaumont, Braden MacDonald, and Erik Rosolowsky',
       author_email='braden@bradenmacdonald.com',
-      packages=['astrodendro', 'astrodendro.io', 'astrodendro.test'],
-      package_data={'astrodendro.test':['*.npz', 'benchmark_data/*fits']},
+      packages=['astrodendro', 'astrodendro.io', 'astrodendro.test',
+                'astrodendro.datasets'],
+      package_data={'astrodendro.test':['*.npz', 'benchmark_data/*fits'],
+                    'astrodendro.datasets': ['*fits']},
       provides=['astrodendro'],
       requires=['numpy'],
       cmdclass={'build_py': build_py, 'test': DendroTest},