Add feature to download global DEM dataset

[adehecq]

We could have a feature similar to this directly implemented in xDEM to download global DEM dataset for any region.

[rhugonnet]

Agreed! OpenTopo is a great resource that does a lot of the work for us.

@ShashankBice Where do you think this code should live?

I see three options:

• In this separate repo you set up (maybe not the easiest for long-term maintenance: need to add its own doc/tests/CI, etc),
• In a new "DEM products" repo we could start, not necessarily limited to OpenTopo,
• In a new xDEM module (e.g., xdem/products.py)?

[erikmannerfelt]

That would be an awesome feature!

[ShashankBice]

Hi Romain, Thanks for the ping! @dshean also initiated a thread about this on slack when I was on vacation. The idea there was similar, to have something more generic, and include both opentopography and @friedrichknuth's geodata (DEM downloader from Microsoft Planetary Computer) functionalities. I would be open to options 2/3, with preference for option 2. It might be easier for new folks to query or learn about a "reference_dem package" or whatever the name we come up with, than for it to be embedded in the larger xdem repo. I am thinking of what happened with the centerline tool embedded within OGGM stack, and I think for users wanting to have reference DEM for their study, the standalone repo would make more sense. We can inherit the xdem features and tests in the repository, and I would love to learn how to contribute to a repository properly with tests and everything. Codeblocks in xdem could then get the DEMs where needed from that package we come up with. We could also rename the current opentopo DEM repository to something else and write everything there. Really open to what all of you suggest! Let me know, in the meantime, I will work on implementing a geoid<-->ellipsoid bridge. Based on what we all decide, we could properly build the repo, maybe when you get back to Wilcox next time :wink:

Have a great weekend everyone, Cheers, Shashank

[rhugonnet]

Sounds good @ShashankBice! I agree option 2 makes the most sense, it can be its own package with very few dependencies, and be wrapped in xDEM (and other packages) for users who want to access things from there! :slightly_smiling_face:

I agree we'd have to change the name. For the location: doesn't matter, let's see what @friedrichknuth & @dshean think! Building on top of opentopo_dem or geodata are good options, but I think this choice depends on who'd be motivated to be a long-term maintainer (one person very motivated: personal repo; group of people: organization like GlacioHack or UWCryo, but institutional groups are a bit more temporary)

On the 3D CRS transformations: in case you're not aware (to avoid duplicating efforts!), since May xDEM supports transformations from/to any 3D CRS: https://xdem.readthedocs.io/en/stable/vertical_ref.html. There's a big test suite: https://github.com/GlacioHack/xdem/blob/main/xdem/vcrs.py and I had to open some PRs in pyproj to make this consistent and operational for any possible 2D/3D CRS combination.

For managing the 3D CRSs: I think it's important to properly define (even overwrite) the metadata if necessary. That could be done in the new repo. But transformations should probably not live there as they require so many more geospatial dependencies, we can leave them to the user (and point towards pyproj, xDEM, etc...).

Have a great weekend :smile:

[friedrichknuth]

+1 for option 2: a simple package with few dependancies that is solely dedicated to public DEM data retrieval

+1 host in organization like GlacioHack or UWCryo: ideally unaffiliated with a single institution to maximize likelihood for collective contributions and maintenance

Happy to proceed with opentopo_dem, transfer geodata, or scrap and start fresh!

[ShashankBice]

Thanks Romain, The documentation on xdem 3D suite is awesome, excited to try it out :D On your second point, I do think however providing outputs with ellipsoidal height might be valuable, consider this situation: looking at myself 7 to 8 years ago, i wanted to convert from geoid to ellipsoid for a project, and i did it somehow, but I did not know at that point if I was using the correct realizations etc.. This is a complex topic, and to make publicly available DEMs more accessible, it would be good if we took care of the the mapping properly (atleast the vertical shift) on our end, using the new xDEM functionality. The average user might just be happy with the standard ellipsoidal height we return with respect to the latest/most widely used realization using the dictionary mapping we have. We could print out a note or something for more advanced users to take a look at the 3D crs documentation on xdem and use xdem directly for more advanced/customized changes. How do you feel about this?

Cheers, Shashank

[rhugonnet]

Thanks! I'm not sure I completely understand your point, but I think I agree!

