VHR example¶

Let’s use EOReader with Very High Resolution data.

Warning:

We do not provide Pleiades data

You will need matplotlib to complete this tutorial

import os
import glob

# First of all, we need some VHR data, let's use Pleiades data
path = os.path.abspath(glob.glob(os.path.join("/", "*", "DATA", "PRODS", "PLEIADES", "5547047101", "IMG_PHR1A_PMS_001"))[0])

# Create logger
import logging

logger = logging.getLogger("eoreader")
logger.setLevel(logging.INFO)

# create console handler and set level to debug
ch = logging.StreamHandler()
ch.setLevel(logging.INFO)

# create formatter
formatter = logging.Formatter('%(message)s')

# add formatter to ch
ch.setFormatter(formatter)

# add ch to logger
logger.addHandler(ch)

from eoreader.reader import Reader

# Create the reader
eoreader = Reader()

from eoreader.bands.alias import *

# Open your product
prod = eoreader.open(path)
print(f"Acquisition datetime: {prod.datetime}")
print(f"Condensed name: {prod.condensed_name}")

# 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

Acquisition datetime: 2020-05-11 02:31:58
Condensed name: 20200511T023158_PLD_ORT_PMS

# Here, if you want to orthorectify or pansharpen your data manually, you can set your stack here.
# If you do not provide this stack but you give a non-orthorectified product to EOReader 
# (ie. SEN or PRJ products for Pleiades), you must provide a DEM to orthorectify correctly the data
# prod.ortho_stack = ""
from eoreader.env_vars import DEM_PATH
os.environ[DEM_PATH] = os.path.abspath(glob.glob(os.path.join("/", "*", "DS2", "BASES_DE_DONNEES", "GLOBAL", "MERIT_Hydrologically_Adjusted_Elevations", "MERIT_DEM.vrt"))[0])

# Open here some more interesting geographical data: extent
extent = prod.extent()
extent.geometry.to_crs("EPSG:4326").iat[0]  # Display

# Open here some more interesting geographical data: footprint
footprint = prod.footprint()
footprint.geometry.to_crs("EPSG:4326").iat[0]  # Display

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

# The nan corresponds to the nodata you see on the footprint
%matplotlib inline

# Plot a subsampled version
band_dict[GREEN][:, ::10, ::10].plot()

<matplotlib.collections.QuadMesh at 0x220c7ce6240>

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

<matplotlib.collections.QuadMesh at 0x220c7aa2f98>

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

<matplotlib.collections.QuadMesh at 0x220c7b76c50>

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

<matplotlib.collections.QuadMesh at 0x220c7c38f60>

# Error in plotting with a list
if "long_name" in stack.attrs:
    stack.attrs.pop("long_name")

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