[PROJ] PROJ 6.1.0RC1

[Even Rouault]

...OK, I didn't realize that the failures where not necessarily at the end of the log file...

1) So one failure is:
""""
Test:   01 Texel
Exec:   cs2cs EPSG:4937 +to EPSG:28992+5709 -f %.4f
Input:  53.160753042 4.824761912 42.8614

Output: 117380.1203	575040.3399 0.9985
Expect: 117380.1200	575040.3400 1.0000
NAP exceeds 0.001 meters: 0.00149999999999995
Test FAILED: From ETRS89 to RD/NAP - 01 Texel
"""

I will save you the details and headaches, but the gist of the issue was that the vertical transformation
and horizontal transformation where applied in the wrong order.
There were 2 separate causes:
- the EPSG dataset doesn't register the interpolationCRS for "ETRS89 to NAP height (1)" vertical transformation (naptrans2008.gtx).
  According to the PDF you provided, this must be the Amersfoort datum. As this information is missing, the code assumes that
  the horizontal CRS to use for vertical transformation is ETRS89.
- The code didn't take properly into account the interpolationCRS for such geog3D <--> compoundCRS transformations

I've addressed both issues by the followup commit
https://github.com/OSGeo/proj.4/pull/1454/commits/ae16cff6877a39903cd049491fbb1d96a4149a96
of
https://github.com/OSGeo/proj.4/pull/1454

With that I now get

echo "53.160753042 4.824761912 42.8614" | PROJ_LIB=$PWD:data src/cs2cs -f %.15f EPSG:4937 +to EPSG:7415
117380.120258881375776	575040.339879254694097 1.000069467784584

(instead of the custom compoundCRS "EPSG:28992+5709", you can equivalently use the registered compoundCRS of the EPSG database "EPSG:7415")

2) For
"""
Test:   07 No_rd&geoid
Exec:   cs2cs EPSG:28992+5709 +to EPSG:4937 -f %.9f
Input:  100000.6700 300000.8900

Output: 50.687420379	4.608971888 0.000000000
Expect: 50.687420405	4.608971812
Latitude exceeds 1e-08 degrees: 2.59999950458223e-08
Longitude exceeds 1e-08 degrees: 7.60000000710193e-08
Test FAILED: From RD/NAP to ETRS89 - 07 No_rd&geoid (Expected output overriden: *	* ^-?(\d+\.\d+|inf)$)
"""

PROJ will indeed detect that the point is outside of the RD grid, and then it will use
the EPSG:15739 "Amersfoort to ETRS89 (3)" Helmert-based transformation. 

$ echo "100000.6700 300000.8900 0" | src/cct -d 9 +proj=pipeline +step +inv +proj=sterea +lat_0=52.1561605555556 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +step +proj=push +v_3 +step +proj=cart +ellps=bessel +step +proj=helmert +x=565.2369 +y=50.0087 +z=465.658 +rx=0.406857330322398 +ry=-0.350732676542563 +rz=1.8703473836068 +s=4.0812 +convention=coordinate_frame +step +inv +proj=cart +ellps=GRS80 +step +proj=pop +v_3 +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1
 50.687420379    4.608971888   0.000000000           inf

I'm not sure how they get their expected result. It should be noted that proj 5.2.0 fails to transform in that case:
$ echo "100000.6700 300000.8900" |  ~/proj/install-proj-5.2.0/bin/cs2cs -f %.15f +init=rdnap:rdnap +to +init=epsg:4258 
Rel. 5.2.0, September 15th, 2018
<cs2cs>: while processing file: <stdin>, line 1
pj_transform(): point not within available datum shift grids
*	* inf

Even

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