Using an external DEM

[EJFielding]

We (Tonian Robinson, @tonianr, and I) were having a discussion under a closed Issue (#329 (comment)) about how to prepare an external DEM for use in ISCE2. It was suggested that we have the discussion here in the newly enabled Discussion section.

I suggested that there are ways to use a GeoTIFF file as a DEM in ISCE2 if you create the necessary .xml and .vrt files, but you are better off converting it to a standard ISCE raster format with the GDAL gdal_translate -of ISCE command.

The other big consideration is getting an external DEM into the horizontal and vertical datum that ISCE2 expects: WGS84 horizontal and WGS84 ellipsoid vertical datums.

[tonianr]

Steps I used to transform data from NAVD88 vertical datum to WGS84 Ellipsoid vertical datum (necessary for ISCE2):

From NAVD88 to WGS84 Ellipsoid
gdalwarp -overwrite Hernando_2007_FDEM_WGS84_NAVD88.tif Hernando_2007_FDEM_WGS84_WGS84.tif -s_srs EPSG:3857+5703 -t_srs EPSG:4979

Since my input is projected to WGS84 Ancillary (EPSG 3857) with vertical datum NAVD88 (EPSG 5703) and the desired output for ISCE2 is EPSG 4979 for WGS84 horizontal and vertical height

Tif to ISCE format using GDAL
gdal_translate -of ISCE Hernando_2007_FDEM_WGS84_WGS84.tif Hernando_07_FDEM_WGS84_WGS84.dem.wgs84

Create ISCE format .vrt
import isce
from applications.gdal2isce_xml import gdal2isce_xml
fname = 'Hernando_07_FDEM_WGS84_WGS84.dem.wgs84'
xml_file = gdal2isce_xml(fname)

[EJFielding]

The gdalwarp -s_srs EPSG:3857+5703 -t_srs EPSG:4979 will only transform the horizontal position of the DEM. It won't change the vertical datum. I thought this would work, but I tried it a while ago and found that the vertical datum was not changed.
You need to use one of those Geoid correction files like the one you used in the posting in the "Issues" thread to do the vertical datum conversion.

[vincentschut]

Indeed. Chiming in here, because I recently had to do something like that: converting Copernicus dem (vertical datum: egm08) to wgs84 ellipsoidal. What appears to work for me is this:

gdalwarp -s_srs "+proj=longlat +datum=WGS84 +no_defs +geoidgrids=/usr/share/proj/egm08_25.gtx" -t_srs "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs" <dem_egm08.tif> <dem_wgs84.tif>

When comparing the input and output tif in qgis, the dem's values have changed, so apparently this is working.

Maybe you can use something like this, but then with the NAVD88 vertical datum shift file?

[tonianr]

Thank you, I'll definitely try that I stumbled on Proj Datum Grid which stated that the EGM2008 model " transforms physical heights to WGS84 ellipsoidal heights" and I was confused as to how to use it for the conversion.

[EJFielding]

For the Copernicus GLO-30 DEM Geoid EGM2008 correction, there is also this software developed by the JPL ARIA project:
https://github.com/aria-jpl/dem_stitcher

[tonianr]

Here's how I finally got to convert my NAD83 NADV88 datum LiDAR data to a WGS84 ref ellipsoid GeoTIFF file for ISCE stack processing.
Note I am using the LiDAR to create a DSM.

1. Downloaded and used NOAA vdatum to convert the .las files to WGS84 UTM. (UTM since I discovered that Arcgis pro would not interpolate my LiDAR to a raster tif unless it's in UTM, otherwise you will receive the error "000210: Cannot create output " ) This creates a *_NEW.las file.:

NOTE: the conversion took some weeks for me since I had about 200 Lidar tiles (0.5 - 1 GB in size)
NOTE: You typically receive lidar data as .laz files; my .laz were unzipped in GIS before using vdatum

2. In ArcGIS pro:

a) Define Projection for each tile, set the coordinate system to WGS_1984_UTM_Zone_17N
b) Create LAS Dataset, this creates a .lasd file, use about 5 tiles to start then add the other .las files to the .lasd using Add Files to LAS Dataset
c) In Catalog window right-click on the compiled .lasd dataset got to >> Statistics >>click update
d) For DSM creation >> add the .lasd dataset to the map>> under LAS Dataset Layer >> Appearancec >> LAS Points >> First Return Points
e) Uncheck "noise": right-click on .lasd layer >> LAS Filter >> uncheck Noise in Classification Codes
f) LAS Dataset to Raster:

g) Project Raster to project from WGS84 UTM to GCS_WGS_1984:

h) Export raster only defining GCS_WGS_1984 coordinate system do not specify a vertical coordinate system, else GDAL will give an error: "error PROJ: proj_create_from_database: crs not found in gdal"

Plotting tif with gdal:

import elevation
filename = "HERNANDO_2018_LiDAR_GCS_VCS_WGS84m.tif"
gdal_data = gdal.Open(filename)
gdal_band = gdal_data.GetRasterBand(1)
nodataval = gdal_band.GetNoDataValue()

convert to a numpy array:
data_array = gdal_data.ReadAsArray().astype(np.float64)
data_array

replace missing values if necessary:
if np.any(data_array == nodataval):
data_array[data_array == nodataval] = np.nan

Plot our data w/ 'contourf':
fig = plt.figure(figsize = (12, 8))
ax = fig.add_subplot(111)
plt.contourf(data_array, cmap = "viridis",
levels = list(range(-30, 70, 5)))
plt.title("HERNANDO LIDAR ")
cbar = plt.colorbar()
plt.gca().set_aspect('equal', adjustable='box')
plt.show()

Convert tif DSM to ISCE format

!gdal_translate -of ISCE HERNANDO_2018_LiDAR_GCS_VCS_WGS84m.tif HERNANDO_2018_LiDAR_GCS_VCS_WGS84m.dem.wgs84

after using gdal_translate to generate .vrt for ISCE:
import isce
from applications.gdal2isce_xml import gdal2isce_xml
fname = 'HERNANDO_2018_LiDAR_GCS_VCS_WGS84m.dem.wgs84'
xml_file = gdal2isce_xml(fname)

After running stackStripMap.py, so far I have run:
run_01_crop, run_02_reference, run_03_focus_split and currently run_04_geo2rdr_coarseResamp with out any errors.
I'll update if any; also not sure if I'd know if they are DEM related, but I have run this same processing before with the SRTM DEM, and I am repeating the exact steps.

[alejobeap]

Dear All
I would like to change my Copernicus DEM into ISCE format or to radar coordinates. Can you help me? Thank you.
Alejandro

[pizzawoo]

Did you find the solution ? I've met the same situation now. I can not figure out how to create the necessary .xml and .vrt files.

[EJFielding]

The easiest way to get Copernicus DEM products into the ISCE2 format is to use the sardem package, which will do the geoid correction and write out the files you need when you use the -isce option:
https://github.com/scottstanie/sardem

[pizzawoo]

Thanks for the quick reply. I should try this later.