Why use xlandsat?#
One of the main features of xlandsat is the ability to easily read scenes
downloaded from USGS EarthExplorer into
xarray.Dataset along with all of the available metadata, which is very
useful for processing and plotting multidimensional array data.
We offer a simple and easy-to-use tool for smaller scale analysis and
visualization, which we can do thanks to the power of xarray.
When reading the data from the TIF files files using other tools like
rioxarray, some of
the rich metadata can be missing since it’s not always present in the TIF files
Things like conversion factors, units, data provenance, WRS path/row numbers,
We take care of fetching that information from the
*_MTL.txt files provided
by EarthExplorer so that xarray can use it, for example when annotating plots.
All functionality in xlandsat is available from the base namespace of the
xlandsat package. This means that you can access all of them with a
# xlandsat is usually imported as xls import xlandsat as xls
Download a sample scene#
xlandsat.datasets module includes functions for downloading some
sample scenes that we can use. These are cropped to smaller regions of interest
in order to save download time. But everything we do here with these sample
scenes is exactly the same as you would do with a full scene from
As an example, lets download a
.tar archive of a Landsat 9 scene of the
city of Manaus, in the Brazilian Amazon:
path_to_archive = xls.datasets.fetch_manaus() print(path_to_archive)
fetch_manaus function downloads the data file
and returns the path to the archive on your machine as an
The rest of this tutorial can be executed with your own data by changing the
path_to_archive to point to your data file instead.
The path can also point to a folder with the
.TIF and the
file instead of a
Running the code above will only download the data once. We use Pooch to handle the downloads and it’s smart enough to check if the file already exists on your computer.
Load the scene#
Now that we have the path to the tar archive of the scene, we can use
xlandsat.load_scene to read the bands and metadata directly from the
scene = xls.load_scene(path_to_archive) scene
<xarray.Dataset> Dimensions: (easting: 851, northing: 468) Coordinates: * easting (easting) float64 8.325e+05 8.325e+05 ... 8.58e+05 8.58e+05 * northing (northing) float64 -3.55e+05 -3.55e+05 ... -3.41e+05 -3.41e+05 Data variables: blue (northing, easting) float16 0.05908 0.05896 ... 0.06702 0.06665 green (northing, easting) float16 0.07446 0.0752 ... 0.09131 0.09082 red (northing, easting) float16 0.06152 0.06201 ... 0.1055 0.1046 nir (northing, easting) float16 0.2915 0.293 ... 0.06409 0.06445 swir1 (northing, easting) float16 0.1411 0.1434 ... 0.0509 0.05042 swir2 (northing, easting) float16 0.07996 0.08093 ... 0.04993 0.04968 Attributes: (12/19) Conventions: CF-1.8 title: Landsat 9 scene from 2023-07-23 (path/row=231... digital_object_identifier: https://doi.org/10.5066/P9OGBGM6 origin: Image courtesy of the U.S. Geological Survey landsat_product_id: LC09_L2SP_231062_20230723_20230802_02_T1 processing_level: L2SP ... ... ellipsoid: WGS84 date_acquired: 2023-07-23 scene_center_time: 14:12:31.2799050Z wrs_path: 231 wrs_row: 62 mtl_file: GROUP = LANDSAT_METADATA_FILE\n GROUP = PROD...
scene variable at the end of a code cell in a Jupyter
notebook will display a nice preview of the data. This is very useful for
looking up metadata and seeing which bands were loaded.
The scene is an
xarray.Dataset. It contains general metadata for the
scene and all of the bands available in the archive as
The bands each have their own set of metadata as well and can be accessed by
<xarray.DataArray 'nir' (northing: 468, easting: 851)> array([[0.2915 , 0.293 , 0.2954 , ..., 0.2905 , 0.2905 , 0.2852 ], [0.2744 , 0.2788 , 0.2876 , ..., 0.2905 , 0.2837 , 0.2866 ], [0.2793 , 0.2769 , 0.2632 , ..., 0.3145 , 0.3247 , 0.3306 ], ..., [0.1864 , 0.1781 , 0.1835 , ..., 0.063 , 0.0641 , 0.06494], [0.1761 , 0.1761 , 0.2042 , ..., 0.06323, 0.0635 , 0.0642 ], [0.1791 , 0.1818 , 0.2052 , ..., 0.0646 , 0.0641 , 0.06445]], dtype=float16) Coordinates: * easting (easting) float64 8.325e+05 8.325e+05 ... 8.58e+05 8.58e+05 * northing (northing) float64 -3.55e+05 -3.55e+05 ... -3.41e+05 -3.41e+05 Attributes: long_name: near-infrared units: reflectance number: 5 filename: LC09_L2SP_231062_20230723_20230802_02_T1_SR_B5.TIF scaling_mult: 2e-05 scaling_add: -0.1
Plot some reflectance bands#
Now we can use the
xarray.DataArray.plot method to make plots of
individual bands with
matplotlib. A bonus is that
xarray uses the
xlandsat.load_scene inserts into the scene to
automatically add labels and annotations to the plot:
import matplotlib.pyplot as plt band_names = list(scene.data_vars.keys()) fig, axes = plt.subplots( len(band_names), 1, figsize=(8, 16), layout="compressed", ) # Set the title using metadata from each scene fig.suptitle(scene.attrs["title"]) for band, ax in zip(band_names, axes.ravel()): # Make a pseudocolor plot of the band scene[band].plot(ax=ax) # Set the aspect to equal so that pixels are squares, not rectangles ax.set_aspect("equal") plt.show()
Checkout some of the other things that you can do with xlandsat:
Plus, by getting the data into an
xarray.Dataset, xlandsat opens the
door for a huge range of operations. You now have access to everything that
xarray can do: reduction, slicing, grouping, saving to cloud-optimized
formats, and much more. So go off and do something cool!