Contents

VHR example

VHR example

Let’s use EOReader with Very High Resolution data.

Imports

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

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)

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
base = prod.extent.plot(color='cyan', edgecolor='black')
prod.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.
  pd.Int64Index,

<AxesSubplot:>
../_images/VHR_7_2.png

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
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 dict of xarray.DataArray
band_dict = prod.load(ok_bands)
band_dict[GREEN]

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)
Coordinates:
  * 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
Attributes:
    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 0x7fcad4628820>
../_images/VHR_11_1.png 
# Plot a subsampled version
band_dict[NDVI][:, ::10, ::10].plot()

<matplotlib.collections.QuadMesh at 0x7fcad45b64f0>
../_images/VHR_12_1.png 
# Plot a subsampled version
band_dict[CLOUDS][:, ::10, ::10].plot()

<matplotlib.collections.QuadMesh at 0x7fcad44e8610>
../_images/VHR_13_1.png 
# Plot a subsampled version
band_dict[HILLSHADE][:, ::10, ::10].plot()

<matplotlib.collections.QuadMesh at 0x7fcad473be20>
../_images/VHR_14_1.png

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])
../_images/VHR_17_0.png

previous

SAR example

next

Advanced