[PROJ] Strange values with NAVD88

[Even Rouault]

Javier,

> Doing some tests with NAVD88 (EPSG:5703), I got some unexpected results.
> See the next 5 commands
> 
> a) echo 40.69 -75.04 0 | PROJ_NETWORK=ON cs2cs EPSG:4979 EPSG:4269+5703 -d 9
> 40.689991531 -75.039999022 35.046181369
> b) echo 40.69 -75.04 0 | PROJ_NETWORK=ON cs2cs EPSG:4979 EPSG:4152+5703 -d 9
> 40.689991800 -75.039998369 35.046181369
> c) echo 40.69 -75.04 0 | PROJ_NETWORK=ON cs2cs EPSG:4979 EPSG:6318+5703 -d 9
> 40.690000000 -75.040000000 33.784496002
> d) echo 40.69 -75.04 0 | PROJ_NETWORK=ON cs2cs EPSG:4979 EPSG:4152 -d 9
>      40.690000000 -75.040000000 0.000000000
> e) echo 40.69 -75.04 0 | PROJ_NETWORK=ON cs2cs EPSG:4957 EPSG:4152+5703 -d 9
> 40.690000000 -75.040000000 33.786615753

I initially wrote """using "projinfo -s X -t Y --spatial-test intersects --
grid-check none" should answer most of your questions.""", which is not false, 
but here you fall into situations that exercice the oracle in PROJ in very 
subtle ways.

Adding --bbox -76,40,-75,41 helps restricting the set of transformations, 
except that it has some very subtle differences in the sorting order, which 
don't always exactly lead to the same result than when using cs2cs, which has 
no prior bbox knowledge... (cs2cs will basically use the first transformation 
listed by projinfo --spatial-test intersects --grid-check none for which the 
point of interest falls into, and that has the best accuracy)

The issue here is with "NAD83(HARN) to WGS 84 (1)", a noop, with 1 m accuracy 
that has a wider extent of use than "NAD83(HARN) to WGS 84 (3)", a Helmert 
based one, with 1m accuracy too, but a lesser extent. So when using --bbox 
-76,40,-75,41, due to the intersection of the area of use of both 
transformation with the specified bbox being the same, "NAD83(HARN) to WGS 84 
(3)" will be prefered over "NAD83(HARN) to WGS 84 (1)", due to it's (3) higher 
version number... Without --bbox, since  "NAD83(HARN) to WGS 84 (1)" has a 
larger extent, it is prefered. Subtle...
 
> Comparing b) and d), why are XY coordinates changed when NAVD88 is
> included? It is the same using EPSG:4957 instead of EPSG:4152 in d). 

In b), the restriction of the area of use due to the NAVD88 transform makes 
"NAD83(HARN) to WGS 84 (3)" (Helmert transformation, done in 3D) to be used 
(since then the transformations using "NAD83(HARN) to WGS 84 (3)" and 
"NAD83(HARN) to WGS 84 (1)" have the same area of use, and thus (3) is 
prefered)

In d), there is no such restriction to the area of use of the transformation, 
so "NAD83(HARN) to WGS 84 (1)" is prefered due to its wider are of use. 
Perhaps the selection logic in proj_trans() should have an extra case for 
"when there are 2 transformations with same accuracy and for which the point 
to transform falls into, select the one with the smallest area of use since it 
is likely to be more accurate". I don't know. Adds yet another fragile 
heuristics to many existing ones.

What sucks here is that both transformations in EPSG are advertized with the 
same accuracy, and that cause PROJ to be ultra sensitive. If the Helmert one 
was advertized with a slightly better accuracy, things would be more 
predictable.

> I thought that NAD83(HARN) and WGS84 were considered the same.

See above. There is a transformation in EPSG for which this is true, and 
another one, not.

> In c), the altitude has more than 1 meter difference. Why? Am I doing or
> assuming something wrong? Checking all the american geoids in this area in
> cdn.proj.org, there are no more than a few centimeters difference.

b) has a Helmert transformation done in 3D for WGS 84 to NAD83(HARN), whereas 
there is no registered transformation between WGS 84 and NAD83(2011), so c) 
only applies the geoid transformation.

> I do not know if it is related, but I have only seen ballpark (horizontal)
> transformations between NAD83(HARN) and NAD83(2011). Is that correct? I was
> expecting some transformation, either direct or a grid.

We probably lack NADCON5 grids here:
https://github.com/OSGeo/PROJ/issues/2194

(Anyway using generic WGS84, you shouldn't expect sub-metre accuracy...)

Even

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