Account for vertical datum offsets

[dshean]

Can create a dictionary for each available global DEM dataset and the required vertical datum transformation. For example COP30 requires EGM2008 to WGS84 offset grid: https://github.com/dshean/demcoreg/blob/master/demcoreg/get_refdem.sh#L47

Some of the available Global DEM datasets are already corrected for ellipsoidal heights, specified with the '_E' suffix (e.g. 'SRTMGL1_E')

[ShashankBice]

I agree with this idea. I think we can create a dict as you mention, and then let the user specify all DEM types with *_E extensions. We can then have a if else statement to check if the _E file for that is not provided by opentopography by default, we will do the transformations ourselves, as in the notebook here. This step should also come up with a bold message on the console saying that the vertical transformation is being done downstream of opentopography, on our end!

Let me know if this sounds like a plan,
Shashank

[dshean]

Sounds good!

Would also be good to verify that we can reproduce the OpenTopography _E products using our internal PROJ geoid offset corrections. There may be some discrepancies between how OT prepared the _E heights vs. how we're doing it.

In general, this process of overriding the default CRS in the geotiff header for DEMs lacking 3D CRS is something we are going to have to do over and over again. There are some examples in SlideRule repo for 3DEP and PGC ArcticDEM/REMA products. So we can prototype some frameworks here for OT examples using a dict, but good to think about options for a more general implementation that can be used by other packages. Can discuss.

[ShashankBice]

On your first point, they are not exactly one to one, see the left plot on code line 52 of this notebook.

From the products currently available on open topography, the following have ellipsoidal heights available directly from API itself:

• SRTMGL1 <--> SRTMGL1_E
• AW3D30 <--> AW3D30_E

These products are referenced with respect to the following geoid model and do not have existing ellipsoidal versions:

• SRTMGL3, EGM96 GEOID
• SRTM15Plus, includes heights over water surfaces, so I am less inclined to provide an ellipsoidal version for this
• NASADEM, EGM96 GEOID, EPSG: 5773
• COP30, EGM2008, EPSG: 3855
• COP90 EGM2008, EPSG: 3855
• EU_DTM, EGM2008, EPSG: 3855 Note this one also has a non-WGS84 horizontal CRS (ETRS89-extended / LAEA Europe [EPSG: 3035]

The GEDI_L3 product is already with respect to WGS84 vertical ellipsoid. I suggest we do not touch it for now.

For now we can provide the option to return ellipsoidal versions of the DEMs listed in second block. I am thinking we leave it to WGS84 ellipsoid for now, and if the users want a particular version of horizontal CRS, they can provide EPSG code for the horizontal part. This should be a good start, as I think epoch-aware CRS transformations will take some time for the community to fully digest, and we can revisit that once documentation is solid, with proper WKT2 strings. What do you think @dshean? I will work on adding these functions then today 😊

[dshean]

Seems like a reasonable plan.

I would not exclude SRTM15Plus - while it may include water, it has land elevations as well. Should be the same CRS as SRTMGL3 products.

I think it could be good to explicitly set compound CRS with vertical EPSG:4979 for the WGS84 ellipsoid height output datasets (and definitions for input datasets that do not require vertical offset). I think most of the output DEMs are in local UTM zone?

I think most of these inputs are using WGS84 (or GRS80 for EPSG: 3035), so don't need to worry as much about different ellipsoid models and the need to define a specific WGS84 realization to enable proper 3D transformations.

[ShashankBice]

The latest code has three major improvements:

• Ellipsoidal heights with respect to WGS84 datum can now be accessed for all data sets except EU_DTM
• The user can specify output horizontal CRS as EPSG codes using -out_proj key
• A final gdal_edit.py call promotes the DEM file to 3D CRS based on final CRS combination selected by user.

[ShashankBice]

Back from my walk, I have a question @dshean, how do we correctly specify input_crs for EU_DTM?
The horizontal CRS uses GRS80 Ellipsoid as datum, while heights are with respect to EGM:2008 geoid model, which I think is defined with respect to WGS84 ellipsoid. When trying to use a compound CRS, gdal complains ERROR 1: PROJ: proj_create: The 'vertical' geographic CRS is not equivalent to the geographic CRS of the horizontal part when trying to run the following command: gdalwarp -r cubic -co COMPRESS=LZW -co TILED=YES -co BIGTIFF=IF_SAFER -s_srs 'EPSG:3035+EPSG:3855' -t_srs 'EPSG:3035+EPSG:7912' EU_DTM_E_temp.tif EU_DTM_E.tif

Maybe I am missing something simple, so all input is appreciated!

[dshean]

Yeah, the compound CRS need to use the same ellipsoid. Even though latest definitions of GRS80 and WGS84 are very similar (within mm), they're not exactly the same.

The issue here seems to be the t_srs EPSG:3035+EPSG:7912.

projinfo EPSG:3035+EPSG:3855 works fine.
projinfo EPSG:3035+EPSG:7912 gives the error you mention.

I tried combinations of other EPSG for ellipsoid heights, but non worked out of the box.

I think -t_srs "$(projinfo EPSG:3035 --3d -o WKT2:2019)" will give you the appropriate WKT2 that can be used for output CRS in the gdalwarp command.

[dshean]

This is my first time working with the ETRF ensemble datum. Looks like similar issues for WGS84 Ensemble, only the ETRF has much better ensemble accuracy (0.1 m as opposed to 2.0 m)