proj_create_crs_to_crs: better deal with coordinates outside of bbox (fixes #1329)

In case several coordinate operations are returned for a CRS to CRS transformation,
we currently determine the one to use by selecting the first operation whose
bounding box contains the input point.

This commit adds a fallback case where after doing that first iteration and finding
no appropriate candidate, we try again by selecting the first operation available
that does not involve grid based transformations.

2 files changed, 55 insertions, 0 deletions

diff --git a/src/4D_api.cpp b/src/4D_api.cpp
index 4f13f238..d0b2748e 100644
--- a/src/4D_api.cpp
+++ b/src/4D_api.cpp
@@ -51,6 +51,7 @@
 #include "proj/common.hpp"
 #include "proj/coordinateoperation.hpp"
 #include "proj/internal/internal.hpp"
+#include "proj/internal/io_internal.hpp"
 
 using namespace NS_PROJ::internal;
 
@@ -192,6 +193,8 @@ similarly, but prefers the 2D resp. 3D interfaces if available.
         direction = opposite_direction(direction);
 
     if( !P->alternativeCoordinateOperations.empty() ) {
+        // Do a first pass and select the first coordinate operation whose area
+        // of use is compatible with the input coordinate
         int i = 0;
         for( const auto &alt: P->alternativeCoordinateOperations ) {
             if( direction == PJ_FWD ) {
@@ -223,6 +226,35 @@ similarly, but prefers the 2D resp. 3D interfaces if available.
             }
             i ++;
         }
+
+        // In case we did not find an operation whose area of use is compatible
+        // with the input coordinate, then goes through again the list, and
+        // use the first operation that does not require grids.
+        i = 0;
+        for( const auto &alt: P->alternativeCoordinateOperations ) {
+            auto coordOperation = dynamic_cast<
+            NS_PROJ::operation::CoordinateOperation*>(alt.pj->iso_obj.get());
+            if( coordOperation ) {
+                if( coordOperation->gridsNeeded(P->ctx->cpp_context ?
+                    P->ctx->cpp_context->databaseContext.as_nullable() :
+                    nullptr).empty() ) {
+                    if( P->iCurCoordOp != i ) {
+                        std::string msg("Using coordinate operation ");
+                        msg += alt.name;
+                        pj_log(P->ctx, PJ_LOG_TRACE, msg.c_str());
+                        P->iCurCoordOp = i;
+                    }
+                    if( direction == PJ_FWD ) {
+                        return pj_fwd4d( coord, alt.pj );
+                    }
+                    else {
+                        return pj_inv4d( coord, alt.pj );
+                    }
+                }
+            }
+            i++;
+        }
+
         proj_errno_set (P, EINVAL);
         return proj_coord_error ();
     }
diff --git a/test/unit/gie_self_tests.cpp b/test/unit/gie_self_tests.cpp
index 3f99b1b3..399f51e5 100644
--- a/test/unit/gie_self_tests.cpp
+++ b/test/unit/gie_self_tests.cpp
@@ -748,4 +748,27 @@ TEST(gie, proj_create_crs_to_crs_PULKOVO42_ETRS89) {
     proj_destroy(P);
 }
 
+// ---------------------------------------------------------------------------
+
+TEST(gie, proj_create_crs_to_crs_outside_area_of_use) {
+
+    // See https://github.com/OSGeo/proj.4/issues/1329
+    auto P = proj_create_crs_to_crs(PJ_DEFAULT_CTX, "EPSG:4275", "EPSG:4807",
+                                    nullptr);
+    ASSERT_TRUE(P != nullptr);
+    PJ_COORD c;
+
+    EXPECT_EQ(P->fwd, nullptr);
+
+    // Test point outside area of use of both candidate coordinate operations
+    c.xyz.x = 58; // Lat in deg
+    c.xyz.y = 5;  // Long in deg
+    c.xyz.z = 0;
+    c = proj_trans(P, PJ_FWD, c);
+    EXPECT_NEAR(c.xy.x, 64.44444444444444, 1e-9); // Lat in grad
+    EXPECT_NEAR(c.xy.y, 2.958634259259258, 1e-9); // Long in grad
+
+    proj_destroy(P);
+}
+
 } // namespace