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 logging
import tempenv

import matplotlib.pyplot as plt

from sertit import logs

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 

# Create logger
logger = logging.getLogger("eoreader")
logs.init_logger(logger)
# Open the product
path = os.path.join("/home", "data", "DATA", "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", "data", "DS2", "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)     
2022-10-10 14:20:10,895 - [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]
2022-10-10 14:20:10,903 - [DEBUG] - Loading DEM bands ['DEM']
2022-10-10 14:20:10,903 - [DEBUG] - Warping DEM for 20200114T065229_S2_T40REQ_L2A_094749
2022-10-10 14:20:10,906 - [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,        
    }
)
2022-10-10 14:20:16,154 - [DEBUG] - Loading DEM bands ['SLOPE', 'DEM', 'HILLSHADE']
2022-10-10 14:20:16,155 - [DEBUG] - Warping DEM for 20200114T065229_S2_T40REQ_L2A_094749
2022-10-10 14:20:16,158 - [DEBUG] - Using DEM: /home/data/DS2/BASES_DE_DONNEES/GLOBAL/MERIT_Hydrologically_Adjusted_Elevations/MERIT_DEM.vrt
2022-10-10 14:20:17,640 - [DEBUG] - Computing slope for 20200114T065229_S2_T40REQ_L2A_094749
2022-10-10 14:20:23,971 - [DEBUG] - Already existing DEM for 20200114T065229_S2_T40REQ_L2A_094749. Skipping process.
2022-10-10 14:20:24,276 - [DEBUG] - Already existing DEM for 20200114T065229_S2_T40REQ_L2A_094749. Skipping process.
2022-10-10 14:20:24,276 - [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
../_images/4d1e9c26293d8a9500596f2066cdec7bf6085f44489e2820f36e2ed784e5e453.png

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", "data", "DATA", "PRODS", "S1", "S1B_IW_GRDH_1SDV_20191215T180300_20191215T180325_019379_0249B2_C99C.SAFE")
prod = Reader().open(path)
# Orthorectifying with COPDEM-30
vv = prod.load(VV, resolution=prod.resolution*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: 2048 x 2048 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
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
 done.
2022-10-10 14:20:35,715 - [DEBUG] - Pre-process SAR image
2022-10-10 14:22:16,166 - [DEBUG] - Converting DIMAP to GeoTiff
<matplotlib.collections.QuadMesh at 0x7f3ff419fa00>
../_images/162da9f06da2e417b80ffd6e5d3ac8d9eff831c822b95e02772091d2782a99e0.png
# Orthorectifying with GETASSE
prod.clean_tmp()
with tempenv.TemporaryEnvironment(
    {
        SNAP_DEM_NAME: SnapDems.GETASSE30.value
    }
):
    vv2 = prod.load(VV, resolution=prod.resolution*100)[VV]
    vv2.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: 2048 x 2048 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
version = 3.1
10%20%30%40%50%60%70%.80%90% done.
2022-10-10 14:22:18,784 - [DEBUG] - Pre-process SAR image
2022-10-10 14:23:26,052 - [DEBUG] - Converting DIMAP to GeoTiff
../_images/1a1e51d87c8333e699604dffb9cfb4fe50e8864969fab2db9beaafa228bfb114.png
# Orthorectifying with custom TIF DEM (EUDEM)
with tempenv.TemporaryEnvironment(
    {
        SNAP_DEM_NAME: SnapDems.EXT_DEM.value,
        DEM_PATH: dem_tif
    }
):
    vv = prod.load(VV, resolution=prod.resolution*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", "data", "DATA", "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)
2022-10-10 14:23:27,082 - [DEBUG] - Loading bands ['RED']
2022-10-10 14:23:27,083 - [INFO] - Manually orthorectified stack not given by the user. Reprojecting whole stack, this may take a while. (May be inaccurate on steep terrain, depending on the DEM resolution)
2022-10-10 14:23:27,083 - [ERROR] - As you are using a non orthorectified VHR product (/home/data/DATA/PRODS/PLEIADES/3302499201/IMG_PHR1A_MS_004), you must provide a valid DEM through the EOREADER_DEM_PATH environment variable
# With environment variable
with tempenv.TemporaryEnvironment({DEM_PATH: dem}):
    red = prod.load(RED)[RED]

red.plot()
2022-10-10 14:23:27,089 - [DEBUG] - Loading bands ['RED']
2022-10-10 14:23:27,090 - [INFO] - Manually orthorectified stack not given by the user. Reprojecting whole stack, this may take a while. (May be inaccurate on steep terrain, depending on the DEM resolution)
2022-10-10 14:23:28,232 - [DEBUG] - Orthorectifying data with /home/data/DS2/BASES_DE_DONNEES/GLOBAL/COPDEM_30m/COPDEM_30m.vrt
2022-10-10 14:25:05,877 - [DEBUG] - Read RED
2022-10-10 14:25:06,890 - [DEBUG] - Converting RED to reflectance
2022-10-10 14:25:06,984 - [DEBUG] - Manage nodata for band RED
ERROR 1: Cannot import urn:ogc:def:derivedCRSType:OGC:1.0:image due to ALLOW_FILE_ACCESS=NO
2022-10-10 14:25:07,103 - [INFO] - Orthorectifying ROI
2022-10-10 14:25:07,164 - [DEBUG] - 	Rasterizing ROI
2022-10-10 14:25:07,204 - [DEBUG] - 	Reprojecting ROI
2022-10-10 14:25:07,205 - [DEBUG] - Orthorectifying data with /home/data/DS2/BASES_DE_DONNEES/GLOBAL/COPDEM_30m/COPDEM_30m.vrt
2022-10-10 14:26:33,052 - [DEBUG] - 	Revectorizing ROI
<matplotlib.collections.QuadMesh at 0x7f3fdc762b90>
../_images/b0bbe6b1644716c27f1f13805a0e3e05b571997fb930161578cede3b3756f8c5.png
# With keyword
red = prod.load(RED, **{DEM_KW: dem})[RED]