download_dem() via py3dep (USGS 3DEP) incorrectly assumes ellipsoid

[brentwilder]

Currently download_dem() function within the pipeline automatically gets USGS 3DEP 1m elevation product IF the user does not provide a DEM.

It incorrectly assumes USGS DEM are in ellipsoid and/or that performing projection in rasterio takes care of vertical datum change.

Options @Surfix and I discussed to fix this:

1. We could remove this entirely and perhaps include as a separate utility script , prior to the new geoid-tool?
2. We could add an if-statement to the download_dem() function to change the vertical datum to ellipsoid if the input lidar point cloud is also in ellipsoid (using Ames Stereo Pipeline method)
3. We could remove this component entirely, and just point users to get their own DEM either from USGS 3DEP, https://www.usgs.gov/tools/download-data-maps-national-map, or other. This may be best use case, because then we're not constrained to define what resolution they want (1, 10, 30, 90 meters). Also, there were some serious issues I was having early on when we tested because py3dep uses a dynamic map over the static maps data products, which introduced a lot of tiling artifacts in testing, see closed issue here: Edge artifacts ~60m present in USGS dem.tif #2 . Reference DEMs should be of the highest quality for this tool to work best, and so, it would make sense for users to have more control with how/where getting DEM data.

[ZachHoppinen]

I think removing this makes sense.

The py3dep module we are using is deprecated and likely to start breaking even more going forward.
from their github:

Users are encouraged to using the new HyRiver package called Seamless3DEP which is a lightweight and efficient package providing access to both static DEMs and dynamic 3DEP service. This package is a part of the HyRiver software stack and is the recommended package for accessing 3DEP data. For more information please visit its documentation. For the time being, Py3DEP will be maintained for bug fixes and minor updates.

Also, from version 0.19, the default CRS for get_map was changed from EPSG:4326 to EPSG:5070. This is due to a recent issue with the 3DEP service where the it returns invalid results when the requested CRS in not in a projected coordinate system.

So we either need to migrate to the new Seamless3DEP (which I have code for but would take some work) or better would be just to tell the user to provide their dem and add some checks.

Also, this kind of functionality isn't really the point of this package and we should focus on the core competencies of this.

[brentwilder]

@ZachHoppinen I totally agree with you. This was reeeally only built into it in the early days for testing.

Further, if the user provided their own DEM we have this new geoid-tool that can do the checks and/or transform the vertical datum if needed… Pending pushing the latest pull request 🙃 …

[jomey]

I agree with @ZachHoppinen as well.

With the new geoid tool, the step of having a proper reference DEM is also moved outside the pipeline.

There are too many options to get a reference DEM and we should not dictate one with this tool.

[brentwilder]

Sounds good all, I will link this with the current PR and update the code to remove this.