Custom stacks#
Let’s use EOReader with custom stacks.
# EOReader Imports
import os
import xarray as xr
from eoreader.reader import Reader
from eoreader.products import SensorType
from eoreader.bands import BLUE, GREEN, RED, NIR, SWIR_1, VV, VV_DSPK, SLOPE, HILLSHADE
from sertit import display
reader = Reader()
# Create logger
import logging
from sertit import logs
logger = logging.getLogger("eoreader")
logs.init_logger(logger)
# Set a DEM
from eoreader.env_vars import DEM_PATH
os.environ[DEM_PATH] = os.path.join("/home", "data", "DS2", "BASES_DE_DONNEES", "GLOBAL", "COPDEM_30m",
"COPDEM_30m.vrt")
Custom stack with minimum data#
For both SAR and optical stacks, the two minimum keywords to provide are:
band_map
: a dictionary mapping the satellite band to the band number (starting to 1, in GDAL style)sensor_type
: EitherSAR
orOPTICAL
(a string or a SensorType Enum)
# Paths
stack_folder = os.path.join("/home", "data", "DS3", "CI", "eoreader", "others")
opt_path = os.path.join(stack_folder, "20200310T030415_WV02_Ortho_BGRN_STK.tif")
sar_path = os.path.join(stack_folder, "20210827T162210_ICEYE_SC_GRD_STK.tif")
# Optical minimum example
opt_prod = reader.open(opt_path,
custom=True,
sensor_type="OPTICAL", # With a string
band_map={BLUE: 1, GREEN: 2, RED: 3, NIR: 4, SWIR_1: 5})
opt_prod
eoreader.CustomProduct '20200310T030415_WV02_Ortho_BGRN_STK'
Attributes:
condensed_name: 20230531T121039_CUSTOM_CUSTOM
path: /home/data/DS3/CI/eoreader/others/20200310T030415_WV02_Ortho_BGRN_STK.tif
constellation: CUSTOM
sensor type: Optical
product type: CUSTOM
default pixel size: 8.0
default resolution: None
acquisition datetime: 2023-05-31T12:10:39.384024
band mapping:
BLUE: 1
GREEN: 2
RED: 3
NIR: 4
SWIR_1: 5
needs extraction: False
opt_stack = opt_prod.stack([BLUE, GREEN, RED])
2023-05-31 12:10:39,401 - [DEBUG] - Loading bands ['BLUE', 'GREEN', 'RED']
2023-05-31 12:10:41,648 - [DEBUG] - Stacking
xr.plot.imshow(opt_stack.copy(data=display.scale(opt_stack.data)))
<matplotlib.image.AxesImage at 0x7f21d8d2ec50>

