Add explicit NaN handling to proj_trans and gie (#3603) (fixes #3596)

docs/source/development/reference/functions.rst

@@ -282,6 +282,13 @@ Coordinate transformation
 
     Transform a single :c:type:`PJ_COORD` coordinate.
 
+    If the input coordinate contains any NaNs you are guaranteed to get a
+    coordinate with all NaNs as a result.
+
+    .. versionchanged:: 9.2.0
+
+        Define NaN handling. Prior NaN handling behavior was undefined.
+
     :param P: Transformation object
     :type P: :c:type:`PJ` *
     :param `direction`: Transformation direction.
@@ -780,6 +787,13 @@ Various
     distance of the starting point :c:data:`coo` and the resulting
     coordinate after :c:data:`n` iterations back and forth.
 
+    If the input coordinate has any NaNs and the expected output of all NaNs
+    is returned then the final distance will be 0.
+
+    .. versionchanged:: 9.2.0
+
+        Define expected NaN distance of 0.
+
     :param P: Transformation object
     :type P: :c:type:`PJ` *
     :param `direction`: Starting direction of transformation

include/proj/internal/internal.hpp

@@ -128,7 +128,7 @@ inline bool starts_with(const std::string &str, const char *prefix) noexcept {
 
 bool ci_less(const std::string &a, const std::string &b) noexcept;
 
-bool ci_starts_with(const char *str, const char *prefix) noexcept;
+PROJ_DLL bool ci_starts_with(const char *str, const char *prefix) noexcept;
 
 bool ci_starts_with(const std::string &str, const std::string &prefix) noexcept;
 

scripts/reference_exported_symbols.txt

@@ -335,6 +335,7 @@ osgeo::proj::HorizontalShiftGridSet::reopen(pj_ctx*)
 osgeo::proj::internal::ci_equal(std::string const&, char const*)
 osgeo::proj::internal::ci_equal(std::string const&, std::string const&)
 osgeo::proj::internal::ci_find(std::string const&, char const*)
+osgeo::proj::internal::ci_starts_with(char const*, char const*)
 osgeo::proj::internal::c_locale_stod(std::string const&)
 osgeo::proj::internal::replaceAll(std::string const&, std::string const&, std::string const&)
 osgeo::proj::internal::split(std::string const&, char)

src/4D_api.cpp

@@ -48,6 +48,7 @@
 #include "proj.h"
 #include "proj_experimental.h"
 #include "proj_internal.h"
+#include <cmath> /* for isnan */
 #include <math.h>
 
 #include "proj/common.hpp"
@@ -161,6 +162,16 @@ double proj_xyz_dist(PJ_COORD a, PJ_COORD b) {
     return hypot(proj_xy_dist(a, b), a.xyz.z - b.xyz.z);
 }
 
+static bool inline coord_is_all_nans(PJ_COORD coo) {
+    return std::isnan(coo.v[0]) && std::isnan(coo.v[1]) &&
+           std::isnan(coo.v[2]) && std::isnan(coo.v[3]);
+}
+
+static bool inline coord_has_nans(PJ_COORD coo) {
+    return std::isnan(coo.v[0]) || std::isnan(coo.v[1]) ||
+           std::isnan(coo.v[2]) || std::isnan(coo.v[3]);
+}
+
 /* Measure numerical deviation after n roundtrips fwd-inv (or inv-fwd) */
 double proj_roundtrip(PJ *P, PJ_DIRECTION direction, int n, PJ_COORD *coord) {
     int i;
@@ -189,6 +200,11 @@ double proj_roundtrip(PJ *P, PJ_DIRECTION direction, int n, PJ_COORD *coord) {
     /* finally, we take the last half-step */
     t = proj_trans(P, opposite_direction(direction), t);
 
+    /* if we start with any NaN, we expect all NaN as output */
+    if (coord_has_nans(org) && coord_is_all_nans(t)) {
+        return 0.0;
+    }
+
     /* checking for angular *input* since we do a roundtrip, and end where we
      * begin */
     if (proj_angular_input(P, direction))
@@ -481,7 +497,10 @@ PJ_COORD proj_trans(PJ *P, PJ_DIRECTION direction, PJ_COORD coord) {
         0; // dummy value, to be used by proj_trans_get_last_used_operation()
     if (P->hasCoordinateEpoch)
         coord.xyzt.t = P->coordinateEpoch;
-    if (direction == PJ_FWD)
+    if (coord_has_nans(coord))
+        coord.v[0] = coord.v[1] = coord.v[2] = coord.v[3] =
+            std::numeric_limits<double>::quiet_NaN();
+    else if (direction == PJ_FWD)
         pj_fwd4d(coord, P);
     else
         pj_inv4d(coord, P);

src/apps/gie.cpp

@@ -116,6 +116,7 @@ Thomas Knudsen, thokn@sdfe.dk, 2017-10-01/2017-10-08
 #include "proj.h"
 #include "proj_internal.h"
 #include "proj_strtod.h"
+#include <cmath> /* for isnan */
 #include <math.h>
 
 #include "optargpm.h"
@@ -753,7 +754,7 @@ static PJ_COORD parse_coord(const char *args) {
         // test failures when testing points at edge of grids.
         // For example 1501000.0 becomes 1501000.000000000233
         double d = proj_strtod(prev, &endp);
-        if (*endp != '\0' && !isspace(*endp)) {
+        if (!std::isnan(d) && *endp != '\0' && !isspace(*endp)) {
             double dms = PJ_TODEG(proj_dmstor(prev, &dmsendp));
             /* TODO: When projects.h is removed, call proj_dmstor() in all cases
              */
@@ -843,7 +844,7 @@ static int roundtrip(const char *args) {
     coo = proj_angular_input(T.P, T.dir) ? torad_coord(T.P, T.dir, T.a) : T.a;
 
     r = proj_roundtrip(T.P, T.dir, ntrips, &coo);
-    if (r <= d)
+    if ((std::isnan(r) && std::isnan(d)) || r <= d)
         return another_succeeding_roundtrip();
 
     if (T.verbosity > -1) {
@@ -1052,10 +1053,13 @@ static int expect(const char *args) {
         co = proj_trans (T.P->axisswap, T.dir, co);
     }
 #endif
-    if (proj_angular_output(T.P, T.dir))
+    if (std::isnan(co.v[0]) && std::isnan(ce.v[0])) {
+        d = 0.0;
+    } else if (proj_angular_output(T.P, T.dir)) {
         d = proj_lpz_dist(T.P, ce, co);
-    else
+    } else {
         d = proj_xyz_dist(co, ce);
+    }
 
     // Test written like that to handle NaN
     if (!(d <= T.tolerance))

src/apps/proj_strtod.cpp

@@ -85,14 +85,22 @@ Thomas Knudsen, thokn@sdfe.dk, 2017-01-17/2017-09-18
 
 ***********************************************************************/
 
+#define FROM_PROJ_CPP
+
 #include "proj_strtod.h"
 
 #include <ctype.h>
 #include <errno.h>
 #include <float.h>  /* for HUGE_VAL */
 #include <math.h>   /* for pow() */
 #include <stdlib.h> /* for abs */
-#include <string.h> /* for strchr */
+#include <string.h> /* for strchr, strncmp */
+
+#include <limits>
+
+#include "proj/internal/internal.hpp"
+
+using namespace NS_PROJ::internal;
 
 double proj_strtod(const char *str, char **endptr) {
     double number = 0, integral_part = 0;
@@ -123,6 +131,13 @@ double proj_strtod(const char *str, char **endptr) {
         return 0;
     }
 
+    /* NaN */
+    if (ci_starts_with(p, "NaN")) {
+        if (endptr)
+            *endptr = const_cast<char *>(p + 3);
+        return std::numeric_limits<double>::quiet_NaN();
+    }
+
     /* non-numeric? */
     if (nullptr == strchr("0123456789+-._", *p)) {
         if (endptr)

test/gie/more_builtins.gie

@@ -787,6 +787,19 @@ accept    9.666666666666666 47.333333333333336 473.000
 expect    9.666666666666666 47.333333333333336 472.690
 roundtrip  1
 
+-------------------------------------------------------------------------------
+# Test NaN handling
+# When given NaNs, return NaNs
+-------------------------------------------------------------------------------
+-------------------------------------------------------------------------------
+operation +proj=laea +lat_0=90 +lon_0=-150 +datum=WGS84 +units=m
+-------------------------------------------------------------------------------
+direction forward
+tolerance 0
+
+accept    NaN NaN NaN NaN
+expect    NaN NaN NaN NaN
+roundtrip  1
 
 -------------------------------------------------------------------------------
 # No-op