[PROJ] Strange values with NAVD88

[Even Rouault]

> First, I don't mean to contradict Even, but I will caution you that Even
> is speaking from the viewpoint of proj and the EPSG databse,which is a
> highly sensible place to be speaking from, especially on this lst.  I am
> coming from a farther-from-proj perspective that might be described as
> NGS-like.  (I do not speak for NGS; I'm just a reader of their
> publications.  But I think if you asked them, they might say similar
> things.)

Exchanges and brainstorming are welcome. I'm perfectly aware that there are 
many limitations to the approach taken and that humans can come up with more 
relevant transforms. That's why projinfo output should be examined to look at 
all alternatives that PROJ could find and potentially pick up something 
different than the "best" / default transformation.
Using cs2cs is a bit like searching an information in your favorite search 
engine and always trusting the top result to be the one you need.
I'm not sure how to overcome those limitations while keeping a data-driven 
approach with as generic as possible heuristics to make the best use of it. 
Perhaps adding some auxiliary table(s) in the database that would influence 
how PROJ ranks its results to override its default heuristics, but maintaining 
that could become complicated.

> I therefore expect a null transform from 4979 to 4269, both in lat-lon
> and in HAE.
> I cannot explain the horizontal difference.  As I understand proj
> doctrine, no transform is supportable here.

Here's what PROJ come from for that case:

Inverse of NAD83(HARN) to WGS 84 (3) + NAD83(HARN) to NAVD88 height (8) + 
Inverse of NAD83 to NAD83(HARN) (46)

+proj=pipeline
  +step +proj=axisswap +order=2,1
  +step +proj=unitconvert +xy_in=deg +xy_out=rad
  +step +proj=cart +ellps=WGS84
  +step +inv +proj=helmert +x=-0.991 +y=1.9072 +z=0.5129 
+rx=-0.0257899075194932
        +ry=-0.0096500989602704 +rz=-0.0116599432323421 +s=0
        +convention=coordinate_frame
  +step +inv +proj=cart +ellps=GRS80
  +step +inv +proj=vgridshift +grids=us_noaa_g1999u08.tif +multiplier=1
  +step +inv +proj=hgridshift +grids=us_noaa_pahpgn.tif
  +step +proj=unitconvert +xy_in=rad +xy_out=deg
  +step +proj=axisswap +order=2,1

So chaining:
- Helmert between WGS 84 and NAD83(HARN) done in 3D. The remarks of this 
transformation mention that it comes from 'ITRF96 to NAD83(CORS96) (1)'  
"ignoring time-dependent parameters and assuming ITRF96(1997.0) and WGS 84, 
plus NAD83(CORS96) and NAD83(HARN), can be considered the same within the 
accuracy of the transformation"
- NAD83(HARN) to NAVD88 height (8) using Geoid1999
- NAD83 to NAD83(HARN) (46) using NADCON4 grid for Pennsylvania

This is convoluted and questionable. The irony is that we use Geoid1999, 
because the "NAD83 to NAVD88 height" transformations using Geoid2003 have been 
deprecated (probably because NAD83(86) ellipsoidal height is not defined. 
Those transformations using Geoid2003 have been superseded by one between 
NAD83(FBN) and NAVD88), so PROJ needs to find another datum that has 
transformations to NAVD88 

projinfo also proposes alternatives using NAD83(NSRS2007) instead of 
NAD83(HARN), and NAD83(NSRS2007) to NAVD88 height uses Geoid2009. But when 
PROJ "computes" the overall accuracy of that alternative, it is slightly worse 
than the one using NAD83(HARN) : 1.2m vs 1.1m ...

But as we've underlined, anything involving "WGS 84" datum ensemble to NAD83 
(something) is going to be fuzzily defined, since that WGS 84 is always 
approximated to another datum, and which datum it is approximated to isn't 
necessarily the same depending on the transformation.

> 
> > b) echo 40.69 -75.04 0 | PROJ_NETWORK=ON cs2cs EPSG:4979 EPSG:4152+5703 -d
> > 9> 
> >      40.689991800 -75.039998369 35.046181369
> 
> 4152 is NAD83(HARN). 
> It's a really good question what's going on here about different
> horizontal coords.

See above and my previous answer. There are 2 candidate transformations in 
EPSG between WGS 84 and NAD83(HARN). One that is a null transformation ("For 
many purposes NAD83(HARN) can be considered to be coincident with WGS 84.") 
and another one, used in the first step of a) that is a Helmert transformation 
coming from ITRF96 to NAD83(CORS96)

> > transformations between NAD83(HARN) and NAD83(2011)

> There should be a grid transform. See NCAT that I linked above.  It is
> entirely possible that this precision of NGS transform is not yet
> incorporated in proj. 

Yep, as mentionned in my previous email, NADCON5 grids which offer 
transformations between all successive "generations" of NAD83 are not yet 
used.

-- 
Spatialys - Geospatial professional services
http://www.spatialys.com