VHR example

Let’s use EOReader with Very High Resolution data.


import os
import glob
import logging

# EOReader
from eoreader.reader import Reader
from eoreader.bands import *
from eoreader.env_vars import DEM_PATH

Create the logger

# Create logger

logger = logging.getLogger("eoreader")

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

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

# add formatter to ch

# add ch to 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 = glob.glob(os.path.join("/home", "data", "DATA", "PRODS", "PLEIADES", "5547047101", "IMG_PHR1A_PMS_001"))[0]
eoreader = Reader()
prod = eoreader.open(path, remove_tmp=True)
EOReader PldProduct
	condensed_name: 20200511T023158_PLD_ORT_PMS
	name: PHR1A_PMS_202005110231585_ORT_5547047101
	path: /home/data/DATA/PRODS/PLEIADES/5547047101/IMG_PHR1A_PMS_001
	platform: Pleiades
	sensor type: Optical
	product type: Ortho Single Image
	default resolution: 0.5
	acquisition datetime: 2020-05-11T02:31:58
	band mapping:
		BLUE: 3
		GREEN: 2
		RED: 1
		NIR: 4
	tile name: N/A
	needs_extraction: False
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
# 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)
/opt/conda/lib/python3.9/site-packages/geopandas/io/file.py:362: FutureWarning: pandas.Int64Index is deprecated and will be removed from pandas in a future version. Use pandas.Index with the appropriate dtype instead.

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 = "/path/to/ortho_stack.tif"

Load some bands

# Select the bands you want to load

# 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
# Load those bands as a dict of xarray.DataArray
band_dict = prod.load(ok_bands)
Reprojecting band RED to UTM with a 0.5 m resolution.
Reprojecting band GREEN to UTM with a 0.5 m resolution.
Reprojecting band NIR to UTM with a 0.5 m resolution.
<xarray.DataArray '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)
  * band         (band) int64 1
  * 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
    spatial_ref  int64 0
    scale_factor:      1.0
    add_offset:        0.0
    long_name:         GREEN
    sensor:            Pleiades
    sensor_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
# 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 0x7f77782a5670>
# Plot a subsampled version
band_dict[NDVI][:, ::10, ::10].plot()
<matplotlib.collections.QuadMesh at 0x7f77781821f0>
# Plot a subsampled version
band_dict[CLOUDS][:, ::10, ::10].plot()
<matplotlib.collections.QuadMesh at 0x7f77780b4310>
# Plot a subsampled version
band_dict[HILLSHADE][:, ::10, ::10].plot()
<matplotlib.collections.QuadMesh at 0x7f7778298c40>

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