Custom stacks
Contents
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
logs.init_logger(logging.getLogger("eoreader"))
# 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
---------------------------------------------------------------------------
CPLE_OpenFailedError Traceback (most recent call last)
File rasterio/_base.pyx:307, in rasterio._base.DatasetBase.__init__()
File rasterio/_base.pyx:218, in rasterio._base.open_dataset()
File rasterio/_err.pyx:221, in rasterio._err.exc_wrap_pointer()
CPLE_OpenFailedError: /home/data/DS3/CI/eoreader/others/20200310T030415_WV02_Ortho_BGRN_STK.tif: No such file or directory
During handling of the above exception, another exception occurred:
RasterioIOError Traceback (most recent call last)
Cell In[5], line 2
1 # Optical minimum example
----> 2 opt_prod = reader.open(opt_path,
3 custom=True,
4 sensor_type="OPTICAL", # With a string
5 band_map={BLUE: 1, GREEN: 2, RED: 3, NIR: 4, SWIR_1: 5})
6 opt_prod
File ~/checkouts/readthedocs.org/user_builds/eoreader/checkouts/latest/eoreader/reader.py:467, in Reader.open(self, product_path, archive_path, output_path, method, remove_tmp, custom, constellation, **kwargs)
464 if custom:
465 from eoreader.products import CustomProduct
--> 467 prod = CustomProduct(
468 product_path=product_path,
469 archive_path=archive_path,
470 output_path=output_path,
471 remove_tmp=remove_tmp,
472 constellation=constellation,
473 **kwargs,
474 )
475 else:
476 prod = None
File ~/checkouts/readthedocs.org/user_builds/eoreader/checkouts/latest/eoreader/products/custom_product.py:107, in CustomProduct.__init__(self, product_path, archive_path, output_path, remove_tmp, **kwargs)
105 super_kwargs = kwargs.copy()
106 super_kwargs.pop("constellation", None)
--> 107 super().__init__(
108 product_path, archive_path, output_path, remove_tmp, **super_kwargs
109 )
File ~/checkouts/readthedocs.org/user_builds/eoreader/checkouts/latest/eoreader/products/product.py:244, in Product.__init__(self, product_path, archive_path, output_path, remove_tmp, **kwargs)
241 self._set_instrument()
243 # Post initialization
--> 244 self._post_init(**kwargs)
246 # Set product type, needs to be done after the post-initialization
247 self._set_product_type()
File ~/checkouts/readthedocs.org/user_builds/eoreader/checkouts/latest/eoreader/products/custom_product.py:177, in CustomProduct._post_init(self, **kwargs)
175 # Check CRS
176 try:
--> 177 crs = self.crs() # noqa
178 except InvalidProductError as msg:
179 LOGGER.warning(msg)
File ~/checkouts/readthedocs.org/user_builds/eoreader/envs/latest/lib/python3.8/site-packages/methodtools.py:72, in _LruCacheWire.__call__(self, *args, **kwargs)
70 def __call__(self, *args, **kwargs):
71 # descriptor detection support - never called
---> 72 return self.__call__(*args, **kwargs)
File ~/checkouts/readthedocs.org/user_builds/eoreader/checkouts/latest/eoreader/__init__.py:34, in cache.<locals>.wrapper(*args, **kwargs)
31 @lru_cache()
32 @wraps(func)
33 def wrapper(*args, **kwargs):
---> 34 return func(*args, **kwargs)
File ~/checkouts/readthedocs.org/user_builds/eoreader/checkouts/latest/eoreader/products/custom_product.py:310, in CustomProduct.crs(self)
302 @cache
303 def crs(self) -> crs.CRS:
304 """
305 Get UTM projection of stack.
306
307 Returns:
308 crs.CRS: CRS object
309 """
--> 310 with rasterio.open(str(self.path)) as ds:
311 def_crs = ds.crs
313 if def_crs.is_projected:
File ~/checkouts/readthedocs.org/user_builds/eoreader/envs/latest/lib/python3.8/site-packages/rasterio/env.py:451, in ensure_env_with_credentials.<locals>.wrapper(*args, **kwds)
448 session = DummySession()
450 with env_ctor(session=session):
--> 451 return f(*args, **kwds)
File ~/checkouts/readthedocs.org/user_builds/eoreader/envs/latest/lib/python3.8/site-packages/rasterio/__init__.py:304, in open(fp, mode, driver, width, height, count, crs, transform, dtype, nodata, sharing, **kwargs)
301 path = _parse_path(raw_dataset_path)
303 if mode == "r":
--> 304 dataset = DatasetReader(path, driver=driver, sharing=sharing, **kwargs)
305 elif mode == "r+":
306 dataset = get_writer_for_path(path, driver=driver)(
307 path, mode, driver=driver, sharing=sharing, **kwargs
308 )
File rasterio/_base.pyx:309, in rasterio._base.DatasetBase.__init__()
RasterioIOError: /home/data/DS3/CI/eoreader/others/20200310T030415_WV02_Ortho_BGRN_STK.tif: No such file or directory
opt_stack = opt_prod.stack([BLUE, GREEN, RED])
xr.plot.imshow(opt_stack.copy(data=display.scale(opt_stack.data)))
<matplotlib.image.AxesImage at 0x7f3f735a6be0>

