3D Geospatial Transformations (Geoid <-> Ellipsoid heights)

[scottyhq]

I'm interested in using rioxarray directly for geodetic transforms. For example, it's common to need to apply a reprojection pipeline to go from Geoid to Ellipsoid heights with datasets like Copernicus DEM

I think this might require opening an issue or discussion on rasterio since it's rasterio.warp.reproject() under the hood, but figured I'd start here!

projinfo -s EPSG:4326+3855 -t EPSG:7661 -o PROJ --hide-ballpark --spatial-test intersects

PROJ string:
+proj=pipeline
  +step +proj=axisswap +order=2,1
  +step +proj=unitconvert +xy_in=deg +xy_out=rad
  +step +proj=vgridshift +grids=us_nga_egm08_25.tif +multiplier=1
  +step +proj=unitconvert +xy_in=rad +xy_out=deg
  +step +proj=axisswap +order=2,1

☝️GDAL figures out how to do this automatically:
gdalwarp -s_srs EPSG:4326+3855 -t_srs EPSG:7661 $INPUT $OUTPUT

Similarly, I'd naively expect daEllipsoid = daGeoid.rio.reproject(7661) to work. But right now rioxarray does not throw any errors or warnings, but also does not do the correct thing which is retrieving us_nga_egm08_25.tif to apply a vertical shift.

Full notebook with logging outputs captured here https://gist.github.com/scottyhq/2153cb01cf951680a4dcb9e1e84b69d7

[snowman2]

Related rasterio issue: rasterio/rasterio#1990

[snowman2]

But right now rioxarray does not throw any errors or warnings, but also does not do the correct thing which is retrieving us_nga_egm08_25.tif to apply a vertical shift.

I got PROJ_DEBUG to work:
https://gist.github.com/snowman2/91a2c99fc9e6f958cbc8b91970879c72

It appears to be retrieving the correct file in the transformation. This leads me to believe that this difference is likely related to this rasterio issue: rasterio/rasterio#1990

[scottyhq]

Interesting! thanks @snowman2. Any idea why PROJ_DEBUG is picked up by os.environ[] but not by `rasterio.Env()'? It's definitely useful in my experience to see what transform is being used. Especially in cases where multiple candidate transforms are found (NAVD88->WGS84, etc). Not being familiar with the GDAL internals it's somewhat unclear what the scope is of that linked rasterio issue. Should I open a separate issue on the rasterio tracker and link to this discussion?

[snowman2]

PROJ_DEBUG is a global environment variable not specific to GDAL: https://proj.org/usage/environmentvars.html#envvar-PROJ_DEBUG

[snowman2]

Not being familiar with the GDAL internals it's somewhat unclear what the scope is of that linked rasterio issue.

rasterio reproject uses a different method under the hood than gdalwarp. The issue is to update the internals of the reproject method so they are the same.

I am thinking that once it is updated, the results will match.

[snowman2]

Should I open a separate issue on the rasterio tracker and link to this discussion?

It wouldn't hurt to open a discussion in the rasterio repo.

[dshean]

I remember running into issues like this around network grid functionality: https://github.com/OSGeo/PROJ/blob/master/docs/source/usage/network.rst#how-to-enable-network-capabilities-

I quickly looked at the gist, and see you set several env, but nothing around network options. Not sure if this is the problem. The cludgy workaround is to run projsync --all and store all grids locally, where GDAL/PROJ knows to find them by default.

[scottyhq]

Yes - looks like from the above, the network grids are retrieved without any extra settings, but to get them all locally with https://github.com/conda-forge/proj-data-feedstock (conda install proj-data) ~600MB.