[PROJ] NTv2 grid file from Bessel ellipsoid to Bessel ellipsoid

[Even Rouault]

Jochem,

> 
> I noticed something odd in cs2cs results, which has nothing to do with the
> new Proj. version, but could be a old bug. I'm using the OSGeo4W Shell of
> QGIS 3.0.1 (Proj.493).
> 
> BACKGROUND
> I am testing with different horizontal grid shift files, since we are
> defining a new official coordinate tranformation between ETRS89
> (ETRF2000(R14)) and the Dutch RD (new) grid and NAP heights. We like to
> provide two variants. First, a easy implementation variant with a grid shift
> file that handles both the datum transformation (7 parameter Helmert
> transformation from Bessel to GRS80) and the datum distortions (correction
> grid on Bessel ellipsoid) in a single NTv2 (.gsb) grid file. Since this
> method can only be used within the bounding box of the grid file, we also
> define a grid without the datum transformation and publish the 7 paramaters
> for the datum transformation separately. This 7 parameter transformation
> together with the second grid should give the same results as the first
> grid within the bounding box of the grid file, but this second method will
> also provide useful results outside the bounding box of the grid.
> 
> PROBLEM
> I do not get the results I expected. This is caused by something odd in
> cs2cs. To demonstrate this, I created a NTv2 grid file with only zero
> corrections (see ascii grid file below). If I transform coordinates on the
> Bessel ellipsoid with this grid (command 1), I do not get unchanged results
> (command 2), but I get results which are the same as I get when changing the
> ellipsoid (command 3).

Yes this is somewhat logical given what the code does (which of course is true 
of any bug :-)). But here I believe this is more a feature than a bug. This 
reflects the somewhat opacity/magic of the historical "early-binding" 
approach.

Command 1:
You use +ellps=bessel +nadgrids=zeros.gsb in the left term, so there is first 
a geographic shift using the grid values to go from bessel to WGS84 geographic 
coordinates (but if you change bessel by anything else like R=1 this will give 
the same result), which is null here given your grid, and then a 
transformation from WGS84 geographic coordinates to WGS84 earth centric 
cartesian coordinates. The right term involves an ellipsoid change to Bessel 
when doing the cartesian->geographic step, hence the modified value on the 
latitude

Note: you could have also use the provided 'null' grid instead of zeros.gsb

Command 3:
+ellps=GRS80 +towgs84=0,0,0,0,0,0,0 ==> transformation from GRS80 geographic 
coordinates to GRS80 earth centric cartesian coordinates, followed by a null 
shift to get WGS84 earth centric cartesian coordinates. As GRS80 ~= WGS84, the 
effect of this is identical to command 1

Command 2:
+ellps=bessel +towgs84=0,0,0,0,0,0,0 is found in the left and right terms, so 
this is a no-operation.


For what you describe in your needs, you'd better use PROJ 5.2 or the nearly 
released PROJ 6.0 and use PROJ pipelines with the 'cart', 'helmert' and 
'hgridshift' operations that will avoid using the WGS84 pivot implied by 
towgs84&nadgrids, which is likely not appropriate for the tranformations you 
want to accomplish.
PROJ 6 might be even more appropriate than PROJ 5.2 since we have now the new 
push/pop operators that can be used to avoid modifying the height when doing a 
cartesian -> helmert -> inverse cartesian chain.

Example of a pipeline you may use with the cct program:
+proj=pipeline
	+step +proj=push +v_3
	+step +proj=cart +ellps=GRS80
	+step +proj=helmert +x=1 +y=1 +z=-1
	+step +inv +proj=cart +ellps=WGS84
	+step +proj=pop +v_3 

And the pipeline approach is also more natural to chain together a grid shift 
and a Helmert tranformation, as you explain in the second method you want to 
create.

Even

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