opt_stack
<xarray.DataArray 'BLUE GREEN RED' (z: 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: spatial_ref int64 0 * 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 * z (z) MultiIndex - variable (z) object 'BLUE' 'GREEN' 'RED' - band (z) int64 1 1 1 Attributes: long_name: BLUE GREEN RED sensor: CUSTOM sensor_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 product_type: CUSTOM acquisition_date: 20220111T142955 condensed_name: 20220111T142955_CUSTOM_CUSTOM
# 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
2022-01-11 14:29:57,786 - [DEBUG] - Warping DEM for 20220111T142956_CUSTOM_CUSTOM
2022-01-11 14:29:57,791 - [DEBUG] - Using DEM: /home/data/DS2/BASES_DE_DONNEES/GLOBAL/COPDEM_30m/COPDEM_30m.vrt
2022-01-11 14:30:00,399 - [DEBUG] - Computing slope for 20220111T142956_CUSTOM_CUSTOM
sar_stack = sar_prod.stack([SLOPE, VV, VV_DSPK])
xr.plot.imshow(sar_stack.copy(data=display.scale(sar_stack.data)))
<matplotlib.image.AxesImage at 0x7f3f6cbd9940>

sar_stack
<xarray.DataArray 'SLOPE VV VV_DSPK' (z: 3, y: 2748, x: 2967)> array([[[1.1417845 , 0.9661645 , 0.88848215, ..., 0. , 0. , 0. ], [0.91908467, 0.8988768 , 0.9166924 , ..., 0. , 0. , 0. ], [1.0019214 , 0.84933126, 0.86957526, ..., 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: spatial_ref int64 0 * 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 * z (z) MultiIndex - variable (z) object 'SLOPE' 'VV' 'VV_DSPK' - band (z) int64 1 1 1 Attributes: long_name: SLOPE VV VV_DSPK sensor: CUSTOM sensor_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 product_type: CUSTOM acquisition_date: 20220111T142956 condensed_name: 20220111T142956_CUSTOM_CUSTOM
# 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)
<AxesSubplot:>

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)
<AxesSubplot:>

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 usedacquisition_datetime
: product acquisition datetime. If not provided, the datetime of the creation of the object will be usedplatform
: product platform. If not provided,CUSTOM
will be set. Either a string of a Platform enum.product_type
: product type. If not provided,CUSTOM
will be set.default_resolution
: product default resolution. If not provided, the stack resolution 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",
acquisition_datetime="20200310T030415",
sensor_type=SensorType.OPTICAL,
platform="WV02",
product_type="Ortho",
default_resolution=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]
2022-01-11 14:30:03,394 - [DEBUG] - Warping DEM for 20200310T030415_WV02_Ortho
2022-01-11 14:30:03,397 - [DEBUG] - Using DEM: /home/data/DS2/BASES_DE_DONNEES/GLOBAL/COPDEM_30m/COPDEM_30m.vrt
2022-01-11 14:30:03,953 - [DEBUG] - Computing hillshade DEM for 20200310T030415_WV02_Ortho
hillshade.plot()
<matplotlib.collections.QuadMesh at 0x7f3f507e9640>

hillshade
<xarray.DataArray 'HILLSHADE' (band: 1, y: 2237, x: 1244)> array([[[243.5148 , 244.15515, 244.76294, ..., 239.26723, 239.37088, 239.47406], [241.75117, 242.62433, 243.46106, ..., 239.39891, 239.43938, 239.48117], [239.7554 , 240.88208, 241.96498, ..., 239.66376, 239.59615, 239.53004], ..., [247.74507, 248.3067 , 246.9472 , ..., 239.20195, 239.27254, 239.39708], [247.93413, 248.49562, 246.84186, ..., 239.26035, 239.3623 , 239.51692], [248.09465, 248.6562 , 246.72041, ..., 239.2324 , 239.33617, 239.4908 ]]], dtype=float32) Coordinates: * band (band) int64 1 * 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 Attributes: scale_factor: 1.0 add_offset: 0.0 long_name: HILLSHADE sensor: WorldView-2 sensor_id: WV02 product_path: /home/data/DS3/CI/eoreader/others/20200310T030415_WV02... product_name: 20200310T030415_WV02_Ortho product_filename: 20200310T030415_WV02_Ortho_BGRN_STK product_type: Ortho acquisition_date: 20200310T030415 condensed_name: 20200310T030415_WV02_Ortho
# SAR
sar_prod = reader.open(
sar_path,
custom=True,
sensor_type=SensorType.SAR,
name="20210827T162210_ICEYE_SC_GRD",
acquisition_datetime="20210827T162210",
platform="ICEYE",
product_type="GRD",
default_resolution=6.0,
band_map={VV: 1, VV_DSPK: 2},
)
{<SarBandNames.VV: 'VV'>: <xarray.DataArray 'VV' (band: 1, y: 2748, x: 2967)>
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 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
Attributes:
scale_factor: 1.0
add_offset: 0.0
long_name: VV
sensor: ICEYE
sensor_id: ICEYE
product_path: /home/data/DS3/CI/eoreader/others/20210827T162210_ICEY...
product_name: 20210827T162210_ICEYE_SC_GRD
product_filename: 20210827T162210_ICEYE_SC_GRD_STK
product_type: GRD
acquisition_date: 20210827T162210
condensed_name: 20210827T162210_ICEYE_GRD}
from pprint import pprint
from eoreader import utils
# Read and display metadata
mtd, _ = sar_prod.read_mtd()
pprint(utils.quick_xml_to_dict(mtd))