Overview
Contents
Overview#
The library#
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
single import:
# xlandsat is usually imported as xls
import xlandsat as xls
Download a sample scene#
As an example, lets download a tar archive of a Landsat 8 scene of the Brumadinho tailings dam disaster that happened in 2019 in Brazil. The archive is available on figshare at https://doi.org/10.6084/m9.figshare.21665630.v1 and includes scenes from before and after the disaster as both the full scene and a cropped version.
We’ll use pooch
to download the scenes from before and after the
disaster to our computer.
To save space and bandwidth, we’ll use the cropped version here.
import pooch
path_before = pooch.retrieve(
"doi:10.6084/m9.figshare.21665630.v1/cropped-before.tar.gz",
known_hash="md5:d2a503c944bb7ef3b41294d44b77e98c",
)
path_after = pooch.retrieve(
"doi:10.6084/m9.figshare.21665630.v1/cropped-after.tar.gz",
known_hash="md5:4ae61a2d7a8b853c727c0c433680cece",
)
print(path_before)
print(path_after)
/home/runner/.cache/pooch/cc3e48c8ee8ca9c9ccb64e2a0e798286-cropped-before.tar.gz
/home/runner/.cache/pooch/7450bb0a799deea586280f7c3b539898-cropped-after.tar.gz
Tip
Running the code above will only download the data once since Pooch is
smart enough to check if the file already exists on your computer.
See pooch.retrieve
for more information.
See also
If you want to use the full scene instead of the cropped version,
substitute the file name to full-after.tar.gz
and
full-before.tar.gz
. Don’t forget to also update the MD5 hashes
accordingly, which you can find on the figshare archive.
Load the scenes#
Now that we have paths to the tar archives of the scenes, we can use
xlandsat.load_scene
to read the bands and metadata directly from the
archives (without unpacking):
before = xls.load_scene(path_before)
before
<xarray.Dataset> Dimensions: (easting: 400, northing: 300) Coordinates: * easting (easting) float64 5.835e+05 5.835e+05 ... 5.954e+05 5.955e+05 * northing (northing) float64 -2.232e+06 -2.232e+06 ... -2.223e+06 -2.223e+06 Data variables: blue (northing, easting) float16 0.07288 0.07373 ... 0.06506 0.06653 green (northing, easting) float16 0.09851 0.099 ... 0.0835 0.08716 red (northing, easting) float16 0.1035 0.1041 ... 0.07959 0.08423 nir (northing, easting) float16 0.2803 0.2749 0.3203 ... 0.2118 0.2267 swir1 (northing, easting) float16 0.2467 0.2474 0.2147 ... 0.172 0.1769 swir2 (northing, easting) float16 0.1571 0.1571 0.1281 ... 0.1171 0.1206 Attributes: (12/62) Conventions: CF-1.8 title: Landsat 8 scene from 2019-01-14 (path/ro... origin: Image courtesy of the U.S. Geological Su... digital_object_identifier: https://doi.org/10.5066/P9OGBGM6 landsat_product_id: LC08_L2SP_218074_20190114_20200829_02_T1 processing_level: L2SP ... ... temperature_mult_band_st_b10: 0.00341802 temperature_add_band_st_b10: 149.0 reflectance_mult_band_8: 2e-05 reflectance_mult_band_9: 2e-05 reflectance_add_band_8: -0.1 reflectance_add_band_9: -0.1
And the after scene:
after = xls.load_scene(path_after)
after
<xarray.Dataset> Dimensions: (easting: 400, northing: 300) Coordinates: * easting (easting) float64 5.844e+05 5.844e+05 ... 5.963e+05 5.964e+05 * northing (northing) float64 -2.232e+06 -2.232e+06 ... -2.223e+06 -2.223e+06 Data variables: blue (northing, easting) float16 0.0686 0.07043 ... 0.05823 0.0564 green (northing, easting) float16 0.1027 0.09839 ... 0.07593 0.07043 red (northing, easting) float16 0.09778 0.09778 ... 0.06799 0.06177 nir (northing, easting) float16 0.2988 0.2715 0.2881 ... 0.2637 0.251 swir1 (northing, easting) float16 0.2311 0.2274 0.2316 ... 0.1608 0.142 swir2 (northing, easting) float16 0.145 0.1442 0.144 ... 0.09961 0.08655 Attributes: (12/62) Conventions: CF-1.8 title: Landsat 8 scene from 2019-01-30 (path/ro... origin: Image courtesy of the U.S. Geological Su... digital_object_identifier: https://doi.org/10.5066/P9OGBGM6 landsat_product_id: LC08_L2SP_218074_20190130_20200829_02_T1 processing_level: L2SP ... ... temperature_mult_band_st_b10: 0.00341802 temperature_add_band_st_b10: 149.0 reflectance_mult_band_8: 2e-05 reflectance_mult_band_9: 2e-05 reflectance_add_band_8: -0.1 reflectance_add_band_9: -0.1
Did you notice?
If you look carefully at the coordinates for each scene, you may notice
that they don’t exactly coincide in area. That’s OK since xarray
knows how to take the pixel coordinates into account when doing
mathematical operations like calculating indices and differences between
scenes.
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
metadata that xlandsat.load_scene
inserts into the scene to
automatically add labels and annotations to the plot.
import matplotlib.pyplot as plt
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 12))
# Make the pseudocolor plots of the near infrared band
before.nir.plot(ax=ax1)
after.nir.plot(ax=ax2)
# Set the title using metadata from each scene
ax1.set_title(f"Before: {before.attrs['title']}")
ax2.set_title(f"After: {after.attrs['title']}")
# Set the aspect to equal so that pixels are squares, not rectangles
ax1.set_aspect("equal")
ax2.set_aspect("equal")
plt.show()
![_images/overview_4_0.png](_images/overview_4_0.png)
Make composites#
Plotting individual bands is good but we usually want to make some composite
images to visualize information from multiple bands at once.
Let’s make both RGB (true color) and CIR (color infrared) composites for both
of our scenes using xlandsat.composite
:
# Make the composite and add it as a variable to the scene
before = before.assign(rgb=xls.composite(before, rescale_to=[0, 0.2]))
cir_bands = ("nir", "red", "green")
before = before.assign(
cir=xls.composite(before, bands=cir_bands, rescale_to=[0, 0.4]),
)
before
<xarray.Dataset> Dimensions: (easting: 400, northing: 300, channel: 4) Coordinates: * easting (easting) float64 5.835e+05 5.835e+05 ... 5.954e+05 5.955e+05 * northing (northing) float64 -2.232e+06 -2.232e+06 ... -2.223e+06 -2.223e+06 * channel (channel) <U5 'red' 'green' 'blue' 'alpha' Data variables: blue (northing, easting) float16 0.07288 0.07373 ... 0.06506 0.06653 green (northing, easting) float16 0.09851 0.099 ... 0.0835 0.08716 red (northing, easting) float16 0.1035 0.1041 ... 0.07959 0.08423 nir (northing, easting) float16 0.2803 0.2749 0.3203 ... 0.2118 0.2267 swir1 (northing, easting) float16 0.2467 0.2474 0.2147 ... 0.172 0.1769 swir2 (northing, easting) float16 0.1571 0.1571 0.1281 ... 0.1171 0.1206 rgb (northing, easting, channel) uint8 132 125 92 255 ... 111 84 255 cir (northing, easting, channel) uint8 178 66 62 255 ... 144 53 55 255 Attributes: (12/62) Conventions: CF-1.8 title: Landsat 8 scene from 2019-01-14 (path/ro... origin: Image courtesy of the U.S. Geological Su... digital_object_identifier: https://doi.org/10.5066/P9OGBGM6 landsat_product_id: LC08_L2SP_218074_20190114_20200829_02_T1 processing_level: L2SP ... ... temperature_mult_band_st_b10: 0.00341802 temperature_add_band_st_b10: 149.0 reflectance_mult_band_8: 2e-05 reflectance_mult_band_9: 2e-05 reflectance_add_band_8: -0.1 reflectance_add_band_9: -0.1
The composites have a similar layout as the bands but with an extra
"channel"
dimension corresponding to red, green, blue, and
alpha/transparency. The values are scaled to the [0, 255] range and the
composite is an array of unsigned 8-bit integers.
Transparency
If any of the bands used for the composite have NaNs, those pixels will have their transparency set to the maximum value of 255. If there are no NaNs in any band, then the composite will only have 3 channels (red, green, blue).
Now do the same for the after scene:
after = after.assign(rgb=xls.composite(after, rescale_to=[0, 0.2]))
after = after.assign(
cir=xls.composite(after, bands=cir_bands, rescale_to=[0, 0.4]),
)
after
<xarray.Dataset> Dimensions: (easting: 400, northing: 300, channel: 4) Coordinates: * easting (easting) float64 5.844e+05 5.844e+05 ... 5.963e+05 5.964e+05 * northing (northing) float64 -2.232e+06 -2.232e+06 ... -2.223e+06 -2.223e+06 * channel (channel) <U5 'red' 'green' 'blue' 'alpha' Data variables: blue (northing, easting) float16 0.0686 0.07043 ... 0.05823 0.0564 green (northing, easting) float16 0.1027 0.09839 ... 0.07593 0.07043 red (northing, easting) float16 0.09778 0.09778 ... 0.06799 0.06177 nir (northing, easting) float16 0.2988 0.2715 0.2881 ... 0.2637 0.251 swir1 (northing, easting) float16 0.2311 0.2274 0.2316 ... 0.1608 0.142 swir2 (northing, easting) float16 0.145 0.1442 0.144 ... 0.09961 0.08655 rgb (northing, easting, channel) uint8 124 131 87 255 ... 78 89 71 255 cir (northing, easting, channel) uint8 190 62 65 255 ... 160 39 44 255 Attributes: (12/62) Conventions: CF-1.8 title: Landsat 8 scene from 2019-01-30 (path/ro... origin: Image courtesy of the U.S. Geological Su... digital_object_identifier: https://doi.org/10.5066/P9OGBGM6 landsat_product_id: LC08_L2SP_218074_20190130_20200829_02_T1 processing_level: L2SP ... ... temperature_mult_band_st_b10: 0.00341802 temperature_add_band_st_b10: 149.0 reflectance_mult_band_8: 2e-05 reflectance_mult_band_9: 2e-05 reflectance_add_band_8: -0.1 reflectance_add_band_9: -0.1
Plot the composites#
Composites can be plotted using xarray.DataArray.plot.imshow
(using
plot
won’t work and will display histograms instead).
Let’s make the before and after figures again for each of the composites we
generated.
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 12))
# Plot the composites
before.rgb.plot.imshow(ax=ax1)
after.rgb.plot.imshow(ax=ax2)
# The "long_name" of the composite is the band combination
ax1.set_title(f"Before: {before.rgb.attrs['long_name']}")
ax2.set_title(f"After: {after.rgb.attrs['long_name']}")
ax1.set_aspect("equal")
ax2.set_aspect("equal")
plt.show()
![_images/overview_7_0.png](_images/overview_7_0.png)
And now the CIR composites:
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 12))
before.cir.plot.imshow(ax=ax1)
after.cir.plot.imshow(ax=ax2)
ax1.set_title(f"Before: {before.cir.attrs['long_name']}")
ax2.set_title(f"After: {after.cir.attrs['long_name']}")
ax1.set_aspect("equal")
ax2.set_aspect("equal")
plt.show()
![_images/overview_8_0.png](_images/overview_8_0.png)
Calculating indices#
Producing indices from these scenes is straightforward thanks to
xarray
’s excelled support for coordinate-aware operations.
As an example, let’s calculate the
NDVI:
before = before.assign(
ndvi=(before.nir - before.red) / (before.nir + before.red),
)
after = after.assign(
ndvi=(after.nir - after.red) / (after.nir + after.red),
)
# Set some metadata for xarray to find
before.ndvi.attrs["long_name"] = "normalized difference vegetation index"
before.ndvi.attrs["units"] = "dimensionless"
after.ndvi.attrs["long_name"] = "normalized difference vegetation index"
after.ndvi.attrs["units"] = "dimensionless"
after
<xarray.Dataset> Dimensions: (easting: 400, northing: 300, channel: 4) Coordinates: * easting (easting) float64 5.844e+05 5.844e+05 ... 5.963e+05 5.964e+05 * northing (northing) float64 -2.232e+06 -2.232e+06 ... -2.223e+06 -2.223e+06 * channel (channel) <U5 'red' 'green' 'blue' 'alpha' Data variables: blue (northing, easting) float16 0.0686 0.07043 ... 0.05823 0.0564 green (northing, easting) float16 0.1027 0.09839 ... 0.07593 0.07043 red (northing, easting) float16 0.09778 0.09778 ... 0.06799 0.06177 nir (northing, easting) float16 0.2988 0.2715 0.2881 ... 0.2637 0.251 swir1 (northing, easting) float16 0.2311 0.2274 0.2316 ... 0.1608 0.142 swir2 (northing, easting) float16 0.145 0.1442 0.144 ... 0.09961 0.08655 rgb (northing, easting, channel) uint8 124 131 87 255 ... 78 89 71 255 cir (northing, easting, channel) uint8 190 62 65 255 ... 160 39 44 255 ndvi (northing, easting) float16 0.5073 0.4705 0.4907 ... 0.5903 0.605 Attributes: (12/62) Conventions: CF-1.8 title: Landsat 8 scene from 2019-01-30 (path/ro... origin: Image courtesy of the U.S. Geological Su... digital_object_identifier: https://doi.org/10.5066/P9OGBGM6 landsat_product_id: LC08_L2SP_218074_20190130_20200829_02_T1 processing_level: L2SP ... ... temperature_mult_band_st_b10: 0.00341802 temperature_add_band_st_b10: 149.0 reflectance_mult_band_8: 2e-05 reflectance_mult_band_9: 2e-05 reflectance_add_band_8: -0.1 reflectance_add_band_9: -0.1
And now we can make pseudocolor plots of the NDVI:
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 12))
# Limit the scale to [-1, +1] so the plots are easier to compare
before.ndvi.plot(ax=ax1, vmin=-1, vmax=1, cmap="RdBu_r")
after.ndvi.plot(ax=ax2, vmin=-1, vmax=1, cmap="RdBu_r")
ax1.set_title(f"Before: {before.attrs['title']}")
ax2.set_title(f"After: {after.attrs['title']}")
ax1.set_aspect("equal")
ax2.set_aspect("equal")
plt.show()
![_images/overview_10_0.png](_images/overview_10_0.png)
Finally, we can calculate the change in NDVI from one scene to the other by taking the difference:
ndvi_change = before.ndvi - after.ndvi
ndvi_change.name = "ndvi_change"
ndvi_change.attrs["long_name"] = (
f"NDVI change between {before.attrs['date_acquired']} and "
f"{after.attrs['date_acquired']}"
)
ndvi_change
<xarray.DataArray 'ndvi_change' (northing: 300, easting: 370)> array([[ 0.05908 , 0.06323 , 0.0542 , ..., 0.004395, -0.009766, 0.07324 ], [ 0.0498 , 0.07764 , 0.07495 , ..., 0.0957 , 0.012695, 0.003906], [-0.010254, 0.11743 , 0.0747 , ..., 0.03125 , 0.01807 , 0.04004 ], ..., [-0.000977, 0.01123 , 0.000977, ..., 0.00928 , 0.01367 , 0.00708 ], [ 0.00879 , 0.02344 , 0.01318 , ..., 0.006836, 0.00586 , 0.001221], [-0.01221 , 0.02637 , 0.006836, ..., 0.01343 , 0.01221 , 0.0105 ]], dtype=float16) Coordinates: * easting (easting) float64 5.844e+05 5.844e+05 ... 5.954e+05 5.955e+05 * northing (northing) float64 -2.232e+06 -2.232e+06 ... -2.223e+06 -2.223e+06 Attributes: long_name: NDVI change between 2019-01-14 and 2019-01-30
Did you notice?
The keen-eyed among you may have noticed that the number of points along
the "easting"
dimension has decreased. This is because xarray
only makes the calculations for pixels where the two scenes coincide. In
this case, there was an East-West shift between scenes but our calculations
take that into account.
Now lets plot it:
fig, ax = plt.subplots(1, 1, figsize=(10, 6))
ndvi_change.plot(ax=ax, vmin=-1, vmax=1, cmap="RdBu_r")
ax.set_aspect("equal")
plt.show()
![_images/overview_12_0.png](_images/overview_12_0.png)
There’s some noise in the cloudy areas of both scenes in the northeast but otherwise this plots highlights the area affected by flooding from the dam collapse in bright red at the center.
What now?#
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: interpolation, reduction, slicing, grouping, saving to
cloud-optimized formats, and much more. So go off and do something cool!