Vertical grid shift: do not interpolate node values at nodata value (f…

…ixes #1002)
OSGeo · May 19, 2018 · 66bfdb4 · 66bfdb4
1 parent 39581f3
commit 66bfdb4
Show file tree

Hide file tree

Showing 4 changed files with 64 additions and 15 deletions.
diff --git a/appveyor.yml b/appveyor.yml
@@ -57,6 +57,8 @@ test_script:
   - echo test_script
   - if "%BUILD_TYPE%" == "cmake" set PROJ_LIB=C:\projects\proj-4\nad
   - if "%BUILD_TYPE%" == "nmake" set PROJ_LIB=C:\PROJ\SHARE
+  - if "%BUILD_TYPE%" == "nmake" mkdir %PROJ_LIB%\tests
+  - if "%BUILD_TYPE%" == "nmake" copy nad\tests\*.* %PROJ_LIB%\tests
   - cd %PROJ_LIB%
   - curl -O http://download.osgeo.org/proj/proj-datumgrid-1.7.zip
   - 7z e -aoa -y proj-datumgrid-1.7.zip

diff --git a/nad/tests/test_nodata.gtx b/nad/tests/test_nodata.gtx
diff --git a/src/pj_apply_vgridshift.c b/src/pj_apply_vgridshift.c
@@ -35,6 +35,15 @@
 #include "proj_internal.h"
 #include "projects.h"
 
+static int is_nodata(float value)
+{
+    /* 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 */
+    return value > 1000 || value < -1000 || value == -88.88880f;
+}
+
 static double read_vgrid_value( PJ *defn, LP input, int *gridlist_count_p, PJ_GRIDINFO **tables, struct CTABLE *ct) {
     int  itable = 0;
     double value = HUGE_VAL;
@@ -112,23 +121,49 @@ static double read_vgrid_value( PJ *defn, LP input, int *gridlist_count_p, PJ_GR
             grid_iy2 = ct->lim.phi - 1;
 
         cvs = (float *) ct->cvs;
-        value = cvs[grid_ix + grid_iy * ct->lim.lam]
-            * (1.0-grid_x) * (1.0-grid_y)
-            + cvs[grid_ix2 + grid_iy * ct->lim.lam]
-            * (grid_x) * (1.0-grid_y)
-            + cvs[grid_ix + grid_iy2 * ct->lim.lam]
-            * (1.0-grid_x) * (grid_y)
-            + cvs[grid_ix2 + grid_iy2 * ct->lim.lam]
-            * (grid_x) * (grid_y);
+        {
+            float value_a = cvs[grid_ix + grid_iy * ct->lim.lam];
+            float value_b = cvs[grid_ix2 + grid_iy * ct->lim.lam];
+            float value_c = cvs[grid_ix + grid_iy2 * ct->lim.lam];
+            float value_d = cvs[grid_ix2 + grid_iy2 * ct->lim.lam];
+            double total_weight = 0.0;
+            int n_weights = 0;
+            value = 0.0f;
+            if( !is_nodata(value_a) )
+            {
+                double weight = (1.0-grid_x) * (1.0-grid_y);
+                value += value_a * weight;
+                total_weight += weight;
+                n_weights ++;
+            }
+            if( !is_nodata(value_b) )
+            {
+                double weight = (grid_x) * (1.0-grid_y);
+                value += value_b * weight;
+                total_weight += weight;
+                n_weights ++;
+            }
+            if( !is_nodata(value_c) )
+            {
+                double weight = (1.0-grid_x) * (grid_y);
+                value += value_c * weight;
+                total_weight += weight;
+                n_weights ++;
+            }
+            if( !is_nodata(value_d) )
+            {
+                double weight = (grid_x) * (grid_y);
+                value += value_d * weight;
+                total_weight += weight;
+                n_weights ++;
+            }
+            if( n_weights == 0 )
+                value = HUGE_VAL;
+            else if( n_weights != 4 )
+                value /= total_weight;
+        }
 
     }
-    /* 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 )
-        value = HUGE_VAL;
-
 
     return value;
 }

diff --git a/test/gie/4D-API_cs2cs-style.gie b/test/gie/4D-API_cs2cs-style.gie
@@ -1,3 +1,5 @@
+
+-------------------------------------------------------------------------------
 ===============================================================================
 
 Test the 4D API handling of cs2cs style transformation options.
@@ -278,4 +280,14 @@ expect      1335.8339   7522.963
 -------------------------------------------------------------------------------
 
 
+-------------------------------------------------------------------------------
+Test bugfix of https://github.com/OSGeo/proj.4/issues/1002
+(do not interpolate nodata values)
+-------------------------------------------------------------------------------
+operation   +proj=latlong +ellps=WGS84 +geoidgrids=tests/test_nodata.gtx
+-------------------------------------------------------------------------------
+accept      4.05        52.1        0
+expect      4.05        52.1        -10
+-------------------------------------------------------------------------------
+
 </gie>