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, DEM_VCRS, SNAP_DEM_NAME, SAR_DEF_PIXEL_SIZE
from eoreader.keywords import DEM_KW, SLOPE_KW, HILLSHADE_KW
from eoreader.products import SnapDems
# Create logger
import logging
from sertit import logs
logger = logging.getLogger("eoreader")
logs.init_logger(logger)
# logs.init_logger(logging.getLogger("sertit"))
# 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)
# 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 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: It hasn't been specified to EOReader
try:
prod.load(DEM)[DEM]
except ValueError as msg:
logger.error(msg)
2025-12-23 11:38:29,966 - [ERROR] - DEM path not set, unable to compute DEM bands! Please set the environment variable EOREADER_DEM_PATH or a DEM keyword.
# This works now
with tempenv.TemporaryEnvironment({DEM_PATH: dem}):
prod.load(DEM)[DEM]
2025-12-23 11:38:29,974 - [DEBUG] - Loading DEM bands ['DEM']
2025-12-23 11:38:29,975 - [DEBUG] - Warping DEM for 20200114T065229_S2_T40REQ_L2A_094749
2025-12-23 11:38:29,977 - [DEBUG] - Using DEM: /home/ds2_db2/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 has been loaded before and won't be recomputed
dem_bands = prod.load(
[DEM, SLOPE, HILLSHADE],
pixel_size=60,
**{
DEM_KW: dem,
SLOPE_KW: dtm, # We want a DTM here
HILLSHADE_KW: dem,
}
)
2025-12-23 11:38:35,870 - [DEBUG] - Loading DEM bands ['DEM', 'SLOPE', 'HILLSHADE']
2025-12-23 11:38:35,871 - [INFO] - Already existing DEM for /tmp/tmp24ydc89k/tmp_20200114T065229_S2_T40REQ_L2A_094749/20200114T065229_S2_T40REQ_L2A_094749_DEM_COPDEM_30m.vrt. Skipping process.
2025-12-23 11:38:39,326 - [DEBUG] - Warping DEM for 20200114T065229_S2_T40REQ_L2A_094749
2025-12-23 11:38:39,330 - [DEBUG] - Using DEM: /home/ds2_db2/BASES_DE_DONNEES/GLOBAL/MERIT_Hydrologically_Adjusted_Elevations/MERIT_DEM.vrt
2025-12-23 11:38:39,532 - [DEBUG] - Computing slope for 20200114T065229_S2_T40REQ_L2A_094749
/opt/conda/lib/python3.11/site-packages/rasterio/warp.py:387: NotGeoreferencedWarning: Dataset has no geotransform, gcps, or rpcs. The identity matrix will be returned.
dest = _reproject(
2025-12-23 11:38:42,458 - [INFO] - Already existing DEM for /tmp/tmp24ydc89k/tmp_20200114T065229_S2_T40REQ_L2A_094749/20200114T065229_S2_T40REQ_L2A_094749_DEM_COPDEM_30m.vrt. Skipping process.
2025-12-23 11:38:42,459 - [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, ::5, ::5].plot.imshow(
robust=True,
cmap="mako",
interpolation="gaussian",
interpolation_stage="data"
)
i += 1
Orthorectification#
SAR orthorectification (with SNAP)#
SNAP uses its own DEM mechanics: 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
The default resolution of the orthorectification computed by SNAP can be controlled through the EOREADER_SAR_DEFAULT_PIXEL_SIZE environment variable.
# Open a SAR product
path = os.path.join("/home", "prods", "ICEYE", "SC_124020")
prod = Reader().open(path)
prod
eoreader.IceyeProduct 'SC_124020'
Attributes:
condensed_name: 20210827T162210_ICEYE_VV_SC_GRD
path: /home/prods/ICEYE/SC_124020
constellation: ICEYE
sensor type: SAR
product type: GRD
default pixel size: 6.0
default resolution: 15.0
acquisition datetime: 2021-08-27T16:22:10.612000
band mapping:
VV: VV
VV_DSPK: VV_DSPK
needs extraction: True
orbit direction: DESCENDING
# Orthorectify the data at a coarser resolution to speed up process: Use 30 meters
os.environ[SAR_DEF_PIXEL_SIZE] = "30"
# Orthorectifying with COPDEM-30
ortho = prod.load(VV)[VV]
SNAP Release version 13.0.0
SNAP home: /opt/esa-snap/bin/..
SNAP debug: null
SNAP log level: WARNING
Java home: /opt/esa-snap/jre
Java version: 21.0.6
Processors: 16
Max memory: 29.0 GB
Cache size: 14.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
....11%....22%..
.
.
32%
.
.
.
.43%
.
.
.
.
54%.
.
.
.
65%
.
.
.
.
75%
.
.
.
.86%
.
.
done.
2025-12-23 11:39:36,849 - [DEBUG] - Loading bands ['VV']
SLF4J: Failed to load class "org.slf4j.impl.StaticLoggerBinder".
SLF4J: Defaulting to no-operation (NOP) logger implementation
SLF4J: See http://www.slf4j.org/codes.html#StaticLoggerBinder for further details.
SEVERE: org.esa.snap.core.gpf.operators.tooladapter.ToolAdapterIO: Failed to scan /opt/esa-snap/opttbx [reason: /opt/esa-snap/opttbx]
2025-12-23 11:39:38,185 - [DEBUG] - Pre-process SAR image
SLF4J: Failed to load class "org.slf4j.impl.StaticLoggerBinder".
SLF4J: Defaulting to no-operation (NOP) logger implementation
SLF4J: See http://www.slf4j.org/codes.html#StaticLoggerBinder for further details.
SEVERE: org.esa.snap.core.gpf.operators.tooladapter.ToolAdapterIO: Failed to scan /opt/esa-snap/opttbx [reason: /opt/esa-snap/opttbx]
2025-12-23 11:42:10,344 - [DEBUG] - Converting DIMAP to GeoTiff
2025-12-23 11:42:10,344 - [DEBUG] - Write SAR
2025-12-23 11:42:10,346 - [DEBUG] - Found [PosixPath('/tmp/tmpcuvdls4e/20210827T162210_ICEYE_VV_SC_GRD.data/Sigma0_VV_db.img')] sub-images (for VV).
2025-12-23 11:42:11,592 - [DEBUG] - SAR predictor: 3 (SNAP version: 13.0.0)
ortho[0, ::10, ::10].plot(cmap="mako", robust=True)
<matplotlib.collections.QuadMesh at 0x7f794c20cd10>
Orthorectifying with other DEM than default (GETASSE-30)#
You have to update properly the environment variable with the correct DEM name.
With EOReader’s variable, you can use:
with tempenv.TemporaryEnvironment(
{
SNAP_DEM_NAME: SnapDems.GETASSE30.value
}
):
prod.load(...)
Or directly with strings:
with tempenv.TemporaryEnvironment(
{
"EOREADER_SNAP_DEM_NAME": "GETASSE30"
}
):
prod.load(...)
Orthorectifying with custom DEM (EUDEM)#
You have to update properly the environment variables with the external DEM and its path:
with tempenv.TemporaryEnvironment(
{
SNAP_DEM_NAME: SnapDems.EXT_DEM.value,
DEM_PATH: f"{dem_folder}/EUDEM_v2/eudem_wgs84.tif"
}
):
prod.load(...)
or
with tempenv.TemporaryEnvironment(
{
SNAP_DEM_NAME: "External DEM",
DEM_PATH: f"{dem_folder}/EUDEM_v2/eudem_wgs84.tif"
}
):
prod.load(...)
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 Vantor (GeoEye and WorldView) products.
Warning: The orthorectification may take a while!
Warning: Orthorectification with RPCs should be done with an ellipsoid-based DEM!
EOReader is able to convert your put DEM if you provide the corresponding vertical CRS (passing a ‘vcrs’ keyword or put it in ‘EOREADER_DEM_VCRS’ environment variable).
More insights here: https://xdem.readthedocs.io/en/stable/vertical_ref.html
# Open a VHR product in sensor geometry
path = os.path.join("/home", "prods", "SPOT", "IMG_SPOT7_MS_001_A")
prod = Reader().open(path)
# Error when orthorectifying the RED band: the DEM hasn't been specified
try:
prod.load(RED)[RED]
except ValueError as msg:
logger.error(msg)
2025-12-23 11:42:53,677 - [DEBUG] - Loading bands ['RED']
2025-12-23 11:42:53,678 - [INFO] - Manually orthorectified stack not given by the user. Reprojecting whole stack, this may take a while. (Might be inaccurate on steep terrain, depending on the DEM pixel size).
2025-12-23 11:42:53,678 - [ERROR] - As you are using a non orthorectified product (/home/prods/SPOT/IMG_SPOT7_MS_001_A), you must provide a valid DEM through the EOREADER_DEM_PATH environment variable
# With environment variable
# COPDEM-30 is based on EGM08
# (however, EOReader is able to detect COPDEM vertical CRS thanks to the filename:
# it's looking for 'COPDEM' or 'Copernicus' in it)
with tempenv.TemporaryEnvironment({DEM_PATH: dem, DEM_VCRS: "EGM08"}):
red = prod.load(RED, chunks=None)[RED]
2025-12-23 11:42:53,684 - [DEBUG] - Loading bands ['RED']
2025-12-23 11:42:53,685 - [INFO] - Manually orthorectified stack not given by the user. Reprojecting whole stack, this may take a while. (Might be inaccurate on steep terrain, depending on the DEM pixel size).
Reprojection with RPCs doesn't work with Dask. Computing the raster.
2025-12-23 11:43:11,668 - [DEBUG] - Read RED
2025-12-23 11:43:14,183 - [DEBUG] - Manage nodata for band RED
2025-12-23 11:43:14,183 - [DEBUG] - Load nodata
2025-12-23 11:43:14,244 - [INFO] - Orthorectifying ROI
2025-12-23 11:43:14,244 - [DEBUG] - Rasterizing ROI
2025-12-23 11:43:36,392 - [DEBUG] - Reprojecting ROI
Reprojection with RPCs doesn't work with Dask. Computing the raster.
2025-12-23 11:44:01,180 - [DEBUG] - Revectorizing ROI
2025-12-23 11:44:05,162 - [DEBUG] - Rasterizing ROI mask
2025-12-23 11:44:06,132 - [DEBUG] - Set nodata mask
2025-12-23 11:44:06,221 - [DEBUG] - Converting RED to reflectance (if needed)
2025-12-23 11:44:06,370 - [DEBUG] - Clip the reflectance array to 0 as minimum value (in some cases, reflectance can have higher value than 1)
# With keyword
red = prod.load(RED, **{DEM_KW: dem})[RED]
red[0, ::10, ::10].plot(cmap="mako", robust=True)
<matplotlib.collections.QuadMesh at 0x7f7914e49810>
