VHR example#

Let’s use EOReader with Very High Resolution data.

Imports#

import os

# EOReader
from eoreader.reader import Reader
from eoreader.bands import GREEN, NDVI, TIR_1, CLOUDS, HILLSHADE, to_str
from eoreader.env_vars import DEM_PATH

Create the logger#

# Create logger
import logging
from sertit import logs

logger = logging.getLogger("eoreader")
logs.init_logger(logger)

Open the VHR product#

Please be aware that EOReader will always work in UTM projection.
So if you give WGS84 data, EOReader will reproject the stacks and this can be time-consuming

# Set a DEM
os.environ[DEM_PATH] = os.path.join(
    "/home", "data", "DS2", "BASES_DE_DONNEES", "GLOBAL",
    "MERIT_Hydrologically_Adjusted_Elevations", "MERIT_DEM.vrt"
)

# Open your product
path = os.path.join("/home", "data", "DATA", "PRODS", "PLEIADES", "5547047101", "IMG_PHR1A_PMS_001")
reader = Reader()
prod = reader.open(path, remove_tmp=True)
prod
eoreader.PldProduct 'PHR1A_PMS_202005110231585_ORT_5547047101'
Attributes:
	condensed_name: 20200511T023158_PLD_ORT_PMS_5547047101
	path: /home/data/DATA/PRODS/PLEIADES/5547047101/IMG_PHR1A_PMS_001
	constellation: Pleiades
	sensor type: Optical
	product type: Ortho Single Image
	default pixel size: 0.5
	default resolution: 0.5
	acquisition datetime: 2020-05-11T02:31:58
	band mapping:
		BLUE: 3
		GREEN: 2
		RED: 1
		NIR: 4
		NARROW_NIR: 4
	needs extraction: False
	cloud cover: 0.0
# Plot the quicklook
prod.plot()
../_images/4400d95d3d2cb809f2e905f2b9207e4ad4ce2e185a5401555f1759290c316ec0.png
# Get the bands information
prod.bands
eoreader.SpectralBand 'RED'
Attributes:
	id: 1
	eoreader_name: RED
	common_name: red
	gsd (m): 0.5
	asset_role: reflectance
	Center wavelength (nm): 650.0
	Bandwidth (nm): 120.0
eoreader.SpectralBand 'GREEN'
Attributes:
	id: 2
	eoreader_name: GREEN
	common_name: green
	gsd (m): 0.5
	asset_role: reflectance
	Center wavelength (nm): 560.0
	Bandwidth (nm): 120.0
eoreader.SpectralBand 'BLUE'
Attributes:
	id: 3
	eoreader_name: BLUE
	common_name: blue
	gsd (m): 0.5
	asset_role: reflectance
	Center wavelength (nm): 495.0
	Bandwidth (nm): 70.0
eoreader.SpectralBand 'NIR'
Attributes:
	id: 4
	eoreader_name: NIR
	common_name: nir
	gsd (m): 0.5
	asset_role: reflectance
	Center wavelength (nm): 840.0
	Bandwidth (nm): 200.0
eoreader.SpectralBand 'NIR'
Attributes:
	id: 4
	eoreader_name: NIR
	common_name: nir
	gsd (m): 0.5
	asset_role: reflectance
	Center wavelength (nm): 840.0
	Bandwidth (nm): 200.0
print(f"Acquisition datetime: {prod.datetime}")
print(f"Condensed name: {prod.condensed_name}")
Acquisition datetime: 2020-05-11 02:31:58
Condensed name: 20200511T023158_PLD_ORT_PMS_5547047101
# 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)
<Axes: >
../_images/4ab09fd92b381295c1efea6b89567e155f41878383bf7bb65b0a18aaf21f4b31.png

Here, if you want to orthorectify or pansharpen your data manually, you can set your stack here.

prod.ortho_stack = "/path/to/ortho_stack.tif"

If you do not provide this stack, but you give a non-orthorectified product to EOReader (i.e. SEN or PRJ products for Pleiades), you must provide a DEM to orthorectify correctly the data.

⚠⚠⚠ DIMAP SEN products are orthorectified using RPCs and not the rigorous sensor model. A shift of several meters can occur. Please refer to this issue.

Load some bands#

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

