Contents

Optical example

Optical example

Let’s use EOReader with optical data.

Imports

import os
from pprint import pformat
import matplotlib.pyplot as plt

# EOReader
from eoreader.reader import Reader
from eoreader.bands import *
from eoreader.env_vars import DEM_PATH

Data

First of all, we need some satellite data. Let’s open a lightweight Landsat-5 TM collection 2 tile.

path = os.path.join("/home", "data", "DATA", "PRODS", "LANDSATS_COL2", "LT05_L1TP_200030_20111110_20200820_02_T1.tar")

Open the product

First, create the Reader object.
The Reader is a singleton that should be called only once. It can be used several times to open all your satellite products.

No need to extract the product here: Archived Landsat-5 Collection are handled by EOReader

eoreader = Reader()
prod = eoreader.open(path)

Product information

You have opened your product and you have its object in hands You can play a little with it to see what it is inside

print(
    f"Landsat tile: {prod.tile_name}\n\n"
    f"Acquisition datetime: {prod.datetime}\n\n"
    f"Existing bands:\n{pformat([band.value for band in prod.get_existing_bands()])}"
)

Landsat tile: 200030

Acquisition datetime: 2011-11-10 10:36:12

Existing bands:
['BLUE',
 'GREEN',
 'RED',
 'NIR',
 'NARROW_NIR',
 'SWIR_1',
 'SWIR_2',
 'THERMAL_IR_1',
 'THERMAL_IR_2']

# 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_9_1.png

See the difference between footprint and extent hereunder:

Without nodata

With nodata

without_nodata

with_nodata

Load bands

Let’s open some optical bands, cloud bands and maybe DEM bands.

# 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 0x7fd3543f1a60>
../_images/base_13_1.png 
# Plot a subsampled version
band_dict[NDVI][:, ::10, ::10].plot()

<matplotlib.collections.QuadMesh at 0x7fd35e2b13d0>
../_images/base_14_1.png 
# Plot a subsampled version
if HILLSHADE in band_dict:
    band_dict[HILLSHADE][:, ::10, ::10].plot()

Load a stack

You can also stack the bands you requested right before

# 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
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_18_0.png

previous

Basics

next

SAR example