[PROJ] PROJ 6 vertical transformation confusion

[Even Rouault]

Chris,

The answer to your 2 questions is closely linked

> Question 1: is it possible that the vertical offset file transformation 7860
> is not being recognized by PROJ.
> 
> When I run " projinfo -s EPSG:4959 -t EPSG:4959+5759" to transform NZGD2000
> heights to AUCK1946  it only finds one candidate transformation, which is
> via NZVD2009 (4459+4442).  It doesn't seem to find the more accurate
> transformation via NZVD2016 (7840+7860).

The issue is that EPSG:7860 'NZVD2016 height to Auckland 1946 height (1)' uses 
the EPSG:1071 'Vertical Offset by Grid Interpolation (NZLVD)' method which is 
not currently implemented by PROJ (cannot deal with .csv files). If the grid 
was available under a format that PROJ can understand (let's say .gtx), then 
some little code addition (extending the NADCON specific check in https://
github.com/OSGeo/PROJ/blob/master/src/iso19111/coordinateoperation.cpp#L8763) 
and a tweak in the database should make it available.

> 
> Question 2: what determines how hard PROJ looks for transformation paths.

That's a question whose answer is a lengthy document I'm in progress of 
writing :-) It is a subtle mix of logical steps and heuristics.

> Instead fails as below:
> 
> # projinfo -s EPSG:4959+7839 -t EPSG:4959+5759
> createOperations() failed with: Unimplemented

See answer to first point

> It shouldn't have to do this in any case as it should be able to get use
> 7860 operation directly  but even if it can't find that I would have
> thought it would be able to find the indirect path through  NZVD2009.

Well that's one of the many limitations to keep the computation time somewhat 
bounded (otherwise could probably be a travelling salesman-like problem. PROJ 
6 is a sort of routing engine I tend to joke about). As soon as there's a 
direct transformation found in the database, it doesn't try to compute 
concantenated operations through an intermediate. Could/should probably be 
tweaked a bit to check if the direct transformation found is actually 
supported.

Actually I just fixed that. So now it *does* try to use an intermediate when 
the direct path is not supported by PROJ .. but that doesn't still work. As 
far as I can see grepping in the EPSG dataset, one cannot do Auckland 46 -> 
NZVD2016 by using the NZVD2009 intermediate. There is EPSG:4442 "NZVD2009 
height to Auckland 1946 height (1)", but I can't see any direct transformation 
between NZVD2009 and NZVD2016. That would require going through the NZGD2000 
CRS by appling one geoid and one reverse geoid. PROJ only supports researching 
one intermediate CRS as I mentionned above. What could be done however is that 
LINZ would register in the EPSG dataset a NZVD2009 to NZVD2016 concatenated 
operation chaining EPSG:4459 and EPSG:7840 (for the sake of operation 
research, a concatenated operation is seen as a unique operation. It is 
expanded to its component steps afterwards). Of course, you're far better 
placed than me to determine if that's the best route to follow (I could also 
imagine a direct adjustment grid). 

> I am running all these operations on osgeo/proj docker image, proj version
> Rel. 6.3.0, January 1st, 2020.
> 
> BTW I am loving projinfo - it is saving my sanity!

Glad you like it !

If you really want to know what it does, because it just returns results 
filtered and sorted, you can build it with
-DENABLE_TRACING -DTRACE_DATABASE -DTRACE_CREATE_OPERATIONS
You'll get logs of thousands of lines of every step it tries.

Even

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