diff --git a/doc/plot-overlay.rst b/doc/plot-overlay.rst
index a88d2d6..f16b91c 100644
--- a/doc/plot-overlay.rst
+++ b/doc/plot-overlay.rst
@@ -11,7 +11,7 @@ volcano <https://en.wikipedia.org/wiki/Momotombo>`__, Nicaragua.
 We'll overlay the thermal band (only pixels above 320 K) on top of an RGB
 composite to show the ongoing lava flow.
 
-First, we'll import the required packages and load the sample scene:
+First, we'll import the required packages:
 
 .. jupyter-execute::
 
@@ -20,10 +20,18 @@ First, we'll import the required packages and load the sample scene:
     import xarray as xr
     import numpy as np
 
-    path = xls.datasets.fetch_momotombo()
-    scene = xls.load_scene(path)
-    # Fill the missing values due to the volcanic clouds to make it look nicer
-    scene = xls.interpolate_missing(scene)
+Now we can load a Level 1 version of the scene to make the RGB plot (the L2
+scene has a lot of issues from the atmospheric correction) and a Level 2
+version to get the thermal band as surface temperature:
+
+.. jupyter-execute::
+
+    path_l1 = xls.datasets.fetch_momotombo(level=1)
+    scene = xls.load_scene(path_l1)
+
+    path_l2 = xls.datasets.fetch_momotombo(level=2)
+    scene = scene.merge(xls.load_scene(path_l2, bands=["thermal"]))
+
     scene
 
 Now we can plot an RGB composite and thermal band separately to see that they
@@ -32,10 +40,7 @@ have to show:
 .. jupyter-execute::
 
     # Make the composite
-    rgb = xls.composite(scene, rescale_to=(0, 0.6))
-
-    # Histogram equalization for a better looking image
-    rgb = xls.equalize_histogram(rgb, clip_limit=0.02, kernel_size=200)
+    rgb = xls.composite(scene, rescale_to=(0, 0.2))
 
     # Plot the RGB and thermal separately
     fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 12))
diff --git a/xlandsat/datasets.py b/xlandsat/datasets.py
index fd638ee..66e1ba5 100644
--- a/xlandsat/datasets.py
+++ b/xlandsat/datasets.py
@@ -24,8 +24,10 @@
         "LC08_L2SP_204023_20200927_20201006_02_T1-cropped.tar.gz": "md5:3c07e343ccf959be4e5dd5c9aca4e0a4",
         # Liverpool - Panchromatic
         "LC08_L1TP_204023_20200927_20201006_02_T1-cropped.tar.gz": "md5:7d43f8580b8e583d137a93f9ae51a73d",
-        # Momotombo
+        # Momotombo L2
         "LC08_L2SP_017051_20151205_20200908_02_T1-cropped.tar.gz": "md5:8cc2e4c15e65940a7152fc1c8b412aa9",
+        # Momotombo L1
+        "LC08_L1TP_017051_20151205_20200908_02_T1-cropped.tar.gz": "md5:112d42e7adf709ac3a1179bbeedded6d",
         # Roraima
         "LC08_L2SP_232056_20151004_20200908_02_T1-cropped.tar.gz": "md5:f236a8b024aa4a4c62bee294d3bd737f",
         # Manaus
@@ -187,7 +189,7 @@ def fetch_liverpool_panchromatic(untar=False):
     )
 
 
-def fetch_momotombo(untar=False):
+def fetch_momotombo(untar=False, level=2):
     """
     Download a sample scene from the December 2015 Momotombo volcano eruption
 
@@ -209,16 +211,29 @@ def fetch_momotombo(untar=False):
     untar : bool
         If True, unpack the tar archive after downloading and return a path to
         the folder containing the unpacked files instead. Default is False.
+    level : int
+        Which level of data to load. Default is to load Level 2 surface
+        reflectance data. Can also be Level 1 top-of-the-atmosphere
+        reflectance.
 
     Returns
     -------
     path : str
         The path to the downloaded `.tar` file that contains the scene.
     """
-    return _fetch(
-        "LC08_L2SP_017051_20151205_20200908_02_T1-cropped.tar.gz",
-        untar,
-    )
+    if level not in {1, 2}:
+        raise ValueError(f"Invalid data level '{level}'. Must be 1 or 2.")
+    if level == 2:
+        data = _fetch(
+            "LC08_L2SP_017051_20151205_20200908_02_T1-cropped.tar.gz",
+            untar,
+        )
+    else:
+        data = _fetch(
+            "LC08_L1TP_017051_20151205_20200908_02_T1-cropped.tar.gz",
+            untar,
+        )
+    return data
 
 
 def fetch_roraima(untar=False):
diff --git a/xlandsat/tests/test_datasets.py b/xlandsat/tests/test_datasets.py
index 3aa386f..6d1581d 100644
--- a/xlandsat/tests/test_datasets.py
+++ b/xlandsat/tests/test_datasets.py
@@ -20,21 +20,40 @@
 
 
 @pytest.mark.parametrize("untar", [False, True], ids=["archive", "folder"])
-def test_fetching_functions(untar):
-    "Check that the download functions work"
-    functions = [
+@pytest.mark.parametrize(
+    "fetcher",
+    [
         fetch_brumadinho_after,
         fetch_brumadinho_before,
         fetch_liverpool,
         fetch_liverpool_panchromatic,
         fetch_manaus,
-        fetch_momotombo,
         fetch_roraima,
-    ]
-    for func in functions:
-        path = pathlib.Path(func(untar=untar))
-        assert path.exists()
-        if untar:
-            assert path.is_dir()
-        else:
-            assert not path.is_dir()
+    ],
+)
+def test_fetching_functions(untar, fetcher):
+    "Check that the download functions work"
+    path = pathlib.Path(fetcher(untar=untar))
+    assert path.exists()
+    if untar:
+        assert path.is_dir()
+    else:
+        assert not path.is_dir()
+
+
+@pytest.mark.parametrize("untar", [False, True], ids=["archive", "folder"])
+@pytest.mark.parametrize("level", [1, 2], ids=["L1", "L2"])
+@pytest.mark.parametrize(
+    "fetcher",
+    [
+        fetch_momotombo,
+    ],
+)
+def test_fetching_functions_levels(untar, level, fetcher):
+    "Check that the download functions work"
+    path = pathlib.Path(fetcher(untar=untar, level=level))
+    assert path.exists()
+    if untar:
+        assert path.is_dir()
+    else:
+        assert not path.is_dir()