# Be sure they exist for Pleiades sensor
ok_bands = [band for band in bands if prod.has_band(band)]
print(to_str(ok_bands))  # Pleiades doesn't provide TIR and SHADOWS bands
['GREEN', 'NDVI', 'CLOUDS', 'HILLSHADE']
# Load those bands as a xarray.Dataset
band_ds = prod.load(ok_bands)
band_ds[GREEN]
2023-11-02 16:09:39,773 - [DEBUG] - Loading bands ['GREEN', 'NIR', 'RED']
2023-11-02 16:09:39,778 - [DEBUG] - Read GREEN
2023-11-02 16:09:39,853 - [INFO] - Warping DIM_PHR1A_PMS_202005110231585_ORT_5547047101 to UTM with a 0.5 m pixel size.
2023-11-02 16:09:39,931 - [DEBUG] - Reading warped GREEN.
2023-11-02 16:09:40,129 - [DEBUG] - Manage nodata for band GREEN
2023-11-02 16:09:41,010 - [DEBUG] - Converting GREEN to reflectance
2023-11-02 16:10:43,805 - [DEBUG] - Read NIR
2023-11-02 16:10:43,856 - [INFO] - Warping DIM_PHR1A_PMS_202005110231585_ORT_5547047101 to UTM with a 0.5 m pixel size.
2023-11-02 16:10:43,922 - [DEBUG] - Reading warped NIR.
2023-11-02 16:10:43,976 - [DEBUG] - Manage nodata for band NIR
2023-11-02 16:10:44,821 - [DEBUG] - Converting NIR to reflectance
2023-11-02 16:11:47,315 - [DEBUG] - Read RED
2023-11-02 16:11:47,340 - [INFO] - Warping DIM_PHR1A_PMS_202005110231585_ORT_5547047101 to UTM with a 0.5 m pixel size.
2023-11-02 16:11:47,394 - [DEBUG] - Reading warped RED.
2023-11-02 16:11:47,432 - [DEBUG] - Manage nodata for band RED
2023-11-02 16:11:48,265 - [DEBUG] - Converting RED to reflectance
2023-11-02 16:14:57,462 - [DEBUG] - Loading indices ['NDVI']
2023-11-02 16:15:00,112 - [DEBUG] - Loading DEM bands ['HILLSHADE']
2023-11-02 16:15:00,113 - [DEBUG] - Warping DEM for 20200511T023158_PLD_ORT_PMS_5547047101
2023-11-02 16:15:00,115 - [DEBUG] - Using DEM: /home/data/DS2/BASES_DE_DONNEES/GLOBAL/MERIT_Hydrologically_Adjusted_Elevations/MERIT_DEM.vrt
2023-11-02 16:15:00,240 - [DEBUG] - Computing hillshade DEM for PHR1A_PMS_202005110231585_ORT_5547047101
2023-11-02 16:15:23,821 - [DEBUG] - Loading Cloud bands ['CLOUDS']
<xarray.DataArray <SpectralBandNames.GREEN: 'GREEN'> (band: 1, y: 18124,
                                                      x: 16754)>
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
    spatial_ref  int64 0
  * x            (x) float64 7.024e+05 7.024e+05 ... 7.108e+05 7.108e+05
  * y            (y) float64 9.689e+06 9.689e+06 9.689e+06 ... 9.68e+06 9.68e+06
Attributes: (12/13)
    long_name:         GREEN
    constellation:     Pleiades
    constellation_id:  PLD
    product_path:      /home/data/DATA/PRODS/PLEIADES/5547047101/IMG_PHR1A_PM...
    product_name:      PHR1A_PMS_202005110231585_ORT_5547047101
    product_filename:  IMG_PHR1A_PMS_001
    ...                ...
    product_type:      Ortho Single Image
    acquisition_date:  20200511T023158
    condensed_name:    20200511T023158_PLD_ORT_PMS_5547047101
    orbit_direction:   DESCENDING
    radiometry:        reflectance
    cloud_cover:       0.0
# The nan corresponds to the nodata you see on the footprint
# Plot a subsampled version
band_ds[GREEN][:, ::10, ::10].plot()
<matplotlib.collections.QuadMesh at 0x7f0a2c65de10>
../_images/a43def205fabdae01aa3b5fdbb9b682bcb80bd57b9a65626943565d187a61c23.png
# Plot a subsampled version
band_ds[NDVI][:, ::10, ::10].plot()
<matplotlib.collections.QuadMesh at 0x7f0a2c524310>
../_images/25e70f8f6a7020e5b8b77ea33f65c38f0e71324fd169f54b0bd31bbab032a19a.png
# Plot a subsampled version
band_ds[CLOUDS][:, ::10, ::10].plot()
<matplotlib.collections.QuadMesh at 0x7f0a2c74b710>
../_images/5abf0862912c4b5a96c192fa911f9c04324880b5a1959eae6cd612ca2673ce47.png
# Plot a subsampled version
band_ds[HILLSHADE][:, ::10, ::10].plot()
<matplotlib.collections.QuadMesh at 0x7f0a2c5f37d0>
../_images/30a178937e22d5556a27bdca1bfed82920eb1d86fd7c49e0af5f6daaaa35093b.png

Stack some bands#

# You can also stack those bands
stack = prod.stack(ok_bands)
2023-11-02 16:16:56,952 - [DEBUG] - Loading bands ['GREEN', 'NIR', 'RED']
2023-11-02 16:16:56,952 - [DEBUG] - Read GREEN
2023-11-02 16:16:56,971 - [DEBUG] - Read NIR
2023-11-02 16:16:56,989 - [DEBUG] - Read RED
2023-11-02 16:17:19,070 - [DEBUG] - Loading indices ['NDVI']
2023-11-02 16:17:19,088 - [DEBUG] - Loading DEM bands ['HILLSHADE']
2023-11-02 16:17:19,088 - [DEBUG] - Already existing DEM for 20200511T023158_PLD_ORT_PMS_5547047101. Skipping process.
2023-11-02 16:17:19,089 - [DEBUG] - Already existing hillshade DEM for PHR1A_PMS_202005110231585_ORT_5547047101. Skipping process.
2023-11-02 16:17:19,107 - [DEBUG] - Loading Cloud bands ['CLOUDS']
2023-11-02 16:17:49,692 - [DEBUG] - Stacking
/opt/conda/lib/python3.11/site-packages/xarray/core/indexes.py:662: RuntimeWarning: '<' not supported between instances of 'CloudsBandNames' and 'SpectralBandNames', sort order is undefined for incomparable objects.
  new_pd_index = pd_indexes[0].append(pd_indexes[1:])
# 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})

for i in range(nrows):
    stack[i, ::10, ::10].plot(x="x", y="y", ax=axes[i])
../_images/75c9a41e7ce36a40d848da2e0316a9be08315731cad25c949b8e791783c1b9ac.png