Add new option to proj_create_crs_to_crs_from_pj method to force +ove…

…r on transformation operations (#2914)

Fixes #2512

docs/source/development/reference/functions.rst

@@ -201,6 +201,9 @@ paragraph for more details.
     - ALLOW_BALLPARK=YES/NO: can be set to NO to disallow the use of
       :term:`Ballpark transformation` in the candidate coordinate operations.
 
+    - FORCE_OVER=YES/NO: can be set to YES to force the +over flag on the transformation
+      returned by this function.
+
 .. doxygenfunction:: proj_normalize_for_visualization
    :project: doxygen_api
 

src/4D_api.cpp

@@ -1843,6 +1843,7 @@ PJ  *proj_create_crs_to_crs_from_pj (PJ_CONTEXT *ctx, const PJ *source_crs, cons
     const char* authority = nullptr;
     double accuracy = -1;
     bool allowBallparkTransformations = true;
+    bool forceOver = false;
     for (auto iter = options; iter && iter[0]; ++iter) {
         const char *value;
         if ((value = getOptionValue(*iter, "AUTHORITY="))) {
@@ -1859,7 +1860,13 @@ PJ  *proj_create_crs_to_crs_from_pj (PJ_CONTEXT *ctx, const PJ *source_crs, cons
                             "Invalid value for ALLOW_BALLPARK option.");
                 return nullptr;
             }
-        } else {
+        }
+        else if ((value = getOptionValue(*iter, "FORCE_OVER="))) {
+            if (ci_equal(value, "yes")) {
+                forceOver = true;
+            }
+        }
+        else {
             std::string msg("Unknown option :");
             msg += *iter;
             ctx->logger(ctx->logger_app_data, PJ_LOG_ERROR, msg.c_str());
@@ -1913,18 +1920,23 @@ PJ  *proj_create_crs_to_crs_from_pj (PJ_CONTEXT *ctx, const PJ *source_crs, cons
         return nullptr;
     }
 
+    ctx->forceOver = forceOver;
+
     PJ* P = proj_list_get(ctx, op_list, 0);
     assert(P);
 
     if( P == nullptr || op_count == 1 || (area && area->bbox_set) ||
         proj_get_type(source_crs) == PJ_TYPE_GEOCENTRIC_CRS ||
         proj_get_type(target_crs) == PJ_TYPE_GEOCENTRIC_CRS ) {
         proj_list_destroy(op_list);
+        ctx->forceOver = false;
         return P;
     }
 
     auto preparedOpList = pj_create_prepared_operations(ctx, source_crs, target_crs,
                                                    op_list);
+
+    ctx->forceOver = false;
     proj_list_destroy(op_list);
 
     if( preparedOpList.empty() )

src/init.cpp

@@ -589,6 +589,9 @@ pj_init_ctx_with_allow_init_epsg(PJ_CONTEXT *ctx, int argc, char **argv, int all
 
     /* Over-ranging flag */
     PIN->over = pj_param(ctx, start, "bover").i;
+    if (ctx->forceOver) {
+        PIN->over = ctx->forceOver;
+    }
 
     /* Vertical datum geoid grids */
     PIN->has_geoid_vgrids = pj_param(ctx, start, "tgeoidgrids").i;

src/proj_internal.h

@@ -683,6 +683,7 @@ struct pj_ctx{
     void    *logger_app_data = nullptr;
     struct projCppContext* cpp_context = nullptr; /* internal context for C++ code */
     int     use_proj4_init_rules = -1; /* -1 = unknown, 0 = no, 1 = yes */
+    bool     forceOver = false; 
     int     epsg_file_exists = -1; /* -1 = unknown, 0 = no, 1 = yes */
     std::string ca_bundle_path{};
 

test/unit/gie_self_tests.cpp

@@ -890,4 +890,188 @@ TEST(gie, proj_trans_generic) {
     proj_destroy(P);
 }
 
+TEST(gie, proj_create_crs_to_crs_from_pj_force_over) {
+
+    PJ_CONTEXT* ctx;
+
+    ctx = proj_context_create();
+    ASSERT_TRUE(ctx != nullptr);
+
+    auto epsg27700 = proj_create(ctx, "EPSG:27700");
+    ASSERT_TRUE(epsg27700 != nullptr);
+
+    auto epsg4326 = proj_create(ctx, "EPSG:4326");
+    ASSERT_TRUE(epsg4326 != nullptr);
+
+    auto epsg3857 = proj_create(ctx, "EPSG:3857");
+    ASSERT_TRUE(epsg3857 != nullptr);
+
+    {
+        const char* const options[] = { "FORCE_OVER=YES", nullptr };
+        auto P = proj_create_crs_to_crs_from_pj(ctx, epsg4326, epsg3857, nullptr, options);
+        ASSERT_TRUE(P != nullptr);
+        ASSERT_TRUE(P->over);
+        PJ_COORD input;
+        PJ_COORD input_over;
+
+        //Test a point along the equator. 
+        //The same point, but in two different representations.
+        input.xyz.x = 0; // Lat in deg
+        input.xyz.y = 140;  // Long in deg
+        input.xyz.z = 0;
+
+        input_over.xyz.x = 0; // Lat in deg
+        input_over.xyz.y = -220;  // Long in deg
+        input_over.xyz.z = 0;
+
+        auto output = proj_trans(P, PJ_FWD, input);
+        auto output_over = proj_trans(P, PJ_FWD, input_over);
+
+        auto input_inv = proj_trans(P, PJ_INV, output);
+        auto input_over_inv = proj_trans(P, PJ_INV, output_over);
+
+        //Web Mercator x's between 0 and 180 longitude come out positive.
+        //But when forcing the over flag, the -220 calculation makes it flip.
+        EXPECT_GT(output.xyz.x, 0);
+        EXPECT_LT(output_over.xyz.x, 0);
+
+        EXPECT_NEAR(output.xyz.x, 15584728.711058298, 1e-8);
+        EXPECT_NEAR(output_over.xyz.x, -24490287.974520184, 1e-8);
+
+        //The distance from 140 to 180 and -220 to -180 should be pretty much the same.
+        auto dx_o = fabs(output.xyz.x - 20037508.342789244);
+        auto dx_over = fabs(output_over.xyz.x + 20037508.342789244);
+        auto dx = fabs(dx_o - dx_over);
+
+        EXPECT_NEAR(dx, 0, 1e-8);
+
+        //Check the inverse operations get us back close to our original input values.
+        EXPECT_NEAR(input.xyz.x, input_inv.xyz.x, 1e-8);
+        EXPECT_NEAR(input.xyz.y, input_inv.xyz.y, 1e-8);
+        EXPECT_NEAR(input_over.xyz.x, input_over_inv.xyz.x, 1e-8);
+        EXPECT_NEAR(input_over.xyz.y, input_over_inv.xyz.y, 1e-8);
+
+        proj_destroy(P);
+    }
+
+    {
+        //Try again with force over set to anything but YES to verify it didn't do anything.
+        const char* const options[] = { "FORCE_OVER=NO", nullptr };
+        auto P = proj_create_crs_to_crs_from_pj(ctx, epsg4326, epsg3857, nullptr, options);
+        ASSERT_TRUE(P != nullptr);
+        ASSERT_FALSE(P->over);
+        PJ_COORD input;
+        PJ_COORD input_notOver;
+
+        input.xyz.x = 0; // Lat in deg
+        input.xyz.y = 140;  // Long in deg
+        input.xyz.z = 0;
+
+        input_notOver.xyz.x = 0; // Lat in deg
+        input_notOver.xyz.y = -220;  // Long in deg
+        input_notOver.xyz.z = 0;
+
+        auto output = proj_trans(P, PJ_FWD, input);
+        auto output_notOver = proj_trans(P, PJ_FWD, input_notOver);
+
+        EXPECT_GT(output.xyz.x, 0);
+        EXPECT_GT(output_notOver.xyz.x, 0);
+
+        EXPECT_NEAR(output.xyz.x, 15584728.711058298, 1e-8);
+        EXPECT_NEAR(output_notOver.xyz.x, 15584728.711058298, 1e-8);
+
+        proj_destroy(P);
+    }
+
+    {
+        //Try again with no options to verify it didn't do anything.
+        auto P = proj_create_crs_to_crs_from_pj(ctx, epsg4326, epsg3857, nullptr, nullptr);
+        ASSERT_TRUE(P != nullptr);
+        ASSERT_FALSE(P->over);
+        PJ_COORD input;
+        PJ_COORD input_notOver;
+
+        input.xyz.x = 0; // Lat in deg
+        input.xyz.y = 140;  // Long in deg
+        input.xyz.z = 0;
+
+        input_notOver.xyz.x = 0; // Lat in deg
+        input_notOver.xyz.y = -220;  // Long in deg
+        input_notOver.xyz.z = 0;
+
+        auto output = proj_trans(P, PJ_FWD, input);
+        auto output_notOver = proj_trans(P, PJ_FWD, input_notOver);
+
+        EXPECT_GT(output.xyz.x, 0);
+        EXPECT_GT(output_notOver.xyz.x, 0);
+
+        EXPECT_NEAR(output.xyz.x, 15584728.711058298, 1e-8);
+        EXPECT_NEAR(output_notOver.xyz.x, 15584728.711058298, 1e-8);
+
+        proj_destroy(P);
+    }
+
+    {
+        //EPSG:4326 -> EPSG:27700 has more than one coordinate operation candidate.
+        const char* const options[] = { "FORCE_OVER=YES", nullptr };
+        auto P = proj_create_crs_to_crs_from_pj(ctx, epsg4326, epsg27700, nullptr, options);
+        ASSERT_TRUE(P != nullptr);
+        ASSERT_TRUE(P->over);
+        PJ_COORD input;
+        PJ_COORD input_over;
+
+        input.xyz.x = 0; // Lat in deg
+        input.xyz.y = 140;  // Long in deg
+        input.xyz.z = 0;
+
+        input_over.xyz.x = 0; // Lat in deg
+        input_over.xyz.y = -220;  // Long in deg
+        input_over.xyz.z = 0;
+
+        auto output = proj_trans(P, PJ_FWD, input);
+        auto output_over = proj_trans(P, PJ_FWD, input_over);
+
+        //Doesn't actually change the result for this tmerc transformation.
+        EXPECT_NEAR(output.xyz.x, 4980122.749364435, 1e-8);
+        EXPECT_NEAR(output.xyz.y, 14467212.882603768, 1e-8);
+        EXPECT_NEAR(output_over.xyz.x, 4980122.749364435, 1e-8);
+        EXPECT_NEAR(output_over.xyz.y, 14467212.882603768, 1e-8);
+
+        proj_destroy(P);
+    }
+
+    {
+        //Negative test for 27700.
+        const char* const options[] = { "FORCE_OVER=NO", nullptr };
+        auto P = proj_create_crs_to_crs_from_pj(ctx, epsg4326, epsg27700, nullptr, options);
+        ASSERT_TRUE(P != nullptr);
+        ASSERT_FALSE(P->over);
+        PJ_COORD input;
+        PJ_COORD input_over;
+
+        input.xyz.x = 0; // Lat in deg
+        input.xyz.y = 140;  // Long in deg
+        input.xyz.z = 0;
+
+        input_over.xyz.x = 0; // Lat in deg
+        input_over.xyz.y = -220;  // Long in deg
+        input_over.xyz.z = 0;
+
+        auto output = proj_trans(P, PJ_FWD, input);
+        auto output_over = proj_trans(P, PJ_FWD, input_over);
+
+        EXPECT_NEAR(output.xyz.x, 4980122.749364435, 1e-8);
+        EXPECT_NEAR(output.xyz.y, 14467212.882603768, 1e-8);
+        EXPECT_NEAR(output_over.xyz.x, 4980122.749364435, 1e-8);
+        EXPECT_NEAR(output_over.xyz.y, 14467212.882603768, 1e-8);
+
+        proj_destroy(P);
+    }
+
+    proj_destroy(epsg27700);
+    proj_destroy(epsg4326);
+    proj_destroy(epsg3857);
+    proj_context_destroy(ctx);
+}
+
 } // namespace