From 32b04ef99e8e6c747c8242f5050c3b3e321d611d Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Thu, 26 Jan 2023 18:30:37 +0100 Subject: [PATCH] pj_get_suggested_operation(): handle longitudes outside of [-180,180] for coordinate operation selection (fixes #3594) --- src/4D_api.cpp | 96 ++++++++++++++++++++++++++++++++++++ src/proj_internal.h | 23 ++++++--- test/unit/gie_self_tests.cpp | 65 ++++++++++++++++++++++++ 3 files changed, 176 insertions(+), 8 deletions(-) diff --git a/src/4D_api.cpp b/src/4D_api.cpp index d70b1d9fcb..0bb0200483 100644 --- a/src/4D_api.cpp +++ b/src/4D_api.cpp @@ -213,6 +213,19 @@ int pj_get_suggested_operation(PJ_CONTEXT *, PJ_COORD coord) /**************************************************************************************/ { + const auto normalizeLongitude = [](double x) { + if (x > 180.0) { + x -= 360.0; + if (x > 180.0) + x = fmod(x + 180.0, 360.0) - 180.0; + } else if (x < -180.0) { + x += 360.0; + if (x < -180.0) + x = fmod(x + 180.0, 360.0) - 180.0; + } + return x; + }; + // Select the operations that match the area of use // and has the best accuracy. int iBest = -1; @@ -228,11 +241,39 @@ int pj_get_suggested_operation(PJ_CONTEXT *, if (coord.xyzt.x >= alt.minxSrc && coord.xyzt.y >= alt.minySrc && coord.xyzt.x <= alt.maxxSrc && coord.xyzt.y <= alt.maxySrc) { spatialCriterionOK = true; + } else if (alt.srcIsLonLatDegree && coord.xyzt.y >= alt.minySrc && + coord.xyzt.y <= alt.maxySrc) { + const double normalizedLon = normalizeLongitude(coord.xyzt.x); + if (normalizedLon >= alt.minxSrc && + normalizedLon <= alt.maxxSrc) { + spatialCriterionOK = true; + } + } else if (alt.srcIsLatLonDegree && coord.xyzt.x >= alt.minxSrc && + coord.xyzt.x <= alt.maxxSrc) { + const double normalizedLon = normalizeLongitude(coord.xyzt.y); + if (normalizedLon >= alt.minySrc && + normalizedLon <= alt.maxySrc) { + spatialCriterionOK = true; + } } } else { if (coord.xyzt.x >= alt.minxDst && coord.xyzt.y >= alt.minyDst && coord.xyzt.x <= alt.maxxDst && coord.xyzt.y <= alt.maxyDst) { spatialCriterionOK = true; + } else if (alt.dstIsLonLatDegree && coord.xyzt.y >= alt.minyDst && + coord.xyzt.y <= alt.maxyDst) { + const double normalizedLon = normalizeLongitude(coord.xyzt.x); + if (normalizedLon >= alt.minxDst && + normalizedLon <= alt.maxxDst) { + spatialCriterionOK = true; + } + } else if (alt.dstIsLatLonDegree && coord.xyzt.x >= alt.minxDst && + coord.xyzt.x <= alt.maxxDst) { + const double normalizedLon = normalizeLongitude(coord.xyzt.y); + if (normalizedLon >= alt.minyDst && + normalizedLon <= alt.maxyDst) { + spatialCriterionOK = true; + } } } @@ -2589,3 +2630,58 @@ PJ_FACTORS proj_factors(PJ *P, PJ_COORD lp) { return factors; } + +// --------------------------------------------------------------------------- + +//! @cond Doxygen_Suppress +PJCoordOperation::PJCoordOperation(int idxInOriginalListIn, double minxSrcIn, + double minySrcIn, double maxxSrcIn, + double maxySrcIn, double minxDstIn, + double minyDstIn, double maxxDstIn, + double maxyDstIn, PJ *pjIn, + const std::string &nameIn, double accuracyIn, + bool isOffshoreIn) + : idxInOriginalList(idxInOriginalListIn), minxSrc(minxSrcIn), + minySrc(minySrcIn), maxxSrc(maxxSrcIn), maxySrc(maxySrcIn), + minxDst(minxDstIn), minyDst(minyDstIn), maxxDst(maxxDstIn), + maxyDst(maxyDstIn), pj(pjIn), name(nameIn), accuracy(accuracyIn), + isOffshore(isOffshoreIn) { + + const auto IsLonLatOrLatLon = [](const PJ *crs, bool &isLonLatDegreeOut, + bool &isLatLonDegreeOut) { + const auto eType = proj_get_type(crs); + if (eType == PJ_TYPE_GEOGRAPHIC_2D_CRS || + eType == PJ_TYPE_GEOGRAPHIC_3D_CRS) { + const auto cs = proj_crs_get_coordinate_system(crs->ctx, crs); + const char *direction = ""; + double conv_factor = 0; + constexpr double EPS = 1e-14; + if (proj_cs_get_axis_info(crs->ctx, cs, 0, nullptr, nullptr, + &direction, &conv_factor, nullptr, + nullptr, nullptr) && + ci_equal(direction, "East")) { + isLonLatDegreeOut = fabs(conv_factor - M_PI / 180) < EPS; + } else if (proj_cs_get_axis_info(crs->ctx, cs, 1, nullptr, nullptr, + &direction, &conv_factor, nullptr, + nullptr, nullptr) && + ci_equal(direction, "East")) { + isLatLonDegreeOut = fabs(conv_factor - M_PI / 180) < EPS; + } + proj_destroy(cs); + } + }; + + const auto source = proj_get_source_crs(pj->ctx, pj); + if (source) { + IsLonLatOrLatLon(source, srcIsLonLatDegree, srcIsLatLonDegree); + proj_destroy(source); + } + + const auto target = proj_get_target_crs(pj->ctx, pj); + if (target) { + IsLonLatOrLatLon(target, dstIsLonLatDegree, dstIsLatLonDegree); + proj_destroy(target); + } +} + +//! @endcond diff --git a/src/proj_internal.h b/src/proj_internal.h index 112b6876f9..473fa92434 100644 --- a/src/proj_internal.h +++ b/src/proj_internal.h @@ -312,17 +312,16 @@ struct PJCoordOperation { std::string name{}; double accuracy = -1.0; bool isOffshore = false; + bool srcIsLonLatDegree = false; + bool srcIsLatLonDegree = false; + bool dstIsLonLatDegree = false; + bool dstIsLatLonDegree = false; PJCoordOperation(int idxInOriginalListIn, double minxSrcIn, double minySrcIn, double maxxSrcIn, double maxySrcIn, double minxDstIn, double minyDstIn, double maxxDstIn, double maxyDstIn, PJ *pjIn, const std::string &nameIn, - double accuracyIn, bool isOffshoreIn) - : idxInOriginalList(idxInOriginalListIn), minxSrc(minxSrcIn), - minySrc(minySrcIn), maxxSrc(maxxSrcIn), maxySrc(maxySrcIn), - minxDst(minxDstIn), minyDst(minyDstIn), maxxDst(maxxDstIn), - maxyDst(maxyDstIn), pj(pjIn), name(nameIn), accuracy(accuracyIn), - isOffshore(isOffshoreIn) {} + double accuracyIn, bool isOffshoreIn); PJCoordOperation(const PJCoordOperation &) = delete; @@ -333,7 +332,11 @@ struct PJCoordOperation { minyDst(other.minyDst), maxxDst(other.maxxDst), maxyDst(other.maxyDst), pj(proj_clone(ctx, other.pj)), name(std::move(other.name)), accuracy(other.accuracy), - isOffshore(other.isOffshore) {} + isOffshore(other.isOffshore), + srcIsLonLatDegree(other.srcIsLonLatDegree), + srcIsLatLonDegree(other.srcIsLatLonDegree), + dstIsLonLatDegree(other.dstIsLonLatDegree), + dstIsLatLonDegree(other.dstIsLatLonDegree) {} PJCoordOperation(PJCoordOperation &&other) : idxInOriginalList(other.idxInOriginalList), minxSrc(other.minxSrc), @@ -341,7 +344,11 @@ struct PJCoordOperation { maxySrc(other.maxySrc), minxDst(other.minxDst), minyDst(other.minyDst), maxxDst(other.maxxDst), maxyDst(other.maxyDst), name(std::move(other.name)), - accuracy(other.accuracy), isOffshore(other.isOffshore) { + accuracy(other.accuracy), isOffshore(other.isOffshore), + srcIsLonLatDegree(other.srcIsLonLatDegree), + srcIsLatLonDegree(other.srcIsLatLonDegree), + dstIsLonLatDegree(other.dstIsLonLatDegree), + dstIsLatLonDegree(other.dstIsLatLonDegree) { pj = other.pj; other.pj = nullptr; } diff --git a/test/unit/gie_self_tests.cpp b/test/unit/gie_self_tests.cpp index 2a211c3eff..9c164b33a5 100644 --- a/test/unit/gie_self_tests.cpp +++ b/test/unit/gie_self_tests.cpp @@ -1041,6 +1041,71 @@ TEST(gie, proj_create_crs_to_crs_with_area_large) { // --------------------------------------------------------------------------- +TEST(gie, proj_create_crs_to_crs_with_longitude_outside_minus_180_180) { + + // Test bugfix for /~https://github.com/OSGeo/PROJ/issues/3594 + auto P = proj_create_crs_to_crs(PJ_DEFAULT_CTX, "EPSG:4277", "EPSG:4326", + nullptr); + ASSERT_TRUE(P != nullptr); + PJ_COORD c; + + c.xyzt.x = 50; // Lat in deg + c.xyzt.y = -2 + 360; // Long in deg + c.xyzt.z = 0; + c.xyzt.t = HUGE_VAL; + c = proj_trans(P, PJ_FWD, c); + EXPECT_NEAR(c.xy.x, 50.00065628, 1e-8); + EXPECT_NEAR(c.xy.y, -2.00133989, 1e-8); + + c.xyzt.x = 50; // Lat in deg + c.xyzt.y = -2 - 360; // Long in deg + c.xyzt.z = 0; + c.xyzt.t = HUGE_VAL; + c = proj_trans(P, PJ_FWD, c); + EXPECT_NEAR(c.xy.x, 50.00065628, 1e-8); + EXPECT_NEAR(c.xy.y, -2.00133989, 1e-8); + + c.xyzt.x = 50.00065628; // Lat in deg + c.xyzt.y = -2.00133989 + 360; // Long in deg + c.xyzt.z = 0; + c.xyzt.t = HUGE_VAL; + c = proj_trans(P, PJ_INV, c); + EXPECT_NEAR(c.xy.x, 50, 1e-8); + EXPECT_NEAR(c.xy.y, -2, 1e-8); + + c.xyzt.x = 50.00065628; // Lat in deg + c.xyzt.y = -2.00133989 - 360; // Long in deg + c.xyzt.z = 0; + c.xyzt.t = HUGE_VAL; + c = proj_trans(P, PJ_INV, c); + EXPECT_NEAR(c.xy.x, 50, 1e-8); + EXPECT_NEAR(c.xy.y, -2, 1e-8); + + auto Pnormalized = proj_normalize_for_visualization(PJ_DEFAULT_CTX, P); + + c.xyzt.x = -2 + 360; // Long in deg + c.xyzt.y = 50; // Lat in deg + c.xyzt.z = 0; + c.xyzt.t = HUGE_VAL; + c = proj_trans(Pnormalized, PJ_FWD, c); + EXPECT_NEAR(c.xy.x, -2.00133989, 1e-8); + EXPECT_NEAR(c.xy.y, 50.00065628, 1e-8); + + c.xyzt.x = -2.00133989 + 360; // Long in deg + c.xyzt.y = 50.00065628; // Lat in deg + c.xyzt.z = 0; + c.xyzt.t = HUGE_VAL; + c = proj_trans(Pnormalized, PJ_INV, c); + EXPECT_NEAR(c.xy.x, -2, 1e-8); + EXPECT_NEAR(c.xy.y, 50, 1e-8); + + proj_destroy(Pnormalized); + + proj_destroy(P); +} + +// --------------------------------------------------------------------------- + TEST(gie, proj_trans_generic) { // GDA2020 to WGS84 (G1762) auto P = proj_create(