VHR example¶
Let’s use EOReader with Very High Resolution data.
Warning:
import os
import glob
# First of all, we need some VHR data, let's use Pleiades data
path = glob.glob(os.path.join("/home", "data", "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()
# Open your product
prod = eoreader.open(path, remove_tmp=True)
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
from eoreader.bands.alias import *
from eoreader.env_vars import DEM_PATH
# 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 = ""
os.environ[DEM_PATH] = os.path.join("/home", "data", "DS2", "BASES_DE_DONNEES", "GLOBAL", "MERIT_Hydrologically_Adjusted_Elevations", "MERIT_DEM.vrt")
# 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']
# Load those bands as a dict of xarray.DataArray
band_dict = prod.load(ok_bands)
band_dict[GREEN]
2021-06-09 14:11:49,305 - eoreader - DEBUG - Loading bands ['NIR', 'GREEN', 'RED']
2021-06-09 14:11:49,305 - eoreader - DEBUG - Read NIR
2021-06-09 14:11:49,414 - eoreader - INFO - Reprojecting band NIR to UTM with a 0.5 m resolution.
2021-06-09 14:12:16,191 - eoreader - DEBUG - Manage invalid pixels for band NIR
2021-06-09 14:12:25,649 - eoreader - DEBUG - Read GREEN
2021-06-09 14:12:25,760 - eoreader - INFO - Reprojecting band GREEN to UTM with a 0.5 m resolution.
2021-06-09 14:12:53,321 - eoreader - DEBUG - Manage invalid pixels for band GREEN
2021-06-09 14:13:02,613 - eoreader - DEBUG - Read RED
2021-06-09 14:13:02,714 - eoreader - INFO - Reprojecting band RED to UTM with a 0.5 m resolution.
2021-06-09 14:13:29,244 - eoreader - DEBUG - Manage invalid pixels for band RED
2021-06-09 14:13:38,462 - eoreader - DEBUG - Loading index ['NDVI']
2021-06-09 14:13:39,407 - eoreader - DEBUG - Loading DEM bands ['HILLSHADE']
2021-06-09 14:13:39,407 - eoreader - DEBUG - Warping DEM for IMG_PHR1A_PMS_001
2021-06-09 14:13:39,409 - eoreader - DEBUG - Using DEM: D:\_EXTRACTEO\DS2\BASES_DE_DONNEES\GLOBAL\MERIT_Hydrologically_Adjusted_Elevations\MERIT_DEM.vrt
2021-06-09 14:13:43,182 - eoreader - DEBUG - Computing hillshade DEM for IMG_PHR1A_PMS_001
2021-06-09 14:13:54,877 - eoreader - DEBUG - Loading Cloud bands ['CLOUDS']
<xarray.DataArray 'DIM_PHR1A_PMS_202005110231585_ORT_5547047101' (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: * 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 * band (band) int32 2 spatial_ref int32 0 Attributes: (12/19) RADIANCE_BIAS: 0 RADIANCE_CALIBRATION_DATE: 2020-04-01T16:00:00.000Z RADIANCE_GAIN: 8.92 RADIANCE_MEASURE_DESC: Raw radiometric count (DN) to TOA ... RADIANCE_MEASURE_UNCERTAINTY: Specification accuracy value RADIANCE_MEASURE_UNIT: watt/m2/steradians/micrometers ... ... SPECTRAL_RANGE_MEASURE_DESC: Spectral range value of raw radiom... SPECTRAL_RANGE_MEASURE_UNCERTAINTY: Specification accuracy value SPECTRAL_RANGE_MEASURE_UNIT: micrometers SPECTRAL_RANGE_MIN: 0.43 scale_factor: 1.0 add_offset: 0.0
xarray.DataArray
'DIM_PHR1A_PMS_202005110231585_ORT_5547047101'
- band: 1
- y: 18124
- x: 16754
- nan nan nan nan nan nan nan nan ... nan nan nan nan nan nan nan nan
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)
- x(x)float647.024e+05 7.024e+05 ... 7.108e+05
- axis :
- X
- long_name :
- x coordinate of projection
- standard_name :
- projection_x_coordinate
- units :
- metre
array([702448.662765, 702449.162765, 702449.662765, ..., 710824.162765, 710824.662765, 710825.162765])
- y(y)float649.689e+06 9.689e+06 ... 9.68e+06
- axis :
- Y
- long_name :
- y coordinate of projection
- standard_name :
- projection_y_coordinate
- units :
- metre
array([9688620.534414, 9688620.034414, 9688619.534414, ..., 9679560.034414, 9679559.534414, 9679559.034414])
- band(band)int322
array([2])
- spatial_ref()int320
- crs_wkt :
- PROJCS["WGS 84 / UTM zone 50S",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",117],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",10000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32750"]]
- semi_major_axis :
- 6378137.0
- semi_minor_axis :
- 6356752.314245179
- inverse_flattening :
- 298.257223563
- reference_ellipsoid_name :
- WGS 84
- longitude_of_prime_meridian :
- 0.0
- prime_meridian_name :
- Greenwich
- geographic_crs_name :
- WGS 84
- horizontal_datum_name :
- World Geodetic System 1984
- projected_crs_name :
- WGS 84 / UTM zone 50S
- grid_mapping_name :
- transverse_mercator
- latitude_of_projection_origin :
- 0.0
- longitude_of_central_meridian :
- 117.0
- false_easting :
- 500000.0
- false_northing :
- 10000000.0
- scale_factor_at_central_meridian :
- 0.9996
- spatial_ref :
- PROJCS["WGS 84 / UTM zone 50S",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",117],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",10000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32750"]]
- GeoTransform :
- 702448.4127654727 0.5 0.0 9688620.784413863 0.0 -0.5
array(0)
- RADIANCE_BIAS :
- 0
- RADIANCE_CALIBRATION_DATE :
- 2020-04-01T16:00:00.000Z
- RADIANCE_GAIN :
- 8.92
- RADIANCE_MEASURE_DESC :
- Raw radiometric count (DN) to TOA Radiance (L). Formulae L=DN/GAIN+BIAS
- RADIANCE_MEASURE_UNCERTAINTY :
- Specification accuracy value
- RADIANCE_MEASURE_UNIT :
- watt/m2/steradians/micrometers
- SOLAR_IRRADIANCE_CALIBRATION_DATE :
- 2011-12-17T00:00:00.000Z
- SOLAR_IRRADIANCE_MEASURE_DESC :
- Solar irradiance value of raw radiometric Band
- SOLAR_IRRADIANCE_MEASURE_UNCERTAINTY :
- Specification
- SOLAR_IRRADIANCE_MEASURE_UNIT :
- watt/m2/micron
- SOLAR_IRRADIANCE_VALUE :
- 1915.0
- SPECTRAL_RANGE_CALIBRATION_DATE :
- 2011-12-17T00:00:00.000Z
- SPECTRAL_RANGE_MAX :
- 0.55
- SPECTRAL_RANGE_MEASURE_DESC :
- Spectral range value of raw radiometric Band
- SPECTRAL_RANGE_MEASURE_UNCERTAINTY :
- Specification accuracy value
- SPECTRAL_RANGE_MEASURE_UNIT :
- micrometers
- SPECTRAL_RANGE_MIN :
- 0.43
- scale_factor :
- 1.0
- add_offset :
- 0.0
# 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 0x13a3231ee10>
# Plot a subsampled version
band_dict[NDVI][:, ::10, ::10].plot()
<matplotlib.collections.QuadMesh at 0x13a3249e668>
# Plot a subsampled version
band_dict[CLOUDS][:, ::10, ::10].plot()
<matplotlib.collections.QuadMesh at 0x13a8e412470>
# Plot a subsampled version
band_dict[HILLSHADE][:, ::10, ::10].plot()
<matplotlib.collections.QuadMesh at 0x13a8e4d59b0>
# You can also stack those bands
stack = prod.stack(ok_bands)
stack
2021-06-09 14:14:08,764 - eoreader - DEBUG - Loading bands ['NIR', 'GREEN', 'RED']
2021-06-09 14:14:08,764 - eoreader - DEBUG - Read NIR
2021-06-09 14:14:12,928 - eoreader - DEBUG - Read GREEN
2021-06-09 14:14:16,838 - eoreader - DEBUG - Read RED
2021-06-09 14:14:20,723 - eoreader - DEBUG - Loading index ['NDVI']
2021-06-09 14:14:21,655 - eoreader - DEBUG - Loading DEM bands ['HILLSHADE']
2021-06-09 14:14:21,655 - eoreader - DEBUG - Already existing DEM for IMG_PHR1A_PMS_001. Skipping process.
2021-06-09 14:14:21,736 - eoreader - DEBUG - Already existing hillshade DEM for IMG_PHR1A_PMS_001. Skipping process.
2021-06-09 14:14:23,582 - eoreader - DEBUG - Loading Cloud bands ['CLOUDS']
<xarray.DataArray 'NDVI_GREEN_HILLSHADE_CLOUDS' (z: 4, 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]], [[ 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]], [[217., 203., 203., ..., 193., 193., 213.], [203., 203., 203., ..., 193., 193., 193.], [203., 203., 203., ..., 193., 193., 193.], ..., [ nan, nan, nan, ..., 174., 174., 174.], [ nan, nan, nan, ..., 176., 176., 176.], [ nan, nan, nan, ..., 177., 177., 196.]], [[ 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: spatial_ref int32 0 * 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 * z (z) MultiIndex - variable (z) object 'NDVI' 'GREEN' 'HILLSHADE' 'CLOUDS' - band (z) int64 1 1 1 1 Attributes: long_name: ['NDVI', 'GREEN', 'HILLSHADE', 'CLOUDS']
xarray.DataArray
'NDVI_GREEN_HILLSHADE_CLOUDS'
- z: 4
- y: 18124
- x: 16754
- nan nan nan nan nan nan nan nan ... nan nan nan nan nan nan nan nan
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]], [[ 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]], [[217., 203., 203., ..., 193., 193., 213.], [203., 203., 203., ..., 193., 193., 193.], [203., 203., 203., ..., 193., 193., 193.], ..., [ nan, nan, nan, ..., 174., 174., 174.], [ nan, nan, nan, ..., 176., 176., 176.], [ nan, nan, nan, ..., 177., 177., 196.]], [[ 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)
- spatial_ref()int320
- crs_wkt :
- PROJCS["WGS 84 / UTM zone 50S",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",117],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",10000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32750"]]
- semi_major_axis :
- 6378137.0
- semi_minor_axis :
- 6356752.314245179
- inverse_flattening :
- 298.257223563
- reference_ellipsoid_name :
- WGS 84
- longitude_of_prime_meridian :
- 0.0
- prime_meridian_name :
- Greenwich
- geographic_crs_name :
- WGS 84
- horizontal_datum_name :
- World Geodetic System 1984
- projected_crs_name :
- WGS 84 / UTM zone 50S
- grid_mapping_name :
- transverse_mercator
- latitude_of_projection_origin :
- 0.0
- longitude_of_central_meridian :
- 117.0
- false_easting :
- 500000.0
- false_northing :
- 10000000.0
- scale_factor_at_central_meridian :
- 0.9996
- spatial_ref :
- PROJCS["WGS 84 / UTM zone 50S",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",117],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",10000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32750"]]
- GeoTransform :
- 702448.4127654727 0.5 0.0 9688620.784413863 0.0 -0.5
array(0)
- x(x)float647.024e+05 7.024e+05 ... 7.108e+05
array([702448.662765, 702449.162765, 702449.662765, ..., 710824.162765, 710824.662765, 710825.162765])
- y(y)float649.689e+06 9.689e+06 ... 9.68e+06
array([9688620.534414, 9688620.034414, 9688619.534414, ..., 9679560.034414, 9679559.534414, 9679559.034414])
- z(z)MultiIndex(variable, band)
array([('NDVI', 1), ('GREEN', 1), ('HILLSHADE', 1), ('CLOUDS', 1)], dtype=object)
- variable(z)object'NDVI' 'GREEN' 'HILLSHADE' 'CLOUDS'
array(['NDVI', 'GREEN', 'HILLSHADE', 'CLOUDS'], dtype=object)
- band(z)int641 1 1 1
array([1, 1, 1, 1], dtype=int64)
- long_name :
- ['NDVI', 'GREEN', 'HILLSHADE', 'CLOUDS']
# 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])