VHR example
Contents
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]
reader = Reader()
prod = reader.open(path, remove_tmp=True)
prod
EOReader PldProduct
Attributes:
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
NARROW_NIR: 4
needs extraction: False
cloud cover: 0.0
# Plot the quicklook
prod.plot()
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.
pd.Int64Index,
<AxesSubplot:>
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 GREEN to UTM with a 0.5 m resolution.
Reprojecting band NARROW_NIR to UTM with a 0.5 m resolution.
Reprojecting band RED 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:
cleaning_method: nodata
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
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_dict[GREEN][:, ::10, ::10].plot()
<matplotlib.collections.QuadMesh at 0x7f1ea5191be0>
# Plot a subsampled version
band_dict[NDVI][:, ::10, ::10].plot()
<matplotlib.collections.QuadMesh at 0x7f1ea51122b0>
# Plot a subsampled version
band_dict[CLOUDS][:, ::10, ::10].plot()
<matplotlib.collections.QuadMesh at 0x7f1ea50c21f0>
# Plot a subsampled version
band_dict[HILLSHADE][:, ::10, ::10].plot()
<matplotlib.collections.QuadMesh at 0x7f1ea4fe9bb0>
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])
