ensure pj_apply_gridshift() supports points up to epsilon outside mesh (#141)

git-svn-id: http://svn.osgeo.org/metacrs/proj/trunk@2152 4e78687f-474d-0410-85f9-8d5e500ac6b2

4 files changed, 57 insertions, 7 deletions

diff --git a/ChangeLog b/ChangeLog
index 1debfba0..9cd9aa75 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,8 @@
+2012-02-01  Frank Warmerdam  <warmerdam@pobox.com>
+
+	* src/pj_apply_gridshift.c: ensure we try to use grids as long as we
+	are within epsilon of the edge (#141). 
+
 2012-01-31  Frank Warmerdam  <warmerdam@pobox.com>
 
 	* src/nad2bin.c: fix comparison test for -f flag (#139).
diff --git a/nad/td_out.dist b/nad/td_out.dist
index 540d523e..060d14d2 100644
--- a/nad/td_out.dist
+++ b/nad/td_out.dist
@@ -12,3 +12,14 @@ Test MD used where available
 79d58'00.000"W 36d58'00.000"N 0.0	79d57'59.128"W	36d58'0.501"N 0.000
 79d58'00.000"W 37d02'00.000"N 0.0	79d57'59.126"W	37d2'0.501"N 0.000
 79d58'00.000"W 36d58'00.000"N 0.0	79d57'59.128"W	36d58'0.501"N 0.000
+##############################################################
+Test that we use grid shift files even if we are right on the
+edge or even a wee bit outside (#141).
+-5.5 52.0	-5.501106465528	51.999890470284 0.000000000000
+-5.5000000000001 52.0000000000001	-5.501106465529	51.999890470284 0.000000000000
+-5.4999 51.9999	-5.501006458305	51.999790470257 0.000000000000
+-5.5001 52.0	-5.500100000000	52.000000000000 0.000000000000
+-5.5 52.0	-5.498893534472	52.000109529716 0.000000000000
+-5.5000000000001 52.0000000000001	-5.498893534472	52.000109529717 0.000000000000
+-5.4999 51.9999	-5.498793541695	52.000009529743 0.000000000000
+-5.5001 52.0	-5.500100000000	52.000000000000 0.000000000000
diff --git a/nad/testdatumfile b/nad/testdatumfile
index 0285e4cf..e03d1d8b 100755
--- a/nad/testdatumfile
+++ b/nad/testdatumfile
@@ -68,6 +68,31 @@ $EXE +proj=latlong +ellps=clrk66 +nadgrids=conus \
 79d58'00.000"W 36d58'00.000"N 0.0
 EOF
 #
+echo "##############################################################" >> ${OUT}
+echo "Test that we use grid shift files even if we are right on the" >> ${OUT}
+echo "edge or even a wee bit outside (#141)."  >> ${OUT}
+#
+# Our test points are (1) right on mesh corner, (2) outside but within 
+# epsilon (3) inside a bit (4) outside by more than epsilon
+#
+$EXE +proj=latlong +ellps=WGS84 +nadgrids=ntf_r93.gsb \
+ +to +proj=latlong +datum=WGS84 \
+ -E -f "%.12f" >>${OUT} <<EOF
+-5.5 52.0
+-5.5000000000001 52.0000000000001
+-5.4999 51.9999
+-5.5001 52.0
+EOF
+#
+$EXE +proj=latlong +datum=WGS84 \
+ +to +proj=latlong +ellps=WGS84 +nadgrids=ntf_r93.gsb \
+ -E -f "%.12f" >>${OUT} <<EOF
+-5.5 52.0
+-5.5000000000001 52.0000000000001
+-5.4999 51.9999
+-5.5001 52.0
+EOF
+#
 ##############################################################################
 # Done! 
 # do 'diff' with distribution results
diff --git a/src/pj_apply_gridshift.c b/src/pj_apply_gridshift.c
index 16f59767..366026db 100644
--- a/src/pj_apply_gridshift.c
+++ b/src/pj_apply_gridshift.c
@@ -140,11 +140,15 @@ int pj_apply_gridshift_3( projCtx ctx, PJ_GRIDINFO **tables, int grid_count,
         {
             PJ_GRIDINFO *gi = tables[itable];
             struct CTABLE *ct = gi->ct;
+            double epsilon = (fabs(ct->del.phi)+fabs(ct->del.lam))/10000.0;
 
             /* 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 )
+            if( ct->ll.phi - epsilon > input.phi 
+                || ct->ll.lam - epsilon > input.lam
+                || (ct->ll.phi + (ct->lim.phi-1) * ct->del.phi + epsilon 
+                    < input.phi)
+                || (ct->ll.lam + (ct->lim.lam-1) * ct->del.lam + epsilon 
+                    < input.lam) )
                 continue;
 
             /* If we have child nodes, check to see if any of them apply. */
@@ -155,10 +159,15 @@ int pj_apply_gridshift_3( projCtx ctx, PJ_GRIDINFO **tables, int grid_count,
                 for( child = gi->child; child != NULL; child = child->next )
                 {
                     struct CTABLE *ct1 = child->ct;
-
-                    if( ct1->ll.phi > input.phi || ct1->ll.lam > input.lam
-                      || ct1->ll.phi+(ct1->lim.phi-1)*ct1->del.phi < input.phi
-                      || ct1->ll.lam+(ct1->lim.lam-1)*ct1->del.lam < input.lam)
+                    double epsilon = 
+                        (fabs(ct1->del.phi)+fabs(ct1->del.lam))/10000.0;
+
+                    if( ct1->ll.phi - epsilon > input.phi 
+                        || ct1->ll.lam - epsilon > input.lam
+                        || (ct1->ll.phi+(ct1->lim.phi-1)*ct1->del.phi + epsilon 
+                            < input.phi)
+                        || (ct1->ll.lam+(ct1->lim.lam-1)*ct1->del.lam + epsilon 
+                            < input.lam) )
                         continue;
 
                     break;