DEM management
Contents
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.alias import *
from eoreader.products import SnapDems
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.00", "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-02-01 09:55:16,672 - [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-02-01 09:55:16,718 - [DEBUG] - Loading DEM bands ['DEM']
2022-02-01 09:55:16,718 - [DEBUG] - Warping DEM for 20200114T065229_S2_T40REQ_L2A_094749
2022-02-01 09:55:16,721 - [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-02-01 09:55:22,001 - [DEBUG] - Loading DEM bands ['SLOPE', 'DEM', 'HILLSHADE']
2022-02-01 09:55:22,001 - [DEBUG] - Warping DEM for 20200114T065229_S2_T40REQ_L2A_094749
2022-02-01 09:55:22,005 - [DEBUG] - Using DEM: /home/data/DS2/BASES_DE_DONNEES/GLOBAL/MERIT_Hydrologically_Adjusted_Elevations/MERIT_DEM.vrt
2022-02-01 09:55:22,985 - [DEBUG] - Computing slope for 20200114T065229_S2_T40REQ_L2A_094749
2022-02-01 09:55:24,639 - [DEBUG] - Already existing DEM for 20200114T065229_S2_T40REQ_L2A_094749. Skipping process.
2022-02-01 09:55:24,754 - [DEBUG] - Already existing DEM for 20200114T065229_S2_T40REQ_L2A_094749. Skipping process.
2022-02-01 09:55:24,754 - [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", "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()
2022-02-01 09:55:30,091 - [DEBUG] - Pre-process SAR image
SNAP Release version 8.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
Copernicus_DSM_COG_10_N44_00_W002_00_DEM.tif
Copernicus_DSM_COG_10_N43_00_W002_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_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
2022-02-01 09:57:25,072 - [DEBUG] - Converting DIMAP to GeoTiff
done.
<matplotlib.collections.QuadMesh at 0x7fefd7f402b0>
# 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()
2022-02-01 09:57:30,718 - [DEBUG] - Pre-process SAR image
SNAP Release version 8.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%
2022-02-01 09:58:37,867 - [DEBUG] - Converting DIMAP to GeoTiff
done.
# 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-02-01 09:58:38,891 - [DEBUG] - Loading bands ['RED']
2022-02-01 09:58:38,892 - [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-02-01 09:58:38,893 - [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-02-01 09:58:38,901 - [DEBUG] - Loading bands ['RED']
2022-02-01 09:58:38,902 - [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-02-01 09:58:40,300 - [DEBUG] - Orthorectifying data with /home/data/DS2/BASES_DE_DONNEES/GLOBAL/COPDEM_30m/COPDEM_30m.vrt
2022-02-01 10:00:22,297 - [DEBUG] - Read RED
2022-02-01 10:00:23,416 - [DEBUG] - Manage nodata for band RED
2022-02-01 10:00:23,516 - [INFO] - Orthorectifying ROI
2022-02-01 10:00:23,557 - [DEBUG] - Rasterizing ROI
2022-02-01 10:00:23,592 - [DEBUG] - Reprojecting ROI
2022-02-01 10:00:23,593 - [DEBUG] - Orthorectifying data with /home/data/DS2/BASES_DE_DONNEES/GLOBAL/COPDEM_30m/COPDEM_30m.vrt
2022-02-01 10:01:56,981 - [DEBUG] - Revectorizing ROI
/opt/conda/lib/python3.9/site-packages/geopandas/io/file.py:362: FutureWarning: pandas.Int64Index is deprecated and will be removed from pandas in a future version. Use pandas.Index with the appropriate dtype instead.
pd.Int64Index,
<matplotlib.collections.QuadMesh at 0x7fefd76a18b0>
# With keyword
red = prod.load(RED, **{DEM_KW: dem})[RED]