opt_stack
<xarray.DataArray 'BLUE_GREEN_RED' (bands: 3, y: 2237, x: 1244)> array([[[ nan, nan, nan, ..., 0.02729181, 0.03021449, 0.0321508 ], [ nan, nan, nan, ..., 0.03289769, 0.03252383, 0.03231718], [ nan, nan, nan, ..., 0.03253607, 0.03250813, 0.03260763], ..., [ nan, nan, nan, ..., nan, nan, nan], [ nan, nan, nan, ..., nan, nan, nan], [ nan, nan, nan, ..., nan, nan, nan]], [[ nan, nan, nan, ..., 0.0325688 , 0.03575394, 0.03786882], [ nan, nan, nan, ..., 0.03874811, 0.0377332 , 0.0372853 ], [ nan, nan, nan, ..., 0.03795209, 0.03785328, 0.03810363], ... [ nan, nan, nan, ..., nan, nan, nan], [ nan, nan, nan, ..., nan, nan, nan], [ nan, nan, nan, ..., nan, nan, nan]], [[ nan, nan, nan, ..., 0.02202989, 0.02403895, 0.02508134], [ nan, nan, nan, ..., 0.02564428, 0.02424301, 0.02346394], [ nan, nan, nan, ..., 0.0244639 , 0.02421321, 0.02448287], ..., [ nan, nan, nan, ..., nan, nan, nan], [ nan, nan, nan, ..., nan, nan, nan], [ nan, nan, nan, ..., nan, nan, nan]]], dtype=float32) Coordinates: * x (x) float64 3.044e+05 3.044e+05 ... 3.143e+05 3.143e+05 * y (y) float64 1.459e+06 1.459e+06 ... 1.441e+06 1.441e+06 spatial_ref int64 0 * bands (bands) object MultiIndex * variable (bands) object SpectralBandNames.BLUE ... SpectralBandNames.RED * band (bands) int64 1 1 1 Attributes: long_name: BLUE GREEN RED constellation: CUSTOM constellation_id: CUSTOM product_path: /home/data/DS3/CI/eoreader/others/20200310T030415_WV02... product_name: 20200310T030415_WV02_Ortho_BGRN_STK product_filename: 20200310T030415_WV02_Ortho_BGRN_STK instrument: CUSTOM product_type: CUSTOM acquisition_date: 20230531T121041 condensed_name: 20230531T121039_CUSTOM_CUSTOM orbit_direction: None
# SAR minimum example
sar_prod = reader.open(sar_path,
custom=True,
sensor_type=SensorType.SAR, # With the Enum
band_map={VV: 1, VV_DSPK: 2})
sar_prod
eoreader.CustomProduct '20210827T162210_ICEYE_SC_GRD_STK'
Attributes:
condensed_name: 20230531T121042_CUSTOM_CUSTOM
path: /home/data/DS3/CI/eoreader/others/20210827T162210_ICEYE_SC_GRD_STK.tif
constellation: CUSTOM
sensor type: SAR
product type: CUSTOM
default pixel size: 48.0
default resolution: None
acquisition datetime: 2023-05-31T12:10:42.653775
band mapping:
VV: 1
VV_DSPK: 2
needs extraction: False
sar_stack = sar_prod.stack([SLOPE, VV, VV_DSPK])
2023-05-31 12:10:42,664 - [DEBUG] - Loading bands ['VV', 'VV_DSPK']
2023-05-31 12:10:42,710 - [DEBUG] - Loading DEM bands ['SLOPE']
2023-05-31 12:10:42,711 - [DEBUG] - Warping DEM for 20230531T121042_CUSTOM_CUSTOM
2023-05-31 12:10:42,714 - [DEBUG] - Using DEM: /home/data/DS2/BASES_DE_DONNEES/GLOBAL/COPDEM_30m/COPDEM_30m.vrt
2023-05-31 12:10:43,541 - [DEBUG] - Computing slope for 20230531T121042_CUSTOM_CUSTOM
2023-05-31 12:10:53,534 - [DEBUG] - Stacking
xr.plot.imshow(sar_stack.copy(data=display.scale(sar_stack.data)))
<matplotlib.image.AxesImage at 0x7f21d8f34940>

sar_stack
<xarray.DataArray 'SLOPE_VV_VV_DSPK' (bands: 3, y: 2748, x: 2967)> array([[[1.1416357 , 0.96605515, 0.88842905, ..., 0. , 0. , 0. ], [0.91910994, 0.8988673 , 0.9167366 , ..., 0. , 0. , 0. ], [1.001919 , 0.84930086, 0.8696762 , ..., 0. , 0. , 0. ], ..., [0. , 0. , 0. , ..., 0. , 0. , 0. ], [0. , 0. , 0. , ..., 0. , 0. , 0. ], [0. , 0. , 0. , ..., 0. , 0. , 0. ]], [[ 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]]], dtype=float32) Coordinates: * x (x) float64 6.7e+05 6.701e+05 6.701e+05 ... 8.124e+05 8.124e+05 * y (y) float64 1.113e+04 1.109e+04 ... -1.206e+05 -1.207e+05 spatial_ref int64 0 * bands (bands) object MultiIndex * variable (bands) object DemBandNames.SLOPE ... SarBandNames.VV_DSPK * band (bands) int64 1 1 1 Attributes: long_name: SLOPE VV VV_DSPK constellation: CUSTOM constellation_id: CUSTOM product_path: /home/data/DS3/CI/eoreader/others/20210827T162210_ICEY... product_name: 20210827T162210_ICEYE_SC_GRD_STK product_filename: 20210827T162210_ICEYE_SC_GRD_STK instrument: CUSTOM product_type: CUSTOM acquisition_date: 20230531T121053 condensed_name: 20230531T121042_CUSTOM_CUSTOM orbit_direction: None
# You can compute the footprint and the extent
extent = opt_prod.extent()
footprint = opt_prod.footprint()
base = extent.plot(color='cyan', edgecolor='black')
footprint.plot(ax=base, color='blue', edgecolor='black', alpha=0.5)
<Axes: >

extent = sar_prod.extent()
footprint = sar_prod.footprint()
base = extent.plot(color='cyan', edgecolor='black')
footprint.plot(ax=base, color='blue', edgecolor='black', alpha=0.5)
<Axes: >

