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

Open the product#

First, open a Landsat-9 OLCI collection 2 product.

path = os.path.join("/home", "data", "DATA", "PRODS", "LANDSATS_COL2", "LC09_L1TP_200030_20220201_20220201_02_T1.tar")
reader = Reader()
prod = reader.open(path)
prod
EOReader L9Product
Attributes:
	condensed_name: 20220201T104852_L9_200030_OLCI
	name: LC09_L1TP_200030_20220201_20220201_02_T1
	path: /home/data/DATA/PRODS/LANDSATS_COL2/LC09_L1TP_200030_20220201_20220201_02_T1.tar
	platform: Landsat-9
	sensor type: Optical
	product type: OLCI
	default resolution: 30.0
	acquisition datetime: 2022-02-01T10:48:52
	band mapping:
		COASTAL_AEROSOL: 1
		BLUE: 2
		GREEN: 3
		RED: 4
		NIR: 5
		NARROW_NIR: 5
		CIRRUS: 9
		SWIR_1: 6
		SWIR_2: 7
		THERMAL_IR_1: 10
		THERMAL_IR_2: 11
		PANCHROMATIC: 8
	needs extraction: False
	cloud cover: 49.31
	tile name: 200030

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

# Plot the quicklook
prod.plot()
../_images/optical_5_0.png
# Some other information
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: 2022-02-01 10:48:52

Existing bands:
['COASTAL_AEROSOL',
 'BLUE',
 'GREEN',
 'RED',
 'NIR',
 'NARROW_NIR',
 'CIRRUS',
 'SWIR_1',
 'SWIR_2',
 'THERMAL_IR_1',
 'THERMAL_IR_2',
 'PANCHROMATIC']
# Retrieve the UTM CRS of the tile
prod.crs()
CRS.from_epsg(32630)
# Open here some more interesting geographical data: extent and footprint
extent = prod.extent()
footprint = prod.footprint()

base = extent.plot(color='cyan', edgecolor='black')
footprint.plot(ax=base, color='blue', edgecolor='black', alpha=0.5)
<AxesSubplot:>
../_images/optical_8_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 some 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: 7791, x: 7681)>
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.193e+05 5.193e+05 ... 7.497e+05 7.497e+05
  * y            (y) float64 4.899e+06 4.899e+06 ... 4.665e+06 4.665e+06
    spatial_ref  int64 0
Attributes:
    cleaning_method:   nodata
    long_name:         GREEN
    sensor:            Landsat-9
    sensor_id:         L9
    product_path:      /home/data/DATA/PRODS/LANDSATS_COL2/LC09_L1TP_200030_2...
    product_name:      LC09_L1TP_200030_20220201_20220201_02_T1
    product_filename:  LC09_L1TP_200030_20220201_20220201_02_T1
    product_type:      OLCI
    acquisition_date:  20220201T104852
    condensed_name:    20220201T104852_L9_200030_OLCI
    orbit_direction:   DESCENDING
    radiometry:        reflectance
    cloud_cover:       49.31
# 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 0x7ff670813d00>
../_images/optical_12_1.png
# Plot a subsampled version
band_dict[NDVI][:, ::10, ::10].plot()
<matplotlib.collections.QuadMesh at 0x7ff6707bef40>
../_images/optical_13_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: 7791, x: 7681)>
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.193e+05 5.193e+05 ... 7.497e+05 7.497e+05
  * y            (y) float64 4.899e+06 4.899e+06 ... 4.665e+06 4.665e+06
  * z            (z) MultiIndex
  - variable     (z) object 'GREEN' 'NDVI' 'CLOUDS'
  - band         (z) int64 1 1 1
Attributes:
    long_name:         GREEN NDVI CLOUDS
    sensor:            Landsat-9
    sensor_id:         L9
    product_path:      /home/data/DATA/PRODS/LANDSATS_COL2/LC09_L1TP_200030_2...
    product_name:      LC09_L1TP_200030_20220201_20220201_02_T1
    product_filename:  LC09_L1TP_200030_20220201_20220201_02_T1
    product_type:      OLCI
    acquisition_date:  20220201T104852
    condensed_name:    20220201T104852_L9_200030_OLCI
    orbit_direction:   DESCENDING
    radiometry:        reflectance
    cloud_cover:       49.31
# 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/optical_17_0.png

Radiometric processing#

EOReader allows you to load the band array as provided (either in DN, scaled radiance or reflectance). However, by default, EOReader converts all optical bands in reflectance (except for the thermal bands that are left if brilliance temperature)

The radiometric processing is described in the band’s attribute, either reflectance or as is

# Reflectance band
from eoreader.keywords import TO_REFLECTANCE
prod.load(RED, **{TO_REFLECTANCE: True})
{<OpticalBandNames.RED: 'RED'>: <xarray.DataArray 'RED' (band: 1, y: 7791, x: 7681)>
 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.193e+05 5.193e+05 ... 7.497e+05 7.497e+05
   * y            (y) float64 4.899e+06 4.899e+06 ... 4.665e+06 4.665e+06
     spatial_ref  int64 0
 Attributes:
     scale_factor:      1.0
     add_offset:        0.0
     long_name:         RED
     sensor:            Landsat-9
     sensor_id:         L9
     product_path:      /home/data/DATA/PRODS/LANDSATS_COL2/LC09_L1TP_200030_2...
     product_name:      LC09_L1TP_200030_20220201_20220201_02_T1
     product_filename:  LC09_L1TP_200030_20220201_20220201_02_T1
     product_type:      OLCI
     acquisition_date:  20220201T104852
     condensed_name:    20220201T104852_L9_200030_OLCI
     orbit_direction:   DESCENDING
     radiometry:        reflectance
     cloud_cover:       49.31}
# As is band
prod.load(RED, **{TO_REFLECTANCE: False})
{<OpticalBandNames.RED: 'RED'>: <xarray.DataArray 'RED' (band: 1, y: 7791, x: 7681)>
 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.193e+05 5.193e+05 ... 7.497e+05 7.497e+05
   * y            (y) float64 4.899e+06 4.899e+06 ... 4.665e+06 4.665e+06
     spatial_ref  int64 0
 Attributes:
     scale_factor:      1.0
     add_offset:        0.0
     cleaning_method:   nodata
     long_name:         RED
     sensor:            Landsat-9
     sensor_id:         L9
     product_path:      /home/data/DATA/PRODS/LANDSATS_COL2/LC09_L1TP_200030_2...
     product_name:      LC09_L1TP_200030_20220201_20220201_02_T1
     product_filename:  LC09_L1TP_200030_20220201_20220201_02_T1
     product_type:      OLCI
     acquisition_date:  20220201T104852
     condensed_name:    20220201T104852_L9_200030_OLCI
     orbit_direction:   DESCENDING
     radiometry:        as is
     cloud_cover:       49.31}