diff --git a/docs/source/development/reference/functions.rst b/docs/source/development/reference/functions.rst index 564135d21f..bd2f5d545e 100644 --- a/docs/source/development/reference/functions.rst +++ b/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 diff --git a/src/4D_api.cpp b/src/4D_api.cpp index d366aba82e..9445839801 100644 --- a/src/4D_api.cpp +++ b/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,6 +1920,8 @@ 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); @@ -1920,11 +1929,14 @@ PJ *proj_create_crs_to_crs_from_pj (PJ_CONTEXT *ctx, const PJ *source_crs, cons 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() ) diff --git a/src/init.cpp b/src/init.cpp index 61457cb6c5..b45c112857 100644 --- a/src/init.cpp +++ b/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; diff --git a/src/proj_internal.h b/src/proj_internal.h index 6edb6aec04..330869e002 100644 --- a/src/proj_internal.h +++ b/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{}; diff --git a/test/unit/gie_self_tests.cpp b/test/unit/gie_self_tests.cpp index 5f96a93a92..23684e68bd 100644 --- a/test/unit/gie_self_tests.cpp +++ b/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