VHR example

Contents

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

/home/docs/checkouts/readthedocs.org/user_builds/eoreader/envs/latest/lib/python3.9/site-packages/dask/dataframe/__init__.py:42: FutureWarning: 
Dask dataframe query planning is disabled because dask-expr is not installed.

You can install it with `pip install dask[dataframe]` or `conda install dask`.
This will raise in a future version.

  warnings.warn(msg, FutureWarning)

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", "ds2_db2", "BASES_DE_DONNEES", "GLOBAL", "COPDEM_30m", "COPDEM_30m.vrt"
)

# Open your product
path = os.path.join("/home", "prods", "PLEIADES", "5547047101", "IMG_PHR1A_PMS_001")
reader = Reader()
prod = reader.open(path, remove_tmp=True)
prod

2025-05-05 11:13:38,413 - [WARNING] - There is no existing products in EOReader corresponding to /home/prods/PLEIADES/5547047101/IMG_PHR1A_PMS_001.
2025-05-05 11:13:38,414 - [INFO] - Your given path may not be a satellite image. If it is, maybe the product isn't handled by EOReader. If you are sure this product is handled, it is either corrupted or you may need to go deeper in the filetree to find the correct path to give.
2025-05-05 11:13:38,415 - [DEBUG] - Please look at what folder you should give to EOReader by accessing the documentation: https://eoreader.readthedocs.io/latest/main_features.html#recognized-paths

# Plot the quicklook
prod.plot()

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In[4], line 2
      1 # Plot the quicklook
----> 2 prod.plot()

AttributeError: 'NoneType' object has no attribute 'plot'

# 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/f095a4abc9d548967541a94601c5146822d0f7c90a59df68072f286e57304786.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, pixel_size=2.0)
band_ds[GREEN]

2023-05-31 11:47:20,447 - [DEBUG] - Loading bands ['GREEN', 'NIR', 'RED']
2023-05-31 11:47:20,453 - [DEBUG] - Read GREEN
2023-05-31 11:47:20,573 - [INFO] - Warping DIM_PHR1A_PMS_202005110231585_ORT_5547047101 to UTM with a 0.5 m pixel size.
2023-05-31 11:47:20,844 - [DEBUG] - Reading warped GREEN.
2023-05-31 11:47:21,261 - [DEBUG] - Manage nodata for band GREEN
2023-05-31 11:47:28,021 - [DEBUG] - Converting GREEN to reflectance
2023-05-31 11:58:09,959 - [DEBUG] - Read NIR
2023-05-31 11:58:10,429 - [INFO] - Warping DIM_PHR1A_PMS_202005110231585_ORT_5547047101 to UTM with a 0.5 m pixel size.
2023-05-31 11:58:10,550 - [DEBUG] - Reading warped NIR.
2023-05-31 11:58:11,231 - [DEBUG] - Manage nodata for band NIR
2023-05-31 11:58:12,656 - [DEBUG] - Converting NIR to reflectance
2023-05-31 11:59:40,816 - [DEBUG] - Read RED
2023-05-31 11:59:40,902 - [INFO] - Warping DIM_PHR1A_PMS_202005110231585_ORT_5547047101 to UTM with a 0.5 m pixel size.
2023-05-31 11:59:40,964 - [DEBUG] - Reading warped RED.
2023-05-31 11:59:41,016 - [DEBUG] - Manage nodata for band RED
2023-05-31 11:59:42,395 - [DEBUG] - Converting RED to reflectance

# The nan corresponds to the nodata you see on the footprint
# Plot a subsampled version
band_ds[GREEN][:, ::10, ::10].plot()

# Plot a subsampled version
band_ds[NDVI][:, ::10, ::10].plot()

# Plot a subsampled version
band_ds[CLOUDS][:, ::10, ::10].plot()

# Plot a subsampled version
band_ds[HILLSHADE][:, ::10, ::10].plot()

Stack some bands#

# You can also stack those bands
stack = prod.stack(ok_bands)

# 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])