Custom stack with full data#
If you know them, it is best to give EOReader all the data you know about your stack:
name
: product name. If not provided, the filename will be useddatetime
: product acquisition datetime. If not provided, the datetime of the creation of the object will be usedconstellation
: product constellation. If not provided,CUSTOM
will be set. Either a string of a Constellation enum.product_type
: product type. If not provided,CUSTOM
will be set.pixel_size
: product default pixel size. If not provided, the stack pixel size will be used.
For optical products, two additional keyword can be set to compute the hillshade band:
sun_azimuth
sun_zenith
# Optical
opt_prod = reader.open(
opt_path,
custom=True,
name="20200310T030415_WV02_Ortho",
datetime="20200310T030415",
sensor_type=SensorType.OPTICAL,
constellation="WV02",
product_type="Ortho",
pixel_size=2.0,
sun_azimuth=10.0,
sun_zenith=20.0,
band_map={BLUE: 1, GREEN: 2, RED: 3, NIR: 4, SWIR_1: 5},
)
hillshade = opt_prod.load(HILLSHADE)[HILLSHADE]
2023-05-31 12:10:57,775 - [DEBUG] - Loading DEM bands ['HILLSHADE']
2023-05-31 12:10:57,776 - [DEBUG] - Warping DEM for 20200310T030415_WV02_Ortho
2023-05-31 12:10:57,778 - [DEBUG] - Using DEM: /home/data/DS2/BASES_DE_DONNEES/GLOBAL/COPDEM_30m/COPDEM_30m.vrt
2023-05-31 12:10:58,429 - [DEBUG] - Computing hillshade DEM for 20200310T030415_WV02_Ortho
hillshade.plot()
<matplotlib.collections.QuadMesh at 0x7f21d8fa1690>

hillshade
<xarray.DataArray <DemBandNames.HILLSHADE: 'HILLSHADE'> (band: 1, y: 8948, x: 4976)> array([[[247.08582, 247.24919, 247.46602, ..., 239.47937, 239.50517, 239.53088], [247.11534, 247.26117, 247.46373, ..., 239.47493, 239.50067, 239.52643], [247.14279, 247.27344, 247.46161, ..., 239.47063, 239.4964 , 239.52213], ..., [247.89543, 248.04077, 248.18419, ..., 239.4887 , 239.52728, 239.5658 ], [247.07147, 247.2222 , 247.37251, ..., 239.4823 , 239.52089, 239.55937], [246.18985, 246.34525, 246.5031 , ..., 239.47559, 239.51416, 239.55264]]], dtype=float32) Coordinates: * x (x) float64 3.044e+05 3.044e+05 ... 3.143e+05 3.143e+05 * y (y) float64 1.459e+06 1.459e+06 ... 1.441e+06 1.441e+06 * band (band) int64 1 spatial_ref int64 0 Attributes: long_name: HILLSHADE constellation: WorldView-2 constellation_id: WV02 product_path: /home/data/DS3/CI/eoreader/others/20200310T030415_WV02... product_name: 20200310T030415_WV02_Ortho product_filename: 20200310T030415_WV02_Ortho_BGRN_STK instrument: CUSTOM product_type: Ortho acquisition_date: 20200310T030415 condensed_name: 20200310T030415_WV02_Ortho orbit_direction: None
# SAR
sar_prod = reader.open(
sar_path,
custom=True,
sensor_type=SensorType.SAR,
name="20210827T162210_ICEYE_SC_GRD",
datetime="20210827T162210",
constellation="ICEYE",
product_type="GRD",
pixel_size=6.0,
band_map={VV: 1, VV_DSPK: 2},
)
from pprint import pprint
from eoreader import utils
# Read and display metadata
mtd, _ = sar_prod.read_mtd()
pprint(utils.quick_xml_to_dict(mtd))
('custom_metadata',
{'band_map': "{'VV': 1, 'VV_DSPK': 2}",
'cloud_cover': 'None',
'constellation': 'ICEYE',
'datetime': '2021-08-27T16:22:10',
'instrument': 'CUSTOM',
'name': '20210827T162210_ICEYE_SC_GRD',
'orbit_direction': 'None',
'pixel_size': '6.0',
'product_type': 'GRD',
'resolution': 'None',
'sensor_type': 'SAR',
'sun_azimuth': 'None',
'sun_zenith': 'None'})