diff --git a/docs/source/development/reference/functions.rst b/docs/source/development/reference/functions.rst index a189551379..a79e1f3c25 100644 --- a/docs/source/development/reference/functions.rst +++ b/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 diff --git a/include/proj/internal/internal.hpp b/include/proj/internal/internal.hpp index dda11cbc14..d7aa0d494e 100644 --- a/include/proj/internal/internal.hpp +++ b/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; diff --git a/scripts/reference_exported_symbols.txt b/scripts/reference_exported_symbols.txt index 6ff75befbf..fc5120376f 100644 --- a/scripts/reference_exported_symbols.txt +++ b/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) diff --git a/src/4D_api.cpp b/src/4D_api.cpp index 44ba788578..58479026eb 100644 --- a/src/4D_api.cpp +++ b/src/4D_api.cpp @@ -48,6 +48,7 @@ #include "proj.h" #include "proj_experimental.h" #include "proj_internal.h" +#include /* for isnan */ #include #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::quiet_NaN(); + else if (direction == PJ_FWD) pj_fwd4d(coord, P); else pj_inv4d(coord, P); diff --git a/src/apps/gie.cpp b/src/apps/gie.cpp index ffb810ad33..2db86b1bbf 100644 --- a/src/apps/gie.cpp +++ b/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 /* for isnan */ #include #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)) diff --git a/src/apps/proj_strtod.cpp b/src/apps/proj_strtod.cpp index 441a78fac2..3a0e81f696 100644 --- a/src/apps/proj_strtod.cpp +++ b/src/apps/proj_strtod.cpp @@ -85,6 +85,8 @@ Thomas Knudsen, thokn@sdfe.dk, 2017-01-17/2017-09-18 ***********************************************************************/ +#define FROM_PROJ_CPP + #include "proj_strtod.h" #include @@ -92,7 +94,13 @@ Thomas Knudsen, thokn@sdfe.dk, 2017-01-17/2017-09-18 #include /* for HUGE_VAL */ #include /* for pow() */ #include /* for abs */ -#include /* for strchr */ +#include /* for strchr, strncmp */ + +#include + +#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(p + 3); + return std::numeric_limits::quiet_NaN(); + } + /* non-numeric? */ if (nullptr == strchr("0123456789+-._", *p)) { if (endptr) diff --git a/test/gie/more_builtins.gie b/test/gie/more_builtins.gie index 847420cef0..0b0f0706e3 100644 --- a/test/gie/more_builtins.gie +++ b/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