Basic example

Basic example

Let’s use EOReader for the first time !

import os

# First of all, we need some satellite data. 
# Let's open a lightweight a Landsat-5 TM collection 2 tile.
path = os.path.join("/home", "data", "DATA", "PRODS", "LANDSATS_COL2", "LT05_L1TP_200030_20111110_20200820_02_T1.tar")
from eoreader.reader import Reader

# Create the reader
# This reader is a singleton can be called once and then open all your data.
eoreader = Reader()
# Open your product
# No need to extract the product here: Archived Landsat-5 Collection are handled by EOReader
prod = eoreader.open(path)
# 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: 200030
Acquisition datetime: 2011-11-10 10:36:12
# Retrieve the UTM CRS of the tile
prod.crs
CRS.from_epsg(32630)
# Open here some more interesting geographical data: extent and footprint
base = prod.extent.plot(color='cyan', edgecolor='black')
prod.footprint.plot(ax=base, color='blue', edgecolor='black', alpha=0.5)
<AxesSubplot:>
../_images/base_6_1.png

See the difference between footprint and extent hereunder:

Without nodata

With nodata

without_nodata

with_nodata

from eoreader.bands import *
from eoreader.env_vars import DEM_PATH

# Select the bands you want to load
bands = [GREEN, NDVI, YELLOW, CLOUDS]

# 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 TM sensor:
ok_bands = [band for band in bands if prod.has_band(band)]
print(to_str(ok_bands))  # Landsat-5 TM doesn't provide YELLOW band
['GREEN', 'NDVI', 'CLOUDS']
# Load those bands as a dict of xarray.DataArray
band_dict = prod.load(ok_bands)
band_dict[GREEN]
<xarray.DataArray 'GREEN' (band: 1, y: 7131, x: 7991)>
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.16e+05 5.16e+05 5.161e+05 ... 7.557e+05 7.557e+05
  * y            (y) float64 4.888e+06 4.888e+06 ... 4.674e+06 4.674e+06
    spatial_ref  int64 0
Attributes:
    long_name:         GREEN
    sensor:            Landsat-5
    sensor_id:         L5
    product_path:      /home/data/DATA/PRODS/LANDSATS_COL2/LT05_L1TP_200030_2...
    product_name:      LT05_L1TP_200030_20111110_20200820_02_T1
    product_filename:  LT05_L1TP_200030_20111110_20200820_02_T1
    product_type:      TM
    acquisition_date:  20111110T103612
    condensed_name:    20111110T103612_L5_200030_TM
# The nan corresponds to the nodata you see on the footprint
# Plot a subsampled version
band_dict[GREEN][:, ::10, ::10].plot()
<matplotlib.collections.QuadMesh at 0x7f506aad02e0>
../_images/base_10_1.png
# Plot a subsampled version
band_dict[NDVI][:, ::10, ::10].plot()
<matplotlib.collections.QuadMesh at 0x7f5073106c10>
../_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 'GREEN NDVI CLOUDS' (z: 3, y: 7131, x: 7991)>
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]],

       [[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:
    spatial_ref  int64 0
  * x            (x) float64 5.16e+05 5.16e+05 5.161e+05 ... 7.557e+05 7.557e+05
  * y            (y) float64 4.888e+06 4.888e+06 ... 4.674e+06 4.674e+06
  * z            (z) MultiIndex
  - variable     (z) object 'GREEN' 'NDVI' 'CLOUDS'
  - band         (z) int64 1 1 1
Attributes:
    long_name:         GREEN NDVI CLOUDS
    sensor:            Landsat-5
    sensor_id:         L5
    product_path:      /home/data/DATA/PRODS/LANDSATS_COL2/LT05_L1TP_200030_2...
    product_name:      LT05_L1TP_200030_20111110_20200820_02_T1
    product_filename:  LT05_L1TP_200030_20111110_20200820_02_T1
    product_type:      TM
    acquisition_date:  20111110T103612
    condensed_name:    20111110T103612_L5_200030_TM
# Plot a subsampled version
import matplotlib.pyplot as plt

nrows = len(stack)
fig, axes = plt.subplots(nrows=nrows, figsize=(2 * nrows, 6 * nrows), subplot_kw={"box_aspect": 1})  # Square plots
for i in range(nrows):
    stack[i, ::10, ::10].plot(x="x", y="y", ax=axes[i])
../_images/base_14_0.png