Points falling on the exact middle of rightmost and uppermost NTv2 cell fail to project

[strk]

When trying to project a point laying exactly on the rightmost column or uppermost raw of an NTv2 gridshift, proj fails the projection complaining that the point falls outside of the grid.

I'm mass-projecting all grid cell centerpoints to produce a new NTv2 from an original one and due to this problem end up with missing informations on the sides.

[warmerdam]

For example:

gdalinfo ntf_r93.gsb
...
Corner Coordinates:
Upper Left  (  -5.5500000,  52.0500000) (  5d33' 0.00"W, 52d 3' 0.00"N)
Lower Left  (  -5.5500000,  40.9500000) (  5d33' 0.00"W, 40d57' 0.00"N)
Upper Right (  10.0500000,  52.0500000) ( 10d 3' 0.00"E, 52d 3' 0.00"N)
Lower Right (  10.0500000,  40.9500000) ( 10d 3' 0.00"E, 40d57' 0.00"N)
Center      (   2.2500000,  46.5000000) (  2d15' 0.00"E, 46d30' 0.00"N)

cs2cs -f "%.10f" +proj=latlong +nadgrids=./ntf_r93.gsb +ellps=GRS80 +to +proj=latlong +datum=WGS84
 
-5.5 52.0
pj_open_lib(./ntf_r93.gsb): call fopen(./ntf_r93.gsb) - succeeded

NTv2 FRANCE   156x111: LL=(-5.5,41) UR=(10,52)

pj_apply_gridshift(): failed to find a grid shift table for
                      location (-5.5000000dW,52.0000000dN)
   tried: ./ntf_r93.gsb
-5.5000000000	52.0000000000 0.0000000000

the problem seems to be that the code in pj_apply_gridshift_3() checks if the point falls within the bounds of the file like this:

           /* skip tables that don't match our point at all.  */
            if( ct->ll.phi > input.phi || ct->ll.lam > input.lam
                || ct->ll.phi + (ct->lim.phi-1) * ct->del.phi < input.phi
                || ct->ll.lam + (ct->lim.lam-1) * ct->del.lam < input.lam )
                continue;

Note that the ct's bounds are in terms of the transformation points (not pixels with areas) and will be in half a pixel from the extents reported by gdalinfo. If if you provide the top left corner (as done in the example) it will still fail the containment test (barring rounding effects in our favor).

Furthermore a case can be made that the grid transformation should be extended slightly outside the mesh of points - perhaps up to half a pixel - to avoid unnecessarily strict edge problems. If not that much, then at least some epsilon out from the mesh.

[warmerdam]

OK, I've made a modest adjustment in the code (r2151) which ensures that points up to some small epsilon outside the mesh will still be transformed, but I do not allow transforming further away from the mesh (ie. I don't do the up-to-half-pixel-away item I proposed initially.)

The patch includes changes to the test suite to test this near edge support.

[warmerdam]

Doh, the patch is actually r2152.

[strk]

It still fails for my usecase, which is projecting the center point of each pixel o the raster extracted by GDAL from the NTv2 file.

ERROR:  transform: couldn't project point (5.55 48.0667 0): point not within available datum shift grids (-48)

The shift grid boundary, according to gdalinfo:

Pixel Size = (0.008333333333333,-0.008333333333333)
Upper Left  (   5.5458333,  48.0708333) (  5d32'45.00"E, 48d 4'15.00"N)
Lower Right (  11.0541667,  45.4625000) ( 11d 3'15.00"E, 45d27'45.00"N)

[strk]

I think the problem in previous comment was fixed by a later commit (r2153 maybe?). Anyway, if it'll represent I'll reopen again.