From 9e2e034665a414b0b96c8f3a07c6846c0a759dfe Mon Sep 17 00:00:00 2001 From: Leonardo Uieda Date: Tue, 9 Apr 2024 17:33:12 -0300 Subject: [PATCH] Add an option to download L1 from Momotombo (#75) The L1 scene looks better and doesn't have the many atmospheric correction issues with the L2 one. --- doc/plot-overlay.rst | 23 +++++++++++------- xlandsat/datasets.py | 27 ++++++++++++++++----- xlandsat/tests/test_datasets.py | 43 ++++++++++++++++++++++++--------- 3 files changed, 66 insertions(+), 27 deletions(-) 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 `__, 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()