[Proj] PROJ.4 4.9.2 and transformation of points outside datum shift grid

[Even Rouault]

Hi Bas,

Reverting 
https://github.com/OSGeo/proj.4/commit/5d2af3a8 on top of current master 
didn't change the behaviour and still resulted in the vertical grid being 
(erroneously) applied as in official 4.9.2. The reason is that cs2cs.c always 
calls pj_transform() with a non-null z pointer (z=0 if not specified), so this 
changeset has no effect on that.

Investigating more, I found the following commit to be the cause :
https://github.com/OSGeo/proj.4/commit/338ea581
" apply patch from #271 to use GTX null value 88.88880f for vertical shifts" 
(message should read -88.88880f by the way to reflect the code)

Looking at the grid (downloaded from 
https://github.com/NaturalGIS/ntv2_transformations_grids_and_sample_data/blob/master/nl/naptrans2008.gtx 
), it appears it doesn't use the official null value -88.88880f  but various big 
negative values like -2147479936, -2147479808 and several others as well.

So I think we need probably a mix of both codes for wider compatibility. 
Something like :

/* nodata?  */
/* GTX official nodata value if  -88.88880f, but some grids also use other */
/* big values for nodata (e.g naptrans2008.gtx has nodata values like
/*  -2147479936), so test them too */
if( value > 1000 || value < -1000 || value == -88.88880f )

Checking >1000 or <-1000 is somewhat arbitrary, but I'm not aware of vgrids 
that have such big offsets (at least on Earth), so it is likely safe. The 
original ticket https://github.com/OSGeo/proj.4/issues/271 doesn't report that 
they can cause problems in practice.

Howard ?

Even


Le mardi 02 février 2016 07:09:26, Sebastiaan Couwenberg a écrit :
> [Resending without attachments to keep size below threshold]
> 
> Debian Bug #813066 [0] revealed a bug in the test handling of the
> conversion between RD/NAP and ETRS89 when PROJ.4 4.9.2 is used.
> 
> Previous versions of PROJ.4 were unable to transform points outside the
> available datum shift grid, as can be seen in the output of the test
> script [1] in the file testrdtrans2008-4.9.1 [2].
> 
> With version 4.9.2 some points can now be calculated as seen in the
> file testrdtrans2008-4.9.2 [3], but for tests 07 to 10 (whose points
> "are outside the region where interpolation between either the NLGEO2004
> geoid or the RD correction grid points is possible") in the direction
> from RD/NAP to ETRS89 the results are entirely incorrect, whereas the
> results for those same tests in opposite direction are mostly correct
> except test 10 (edge_rd).
> 
> The documentation accompanying the rdtrans2008.gsb NTv2-grid and
> naptrans2008.gtx VDatum-grid [4] mentions that the points for test 07 to
> 10 are outside the correction grid and that their "RD coordinates can be
> computed, but the output should be handled with care".
> 
> With earlier versions there was no usable output, 4.9.2 reports
> incorrect values (for some tests).
> 
> It looks like the fix for #292 [5][6] is the cause of this change. I'm
> not sure whether this change is a bug or a feature.
> 
> In de debug logging for 4.9.2 we can see that for tests 07 to 10
> pj_apply_vgridshift is no longer called, instead pj_apply_gridshift is
> used. For tests 07 to 10 no z value is specified, so this looks like the
> intended course of action to not apply the vertical grid shift.
> 
> On the other hand, the previous behaviour caused pj_apply_vgridshift to
> log correctly that it failed to find a grid shift table for the location.
> 
> How should points outside datum shift grids be handled by PROJ.4?
> 
> The previous behaviour seems more correct, but with 4.9.2 some
> transformations can be done that weren't before.
> 
> [0] https://bugs.debian.org/813066
> [1]
> https://sources.debian.net/src/proj-rdnap/2008-4/debian/patches/add-test.pa
> tch/ [2] http://linuxminded.nl/tmp/testrdtrans2008-4.9.1
> [3] http://linuxminded.nl/tmp/testrdtrans2008-4.9.2
> [4]
> https://sources.debian.net/data/non-free/p/proj-rdnap/2008-4/Use%20of%20RDT
> RANS2008%20and%20NAPTRANS2008.pdf [5]
> https://github.com/OSGeo/proj.4/issues/292
> [6]
> https://github.com/OSGeo/proj.4/commit/5d2af3a89be0898458e4f0fe272affc36642
> d304
> 
> Kind Regards,
> 
> Bas

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