From 28f4486331ac272f5e5438950fd4ad95f5d5e828 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Tue, 20 Sep 2022 15:26:45 +0200 Subject: [PATCH 1/6] Database: create a undeprecated copy of EPSG:9123 grid transformation 'NAD83(CSRS) to CGVD28 height (1)' (fixes #3328) --- data/sql/commit.sql | 4 ++-- data/sql/customizations.sql | 44 +++++++++++++++++++++++++++++++++++++ 2 files changed, 46 insertions(+), 2 deletions(-) diff --git a/data/sql/commit.sql b/data/sql/commit.sql index b1f07886b9..e76135e768 100644 --- a/data/sql/commit.sql +++ b/data/sql/commit.sql @@ -73,14 +73,14 @@ FOR EACH ROW BEGIN AND g1.source_crs_code = g2.source_crs_code AND g1.target_crs_auth_name = g2.target_crs_auth_name AND g1.target_crs_code = g2.target_crs_code - WHERE g1.auth_name = 'PROJ' AND g1.code NOT LIKE '%_RESTRICTED_TO_VERTCRS%' AND g2.auth_name = 'EPSG') + WHERE g1.auth_name = 'PROJ' AND g1.code NOT LIKE '%_RESTRICTED_TO_VERTCRS%' AND g2.auth_name = 'EPSG' AND g2.deprecated = 0) OR EXISTS (SELECT 1 FROM grid_transformation g1 JOIN grid_transformation g2 ON g1.source_crs_auth_name = g2.target_crs_auth_name AND g1.source_crs_code = g2.target_crs_code AND g1.target_crs_auth_name = g1.source_crs_auth_name AND g1.target_crs_code = g1.source_crs_code - WHERE g1.auth_name = 'PROJ' AND g1.code NOT LIKE '%_RESTRICTED_TO_VERTCRS%' AND g2.auth_name = 'EPSG'); + WHERE g1.auth_name = 'PROJ' AND g1.code NOT LIKE '%_RESTRICTED_TO_VERTCRS%' AND g2.auth_name = 'EPSG' AND g2.deprecated = 0); SELECT RAISE(ABORT, 'Arg! there is now a EPSG:102100 object. Hack in createFromUserInput() will no longer work') WHERE EXISTS(SELECT 1 FROM crs_view WHERE auth_name = 'EPSG' AND code = '102100'); diff --git a/data/sql/customizations.sql b/data/sql/customizations.sql index 4383f4a028..b209b636a5 100644 --- a/data/sql/customizations.sql +++ b/data/sql/customizations.sql @@ -144,6 +144,50 @@ UPDATE grid_transformation SET accuracy = 2.0 WHERE auth_name = 'EPSG' AND code UPDATE grid_transformation SET accuracy = 2.0 WHERE auth_name = 'EPSG' AND code = '1462'; + +-- Create a PROJ copy of EPSG:9123 "NAD83(CSRS) to CGVD28 height (1)" by +-- removing the deprecation flag. +-- Having a transformation from/to CGVD28 using the "generic" NAD83(CSRS) is much +-- more convenient for low accuracy cases like /~https://github.com/OSGeo/PROJ/issues/3328 +INSERT INTO grid_transformation SELECT + 'PROJ' AS auth_name, + 'EPSG_9123' AS code, + name, + description || ' Imported from EPSG:9123 with deprecation flag removed by PROJ' AS description, + method_auth_name, + method_code, + method_name, + source_crs_auth_name, + source_crs_code, + target_crs_auth_name, + target_crs_code, + accuracy, + grid_param_auth_name, + grid_param_code, + grid_param_name, + grid_name, + grid2_param_auth_name, + grid2_param_code, + grid2_param_name, + grid2_name, + interpolation_crs_auth_name, + interpolation_crs_code, + operation_version, + 0 AS deprecated + FROM grid_transformation WHERE auth_name = 'EPSG' AND code = '9123'; + +INSERT INTO usage SELECT + 'PROJ' AS auth_name, + 'USAGE_PROJ_EPSG_9123' AS code, + object_table_name, + 'PROJ' AS object_auth_name, + 'EPSG_9123' AS object_code, + extent_auth_name, + extent_code, + scope_auth_name, + scope_code + FROM usage WHERE object_table_name = 'grid_transformation' AND object_auth_name = 'EPSG' AND object_code = '9123'; + -- Define the allowed authorities, and their precedence, when researching a -- coordinate operation From 44449575bd793c98ef4af9367acaba2532fded4a Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Tue, 20 Sep 2022 15:27:11 +0200 Subject: [PATCH 2/6] Database: alias HT2_1997.byn to ca_nrc_HT2_2010v70.tif (refs #3328) --- data/sql/grid_alternatives.sql | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/data/sql/grid_alternatives.sql b/data/sql/grid_alternatives.sql index 255ab25329..d93be47373 100644 --- a/data/sql/grid_alternatives.sql +++ b/data/sql/grid_alternatives.sql @@ -47,7 +47,8 @@ VALUES ('CGG2013an83.byn','ca_nrc_CGG2013an83.tif','CGG2013an83.gtx','GTiff','geoid_like',0,NULL,'https://cdn.proj.org/ca_nrc_CGG2013an83.tif',1,1,NULL), ('CGG2013i83.byn','ca_nrc_CGG2013i08.tif','CGG2013i08.gtx','GTiff','geoid_like',0,NULL,'https://cdn.proj.org/ca_nrc_CGG2013i08.tif',1,1,NULL), ('CGG2013n83.byn','ca_nrc_CGG2013n83.tif','CGG2013n83.gtx','GTiff','geoid_like',0,NULL,'https://cdn.proj.org/ca_nrc_CGG2013n83.tif',1,1,NULL), -('HT2_0.byn','ca_nrc_HT2_2010v70.tif','HT2_2010v70.gtx','GTiff','geoid_like',0,NULL,'https://cdn.proj.org/ca_nrc_HT2_2010v70.tif',1,1,NULL), -- deprecated name +('HT2_0.byn','ca_nrc_HT2_2010v70.tif','HT2_2010v70.gtx','GTiff','geoid_like',0,NULL,'https://cdn.proj.org/ca_nrc_HT2_2010v70.tif',1,1,NULL), -- deprecated name. dubious mapping !!! +('HT2_1997.byn','ca_nrc_HT2_2010v70.tif',NULL,'GTiff','geoid_like',0,NULL,'https://cdn.proj.org/ca_nrc_HT2_2010v70.tif',1,1,NULL), -- dubious mapping !!! ('HT2_2010v70.byn','ca_nrc_HT2_2010v70.tif','HT2_2010v70.gtx','GTiff','geoid_like',0,NULL,'https://cdn.proj.org/ca_nrc_HT2_2010v70.tif',1,1,NULL), ('HT2_2010v70_CGG2013a.byn','ca_nrc_HT2_2010v70_CGG2013a.tif',NULL,'GTiff','vgridshift',0,NULL,'https://cdn.proj.org/ca_nrc_HT2_2010v70_CGG2013a.tif',1,1,NULL), -- the PROJ grid is the reverse way of the EPSG one From 67017dbd036082fab4671f337f33db86d51b14a2 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Tue, 20 Sep 2022 15:28:06 +0200 Subject: [PATCH 3/6] createOperations: add missing +proj=push +v_1 +v_2 / +proj=pop +v_1 +v_2 in some vertical transformations (refs #3328) --- .../operation/coordinateoperationfactory.cpp | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/src/iso19111/operation/coordinateoperationfactory.cpp b/src/iso19111/operation/coordinateoperationfactory.cpp index ca88094533..f111f5a0a5 100644 --- a/src/iso19111/operation/coordinateoperationfactory.cpp +++ b/src/iso19111/operation/coordinateoperationfactory.cpp @@ -2099,7 +2099,17 @@ struct MyPROJStringExportableHorizVerticalHorizPROJBased final ->demoteTo2D(std::string(), nullptr) .get(), util::IComparable::Criterion::EQUIVALENT)) { - const int methodEPSGCode = transf->method()->getEPSGCode(); + int methodEPSGCode = transf->method()->getEPSGCode(); + if (methodEPSGCode == 0) { + // If the transformation is actually an inverse transformation, + // we will not get the EPSG code. So get the forward + // transformation. + const auto invTrans = transf->inverse(); + const auto invTransAsTrans = + dynamic_cast(invTrans.get()); + if (invTransAsTrans) + methodEPSGCode = invTransAsTrans->method()->getEPSGCode(); + } const bool bGeocentricTranslation = methodEPSGCode == From 435a78e9977aef0f1331a2c88cd3f5fe04fca9cc Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Tue, 20 Sep 2022 15:28:56 +0200 Subject: [PATCH 4/6] PROJ string pipeline optimizer: add optimization for EPSG:4326+5773 to EPSG:4326+5713 transformation (refs #3328) --- src/iso19111/io.cpp | 120 +++++++++++++++++++++++++++++++++++++++--- test/unit/test_io.cpp | 43 +++++++++++++++ 2 files changed, 157 insertions(+), 6 deletions(-) diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp index d80020059b..09437d3dcf 100644 --- a/src/iso19111/io.cpp +++ b/src/iso19111/io.cpp @@ -8618,7 +8618,8 @@ const std::string &PROJStringFormatter::toString() const { iterNext = std::next(iterNext); } } - if (ok && iterNext != steps.end()) { + ok &= iterNext != steps.end(); + if (ok) { ok = false; auto &nextStep = *iterNext; if (nextStep.name == "axisswap" && @@ -8628,7 +8629,8 @@ const std::string &PROJStringFormatter::toString() const { iterNext = std::next(iterNext); } } - if (ok && iterNext != steps.end()) { + ok &= iterNext != steps.end(); + if (ok) { ok = false; auto &nextStep = *iterNext; if (nextStep.name == "unitconvert" && @@ -8640,7 +8642,8 @@ const std::string &PROJStringFormatter::toString() const { } } auto iterVgridshift = iterNext; - if (ok && iterNext != steps.end()) { + ok &= iterNext != steps.end(); + if (ok) { ok = false; auto &nextStep = *iterNext; if (nextStep.name == "vgridshift") { @@ -8648,7 +8651,8 @@ const std::string &PROJStringFormatter::toString() const { iterNext = std::next(iterNext); } } - if (ok && iterNext != steps.end()) { + ok &= iterNext != steps.end(); + if (ok) { ok = false; auto &nextStep = *iterNext; if (nextStep.name == "unitconvert" && @@ -8659,7 +8663,8 @@ const std::string &PROJStringFormatter::toString() const { iterNext = std::next(iterNext); } } - if (ok && iterNext != steps.end()) { + ok &= iterNext != steps.end(); + if (ok) { ok = false; auto &nextStep = *iterNext; if (nextStep.name == "axisswap" && @@ -8669,7 +8674,8 @@ const std::string &PROJStringFormatter::toString() const { iterNext = std::next(iterNext); } } - if (ok && iterNext != steps.end()) { + ok &= iterNext != steps.end(); + if (ok) { ok = false; auto &nextStep = *iterNext; if (nextStep.name == "pop" && @@ -8677,6 +8683,7 @@ const std::string &PROJStringFormatter::toString() const { nextStep.paramValues[0].keyEquals("v_1") && nextStep.paramValues[1].keyEquals("v_2")) { ok = true; + // iterNext = std::next(iterNext); } } if (ok) { @@ -8688,6 +8695,107 @@ const std::string &PROJStringFormatter::toString() const { } } + // +step +proj=axisswap +order=2,1 + // +step +proj=unitconvert +xy_in=deg +xy_out=rad + // +step +proj=vgridshift ... + // +step +proj=unitconvert +xy_in=rad +xy_out=deg + // +step +proj=axisswap +order=2,1 + // +step +proj=push +v_1 +v_2 + // +step +proj=axisswap +order=2,1 + // +step +proj=unitconvert +xy_in=deg +xy_out=rad + // ==> + // +step +proj=push +v_1 +v_2 + // +step +proj=axisswap +order=2,1 + // +step +proj=unitconvert +xy_in=deg +xy_out=rad + // +step +proj=vgridshift ... + + if (prevStep.name == "axisswap" && prevStepParamCount == 1 && + prevStep.paramValues[0].equals("order", "2,1") && + curStep.name == "unitconvert" && curStepParamCount == 2 && + !curStep.inverted && + curStep.paramValues[0].equals("xy_in", "deg") && + curStep.paramValues[1].equals("xy_out", "rad")) { + auto iterNext = std::next(iterCur); + bool ok = false; + auto iterVgridshift = iterNext; + if (iterNext != steps.end()) { + auto &nextStep = *iterNext; + if (nextStep.name == "vgridshift") { + ok = true; + iterNext = std::next(iterNext); + } + } + ok &= iterNext != steps.end(); + if (ok) { + ok = false; + auto &nextStep = *iterNext; + if (nextStep.name == "unitconvert" && !nextStep.inverted && + nextStep.paramValues.size() == 2 && + nextStep.paramValues[0].equals("xy_in", "rad") && + nextStep.paramValues[1].equals("xy_out", "deg")) { + ok = true; + iterNext = std::next(iterNext); + } + } + ok &= iterNext != steps.end(); + if (ok) { + ok = false; + auto &nextStep = *iterNext; + if (nextStep.name == "axisswap" && + nextStep.paramValues.size() == 1 && + nextStep.paramValues[0].equals("order", "2,1")) { + ok = true; + iterNext = std::next(iterNext); + } + } + auto iterPush = iterNext; + ok &= iterNext != steps.end(); + if (ok) { + ok = false; + auto &nextStep = *iterNext; + if (nextStep.name == "push" && + nextStep.paramValues.size() == 2 && + nextStep.paramValues[0].keyEquals("v_1") && + nextStep.paramValues[1].keyEquals("v_2")) { + ok = true; + iterNext = std::next(iterNext); + } + } + ok &= iterNext != steps.end(); + if (ok) { + ok = false; + auto &nextStep = *iterNext; + if (nextStep.name == "axisswap" && + nextStep.paramValues.size() == 1 && + nextStep.paramValues[0].equals("order", "2,1")) { + ok = true; + iterNext = std::next(iterNext); + } + } + ok &= iterNext != steps.end(); + if (ok) { + ok = false; + auto &nextStep = *iterNext; + if (nextStep.name == "unitconvert" && + nextStep.paramValues.size() == 2 && + !nextStep.inverted && + nextStep.paramValues[0].equals("xy_in", "deg") && + nextStep.paramValues[1].equals("xy_out", "rad")) { + ok = true; + // iterNext = std::next(iterNext); + } + } + + if (ok) { + auto stepVgridshift = *iterVgridshift; + steps.erase(iterPrev, iterPush); + steps.insert(std::next(iterNext), stepVgridshift); + iterPrev = iterPush; + iterCur = std::next(iterPush); + continue; + } + } + ++iterCur; } } diff --git a/test/unit/test_io.cpp b/test/unit/test_io.cpp index 6af58507d2..ac88a81670 100644 --- a/test/unit/test_io.cpp +++ b/test/unit/test_io.cpp @@ -8942,6 +8942,49 @@ TEST(io, projstringformatter_optim_hgridshift_vgridshift_hgridshift_inv) { // --------------------------------------------------------------------------- +TEST(io, projstringformatter_optim_as_uc_vgridshift_uc_as_push_as_uc) { + // Nominal case + { + auto fmt = PROJStringFormatter::create(); + fmt->addStep("axisswap"); + fmt->addParam("order", "2,1"); + + fmt->addStep("unitconvert"); + fmt->addParam("xy_in", "deg"); + fmt->addParam("xy_out", "rad"); + + fmt->addStep("vgridshift"); + fmt->addParam("grids", "foo"); + + fmt->addStep("unitconvert"); + fmt->addParam("xy_in", "rad"); + fmt->addParam("xy_out", "deg"); + + fmt->addStep("axisswap"); + fmt->addParam("order", "2,1"); + + fmt->addStep("push"); + fmt->addParam("v_1"); + fmt->addParam("v_2"); + + fmt->addStep("axisswap"); + fmt->addParam("order", "2,1"); + + fmt->addStep("unitconvert"); + fmt->addParam("xy_in", "deg"); + fmt->addParam("xy_out", "rad"); + + EXPECT_EQ(fmt->toString(), + "+proj=pipeline " + "+step +proj=push +v_1 +v_2 " + "+step +proj=axisswap +order=2,1 " + "+step +proj=unitconvert +xy_in=deg +xy_out=rad " + "+step +proj=vgridshift +grids=foo"); + } +} + +// --------------------------------------------------------------------------- + TEST(io, projparse_longlat) { auto expected = "GEODCRS[\"unknown\",\n" From 966ea2b25de164680770849878a3f2a3f9ed2b84 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Tue, 20 Sep 2022 15:29:15 +0200 Subject: [PATCH 5/6] =?UTF-8?q?Tests:=20test=20EPSG:4326+5773=20to=20EPSG:?= =?UTF-8?q?4326+5713=20(refs=C2=A0#3328)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- test/unit/test_operationfactory.cpp | 57 +++++++++++++++++++++++++++++ 1 file changed, 57 insertions(+) diff --git a/test/unit/test_operationfactory.cpp b/test/unit/test_operationfactory.cpp index 4fa0ba925e..50db58011f 100644 --- a/test/unit/test_operationfactory.cpp +++ b/test/unit/test_operationfactory.cpp @@ -4991,6 +4991,63 @@ TEST(operation, compoundCRS_to_compoundCRS_issue_2720) { // --------------------------------------------------------------------------- +TEST(operation, compoundCRS_to_compoundCRS_issue_3328) { + auto authFactory = + AuthorityFactory::create(DatabaseContext::create(), std::string()); + // "WGS 84 + EGM96 height" + auto srcObj = createFromUserInput("EPSG:4326+5773", + authFactory->databaseContext(), false); + auto src = nn_dynamic_pointer_cast(srcObj); + ASSERT_TRUE(src != nullptr); + + // "WGS 84 + CGVD28 height" + auto dstObj = createFromUserInput("EPSG:4326+5713", + authFactory->databaseContext(), false); + auto dst = nn_dynamic_pointer_cast(dstObj); + ASSERT_TRUE(dst != nullptr); + + auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0); + ctxt->setGridAvailabilityUse( + CoordinateOperationContext::GridAvailabilityUse:: + IGNORE_GRID_AVAILABILITY); + ctxt->setSpatialCriterion( + CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION); + auto list = CoordinateOperationFactory::create()->createOperations( + NN_NO_CHECK(src), NN_NO_CHECK(dst), ctxt); + ASSERT_GE(list.size(), 1U); + EXPECT_EQ(list[0]->nameStr(), "Inverse of WGS 84 to EGM96 height (1) + " + "Inverse of NAD83(CSRS) to WGS 84 (2) + " + "NAD83(CSRS) to CGVD28 height (1) + " + "NAD83(CSRS) to WGS 84 (2)"); + EXPECT_EQ(list[0]->exportToPROJString(PROJStringFormatter::create().get()), + "+proj=pipeline " + "+step +proj=push +v_1 +v_2 " + "+step +proj=axisswap +order=2,1 " + "+step +proj=unitconvert +xy_in=deg +xy_out=rad " + "+step +proj=vgridshift +grids=us_nga_egm96_15.tif +multiplier=1 " + "+step +proj=cart +ellps=WGS84 " + "+step +inv +proj=helmert +x=-0.991 +y=1.9072 +z=0.5129 " + "+rx=-0.0257899075194932 " + "+ry=-0.0096500989602704 +rz=-0.0116599432323421 +s=0 " + "+convention=coordinate_frame " + "+step +inv +proj=cart +ellps=GRS80 " + "+step +inv +proj=vgridshift +grids=ca_nrc_HT2_2010v70.tif " + "+multiplier=1 " + "+step +proj=push +v_3 " + "+step +proj=cart +ellps=GRS80 " + "+step +proj=helmert +x=-0.991 +y=1.9072 +z=0.5129 " + "+rx=-0.0257899075194932 " + "+ry=-0.0096500989602704 +rz=-0.0116599432323421 +s=0 " + "+convention=coordinate_frame " + "+step +inv +proj=cart +ellps=WGS84 " + "+step +proj=pop +v_3 " + "+step +proj=unitconvert +xy_in=rad +xy_out=deg " + "+step +proj=axisswap +order=2,1 " + "+step +proj=pop +v_1 +v_2"); +} + +// --------------------------------------------------------------------------- + TEST( operation, compoundCRS_to_compoundCRS_concatenated_operation_with_two_vert_transformation_and_ballpark_geog) { From ecc13abef3b4cd9d047c0795c80e19c06a167335 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Tue, 20 Sep 2022 15:41:34 +0200 Subject: [PATCH 6/6] createOperations(): tune name and remarks of some vertical transformations --- .../operation/coordinateoperationfactory.cpp | 62 +++++++++++++++---- test/unit/test_c_api.cpp | 5 +- test/unit/test_operationfactory.cpp | 20 +++--- 3 files changed, 61 insertions(+), 26 deletions(-) diff --git a/src/iso19111/operation/coordinateoperationfactory.cpp b/src/iso19111/operation/coordinateoperationfactory.cpp index f111f5a0a5..7eaac8dcc3 100644 --- a/src/iso19111/operation/coordinateoperationfactory.cpp +++ b/src/iso19111/operation/coordinateoperationfactory.cpp @@ -2118,7 +2118,6 @@ struct MyPROJStringExportableHorizVerticalHorizPROJBased final EPSG_CODE_METHOD_GEOCENTRIC_TRANSLATION_GEOGRAPHIC_2D || methodEPSGCode == EPSG_CODE_METHOD_GEOCENTRIC_TRANSLATION_GEOGRAPHIC_3D; - if ((bGeocentricTranslation && !(transf->parameterValueNumericAsSI( EPSG_CODE_PARAMETER_X_AXIS_TRANSLATION) == 0 && @@ -2413,6 +2412,34 @@ static CoordinateOperationNNPtr createHorizVerticalHorizPROJBased( ops.emplace_back(opGeogCRStoDstCRS); } + std::vector opsForRemarks; + std::vector opsForAccuracy; + std::string opName; + if (ops.size() == 3 && opGeogCRStoDstCRS->inverse()->_isEquivalentTo( + opSrcCRSToGeogCRS.get(), + util::IComparable::Criterion::EQUIVALENT)) { + opsForRemarks.emplace_back(opSrcCRSToGeogCRS); + opsForRemarks.emplace_back(verticalTransform); + + // Only taking into account the accuracy of the vertical transform when + // opSrcCRSToGeogCRS and opGeogCRStoDstCRS are reversed and cancel + // themselves would make sense. Unfortunately it causes + // EPSG:4313+5710 (BD72 + Ostend height) to EPSG:9707 + // (WGS 84 + EGM96 height) to use a non-ideal pipeline. + // opsForAccuracy.emplace_back(verticalTransform); + opsForAccuracy = ops; + + opName = verticalTransform->nameStr() + " using "; + if (!starts_with(opSrcCRSToGeogCRS->nameStr(), "Inverse of")) + opName += opSrcCRSToGeogCRS->nameStr(); + else + opName += opGeogCRStoDstCRS->nameStr(); + } else { + opsForRemarks = ops; + opsForAccuracy = ops; + opName = computeConcatenatedName(ops); + } + bool hasBallparkTransformation = false; for (const auto &op : ops) { hasBallparkTransformation |= op->hasBallparkTransformation(); @@ -2426,21 +2453,20 @@ static CoordinateOperationNNPtr createHorizVerticalHorizPROJBased( throw InvalidOperationEmptyIntersection(msg); } auto properties = util::PropertyMap(); - properties.set(common::IdentifiedObject::NAME_KEY, - computeConcatenatedName(ops)); + properties.set(common::IdentifiedObject::NAME_KEY, opName); if (extent) { properties.set(common::ObjectUsage::DOMAIN_OF_VALIDITY_KEY, NN_NO_CHECK(extent)); } - const auto remarks = getRemarks(ops); + const auto remarks = getRemarks(opsForRemarks); if (!remarks.empty()) { properties.set(common::IdentifiedObject::REMARKS_KEY, remarks); } std::vector accuracies; - const double accuracy = getAccuracy(ops); + const double accuracy = getAccuracy(opsForAccuracy); if (accuracy >= 0.0) { accuracies.emplace_back( metadata::PositionalAccuracy::create(toString(accuracy))); @@ -5674,6 +5700,7 @@ void CoordinateOperationFactory::Private::createOperationsCompoundToCompound( const bool hasNonTrivialTargetTransf = hasNonTrivialTransf(opsGeogToTarget); double bestAccuracy = -1; + size_t bestStepCount = 0; if (hasNonTrivialSrcTransf && hasNonTrivialTargetTransf) { const auto opsGeogSrcToGeogDst = createOperations(intermGeogSrc, intermGeogDst, context); @@ -5715,10 +5742,15 @@ void CoordinateOperationFactory::Private::createOperationsCompoundToCompound( op3}, disallowEmptyIntersection)); const double accuracy = getAccuracy(res.back()); + const size_t stepCount = + getStepCount(res.back()); if (accuracy >= 0 && (bestAccuracy < 0 || - accuracy < bestAccuracy)) { + (accuracy < bestAccuracy || + (accuracy == bestAccuracy && + stepCount < bestStepCount)))) { bestAccuracy = accuracy; + bestStepCount = stepCount; } } catch (const std::exception &) { } @@ -5729,11 +5761,12 @@ void CoordinateOperationFactory::Private::createOperationsCompoundToCompound( } const auto createOpsInTwoSteps = - [&res, - bestAccuracy](const std::vector &ops1, - const std::vector &ops2) { + [&res, bestAccuracy, + bestStepCount](const std::vector &ops1, + const std::vector &ops2) { std::vector res2; double bestAccuracy2 = -1; + size_t bestStepCount2 = 0; // In first pass, exclude (horizontal) ballpark operations, but // accept them in second pass. @@ -5768,10 +5801,15 @@ void CoordinateOperationFactory::Private::createOperationsCompoundToCompound( disallowEmptyIntersection)); const double accuracy = getAccuracy(res2.back()); + const size_t stepCount = + getStepCount(res2.back()); if (accuracy >= 0 && (bestAccuracy2 < 0 || - accuracy < bestAccuracy2)) { + (accuracy < bestAccuracy2 || + (accuracy == bestAccuracy2 && + stepCount < bestStepCount2)))) { bestAccuracy2 = accuracy; + bestStepCount2 = stepCount; } } catch (const std::exception &) { } @@ -5782,7 +5820,9 @@ void CoordinateOperationFactory::Private::createOperationsCompoundToCompound( // Keep the results of this new attempt, if there are better // than the previous ones if (bestAccuracy2 >= 0 && - (bestAccuracy < 0 || bestAccuracy2 < bestAccuracy)) { + (bestAccuracy < 0 || (bestAccuracy2 < bestAccuracy || + (bestAccuracy2 == bestAccuracy && + bestStepCount2 < bestStepCount)))) { res = std::move(res2); } }; diff --git a/test/unit/test_c_api.cpp b/test/unit/test_c_api.cpp index a70c092cc6..f17a419523 100644 --- a/test/unit/test_c_api.cpp +++ b/test/unit/test_c_api.cpp @@ -5160,10 +5160,9 @@ TEST_F(CApi, proj_create_vertical_crs_ex_with_geog_crs) { ASSERT_TRUE(name != nullptr); EXPECT_EQ(name, std::string("Inverse of UTM zone 11N + " - "NAD83(2011) to WGS 84 (1) + " "Conversion from myVertCRS to myVertCRS (metre) + " - "Transformation from myVertCRS (metre) to WGS 84 + " - "Inverse of NAD83(2011) to WGS 84 (1)")); + "Transformation from myVertCRS (metre) to WGS 84 " + "using NAD83(2011) to WGS 84 (1)")); auto proj_5 = proj_as_proj_string(m_ctxt, P, PJ_PROJ_5, nullptr); ASSERT_NE(proj_5, nullptr); diff --git a/test/unit/test_operationfactory.cpp b/test/unit/test_operationfactory.cpp index 50db58011f..7f6861ed32 100644 --- a/test/unit/test_operationfactory.cpp +++ b/test/unit/test_operationfactory.cpp @@ -4741,10 +4741,9 @@ TEST(operation, compoundCRS_to_compoundCRS_WGS84_EGM96_to_WGS84_Belfast) { auto list = CoordinateOperationFactory::create()->createOperations( NN_NO_CHECK(srcCrs), NN_NO_CHECK(destCrs), ctxt); ASSERT_GE(list.size(), 1U); - EXPECT_EQ(list[0]->nameStr(), "Inverse of WGS 84 to EGM96 height (1) + " - "Inverse of ETRS89 to WGS 84 (1) + " - "ETRS89 to Belfast height (2) + " - "ETRS89 to WGS 84 (1)"); + EXPECT_EQ(list[0]->nameStr(), + "Inverse of WGS 84 to EGM96 height (1) + " + "ETRS89 to Belfast height (2) using ETRS89 to WGS 84 (1)"); EXPECT_EQ(list[0]->exportToPROJString(PROJStringFormatter::create().get()), "+proj=pipeline +step +proj=axisswap +order=2,1 " "+step +proj=unitconvert +xy_in=deg +xy_out=rad " @@ -5016,9 +5015,8 @@ TEST(operation, compoundCRS_to_compoundCRS_issue_3328) { NN_NO_CHECK(src), NN_NO_CHECK(dst), ctxt); ASSERT_GE(list.size(), 1U); EXPECT_EQ(list[0]->nameStr(), "Inverse of WGS 84 to EGM96 height (1) + " - "Inverse of NAD83(CSRS) to WGS 84 (2) + " - "NAD83(CSRS) to CGVD28 height (1) + " - "NAD83(CSRS) to WGS 84 (2)"); + "NAD83(CSRS) to CGVD28 height (1) " + "using NAD83(CSRS) to WGS 84 (2)"); EXPECT_EQ(list[0]->exportToPROJString(PROJStringFormatter::create().get()), "+proj=pipeline " "+step +proj=push +v_1 +v_2 " @@ -5073,12 +5071,10 @@ TEST( NN_NO_CHECK(src), NN_NO_CHECK(dst), ctxt); ASSERT_GE(list.size(), 1U); EXPECT_EQ(list[0]->nameStr(), - "Ballpark geographic offset from " - "NAD83(CSRS) to NAD83(CSRS)v6 + " "Inverse of NAD83(CSRS)v6 to CGVD28 height (1) + " - "NAD83(CSRS)v6 to CGVD2013(CGG2013) height (1) + " - "Ballpark geographic offset from " - "NAD83(CSRS)v6 to NAD83(CSRS)"); + "NAD83(CSRS)v6 to CGVD2013(CGG2013) height (1) " + "using Ballpark geographic offset " + "from NAD83(CSRS) to NAD83(CSRS)v6"); } #if 0 // Note: below situation is no longer triggered since EPSG v10.066 update