DEM management#
DEM are used in two ways in EOReader:
for computing DEM bands (
DEM,SLOPE,HILLSHADE)for orthorectification (SAR, with SNAP or VHR data)
There are several things you need to know to handle DEMs properly.
# Imports
import os
import tempenv
import matplotlib.pyplot as plt
from eoreader.reader import Reader
from eoreader.bands import DEM, SLOPE, HILLSHADE, VV, RED
from eoreader.env_vars import DEM_PATH, SNAP_DEM_NAME
from eoreader.keywords import DEM_KW, SLOPE_KW, HILLSHADE_KW
from eoreader.products import SnapDems
/home/docs/checkouts/readthedocs.org/user_builds/eoreader/envs/stable/lib/python3.9/site-packages/dask/dataframe/__init__.py:42: FutureWarning:
Dask dataframe query planning is disabled because dask-expr is not installed.
You can install it with `pip install dask[dataframe]` or `conda install dask`.
This will raise in a future version.
warnings.warn(msg, FutureWarning)
# Create logger
import logging
from sertit import logs
logger = logging.getLogger("eoreader")
logs.init_logger(logger)
# Open the product
path = os.path.join("/home", "prods", "S2", "PB 02.07+", "S2B_MSIL2A_20200114T065229_N0213_R020_T40REQ_20200114T094749.SAFE")
prod = Reader().open(path)
2025-05-05 11:29:13,421 - [WARNING] - There is no existing products in EOReader corresponding to /home/prods/S2/PB 02.07+/S2B_MSIL2A_20200114T065229_N0213_R020_T40REQ_20200114T094749.SAFE.
2025-05-05 11:29:13,422 - [INFO] - Your given path may not be a satellite image. If it is, maybe the product isn't handled by EOReader. If you are sure this product is handled, it is either corrupted or you may need to go deeper in the filetree to find the correct path to give.
2025-05-05 11:29:13,423 - [DEBUG] - Please look at what folder you should give to EOReader by accessing the documentation: https://eoreader.readthedocs.io/latest/main_features.html#recognized-paths
# DEM paths
dem_folder = os.path.join("/home", "ds2_db2", "BASES_DE_DONNEES", "GLOBAL")
dtm = os.path.join(dem_folder, "MERIT_Hydrologically_Adjusted_Elevations", "MERIT_DEM.vrt")
dem = os.path.join(dem_folder, "COPDEM_30m", "COPDEM_30m.vrt")
dem_tif = os.path.join(dem_folder, "EUDEM_v2", "eudem_wgs84.tif")
DEM bands#
DEM bands are computed using the DEM stored in the EOREADER_DEM_PATH environment variable.
Without anything specified, DEM bands cannot be computed and an error will be thrown.
Note: SAR products cannot load a HILLSHADE band !
# Error when computing the DEM band: Iy hasn't been specified it
try:
prod.load(DEM)[DEM]
except ValueError as msg:
logger.error(msg)
---------------------------------------------------------------------------
AttributeError Traceback (most recent call last)
Cell In[5], line 3
1 # Error when computing the DEM band: Iy hasn't been specified it
2 try:
----> 3 prod.load(DEM)[DEM]
4 except ValueError as msg:
5 logger.error(msg)
AttributeError: 'NoneType' object has no attribute 'load'
# This works now
with tempenv.TemporaryEnvironment({DEM_PATH: dem}):
prod.load(DEM)[DEM]
2023-05-31 11:46:35,485 - [DEBUG] - Loading DEM bands ['DEM']
2023-05-31 11:46:35,486 - [DEBUG] - Warping DEM for 20200114T065229_S2_T40REQ_L2A_094749
2023-05-31 11:46:35,487 - [DEBUG] - Using DEM: /home/data/DS2/BASES_DE_DONNEES/GLOBAL/COPDEM_30m/COPDEM_30m.vrt
However, in some cases, you may want to specify a specific DEM to load a band from it.
The most common use of this is:
I want to compute the hillshade from a DEM
I want to compute the slope of a DTM as I don’t want the slopes of buildings or trees !
I want those two bands coregistered by default
When specifying a DEM in the load or stack function, you don’t need to set the EOREADER_DEM_PATH environment variable
dem_bands = prod.load(
[DEM, SLOPE, HILLSHADE],
**{
DEM_KW: dem,
SLOPE_KW: dtm, # We want a DTM here
HILLSHADE_KW: dem,
}
)
2023-05-31 11:49:08,486 - [DEBUG] - Loading DEM bands ['DEM', 'SLOPE', 'HILLSHADE']
2023-05-31 11:49:08,489 - [DEBUG] - Already existing DEM for 20200114T065229_S2_T40REQ_L2A_094749. Skipping process.
2023-05-31 11:49:10,335 - [DEBUG] - Warping DEM for 20200114T065229_S2_T40REQ_L2A_094749
2023-05-31 11:49:10,338 - [DEBUG] - Using DEM: /home/data/DS2/BASES_DE_DONNEES/GLOBAL/MERIT_Hydrologically_Adjusted_Elevations/MERIT_DEM.vrt
2023-05-31 11:49:10,599 - [DEBUG] - Computing slope for 20200114T065229_S2_T40REQ_L2A_094749
2023-05-31 11:56:15,503 - [DEBUG] - Already existing DEM for 20200114T065229_S2_T40REQ_L2A_094749. Skipping process.
2023-05-31 11:56:15,513 - [DEBUG] - Computing hillshade DEM for S2B_MSIL2A_20200114T065229_N0213_R020_T40REQ_20200114T094749
# Plot DEM bands
ncols = len(dem_bands)
plt.figure(figsize=(6, 6 * ncols))
i = 0
for key in dem_bands.keys():
axes = plt.subplot(3, 1, i+1)
dem_bands[key][0, ::10, ::10].plot.imshow(robust=True)
i += 1
Orthorectification#
SAR orthorectification (with SNAP)#
SNAP uses its own DEM mecanics: it downloads DEM tiles on its own and everything is managed without any intervention of the user.
By default, EOReader forces SNAP to use the COPDEM-30 (GLO-30) as its a global DTM.
However, this can be modified by setting the environment variable overriding the default DEM name used in SNAP.
The user can choose GETASSE30, SRTM 3Sec, External DEM…
If EOREADER_SNAP_DEM_NAME is set to External DEM, SNAP will use the DEM stored in EOREADER_DEM_PATH as an external DEM.
Warning: It seems that SNAP can only DEM with TIF format and not VRT.
Using such DEM format would result with the following error message: Error: No product reader found for /path/to/dem.vrt
# Open a SAR product
path = os.path.join("/home", "prods", "S1", "S1B_IW_GRDH_1SDV_20191215T180300_20191215T180325_019379_0249B2_C99C.SAFE")
prod = Reader().open(path)
# Orthorectifying with COPDEM-30
vv = prod.load(VV, pixel_size=prod.pixel_size*100)[VV]
vv.plot()
SNAP Release version 9.0.0
SNAP home: /opt/snap/bin/..
SNAP debug: null
SNAP log level: WARNING
Java home: /opt/snap/jre/jre
Java version: 1.8.0_242
Processors: 16
Max memory: 40.9 GB
Cache size: 23.0 GB
Tile parallelism: 14
Tile size: 512 x 512 pixels
To configure your gpt memory usage:
Edit snap/bin/gpt.vmoptions
To configure your gpt cache size and parallelism:
Edit .snap/etc/snap.properties or gpt -c ${cachesize-in-GB}G -q ${parallelism}
Executing processing graph
OpenSearch: https://scihub.copernicus.eu/gnss/search?q=platformname:Sentinel-1 AND platformnumber:B AND producttype:AUX_POEORB AND beginposition:[2019-12-01T00:00:000Z TO 2019-12-31T24:00:000Z]
OpenSearch: 31 total results on 1 pages.
OpenSearch: https://scihub.copernicus.eu/gnss/search?q=platformname:Sentinel-1 AND platformnumber:B AND producttype:AUX_POEORB AND beginposition:[2019-12-01T00:00:000Z TO 2019-12-31T24:00:000Z]
version = 3.1
10%20%30%40%50%60%70%.80%90%Copernicus_DSM_COG_10_N44_00_W005_00_DEM.tif
Copernicus_DSM_COG_10_N44_00_W004_00_DEM.tif
Copernicus_DSM_COG_10_N44_00_W003_00_DEM.tif
Copernicus_DSM_COG_10_N44_00_W002_00_DEM.tif
Copernicus_DSM_COG_10_N44_00_W001_00_DEM.tif
Copernicus_DSM_COG_10_N43_00_W001_00_DEM.tif
Copernicus_DSM_COG_10_N43_00_W005_00_DEM.tif
Copernicus_DSM_COG_10_N43_00_W004_00_DEM.tif
Copernicus_DSM_COG_10_N43_00_W003_00_DEM.tif
Copernicus_DSM_COG_10_N43_00_W002_00_DEM.tif
Copernicus_DSM_COG_10_N42_00_W001_00_DEM.tif
Copernicus_DSM_COG_10_N42_00_W002_00_DEM.tif
Copernicus_DSM_COG_10_N42_00_W005_00_DEM.tif
Copernicus_DSM_COG_10_N42_00_W004_00_DEM.tif
Copernicus_DSM_COG_10_N42_00_W003_00_DEM.tif
2023-05-31 11:58:26,155 - [DEBUG] - Loading bands ['VV']
2023-05-31 11:58:29,995 - [DEBUG] - Pre-process SAR image
Killed
---------------------------------------------------------------------------
FileNotFoundError Traceback (most recent call last)
File /home/data/eoreader/eoreader/products/sar/sar_product.py:432, in SarProduct.get_band_paths(self, band_list, pixel_size, **kwargs)
431 band_id = self.bands[band].id
--> 432 band_paths[band] = files.get_file_in_dir(
433 self._get_band_folder(),
434 f"*{self.condensed_name}_{band_id}.tif",
435 exact_name=True,
436 )
437 except FileNotFoundError:
File /opt/conda/lib/python3.10/site-packages/sertit/files.py:1146, in get_file_in_dir(directory, pattern_str, extension, filename_only, get_list, exact_name)
1145 if len(file_list) == 0:
-> 1146 raise FileNotFoundError(
1147 f"File with pattern {glob_pattern} not found in {directory}"
1148 )
1150 # Return list, file path or file name
FileNotFoundError: File with pattern *20191215T180300_S1_VV_VH_IW_GRD_VV.tif not found in /tmp/tmp4i911y3p/tmp_20191215T180300_S1_VV_VH_IW_GRD
During handling of the above exception, another exception occurred:
RuntimeError Traceback (most recent call last)
File /home/data/eoreader/eoreader/products/sar/sar_product.py:817, in SarProduct._pre_process_sar(self, band, pixel_size, **kwargs)
816 try:
--> 817 misc.run_cli(cmd_list)
818 except RuntimeError as ex:
File /opt/conda/lib/python3.10/site-packages/sertit/misc.py:370, in run_cli(cmd, timeout, check_return_value, in_background, cwd)
369 if check_return_value and retval != 0:
--> 370 raise RuntimeError(f"Exe {cmd[0]} has failed.")
372 return retval, output
RuntimeError: Exe gpt has failed.
The above exception was the direct cause of the following exception:
RuntimeError Traceback (most recent call last)
Cell In[10], line 2
1 # Orthorectifying with COPDEM-30
----> 2 vv = prod.load(VV, pixel_size=prod.pixel_size*100)[VV]
3 vv.plot()
File /home/data/eoreader/eoreader/products/product.py:876, in Product.load(self, bands, pixel_size, size, **kwargs)
874 # Load bands (only once ! and convert the bands to be loaded to correct format)
875 unique_bands = misc.unique(bands)
--> 876 band_xds = self._load(unique_bands, pixel_size, size, **kwargs)
878 # Rename all bands and add attributes
879 for key, val in band_xds.items():
File /home/data/eoreader/eoreader/products/product.py:974, in Product._load(self, bands, pixel_size, size, **kwargs)
972 if unique_bands:
973 LOGGER.debug(f"Loading bands {to_str(unique_bands)}")
--> 974 loaded_bands = self._load_bands(
975 unique_bands, pixel_size=pixel_size, size=size, **kwargs
976 )
978 # Compute index (they conserve the nodata)
979 if index_list:
980 # Collocate bands before indices to ensure the same size to perform operations between bands
File /home/data/eoreader/eoreader/products/sar/sar_product.py:642, in SarProduct._load_bands(self, bands, pixel_size, size, **kwargs)
640 if pixel_size is None and size is not None:
641 pixel_size = self._pixel_size_from_img_size(size)
--> 642 band_paths = self.get_band_paths(bands, pixel_size, **kwargs)
644 # Open bands and get array (resampled if needed)
645 band_arrays = {}
File /home/data/eoreader/eoreader/products/sar/sar_product.py:454, in SarProduct.get_band_paths(self, band_list, pixel_size, **kwargs)
452 band_paths[band] = self._despeckle_sar(speckle_band, **kwargs)
453 else:
--> 454 band_paths[band] = self._pre_process_sar(
455 band, pixel_size, **kwargs
456 )
458 return band_paths
File /home/data/eoreader/eoreader/products/sar/sar_product.py:819, in SarProduct._pre_process_sar(self, band, pixel_size, **kwargs)
817 misc.run_cli(cmd_list)
818 except RuntimeError as ex:
--> 819 raise RuntimeError("Something went wrong with SNAP!") from ex
821 # Convert DIMAP images to GeoTiff
822 LOGGER.debug("Converting DIMAP to GeoTiff")
RuntimeError: Something went wrong with SNAP!
# Orthorectifying with GETASSE
prod.clean_tmp()
with tempenv.TemporaryEnvironment(
{
SNAP_DEM_NAME: SnapDems.GETASSE30.value
}
):
vv2 = prod.load(VV, pixel_size=prod.pixel_size*100)[VV]
vv2.plot()
# Orthorectifying with custom TIF DEM (EUDEM)
with tempenv.TemporaryEnvironment(
{
SNAP_DEM_NAME: SnapDems.EXT_DEM.value,
DEM_PATH: dem_tif
}
):
vv = prod.load(VV, pixel_size=prod.pixel_size*100)[VV]
VHR orthorectification#
Some products are not orthorectified by default (according what product type you bought).
It may be wanted, in order to control yourself the orthorectification.
EOReader can orthorectify VHR products: DIMAP (Pleiades and SPOT) and Maxar (GeoEye and WorldView) products.
Warning: The orthorectification may take a while !
# Open a VHR product in sensor geometry
path = os.path.join("/home", "prods", "PLEIADES", "3302499201", "IMG_PHR1A_MS_004")
prod = Reader().open(path)
# Error when orthorectifying the RED band: Iy hasn't been specified it
try:
prod.load(RED)[RED]
except ValueError as msg:
logger.error(msg)
# With environment variable
with tempenv.TemporaryEnvironment({DEM_PATH: dem}):
red = prod.load(RED)[RED]
red.plot()
# With keyword
red = prod.load(RED, **{DEM_KW: dem})[RED]