Optical example#
Let’s use EOReader with optical data.
Imports#
import os
import matplotlib.pyplot as plt
# EOReader
from eoreader.reader import Reader
from eoreader.bands import RED, GREEN, NDVI, YELLOW, CLOUDS, to_str
Open the product#
First, open a Landsat-9 OLI-TIRS collection 2 product.
path = os.path.join("/home", "prods", "LANDSATS_COL2", "LC09_L1TP_200030_20220201_20220201_02_T1.tar")
reader = Reader()
prod = reader.open(path)
prod
eoreader.LandsatProduct 'LC09_L1TP_200030_20220201_20220201_02_T1'
Attributes:
condensed_name: 20220201T104852_L9_200030_OLI_TIRS
path: /home/prods/LANDSATS_COL2/LC09_L1TP_200030_20220201_20220201_02_T1.tar
constellation: Landsat-9
sensor type: Optical
product type: L1
default pixel size: 30.0
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
SWIR_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(nodata=0)
/opt/conda/lib/python3.11/site-packages/rasterio/__init__.py:356: NotGeoreferencedWarning: Dataset has no geotransform, gcps, or rpcs. The identity matrix will be returned.
dataset = DatasetReader(path, driver=driver, sharing=sharing, **kwargs)
# Get all the bands information
prod.bands
eoreader.SpectralBand 'Coastal aerosol'
Attributes:
id: 1
eoreader_name: COASTAL_AEROSOL
common_name: coastal
gsd (m): 30
asset_role: reflectance
Center wavelength (nm): 440.0
Bandwidth (nm): 20.0
description: Coastal and aerosol studies
eoreader.SpectralBand 'Blue'
Attributes:
id: 2
eoreader_name: BLUE
common_name: blue
gsd (m): 30
asset_role: reflectance
Center wavelength (nm): 480.0
Bandwidth (nm): 60.0
description: Bathymetric mapping, distinguishing soil from vegetation and deciduous from coniferous vegetation
eoreader.SpectralBand 'Green'
Attributes:
id: 3
eoreader_name: GREEN
common_name: green
gsd (m): 30
asset_role: reflectance
Center wavelength (nm): 560.0
Bandwidth (nm): 60.0
description: Emphasizes peak vegetation, which is useful for assessing plant vigor
eoreader.SpectralBand 'Red'
Attributes:
id: 4
eoreader_name: RED
common_name: red
gsd (m): 30
asset_role: reflectance
Center wavelength (nm): 655.0
Bandwidth (nm): 30.0
description: Discriminates vegetation slopes
eoreader.SpectralBand 'Near Infrared (NIR)'
Attributes:
id: 5
eoreader_name: NIR
common_name: nir
gsd (m): 30
asset_role: reflectance
Center wavelength (nm): 865.0
Bandwidth (nm): 30.0
description: Emphasizes biomass content and shorelines
eoreader.SpectralBand 'Near Infrared (NIR)'
Attributes:
id: 5
eoreader_name: NARROW_NIR
common_name: nir08
gsd (m): 30
asset_role: reflectance
Center wavelength (nm): 865.0
Bandwidth (nm): 30.0
description: Emphasizes biomass content and shorelines
eoreader.SpectralBand 'SWIR 1'
Attributes:
id: 6
eoreader_name: SWIR_1
common_name: swir16
gsd (m): 30
asset_role: reflectance
Center wavelength (nm): 1610.0
Bandwidth (nm): 80.0
description: Discriminates moisture content of soil and vegetation; penetrates thin clouds
eoreader.SpectralBand 'SWIR 2'
Attributes:
id: 7
eoreader_name: SWIR_2
common_name: swir22
gsd (m): 30
asset_role: reflectance
Center wavelength (nm): 2200.0
Bandwidth (nm): 180.0
description: Improved moisture content of soil and vegetation; penetrates thin clouds
eoreader.SpectralBand 'Panchromatic'
Attributes:
id: 8
eoreader_name: PANCHROMATIC
common_name: pan
gsd (m): 15
asset_role: reflectance
Center wavelength (nm): 590.0
Bandwidth (nm): 180.0
description: 15 meter resolution, sharper image definition
eoreader.SpectralBand 'Cirrus'
Attributes:
id: 9
eoreader_name: SWIR_CIRRUS
common_name: cirrus
gsd (m): 30
asset_role: reflectance
Center wavelength (nm): 1370.0
Bandwidth (nm): 20.0
description: Improved detection of cirrus cloud contamination
eoreader.SpectralBand 'Thermal Infrared (TIRS) 1'
Attributes:
id: 10
eoreader_name: THERMAL_IR_1
common_name: lwir11
gsd (m): 100
asset_role: brightness_temperature
Center wavelength (nm): 10895.0
Bandwidth (nm): 590.0
description: 100 meter resolution, thermal mapping and estimated soil moisture
eoreader.SpectralBand 'Thermal Infrared (TIRS) 2'
Attributes:
id: 11
eoreader_name: THERMAL_IR_2
common_name: lwir12
gsd (m): 100
asset_role: brightness_temperature
Center wavelength (nm): 12005.0
Bandwidth (nm): 1010.0
description: 100 meter resolution, improved thermal mapping and estimated soil moisture
# Some other information
print(f"Acquisition datetime: {prod.datetime}")
print(f"Condensed name: {prod.condensed_name}")
print(f"Landsat tile: {prod.tile_name}")
Acquisition datetime: 2022-02-01 10:48:52
Condensed name: 20220201T104852_L9_200030_OLI_TIRS
Landsat tile: 200030
# Retrieve the UTM CRS of the tile
prod.crs().to_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)
'_vectorize' function is not lazy yet. Computing the raster.
<Axes: >
See the difference between footprint and extent hereunder:
Without nodata |
With nodata |
|---|---|
|
|
Load bands#
Let’s open some spectral and cloud bands.
# Select some bands you want to load
bands = [GREEN, NDVI, YELLOW, CLOUDS]
# Be sure they exist for Landsat-9 OLI-TIRS sensor:
ok_bands = [band for band in bands if prod.has_band(band)]
print(to_str(ok_bands))
# Landsat-9 OLI-TIRS doesn't provide YELLOW band
['GREEN', 'NDVI', 'CLOUDS']
# Load those bands as a xarray.Dataset
band_ds = prod.load(ok_bands)
band_ds[GREEN]
<xarray.DataArray <SpectralBandNames.GREEN: 'GREEN'> (band: 1, y: 7791, x: 7681)> Size: 239MB
dask.array<clip, shape=(1, 7791, 7681), dtype=float32, chunksize=(1, 1024, 1024), chunktype=numpy.ndarray>
Coordinates:
spatial_ref int64 8B 0
* band (band) int64 8B 1
* y (y) float64 62kB 4.899e+06 4.899e+06 ... 4.665e+06 4.665e+06
* x (x) float64 61kB 5.193e+05 5.193e+05 ... 7.497e+05 7.497e+05
Attributes: (12/14)
path: /tmp/tmp3foe5gt2/tmp_20220201T104852_L9_200030_OLI_TIR...
long_name: GREEN
constellation: Landsat-9
constellation_id: L9
product_path: /home/prods/LANDSATS_COL2/LC09_L1TP_200030_20220201_20...
product_name: LC09_L1TP_200030_20220201_20220201_02_T1
... ...
product_type: L1
acquisition_date: 20220201T104852
condensed_name: 20220201T104852_L9_200030_OLI_TIRS
orbit_direction: DESCENDING
radiometry: reflectance
cloud_cover: 49.31# Plot a subsampled version of the GREEN band
# The nans corresponds to the nodata you see on the footprint
band_ds[GREEN][:, ::10, ::10].plot(robust=True, cmap="mako")
<matplotlib.collections.QuadMesh at 0x7fb0f1300250>
# Plot a subsampled version of the NDVI spectral index
band_ds[NDVI][:, ::10, ::10].plot(robust=True, cmap="mako")
<matplotlib.collections.QuadMesh at 0x7fb0c44f2010>
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' (bands: 3, y: 7792, x: 7696)> Size: 720MB
dask.array<where, shape=(3, 7792, 7696), dtype=float32, chunksize=(1, 1024, 1024), chunktype=numpy.ndarray>
Coordinates:
spatial_ref int64 8B 0
* y (y) float64 62kB 4.899e+06 4.899e+06 ... 4.665e+06 4.665e+06
* x (x) float64 62kB 5.193e+05 5.193e+05 ... 7.501e+05 7.502e+05
* bands (bands) object 24B MultiIndex
* variable (bands) object 24B SpectralBandNames.GREEN ... CloudsBandNam...
* band (bands) int64 24B 1 1 1
Attributes: (12/13)
path: /tmp/tmp3foe5gt2/tmp_20220201T104852_L9_200030_OLI_TIR...
long_name: GREEN NDVI CLOUDS
constellation: Landsat-9
constellation_id: L9
product_path: /home/prods/LANDSATS_COL2/LC09_L1TP_200030_20220201_20...
product_name: LC09_L1TP_200030_20220201_20220201_02_T1
... ...
instrument: OLI-TIRS
product_type: L1
acquisition_date: 20220201T104852
condensed_name: 20220201T104852_L9_200030_OLI_TIRS
orbit_direction: DESCENDING
cloud_cover: 49.31# Plot a subsampled version of the stack
from sertit import display
display_stack = stack[:, ::5, ::5]
display_stack.copy(data=display.scale(display_stack.data)).plot.imshow(robust=True);
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})
<xarray.Dataset> Size: 240MB
Dimensions: (x: 7696, y: 7792, band: 1)
Coordinates:
* x (x) float64 62kB 5.193e+05 5.193e+05 ... 7.502e+05
* y (y) float64 62kB 4.899e+06 4.899e+06 ... 4.665e+06
spatial_ref int64 8B 0
* band (band) int64 8B 1
Data variables:
SpectralBandNames.RED (band, y, x) float32 240MB dask.array<chunksize=(1, 1024, 1024), meta=np.ndarray>
Attributes: (12/13)
long_name: RED
constellation: Landsat-9
constellation_id: L9
product_path: /home/prods/LANDSATS_COL2/LC09_L1TP_200030_20220201_20...
product_name: LC09_L1TP_200030_20220201_20220201_02_T1
product_filename: LC09_L1TP_200030_20220201_20220201_02_T1
... ...
product_type: L1
acquisition_date: 20220201T104852
condensed_name: 20220201T104852_L9_200030_OLI_TIRS
orbit_direction: DESCENDING
radiometry: reflectance
cloud_cover: 49.31# 'As is' band
prod.load(RED, **{TO_REFLECTANCE: False})
<xarray.Dataset> Size: 239MB
Dimensions: (x: 7681, y: 7791, band: 1)
Coordinates:
* x (x) float64 61kB 5.193e+05 5.193e+05 ... 7.497e+05
* y (y) float64 62kB 4.899e+06 4.899e+06 ... 4.665e+06
spatial_ref int64 8B 0
* band (band) int64 8B 1
Data variables:
SpectralBandNames.RED (band, y, x) float32 239MB dask.array<chunksize=(1, 1024, 1024), meta=np.ndarray>
Attributes: (12/13)
long_name: RED
constellation: Landsat-9
constellation_id: L9
product_path: /home/prods/LANDSATS_COL2/LC09_L1TP_200030_20220201_20...
product_name: LC09_L1TP_200030_20220201_20220201_02_T1
product_filename: LC09_L1TP_200030_20220201_20220201_02_T1
... ...
product_type: L1
acquisition_date: 20220201T104852
condensed_name: 20220201T104852_L9_200030_OLI_TIRS
orbit_direction: DESCENDING
radiometry: as is
cloud_cover: 49.31