I was thinking in the same direction in that one of the main problem with vertical CRS (and the steepest learning curve for new DEM users you mention) is not necessarily performing the transformation (even though it's not trivial with current PROJ behaviour), but mostly to define the initial vertical CRS properly (what people would currently have to do with .set_vrs() in xDEM if the input is not recognized).

So, what we could do:

1. In opentopo_dem (or any name we give to the future DEM data package, let's name it demdata), we would properly set the initial vertical CRS in the DEM metadata, so that users don't have to worry about this. This allows the 3D CRS to be read and understood by any GIS software (not just xDEM), but leaves the initial data untouched for reproducibility;
2. We could then point users to xDEM if they want to transform their DEM. Because the initial vertical CRS is already properly defined in 1, it will be read automatically and the transformation is as simple as the one-liner DEM.to_vcrs(), with easy string inputs for the most common vertical CRS of global DEM products ("Ellipsoid", "EGM08", "EGM96").

This also allows to let dependencies live in the most convenient place:

• The demdata package wouldn't need a xDEM dependency (if demdata uses a functionality from xDEM, then xDEM wouldn't be able to call demdata: discouraged circular dependency; or the solution would be to have it live in the same package),
• Even better, demdata wouldn't need a pyproj dependency and stay lightweight: we can write the 3D CRS WKT in full text, or just link to a 3D EPSG code if one exists for that product,
• In xdem, we could link to demdata (for tests, examples, or generally user convenience).

What do you think @ShashankBice?

[dshean]

Appreciate all of the discussion. A few thoughts:

• Potential names: demfetch, fetchdem,demget,getdem
• We need to be careful about "ellipsoid" and definitions of "WGS84 Ellipsoid." Need proper handling of different realizations (e.g., transit vs G1762), ITRF realizations, and associated epoch of data collection. See https://ascelibrary.org/doi/10.1061/%28ASCE%29SU.1943-5428.0000389.
• If we properly define the 3D CRS for a given OT DEM (WKT2 is likely the only way to do this unambiguously until we have better EPSG codes), then we should be able to use standard GDAL/PROJ functionality for all subsequent 3D CRS transformations. This can be done with a gdalwarp command (and a proper 3D CRS WKT2 definition for the t_srs), without having to start Python interpreter, import many modules, load the DEM into memory, and do the transformation. As far as I can tell, the current xDEM vcrs functionality assumes the transformation only requires a simple vertical offset using the same horizontal CRS. This is fine for many cases involving Global DEM products, but won't properly handle cases where the input or output CRS uses a different ellipsoid model (e.g., NAD83(2011)). This is the case for many national or regional elevation datasets in the U.S. and Europe that do not use "WGS84."

[rhugonnet]

Agree!

This can be done with a gdalwarp command (and a proper 3D CRS WKT2 definition for the t_srs), without having to start Python interpreter, import many modules, load the DEM into memory, and do the transformation.

The same will soon become possible using Dask through an Xarray accessor of xDEM's functionalities (relying on rioxarray). The whole point of writing these Python packages is to not rely on GDAL bash that requires I/O at every step, and be able to scale + allow more advanced analysis pipelines.

As far as I can tell, the current xDEM vcrs functionality assumes the transformation only requires a simple vertical offset using the same horizontal CRS.

This a fairly easy fix, PROJ already returns the new x/y in our current implementation if the transform also reprojects horizontally.

[dshean]

The whole point of writing these Python packages is to not rely on GDAL bash that requires I/O at every step, and be able to scale + allow more advanced analysis pipelines.

That may be true for xdem, your preferred analysis workflow, your preference to use Python, and your specific research applications (which require scaleable design). But that is not necessarily true for most anticipated users of this fetchdem tool. The point here is to have a simple utility that can prepare a DEM without a lot of overhead. At the end of the day, we're just talking about storing the output from a single API call: https://github.com/dshean/demcoreg/blob/master/demcoreg/get_refdem.sh#L35.

The last time I talked with the OpenTopography folks, we agreed that it would be best if the GlobalDEM API accepted an optional argument with CRS definition (like an EPSG code), so they could do the transformation/projection on the server, which eliminates the need for a lot of what we're discussing.

Most applications require a single reference DEM for a single study site. This is typically a one-time preprocessing step, with the requirement that a transformed "ready-to-go" DEM is stored on disk. This DEM is then loaded by a range of other processing/analysis software (written in many different languages, including C++, R, Julia, not just Python). And many of these applications rely on command line tools, esp satellite image processing workflows, like SAR processing with ISCE2, orthorectification using GDAL or ASP.

It's fine to include documentation suggesting that xdem can be used for subsequent transformations and analysis. But I don't think it should be presented as the only approach to transform the DEM returned by this tool. We should make it clear that a single gdalwarp command can perform the required 3D CRS transformations.

[dshean]

This a fairly easy fix, PROJ already returns the new x/y in our current implementation if the transform also reprojects horizontally.

Let's talk about this - there are some nuances here, esp for transformations between ellipsoids that have a 3D offset between their origins (requiring a Helmert), as there is a corresponding vertical offset, not just horizontal, and this is not handled with the offset grids.

We've had a lot of good discussion on these topics with SlideRule team, and we had a nice, impromptu group discussion at the ICESat-2 Hackweek about these issues, and started compiling some notes (including more to add in issues): https://github.com/ICESAT-2HackWeek/3D_CRS_Transformation_Resources.

I think defining complete 3D CRS definitions (with epoch, like https://github.com/ICESat2-SlideRule/sliderule/blob/cb8ead40d761c8d637397a28e7b6d53fcf1de3c4/plugins/pgc/plugin/ITRF2014_3031.wkt) using WKT2 is ultimately the safest, simplest way to prepare the complex transformation pipelines required to support the full range of time-dependent elevation measurements out there. It's not going to be easy/fun, and the metadata for existing datasets is incomplete/wrong, but it would be a great service to define these in a centralized location for many of the common elevation data used by multiple communities.

The vcrs implementation is nice, and convenient for many cases, but it will require more effort to wrap all of the latest PROJ support for epoch-aware 3D CRS transformations.

[rhugonnet]

That may be true for xdem, your preferred analysis workflow, your preference to use Python, and your specific research applications (which require scaleable design). But that is not necessarily true for most anticipated users of this fetchdem tool.

Well now I'm a bit lost. Both @ShashankBice's opentopo_dem and @friedrichknuth's geodata are Python tools, so I expected them to at least include a Python API (for synergy with other packages like xDEM)? If the objective is to make the package also usable for a broader audience (with some Python "bin" scripts called in bash), that's completely fine and I think is a good idea.

We've had a lot of good discussion on these topics with SlideRule team, and we had a nice, impromptu group discussion at the ICESat-2 Hackweek about these issues, and started compiling some notes (including more to add in issues): https://github.com/ICESAT-2HackWeek/3D_CRS_Transformation_Resources.

This is a nice effort, thanks for sharing! I had already gone through those notes a couple weeks ago. For how it relates to the DEM class and vcrs, I don't see why this would be a large effort: DEM inherits a datetime attribute from the SatelliteImage class, and already relies on PROJ for transformations to which it can pass anything. So the 3D CRS transformation support will simply evolve as PROJ does (and its pyproj wrappers). If I understand correctly, the core problem is vertical CRS definition, and for this xDEM would rely on this new fetchdata repo we are discussing here (see points further above on how I proposed to split the functionalities). What are you thinking on this @ShashankBice? (I think @dshean goes in the same direction than what I listed, if I understood correctly)

[adehecq]

Hi, just catching up with the now long discussion! Here's my viewpoint:

• option 2 (standalone) module with a minimum of dependency is great. It would of course be great if we can easily link it with xdem, which should not be difficult!
• I also think a broader organization would be better for visibility. It's totally fine if this remains under UW's umbrella, but as pointed out by @friedrichknuth, this may discourage people outside UW to actively contribute. I'm happy to have it within "GlacioHack". About this, GlacioHack could even be renamed and reframed to be less glacier oriented, and maybe be linked to larger organizations like OSGeo...
• saving the VCRS within the saved file would be really great and make it easier to handle vertical references later. This information can sometimes be quite difficult to track down.
• I like "demfetch" or "fetchdem" !

Regarding @dshean's points, I agree that it's not worth for most people to have all the overhead of xdem to just do a vertical datum change once. But I don't know a tool that does that easily, is there any? (I also often uses ASP's dem_geoid for that, but ASP does not come with a smaller overhead 😉 and it's incomplete). At least, on xdem's side we could easily add a CLI tool to run the vertical datum change in a bash script, without starting a Python interpreter. I would actually love to have more of those tools in xdem, as I don't like having to open a Python console for small operations.

Anyway, happy to support this project, whatever direction it takes, as I'm interested in such tools.