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 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 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/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 == 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" 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) {