[PROJ] Converting NAD83/NADV88 heights to HAE

[Even Rouault]

On mercredi 1 juillet 2020 22:45:47 CEST Patrick Young wrote:
> Thanks all for the help!
> 
> Dinking around a little bit more, I found through trial and error that
> eliminating the +vunits seemed to solve it:
> 
> echo 565890.947 4192011.754 *55* | cs2cs +proj=utm +zone=10 +datum=NAD83
> +units=m +no_defs +geoidgrids=us_noaa_g2018u0.tif +to +proj=utm +zone=10
> +datum=WGS84 +units=m +no_defs
> 565890.95       4192011.75 *22.79*
> 
> Somehow, the addition of +vunits=m on the target angers proj:
> 
> echo 565890.947 4192011.754 *55* | cs2cs +proj=utm +zone=10 +datum=NAD83
> +units=m +no_defs +geoidgrids=us_noaa_g2018u0.tif +to +proj=utm +zone=10
> +datum=WGS84 +units=m *+vunits=m *+no_defs
> 565890.91       4192011.91 *55.00*
> 
> Adding +vunits=m to the source seems not to matter.  Wish I could say more
> on why this is the case; is it a bug or is it a malformed request?

When you specify +vunits= in this PROJ.4 legacy syntax, PROJ >= 6 understands it as an 
instruction to create a CompoundCRS with a vertical CRS with a unknown datum, not a 
Projected CRS with ellipsoidal height. And as this is a unknown vertical datum, no 
transformation between it and the WGS84 ellipsoid can be done. This is admitedly a bit 
confusing, non backward compatible and not consistent with the case where you do 
+proj=longlat +vunits=m where it creates a Geographic 3D CRS. I've captured this in https://
github.com/OSGeo/PROJ/issues/2287


Another way of doing this transformation is to use cct and PROJ pipeline syntax:

echo 565890.947 4192011.754 55 | cct \
  +proj=pipeline +step \
      +inv +proj=utm +zone=10 +ellps=GRS80 \
      +step +proj=vgridshift +grids=us_noaa_g2018u0.tif +multiplier=1 \
      +step +proj=utm +zone=10 +ellps=WGS84

The following also works, but you depend here on the fact that geoid2018 is the best 
candidate used currently for the ellipsoidal<-->orthometric height transformation, and 
specifying NAD83(2011) rather than NAD83 is critical to make it sure that geoid2018 is used. 
If using plain NAD83, an older version of the geoid model will be used:

echo 565890.947 4192011.754 55 | cs2cs \
     "NAD83(2011) / UTM zone 10N + NAVD88 height" "WGS84 / UTM zone 10N"

or

echo 565890.947 4192011.754 55 | cs2cs EPSG:6339+5703 EPSG:32610


Even

-- 
Spatialys - Geospatial professional services
http://www.spatialys.com
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/proj/attachments/20200702/f7b4ee2d/attachment.html>