Basic example

Let’s use EOReader for the first time !

Warning: To complete this tutorial you will need matplotlib

import os

# First of all, we need some satellite data. 
# Let's open a lightweight a Landsat-5 MSS collection 2 tile.
path = os.path.abspath("../../CI/DATA/LM05_L1TP_200029_19841014_20200902_02_T2.tar")
from eoreader.reader import Reader

# Create the reader
eoreader = Reader()

# This reader is a singleton can be called once and then open all your data.
# Use it like a logging.getLogger() instance
from eoreader.bands.alias import *

# Open your product
prod = eoreader.open(path) # No need to unzip here
print(prod)
<eoreader.products.optical.l5_product.L5Product object at 0x7f8219d4bfd0>
# Here you have opened your product and you have its object in hands
# You can play a little with it to see what it got inside
print(f"Landsat tile: {prod.tile_name}")
print(f"Acquisition datetime: {prod.datetime}")
Landsat tile: 200029
Acquisition datetime: 1984-10-14 10:18:17
# Open here some more interesting geographical data: extent
extent = prod.extent()
extent.geometry.iat[0]
../_images/base_5_0.svg
# Open here some more interesting geographical data: footprint
footprint = prod.footprint()
footprint.geometry.iat[0]
../_images/base_6_0.svg

See the difference between footprint and extent hereunder:

Without nodata

With nodata

without_nodata

with_nodata

from eoreader.env_vars import DEM_PATH
# Select the bands you want to load
bands = [GREEN, NDVI, TIR_1, SHADOWS]

# Compute DEM band only if you have set a DEM in your environment path
if DEM_PATH in os.environ:
    bands.append(HILLSHADE)

# Be sure they exist for Landsat-5 MSS sensor:
ok_bands = [band for band in bands if prod.has_band(band)]
print(to_str(ok_bands)) # Landsat-5 MSS doesn't provide TIR and SHADOWS bands
['GREEN', 'NDVI']
# Load those bands as a dict of xarray.DataArray
band_dict = prod.load(ok_bands)
band_dict[GREEN]
<xarray.DataArray 'LM05_L1TP_200029_19841014_20200902_02_T2' (band: 1, y: 3473, x: 3909)>
array([[[nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        ...,
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan]]], dtype=float32)
Coordinates:
  * band         (band) int64 1
  * x            (x) float64 5.549e+05 5.55e+05 ... 7.894e+05 7.894e+05
  * y            (y) float64 5.049e+06 5.049e+06 ... 4.841e+06 4.841e+06
    spatial_ref  int64 0
# The nan corresponds to the nodata you see on the footprint
%matplotlib inline

# Plot a subsampled version
band_dict[GREEN][:, ::10, ::10].plot()
<matplotlib.collections.QuadMesh at 0x7f81f2edb890>
../_images/base_10_1.png
# Plot a subsampled version
band_dict[NDVI][:, ::10, ::10].plot()
<matplotlib.collections.QuadMesh at 0x7f81f24c90d0>
../_images/base_11_1.png
# Plot a subsampled version
if HILLSHADE in band_dict:
    band_dict[HILLSHADE][:, ::10, ::10].plot()
# You can also stack those bands
stack = prod.stack(ok_bands)
stack
<xarray.DataArray 'NDVI_GREEN' (z: 2, y: 3473, x: 3909)>
array([[[nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        ...,
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan]],

       [[nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        ...,
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan]]], dtype=float32)
Coordinates:
  * x            (x) float64 5.549e+05 5.55e+05 ... 7.894e+05 7.894e+05
  * y            (y) float64 5.049e+06 5.049e+06 ... 4.841e+06 4.841e+06
    spatial_ref  int64 0
  * z            (z) MultiIndex
  - variable     (z) object 'NDVI' 'GREEN'
  - band         (z) int64 1 1
Attributes:
    long_name:  ['NDVI', 'GREEN']
# Error in plotting with a list
if "long_name" in stack.attrs:
    stack.attrs.pop("long_name")

# Plot a subsampled version
stack[:, ::10, ::10].plot(x="x", y="y", col="z")
<xarray.plot.facetgrid.FacetGrid at 0x7f821311b050>
../_images/base_15_1.png