NTV2 datum shifts fail for points between 10^-5 and 10^-11 pixels outside subgrid

[proj4-bot]

Reported by bpayne on 31 Jan 2013 18:15 UTC
There is a narrow range outside each subgrid in NTV2 in which a point to be transformed is determined to be inside the subgrid, but attempting to calculate the offset within that subgrid fails. This prevents the point from being correctly interpolated within the top level grid and results in a grid shift calculation failure.

To reproduce, attempt to calculate a datum shift (NAD27-Canada-NTV2) for point: 511075.11572697730 5899877.2372716544 (UTM Zone 12 N, METERS). This will result in a failure to calculate the datum shift.

It's a little complicated to explain, so I have included a potential fix here. Essentially, in pj_apply_gridshift_3, a tolerance (1/10,000th of a subgrid pixel)is used to determine that points close to a subgrid (in the example case ALvermil) are actually within the subgrid. However, when the interpolation point is calculated, in nad_intr, only pixels within a different tolerance (1/100000000000 of a subgrid pixel) are assigned a point within the subgrid. This means that an interpolation failure results, when in fact, the higher level grid should be used to calculate the shift. The proposed solution is to adjust the tolerances so that both are approximately 10^-8 (about halfway between the existing tolerance), meaning that the nearest pixels will use the subgrid, but all others will use the main grid.

Example code modifications:

The first half of the fix (accept fewer points outside the subgrid) is in pj_apply_gridshift, approximately line 164,
double epsilon =
(fabs(ct1->del.phi)+fabs(ct1->del.lam))/10000.0;
should be changed to
double epsilon =
(fabs(ct1->del.phi)+fabs(ct1->del.lam))/100000000.0;
Note: there is one other example of the epsilon calculation around line 144. I think this should be left alone as the top level grids do require slightly different handling here (although I would be interested to hear if there are any counter arguments here).

The second half of the fix (allow more points outside the subgrid to be interpolated successfully) is in nad_intr.c, lines 19 and 32.

Change
if (indx.lam == -1 && frct.lam > 0.99999999999) {
and
if (indx.phi == -1 && frct.phi > 0.99999999999) {
to
if (indx.lam == -1 && frct.lam > 0.9999999) {
and
if (indx.phi == -1 && frct.phi > 0.9999999) {

Upon request, I can submit the entire files with updates.

Migrated-From: https://trac.osgeo.org/proj/ticket/209