From ad0c57436c00fdc0af0bb5e7e2c6ff163e0addfc Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Tue, 10 Sep 2019 18:24:10 +0200 Subject: C API: add proj_crs_create_bound_vertical_crs_to_WGS84() --- scripts/reference_exported_symbols.txt | 1 + src/iso19111/c_api.cpp | 46 ++++++++++++++++++++++++++++++++++ src/proj_experimental.h | 4 +++ test/unit/test_c_api.cpp | 31 +++++++++++++++++++++++ 4 files changed, 82 insertions(+) diff --git a/scripts/reference_exported_symbols.txt b/scripts/reference_exported_symbols.txt index 3206c3c8..b12907e5 100644 --- a/scripts/reference_exported_symbols.txt +++ b/scripts/reference_exported_symbols.txt @@ -900,6 +900,7 @@ proj_crs_alter_geodetic_crs proj_crs_alter_parameters_linear_unit proj_crs_create_bound_crs proj_crs_create_bound_crs_to_WGS84 +proj_crs_create_bound_vertical_crs_to_WGS84 proj_crs_create_projected_3D_crs_from_2D proj_crs_get_coordinate_system proj_crs_get_coordoperation diff --git a/src/iso19111/c_api.cpp b/src/iso19111/c_api.cpp index 3a1b66af..fd4090b2 100644 --- a/src/iso19111/c_api.cpp +++ b/src/iso19111/c_api.cpp @@ -1810,6 +1810,52 @@ PJ *proj_crs_create_bound_crs_to_WGS84(PJ_CONTEXT *ctx, const PJ *crs, // --------------------------------------------------------------------------- +/** \brief Returns a BoundCRS, with a transformation to EPSG:4979 using a grid. + * + * The returned object must be unreferenced with proj_destroy() after + * use. + * It should be used by at most one thread at a time. + * + * @param ctx PROJ context, or NULL for default context + * @param vert_crs Object of type VerticalCRS (must not be NULL) + * @param grid_name Grid name (typically a .gtx file) + * @return Object that must be unreferenced with proj_destroy(), or NULL + * in case of error. + * @since 7.0 + */ +PJ *proj_crs_create_bound_vertical_crs_to_WGS84(PJ_CONTEXT *ctx, + const PJ *vert_crs, + const char *grid_name) { + SANITIZE_CTX(ctx); + assert(vert_crs); + assert(grid_name); + auto l_crs = std::dynamic_pointer_cast(vert_crs->iso_obj); + if (!l_crs) { + proj_log_error(ctx, __FUNCTION__, "Object is not a VerticalCRS"); + return nullptr; + } + try { + auto nnCRS = NN_NO_CHECK(l_crs); + auto transformation = + Transformation::createGravityRelatedHeightToGeographic3D( + PropertyMap().set(IdentifiedObject::NAME_KEY, + "unknown to WGS84 ellipsoidal height"), + nnCRS, GeographicCRS::EPSG_4979, nullptr, + std::string(grid_name), std::vector()); + return pj_obj_create( + ctx, + BoundCRS::create(nnCRS, GeographicCRS::EPSG_4979, transformation)); + } catch (const std::exception &e) { + proj_log_error(ctx, __FUNCTION__, e.what()); + if (ctx->cpp_context) { + ctx->cpp_context->autoCloseDbIfNeeded(); + } + return nullptr; + } +} + +// --------------------------------------------------------------------------- + /** \brief Get the ellipsoid from a CRS or a GeodeticReferenceFrame. * * The returned object must be unreferenced with proj_destroy() after diff --git a/src/proj_experimental.h b/src/proj_experimental.h index 5e6bbb13..0452eb4b 100644 --- a/src/proj_experimental.h +++ b/src/proj_experimental.h @@ -313,6 +313,10 @@ PJ PROJ_DLL *proj_crs_create_bound_crs_to_WGS84(PJ_CONTEXT *ctx, const PJ *crs, const char *const *options); +PJ PROJ_DLL *proj_crs_create_bound_vertical_crs_to_WGS84(PJ_CONTEXT *ctx, + const PJ* vert_crs, + const char* grid_name); + /* BEGIN: Generated by scripts/create_c_api_projections.py*/ PJ PROJ_DLL *proj_create_conversion_utm( PJ_CONTEXT *ctx, diff --git a/test/unit/test_c_api.cpp b/test/unit/test_c_api.cpp index b8dde430..3a1c5df5 100644 --- a/test/unit/test_c_api.cpp +++ b/test/unit/test_c_api.cpp @@ -3966,4 +3966,35 @@ TEST_F(CApi, proj_crs_create_projected_3D_crs_from_2D) { } } +// --------------------------------------------------------------------------- + +TEST_F(CApi, proj_crs_create_bound_vertical_crs_to_WGS84) { + + auto vert_crs = proj_create_vertical_crs(m_ctxt, "myVertCRS", "myVertDatum", + nullptr, 0.0); + ObjectKeeper keeper_vert_crs(vert_crs); + ASSERT_NE(vert_crs, nullptr); + + auto bound_crs = proj_crs_create_bound_vertical_crs_to_WGS84( + m_ctxt, vert_crs, "foo.gtx"); + ObjectKeeper keeper_bound_crs(bound_crs); + ASSERT_NE(bound_crs, nullptr); + + auto projCRS = proj_create_from_database(m_ctxt, "EPSG", "32631", + PJ_CATEGORY_CRS, false, nullptr); + ASSERT_NE(projCRS, nullptr); + ObjectKeeper keeper_projCRS(projCRS); + + auto compound_crs = + proj_create_compound_crs(m_ctxt, "myCompoundCRS", projCRS, bound_crs); + ObjectKeeper keeper_compound_crss(compound_crs); + ASSERT_NE(compound_crs, nullptr); + + auto proj_4 = proj_as_proj_string(m_ctxt, compound_crs, PJ_PROJ_4, nullptr); + ASSERT_NE(proj_4, nullptr); + EXPECT_EQ(std::string(proj_4), + "+proj=utm +zone=31 +datum=WGS84 +units=m +geoidgrids=foo.gtx " + "+vunits=m +no_defs +type=crs"); +} + } // namespace -- cgit v1.2.3 From 890c94a730474f057f5237ca07699d6af600ed3f Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Wed, 11 Sep 2019 22:47:32 +0200 Subject: createOperations(): make sure sorting function is transitive (a < b and b < c --> a < c), to get consistent results --- src/iso19111/coordinateoperation.cpp | 50 +++++++++++++++++++++++++-------- test/unit/test_operation.cpp | 54 ++++++++++++++++++++++++++++++++++++ 2 files changed, 92 insertions(+), 12 deletions(-) diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp index 768ef76f..a1612339 100644 --- a/src/iso19111/coordinateoperation.cpp +++ b/src/iso19111/coordinateoperation.cpp @@ -56,7 +56,9 @@ #include #include -#ifdef DEBUG +// #define DEBUG +// #define DEBUG_SORT +#if defined(DEBUG) || defined(DEBUG_SORT) #include #endif @@ -10309,16 +10311,12 @@ struct SortFunction { return false; } - if (iterA->second.hasGrids_ && iterB->second.hasGrids_) { - // Operations where grids are all available go before other - if (iterA->second.gridsAvailable_ && - !iterB->second.gridsAvailable_) { - return true; - } - if (iterB->second.gridsAvailable_ && - !iterA->second.gridsAvailable_) { - return false; - } + // Operations where grids are all available go before other + if (iterA->second.gridsAvailable_ && !iterB->second.gridsAvailable_) { + return true; + } + if (iterB->second.gridsAvailable_ && !iterA->second.gridsAvailable_) { + return false; } // Operations where grids are all known in our DB go before other @@ -10680,7 +10678,35 @@ struct FilterResults { } // Sort ! - std::sort(res.begin(), res.end(), SortFunction(map)); + SortFunction sortFunc(map); + std::sort(res.begin(), res.end(), sortFunc); + +// Debug code to check consistency of the sort function +#ifdef DEBUG_SORT + constexpr bool debugSort = true; +#elif !defined(NDEBUG) + const bool debugSort = getenv("PROJ_DEBUG_SORT_FUNCT") != nullptr; +#endif +#if defined(DEBUG_SORT) || !defined(NDEBUG) + if (debugSort) { + const bool assertIfIssue = + !(getenv("PROJ_DEBUG_SORT_FUNCT_ASSERT") != nullptr); + for (size_t i = 0; i < res.size(); ++i) { + for (size_t j = i + 1; j < res.size(); ++j) { + if (sortFunc(res[j], res[i])) { +#ifdef DEBUG_SORT + std::cerr << "Sorting issue with entry " << i << "(" + << res[i]->nameStr() << ") and " << j << "(" + << res[j]->nameStr() << ")" << std::endl; +#endif + if (assertIfIssue) { + assert(false); + } + } + } + } + } +#endif } // ---------------------------------------------------------------------- diff --git a/test/unit/test_operation.cpp b/test/unit/test_operation.cpp index 21275284..e9dec3f7 100644 --- a/test/unit/test_operation.cpp +++ b/test/unit/test_operation.cpp @@ -4367,6 +4367,9 @@ TEST(operation, geogCRS_to_geogCRS_context_inverse_needed) { { auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0); + ctxt->setGridAvailabilityUse( + CoordinateOperationContext::GridAvailabilityUse:: + IGNORE_GRID_AVAILABILITY); ctxt->setUsePROJAlternativeGridNames(false); auto list = CoordinateOperationFactory::create()->createOperations( authFactory->createCoordinateReferenceSystem("4275"), // NTF @@ -4394,6 +4397,9 @@ TEST(operation, geogCRS_to_geogCRS_context_inverse_needed) { { auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0); + ctxt->setGridAvailabilityUse( + CoordinateOperationContext::GridAvailabilityUse:: + IGNORE_GRID_AVAILABILITY); auto list = CoordinateOperationFactory::create()->createOperations( authFactory->createCoordinateReferenceSystem("4275"), // NTF authFactory->createCoordinateReferenceSystem("4258"), // ETRS89 @@ -4409,6 +4415,9 @@ TEST(operation, geogCRS_to_geogCRS_context_inverse_needed) { { auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0); + ctxt->setGridAvailabilityUse( + CoordinateOperationContext::GridAvailabilityUse:: + IGNORE_GRID_AVAILABILITY); auto list = CoordinateOperationFactory::create()->createOperations( authFactory->createCoordinateReferenceSystem("4258"), // ETRS89 authFactory->createCoordinateReferenceSystem("4275"), // NTF @@ -4459,6 +4468,33 @@ TEST(operation, geogCRS_to_geogCRS_context_ntv1_ntv2_ctable2) { // --------------------------------------------------------------------------- +TEST(operation, geogCRS_to_geogCRS_context_NAD27_to_WGS84) { + auto authFactory = + AuthorityFactory::create(DatabaseContext::create(), "EPSG"); + auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0); + ctxt->setSpatialCriterion( + CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION); + ctxt->setGridAvailabilityUse( + CoordinateOperationContext::GridAvailabilityUse:: + IGNORE_GRID_AVAILABILITY); + + auto list = CoordinateOperationFactory::create()->createOperations( + authFactory->createCoordinateReferenceSystem("4267"), // NAD27 + authFactory->createCoordinateReferenceSystem("4326"), // WGS84 + ctxt); + ASSERT_EQ(list.size(), 78U); + EXPECT_EQ(list[0]->nameStr(), + "NAD27 to WGS 84 (33)"); // 1.0 m, Canada - NAD27 + EXPECT_EQ(list[1]->nameStr(), + "NAD27 to WGS 84 (3)"); // 20.0 m, Canada - NAD27 + EXPECT_EQ(list[2]->nameStr(), + "NAD27 to WGS 84 (79)"); // 5.0 m, USA - CONUS including EEZ + EXPECT_EQ(list[3]->nameStr(), + "NAD27 to WGS 84 (4)"); // 10.0 m, USA - CONUS - onshore +} + +// --------------------------------------------------------------------------- + TEST(operation, vertCRS_to_geogCRS_context) { auto authFactory = AuthorityFactory::create(DatabaseContext::create(), "EPSG"); @@ -4588,6 +4624,9 @@ TEST(operation, geogCRS_to_geogCRS_context_concatenated_operation) { auto authFactory = AuthorityFactory::create(DatabaseContext::create(), "EPSG"); auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0); + ctxt->setGridAvailabilityUse( + CoordinateOperationContext::GridAvailabilityUse:: + IGNORE_GRID_AVAILABILITY); ctxt->setAllowUseIntermediateCRS( CoordinateOperationContext::IntermediateCRSUse::ALWAYS); auto list = CoordinateOperationFactory::create()->createOperations( @@ -4615,6 +4654,9 @@ TEST(operation, geogCRS_to_geogCRS_context_same_grid_name) { auto authFactory = AuthorityFactory::create(DatabaseContext::create(), "EPSG"); auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0); + ctxt->setGridAvailabilityUse( + CoordinateOperationContext::GridAvailabilityUse:: + IGNORE_GRID_AVAILABILITY); auto list = CoordinateOperationFactory::create()->createOperations( authFactory->createCoordinateReferenceSystem("4314"), // DHDN authFactory->createCoordinateReferenceSystem("4258"), // ETRS89 @@ -5253,6 +5295,9 @@ TEST(operation, auto authFactory = AuthorityFactory::create(DatabaseContext::create(), "EPSG"); auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0); + ctxt->setGridAvailabilityUse( + CoordinateOperationContext::GridAvailabilityUse:: + IGNORE_GRID_AVAILABILITY); ctxt->setAllowUseIntermediateCRS( CoordinateOperationContext::IntermediateCRSUse::ALWAYS); auto list = CoordinateOperationFactory::create()->createOperations( @@ -6535,6 +6580,9 @@ TEST(operation, compoundCRS_to_compoundCRS_context) { auto authFactory = AuthorityFactory::create(DatabaseContext::create(), "EPSG"); 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( @@ -6713,6 +6761,9 @@ TEST(operation, compoundCRS_to_geogCRS_3D_context) { { auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0); + ctxt->setGridAvailabilityUse( + CoordinateOperationContext::GridAvailabilityUse:: + IGNORE_GRID_AVAILABILITY); auto list = CoordinateOperationFactory::create()->createOperations( authFactory->createCoordinateReferenceSystem( "7406"), // NAD27 + NGVD29 height (ftUS) @@ -6924,6 +6975,9 @@ TEST(operation, IGNF_LAMB1_TO_EPSG_4326) { auto authFactory = AuthorityFactory::create(DatabaseContext::create(), std::string()); auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0); + ctxt->setGridAvailabilityUse( + CoordinateOperationContext::GridAvailabilityUse:: + IGNORE_GRID_AVAILABILITY); ctxt->setAllowUseIntermediateCRS( CoordinateOperationContext::IntermediateCRSUse::ALWAYS); auto list = CoordinateOperationFactory::create()->createOperations( -- cgit v1.2.3 From 63857c92b271bbcd10df0a032304982011acb2a9 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Tue, 10 Sep 2019 20:21:36 +0200 Subject: Coordinate transformation: improve transformations from/to WGS84 (Gxxxx) Currently very few transformations from/to WGS84 (Gxxxx) are registered in the EPSG database, and there isn't even transformations between WGS84 EPSG:4326 and those ones. Consequently transformations to those realizations often ended up as no-operation, whereas going through WGS84 EPSG:4326 will bring more meaningful results. So register those EPSG:4326<-->WGS 84 (Gxxx) null transformations, and when having WGS 84 (Gxxx) as source/target, consider EPSG:4326 as an intermediate. This change has no effect on the existing direct transformations from/to WGS 84 (Gxxx). --- data/sql/customizations.sql | 17 +++++ data/sql/proj_db_table_defs.sql | 13 ++++ include/proj/io.hpp | 5 ++ src/iso19111/coordinateoperation.cpp | 135 ++++++++++++++++++++++++++++++++++- src/iso19111/factory.cpp | 23 ++++++ test/unit/test_operation.cpp | 67 +++++++++++++++++ 6 files changed, 259 insertions(+), 1 deletion(-) diff --git a/data/sql/customizations.sql b/data/sql/customizations.sql index 5e038de4..98633045 100644 --- a/data/sql/customizations.sql +++ b/data/sql/customizations.sql @@ -65,3 +65,20 @@ INSERT INTO "coordinate_system" VALUES('PROJ','ENh','Cartesian',3); INSERT INTO "axis" VALUES('PROJ','1','Easting','E','east','PROJ','ENh',1,'EPSG','9001'); INSERT INTO "axis" VALUES('PROJ','2','Northing','N','north','PROJ','ENh',2,'EPSG','9001'); INSERT INTO "axis" VALUES('PROJ','3','Ellipsoidal height','h','up','PROJ','ENh',2,'EPSG','9001'); + +---- Preferred hub for geodetic datum ----- + +INSERT INTO "geodetic_datum_preferred_hub" VALUES('EPSG','1152','EPSG','6326'); -- WGS84 (G730) to WGS84 +INSERT INTO "geodetic_datum_preferred_hub" VALUES('EPSG','1153','EPSG','6326'); -- WGS84 (G873) to WGS84 +INSERT INTO "geodetic_datum_preferred_hub" VALUES('EPSG','1154','EPSG','6326'); -- WGS84 (G1150) to WGS84 +INSERT INTO "geodetic_datum_preferred_hub" VALUES('EPSG','1155','EPSG','6326'); -- WGS84 (G1674) to WGS84 +INSERT INTO "geodetic_datum_preferred_hub" VALUES('EPSG','1156','EPSG','6326'); -- WGS84 (G1762) to WGS84 +INSERT INTO "geodetic_datum_preferred_hub" VALUES('EPSG','1166','EPSG','6326'); -- WGS84 (Transit) to WGS84 + +-- Consider all WGS84 related CRS are equivalent with an accuracy of 2m +INSERT INTO "helmert_transformation" VALUES('PROJ','WGS84_TO_WGS84_G730','WGS 84 to WGS 84 (G730)','','Accuracy 2m','EPSG','9603','Geocentric translations (geog2D domain)','EPSG','4326','EPSG','9053','EPSG','1262',2.0,0,0,0,'EPSG','9001',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,'',0); +INSERT INTO "helmert_transformation" VALUES('PROJ','WGS84_TO_WGS84_G873','WGS 84 to WGS 84 (G873)','','Accuracy 2m','EPSG','9603','Geocentric translations (geog2D domain)','EPSG','4326','EPSG','9054','EPSG','1262',2.0,0,0,0,'EPSG','9001',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,'',0); +INSERT INTO "helmert_transformation" VALUES('PROJ','WGS84_TO_WGS84_G1150','WGS 84 to WGS 84 (G1150)','','Accuracy 2m','EPSG','9603','Geocentric translations (geog2D domain)','EPSG','4326','EPSG','9055','EPSG','1262',2.0,0,0,0,'EPSG','9001',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,'',0); +INSERT INTO "helmert_transformation" VALUES('PROJ','WGS84_TO_WGS84_G1674','WGS 84 to WGS 84 (G1674)','','Accuracy 2m','EPSG','9603','Geocentric translations (geog2D domain)','EPSG','4326','EPSG','9056','EPSG','1262',2.0,0,0,0,'EPSG','9001',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,'',0); +INSERT INTO "helmert_transformation" VALUES('PROJ','WGS84_TO_WGS84_G1762','WGS 84 to WGS 84 (G1762)','','Accuracy 2m','EPSG','9603','Geocentric translations (geog2D domain)','EPSG','4326','EPSG','9057','EPSG','1262',2.0,0,0,0,'EPSG','9001',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,'',0); +INSERT INTO "helmert_transformation" VALUES('PROJ','WGS84_TO_WGS84_TRANSIT','WGS 84 to WGS 84 (Transit)','','Accuracy 2m','EPSG','9603','Geocentric translations (geog2D domain)','EPSG','4326','EPSG','8888','EPSG','1262',2.0,0,0,0,'EPSG','9001',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,'',0); diff --git a/data/sql/proj_db_table_defs.sql b/data/sql/proj_db_table_defs.sql index 255de827..6b38e736 100644 --- a/data/sql/proj_db_table_defs.sql +++ b/data/sql/proj_db_table_defs.sql @@ -126,6 +126,19 @@ FOR EACH ROW BEGIN WHERE EXISTS(SELECT 1 FROM area WHERE area.auth_name = NEW.area_of_use_auth_name AND area.code = NEW.area_of_use_code AND area.deprecated != 0) AND NEW.deprecated = 0; END; +-- indicates that if there is no transformation from/into (src_auth_name, src_code), +-- a research going through (hub_auth_name, hub_code) should be made +CREATE TABLE geodetic_datum_preferred_hub( + src_auth_name TEXT NOT NULL CHECK (length(src_auth_name) >= 1), + src_code TEXT NOT NULL CHECK (length(src_code) >= 1), + hub_auth_name TEXT NOT NULL CHECK (length(hub_auth_name) >= 1), + hub_code TEXT NOT NULL CHECK (length(hub_code) >= 1), + + CONSTRAINT unique_geodetic_datum_preferred_hub UNIQUE (src_auth_name, src_code, hub_auth_name, hub_code), + CONSTRAINT fk_geodetic_datum_preferred_hub_src FOREIGN KEY (src_auth_name, src_code) REFERENCES geodetic_datum(auth_name, code), + CONSTRAINT fk_geodetic_datum_preferred_hub_src FOREIGN KEY (hub_auth_name, hub_code) REFERENCES geodetic_datum(auth_name, code) +); + CREATE TABLE vertical_datum ( auth_name TEXT NOT NULL CHECK (length(auth_name) >= 1), code TEXT NOT NULL CHECK (length(code) >= 1), diff --git a/include/proj/io.hpp b/include/proj/io.hpp index e439b9ef..12b3b111 100644 --- a/include/proj/io.hpp +++ b/include/proj/io.hpp @@ -1103,6 +1103,11 @@ class PROJ_GCC_DLL AuthorityFactory { PROJ_INTERNAL crs::CRSNNPtr createCoordinateReferenceSystem(const std::string &code, bool allowCompound) const; + + PROJ_INTERNAL std::list + getPreferredHubGeodeticReferenceFrames( + const std::string &geodeticReferenceFrameCode) const; + //! @endcond protected: diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp index a1612339..f7ea385d 100644 --- a/src/iso19111/coordinateoperation.cpp +++ b/src/iso19111/coordinateoperation.cpp @@ -10171,6 +10171,7 @@ struct CoordinateOperationFactory::Private { const crs::CRSNNPtr &targetCRS; const CoordinateOperationContextNNPtr &context; bool inCreateOperationsWithDatumPivotAntiRecursion = false; + bool inCreateOperationsThroughPreferredHub = false; Context(const crs::CRSNNPtr &sourceCRSIn, const crs::CRSNNPtr &targetCRSIn, @@ -10195,6 +10196,12 @@ struct CoordinateOperationFactory::Private { const crs::GeodeticCRS *geodSrc, const crs::GeodeticCRS *geodDst, Context &context); + static void createOperationsThroughPreferredHub( + std::vector &res, + const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS, + const crs::GeodeticCRS *geodSrc, const crs::GeodeticCRS *geodDst, + Context &context); + static bool hasPerfectAccuracyResult(const std::vector &res, const Context &context); @@ -12016,6 +12023,122 @@ void CoordinateOperationFactory::Private::createOperationsWithDatumPivot( // --------------------------------------------------------------------------- +void CoordinateOperationFactory::Private::createOperationsThroughPreferredHub( + std::vector &res, const crs::CRSNNPtr &sourceCRS, + const crs::CRSNNPtr &targetCRS, const crs::GeodeticCRS *geodSrc, + const crs::GeodeticCRS *geodDst, Private::Context &context) { + + const auto &srcDatum = geodSrc->datum(); + const auto &dstDatum = geodDst->datum(); + + if (!srcDatum || !dstDatum) + return; + const auto &srcDatumIds = srcDatum->identifiers(); + const auto &dstDatumIds = dstDatum->identifiers(); + if (srcDatumIds.empty() || dstDatumIds.empty()) + return; + + const auto &authFactory = context.context->getAuthorityFactory(); + + const auto srcAuthFactory = io::AuthorityFactory::create( + authFactory->databaseContext(), *(srcDatumIds.front()->codeSpace())); + const auto srcPreferredHubs = + srcAuthFactory->getPreferredHubGeodeticReferenceFrames( + srcDatumIds.front()->code()); + + const auto dstAuthFactory = io::AuthorityFactory::create( + authFactory->databaseContext(), *(dstDatumIds.front()->codeSpace())); + const auto dstPreferredHubs = + dstAuthFactory->getPreferredHubGeodeticReferenceFrames( + dstDatumIds.front()->code()); + if (srcPreferredHubs.empty() && dstPreferredHubs.empty()) + return; + + // Currently if we have prefered hubs for both source and target, we + // will use only the one for target, arbitrarily... We could use boths + // but that would complicate things. + if (!srcPreferredHubs.empty() && dstPreferredHubs.empty()) { + std::vector resTmp; + createOperationsThroughPreferredHub(resTmp, targetCRS, sourceCRS, + geodDst, geodSrc, context); + if (!resTmp.empty()) { + resTmp = applyInverse(resTmp); + res.insert(res.end(), resTmp.begin(), resTmp.end()); + } + return; + } + assert(!dstPreferredHubs.empty()); + +#ifdef DEBUG + EnterDebugLevel enterFunction; + debugTrace("createOperationsThroughPreferredHub(" + + objectAsStr(sourceCRS.get()) + "," + + objectAsStr(targetCRS.get()) + ")"); +#endif + + struct AntiRecursionGuard { + Context &context; + + explicit AntiRecursionGuard(Context &contextIn) : context(contextIn) { + assert(!context.inCreateOperationsThroughPreferredHub); + context.inCreateOperationsThroughPreferredHub = true; + } + + ~AntiRecursionGuard() { + context.inCreateOperationsThroughPreferredHub = false; + } + }; + AntiRecursionGuard guard(context); + + std::vector candidatesIntermCRS; + for (const auto &datumHub : dstPreferredHubs) { + auto candidates = + findCandidateGeodCRSForDatum(authFactory, datumHub.get()); + bool addedGeog2D = false; + for (const auto &intermCRS : candidates) { + auto geogCRS = dynamic_cast(intermCRS.get()); + if (geogCRS && + geogCRS->coordinateSystem()->axisList().size() == 2) { + candidatesIntermCRS.emplace_back(intermCRS); + addedGeog2D = true; + break; + } + } + if (!addedGeog2D) { + candidatesIntermCRS.insert(candidatesIntermCRS.end(), + candidates.begin(), candidates.end()); + } + } + + const bool allowEmptyIntersection = true; + for (const auto &intermCRS : candidatesIntermCRS) { +#ifdef DEBUG + EnterDebugLevel loopLevel; + debugTrace("try " + objectAsStr(sourceCRS.get()) + "->" + + objectAsStr(intermCRS.get()) + "->" + + objectAsStr(targetCRS.get()) + ")"); + EnterDebugLevel loopLevel2; +#endif + const auto opsFirst = createOperations(sourceCRS, intermCRS, context); + const auto opsLast = createOperations(intermCRS, targetCRS, context); + for (const auto &opFirst : opsFirst) { + for (const auto &opLast : opsLast) { + if (!opFirst->hasBallparkTransformation() || + !opLast->hasBallparkTransformation()) { + try { + res.emplace_back( + ConcatenatedOperation::createComputeMetadata( + {opFirst, opLast}, !allowEmptyIntersection)); + } catch (const InvalidOperationEmptyIntersection &) { + } + } + } + } + } +} + +// --------------------------------------------------------------------------- + static CoordinateOperationNNPtr createBallparkGeocentricTranslation(const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS) { @@ -12262,7 +12385,17 @@ CoordinateOperationFactory::Private::createOperations( sourceCRS, targetCRS, context.context); res.insert(res.end(), resWithIntermediate.begin(), resWithIntermediate.end()); - doFilterAndCheckPerfectOp = true; + doFilterAndCheckPerfectOp = !res.empty(); + + // If transforming from/to WGS84 (Gxxxx), try through 'neutral' + // WGS84 + if (res.empty() && geodSrc && geodDst && + !context.inCreateOperationsThroughPreferredHub && + !context.inCreateOperationsWithDatumPivotAntiRecursion) { + createOperationsThroughPreferredHub( + res, sourceCRS, targetCRS, geodSrc, geodDst, context); + doFilterAndCheckPerfectOp = !res.empty(); + } } } diff --git a/src/iso19111/factory.cpp b/src/iso19111/factory.cpp index 1f40f1f0..4bf5427d 100644 --- a/src/iso19111/factory.cpp +++ b/src/iso19111/factory.cpp @@ -4981,6 +4981,29 @@ AuthorityFactory::createCompoundCRSFromExisting( // --------------------------------------------------------------------------- +//! @cond Doxygen_Suppress +std::list +AuthorityFactory::getPreferredHubGeodeticReferenceFrames( + const std::string &geodeticReferenceFrameCode) const { + std::list res; + + const std::string sql("SELECT hub_auth_name, hub_code FROM " + "geodetic_datum_preferred_hub WHERE " + "src_auth_name = ? AND src_code = ?"); + auto sqlRes = d->run(sql, {d->authority(), geodeticReferenceFrameCode}); + for (const auto &row : sqlRes) { + const auto &auth_name = row[0]; + const auto &code = row[1]; + res.emplace_back( + d->createFactory(auth_name)->createGeodeticDatum(code)); + } + + return res; +} +//! @endcond + +// --------------------------------------------------------------------------- + //! @cond Doxygen_Suppress FactoryException::FactoryException(const char *message) : Exception(message) {} diff --git a/test/unit/test_operation.cpp b/test/unit/test_operation.cpp index e9dec3f7..a3b49ba9 100644 --- a/test/unit/test_operation.cpp +++ b/test/unit/test_operation.cpp @@ -4495,6 +4495,73 @@ TEST(operation, geogCRS_to_geogCRS_context_NAD27_to_WGS84) { // --------------------------------------------------------------------------- +TEST(operation, geogCRS_to_geogCRS_context_NAD27_to_WGS84_G1762) { + auto authFactory = + AuthorityFactory::create(DatabaseContext::create(), std::string()); + auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0); + ctxt->setSpatialCriterion( + CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION); + ctxt->setGridAvailabilityUse( + CoordinateOperationContext::GridAvailabilityUse:: + IGNORE_GRID_AVAILABILITY); + + auto authFactoryEPSG = + AuthorityFactory::create(DatabaseContext::create(), "EPSG"); + auto list = CoordinateOperationFactory::create()->createOperations( + // NAD27 + authFactoryEPSG->createCoordinateReferenceSystem("4267"), + // WGS84 (G1762) + authFactoryEPSG->createCoordinateReferenceSystem("9057"), ctxt); + ASSERT_GE(list.size(), 78U); + EXPECT_EQ(list[0]->nameStr(), + "NAD27 to WGS 84 (33) + WGS 84 to WGS 84 (G1762)"); + 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 " + "+step +proj=hgridshift +grids=ntv2_0.gsb " + "+step +proj=unitconvert +xy_in=rad +xy_out=deg " + "+step +proj=axisswap +order=2,1"); + EXPECT_EQ(list[1]->nameStr(), + "NAD27 to WGS 84 (3) + WGS 84 to WGS 84 (G1762)"); + EXPECT_EQ(list[2]->nameStr(), + "NAD27 to WGS 84 (79) + WGS 84 to WGS 84 (G1762)"); + EXPECT_EQ(list[3]->nameStr(), + "NAD27 to WGS 84 (4) + WGS 84 to WGS 84 (G1762)"); +} + +// --------------------------------------------------------------------------- + +TEST(operation, geogCRS_to_geogCRS_context_WGS84_G1674_to_WGS84_G1762) { + // Check that particular behaviour with WGS 84 (Gxxx) related to + // 'geodetic_datum_preferred_hub' table and custom no-op transformations + // between WGS 84 and WGS 84 (Gxxx) doesn't affect direct transformations + // to those realizations. + auto authFactory = + AuthorityFactory::create(DatabaseContext::create(), std::string()); + auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0); + + auto authFactoryEPSG = + AuthorityFactory::create(DatabaseContext::create(), "EPSG"); + auto list = CoordinateOperationFactory::create()->createOperations( + // WGS84 (G1674) + authFactoryEPSG->createCoordinateReferenceSystem("9056"), + // WGS84 (G1762) + authFactoryEPSG->createCoordinateReferenceSystem("9057"), ctxt); + ASSERT_EQ(list.size(), 1U); + 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 " + "+step +proj=cart +ellps=WGS84 " + "+step +proj=helmert +x=-0.004 +y=0.003 +z=0.004 +rx=0.00027 " + "+ry=-0.00027 +rz=0.00038 +s=-0.0069 " + "+convention=coordinate_frame " + "+step +inv +proj=cart +ellps=WGS84 " + "+step +proj=unitconvert +xy_in=rad +xy_out=deg " + "+step +proj=axisswap +order=2,1"); +} + +// --------------------------------------------------------------------------- + TEST(operation, vertCRS_to_geogCRS_context) { auto authFactory = AuthorityFactory::create(DatabaseContext::create(), "EPSG"); -- cgit v1.2.3 From eed28e5183579d09e102d1ad72e91fc82005dfe8 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Thu, 12 Sep 2019 15:55:04 +0200 Subject: createOperations(): use more candidates when transforming between a geographic and vertical CRS For example when transforming from NAD83+NAVD88 height to WGS84, there is no transformation between NAVD88 height to WGS84. In that case, use all potential transformations from NAVD88 height to another geographic CRS for the vertical part. --- src/iso19111/coordinateoperation.cpp | 301 +++++++++++++++++++++++++++-------- src/iso19111/factory.cpp | 69 ++++---- test/unit/test_operation.cpp | 15 +- 3 files changed, 287 insertions(+), 98 deletions(-) diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp index f7ea385d..aad86410 100644 --- a/src/iso19111/coordinateoperation.cpp +++ b/src/iso19111/coordinateoperation.cpp @@ -10172,6 +10172,7 @@ struct CoordinateOperationFactory::Private { const CoordinateOperationContextNNPtr &context; bool inCreateOperationsWithDatumPivotAntiRecursion = false; bool inCreateOperationsThroughPreferredHub = false; + bool inCreateOperationsGeogToVertWithIntermediate = false; Context(const crs::CRSNNPtr &sourceCRSIn, const crs::CRSNNPtr &targetCRSIn, @@ -10202,6 +10203,11 @@ struct CoordinateOperationFactory::Private { const crs::GeodeticCRS *geodSrc, const crs::GeodeticCRS *geodDst, Context &context); + static std::vector + createOperationsGeogToVertWithIntermediate(const crs::CRSNNPtr &sourceCRS, + const crs::CRSNNPtr &targetCRS, + Context &context); + static bool hasPerfectAccuracyResult(const std::vector &res, const Context &context); @@ -10950,18 +10956,21 @@ applyInverse(const std::vector &list) { //! @cond Doxygen_Suppress -static void buildSourceAndTargetCRSIds( - const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS, - const CoordinateOperationContextNNPtr &context, - std::list> &sourceIds, - std::list> &targetIds) { +static void buildCRSIds(const crs::CRSNNPtr &crs, + const CoordinateOperationContextNNPtr &context, + std::list> &ids) { const auto &authFactory = context->getAuthorityFactory(); assert(authFactory); const auto &authFactoryName = authFactory->getAuthority(); - const auto findObjectInDB = [&authFactory, &authFactoryName]( - const crs::CRSNNPtr &crs, - std::list> &idList) { + for (const auto &id : crs->identifiers()) { + const auto &authName = *(id->codeSpace()); + const auto &code = id->code(); + if (!authName.empty()) { + ids.emplace_back(authName, code); + } + } + if (ids.empty()) { try { const auto tmpAuthFactory = io::AuthorityFactory::create( authFactory->databaseContext(), @@ -10987,35 +10996,37 @@ static void buildSourceAndTargetCRSIds( matches.front().get(), util::IComparable::Criterion::EQUIVALENT) && !matches.front()->identifiers().empty()) { - const auto &ids = matches.front()->identifiers(); - idList.emplace_back(*(ids[0]->codeSpace()), ids[0]->code()); + const auto &tmpIds = matches.front()->identifiers(); + ids.emplace_back(*(tmpIds[0]->codeSpace()), + tmpIds[0]->code()); } } } catch (const std::exception &) { } - }; - - for (const auto &id : sourceCRS->identifiers()) { - const auto &authName = *(id->codeSpace()); - const auto &code = id->code(); - if (!authName.empty()) { - sourceIds.emplace_back(authName, code); - } - } - if (sourceIds.empty()) { - findObjectInDB(sourceCRS, sourceIds); } +} - for (const auto &id : targetCRS->identifiers()) { - const auto &authName = *(id->codeSpace()); - const auto &code = id->code(); - if (!authName.empty()) { - targetIds.emplace_back(authName, code); - } +// --------------------------------------------------------------------------- + +static std::vector +getCandidateAuthorities(const io::AuthorityFactoryPtr &authFactory, + const std::string &srcAuthName, + const std::string &targetAuthName) { + const auto &authFactoryName = authFactory->getAuthority(); + std::vector authorities; + if (authFactoryName == "any") { + authorities.emplace_back(); } - if (targetIds.empty()) { - findObjectInDB(targetCRS, targetIds); + if (authFactoryName.empty()) { + authorities = authFactory->databaseContext()->getAllowedAuthorities( + srcAuthName, targetAuthName); + if (authorities.empty()) { + authorities.emplace_back(); + } + } else { + authorities.emplace_back(authFactoryName); } + return authorities; } // --------------------------------------------------------------------------- @@ -11027,12 +11038,11 @@ findOpsInRegistryDirect(const crs::CRSNNPtr &sourceCRS, const CoordinateOperationContextNNPtr &context) { const auto &authFactory = context->getAuthorityFactory(); assert(authFactory); - const auto &authFactoryName = authFactory->getAuthority(); std::list> sourceIds; std::list> targetIds; - buildSourceAndTargetCRSIds(sourceCRS, targetCRS, context, sourceIds, - targetIds); + buildCRSIds(sourceCRS, context, sourceIds); + buildCRSIds(targetCRS, context, targetIds); for (const auto &idSrc : sourceIds) { const auto &srcAuthName = idSrc.first; @@ -11041,20 +11051,8 @@ findOpsInRegistryDirect(const crs::CRSNNPtr &sourceCRS, const auto &targetAuthName = idTarget.first; const auto &targetCode = idTarget.second; - std::vector authorities; - if (authFactoryName == "any") { - authorities.emplace_back(); - } - if (authFactoryName.empty()) { - authorities = - authFactory->databaseContext()->getAllowedAuthorities( - srcAuthName, targetAuthName); - if (authorities.empty()) { - authorities.emplace_back(); - } - } else { - authorities.emplace_back(authFactoryName); - } + const auto authorities(getCandidateAuthorities( + authFactory, srcAuthName, targetAuthName)); for (const auto &authority : authorities) { const auto tmpAuthFactory = io::AuthorityFactory::create( authFactory->databaseContext(), @@ -11075,6 +11073,81 @@ findOpsInRegistryDirect(const crs::CRSNNPtr &sourceCRS, } return std::vector(); } + +// --------------------------------------------------------------------------- + +// Look in the authority registry for operations from sourceCRS +static std::vector +findOpsInRegistryDirectFrom(const crs::CRSNNPtr &sourceCRS, + const CoordinateOperationContextNNPtr &context) { + const auto &authFactory = context->getAuthorityFactory(); + assert(authFactory); + + std::list> ids; + buildCRSIds(sourceCRS, context, ids); + + for (const auto &id : ids) { + const auto &srcAuthName = id.first; + const auto &srcCode = id.second; + + const auto authorities( + getCandidateAuthorities(authFactory, srcAuthName, srcAuthName)); + for (const auto &authority : authorities) { + const auto tmpAuthFactory = io::AuthorityFactory::create( + authFactory->databaseContext(), + authority == "any" ? std::string() : authority); + auto res = tmpAuthFactory->createFromCoordinateReferenceSystemCodes( + srcAuthName, srcCode, std::string(), std::string(), + context->getUsePROJAlternativeGridNames(), + context->getGridAvailabilityUse() == + CoordinateOperationContext::GridAvailabilityUse:: + DISCARD_OPERATION_IF_MISSING_GRID, + context->getDiscardSuperseded()); + if (!res.empty()) { + return res; + } + } + } + return std::vector(); +} + +// --------------------------------------------------------------------------- + +// Look in the authority registry for operations to targetCRS +static std::vector +findOpsInRegistryDirectTo(const crs::CRSNNPtr &targetCRS, + const CoordinateOperationContextNNPtr &context) { + const auto &authFactory = context->getAuthorityFactory(); + assert(authFactory); + + std::list> ids; + buildCRSIds(targetCRS, context, ids); + + for (const auto &id : ids) { + const auto &targetAuthName = id.first; + const auto &targetCode = id.second; + + const auto authorities(getCandidateAuthorities( + authFactory, targetAuthName, targetAuthName)); + for (const auto &authority : authorities) { + const auto tmpAuthFactory = io::AuthorityFactory::create( + authFactory->databaseContext(), + authority == "any" ? std::string() : authority); + auto res = tmpAuthFactory->createFromCoordinateReferenceSystemCodes( + std::string(), std::string(), targetAuthName, targetCode, + context->getUsePROJAlternativeGridNames(), + context->getGridAvailabilityUse() == + CoordinateOperationContext::GridAvailabilityUse:: + DISCARD_OPERATION_IF_MISSING_GRID, + context->getDiscardSuperseded()); + if (!res.empty()) { + return res; + } + } + } + return std::vector(); +} + //! @endcond // --------------------------------------------------------------------------- @@ -11089,12 +11162,11 @@ static std::vector findsOpsInRegistryWithIntermediate( const auto &authFactory = context->getAuthorityFactory(); assert(authFactory); - const auto &authFactoryName = authFactory->getAuthority(); std::list> sourceIds; std::list> targetIds; - buildSourceAndTargetCRSIds(sourceCRS, targetCRS, context, sourceIds, - targetIds); + buildCRSIds(sourceCRS, context, sourceIds); + buildCRSIds(targetCRS, context, targetIds); for (const auto &idSrc : sourceIds) { const auto &srcAuthName = idSrc.first; @@ -11103,20 +11175,8 @@ static std::vector findsOpsInRegistryWithIntermediate( const auto &targetAuthName = idTarget.first; const auto &targetCode = idTarget.second; - std::vector authorities; - if (authFactoryName == "any") { - authorities.emplace_back(); - } - if (authFactoryName.empty()) { - authorities = - authFactory->databaseContext()->getAllowedAuthorities( - srcAuthName, targetAuthName); - if (authorities.empty()) { - authorities.emplace_back(); - } - } else { - authorities.emplace_back(authFactoryName); - } + const auto authorities(getCandidateAuthorities( + authFactory, srcAuthName, targetAuthName)); for (const auto &authority : authorities) { const auto tmpAuthFactory = io::AuthorityFactory::create( authFactory->databaseContext(), @@ -12139,6 +12199,51 @@ void CoordinateOperationFactory::Private::createOperationsThroughPreferredHub( // --------------------------------------------------------------------------- +std::vector +CoordinateOperationFactory::Private::createOperationsGeogToVertWithIntermediate( + const crs::CRSNNPtr & /*sourceCRS*/, // geographic CRS + const crs::CRSNNPtr &targetCRS, // vertical CRS + Private::Context &context) { + + std::vector res; + + struct AntiRecursionGuard { + Context &context; + + explicit AntiRecursionGuard(Context &contextIn) : context(contextIn) { + assert(!context.inCreateOperationsGeogToVertWithIntermediate); + context.inCreateOperationsGeogToVertWithIntermediate = true; + } + + ~AntiRecursionGuard() { + context.inCreateOperationsGeogToVertWithIntermediate = false; + } + }; + AntiRecursionGuard guard(context); + + for (int i = 0; i < 2; i++) { + + // Generally EPSG has operations from GeogCrs to VertCRS + auto ops = + i == 0 ? findOpsInRegistryDirectTo(targetCRS, context.context) + : findOpsInRegistryDirectFrom(targetCRS, context.context); + + for (const auto &op : ops) { + const auto tmpCRS = i == 0 ? op->sourceCRS() : op->targetCRS(); + if (tmpCRS && + dynamic_cast(tmpCRS.get())) { + res.emplace_back(i == 0 ? op : op->inverse()); + } + } + if (!res.empty()) + break; + } + + return res; +} + +// --------------------------------------------------------------------------- + static CoordinateOperationNNPtr createBallparkGeocentricTranslation(const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS) { @@ -12350,6 +12455,21 @@ CoordinateOperationFactory::Private::createOperations( } } + // There's no direct transformation from NAVD88 height to WGS84, + // so try to research all transformations from NAVD88 to another + // intermediate GeographicCRS. + if (res.empty() && + !context.inCreateOperationsGeogToVertWithIntermediate && + geogSrc && vertDst) { + res = createOperationsGeogToVertWithIntermediate( + sourceCRS, targetCRS, context); + } else if (res.empty() && + !context.inCreateOperationsGeogToVertWithIntermediate && + geogDst && vertSrc) { + res = applyInverse(createOperationsGeogToVertWithIntermediate( + targetCRS, sourceCRS, context)); + } + if (res.empty() && !sameGeodeticDatum && !context.inCreateOperationsWithDatumPivotAntiRecursion && geodSrc && geodDst) { @@ -12908,6 +13028,61 @@ CoordinateOperationFactory::Private::createOperations( componentsSrc[1]->extractVerticalCRS()) { verticalTransforms = createOperations(componentsSrc[1], targetCRS, context); + bool foundRegisteredTransformWithAllGridsAvailable = false; + for (const auto &op : verticalTransforms) { + if (!op->identifiers().empty() && authFactory) { + bool missingGrid = false; + const auto gridsNeeded = + op->gridsNeeded(authFactory->databaseContext()); + for (const auto &gridDesc : gridsNeeded) { + if (!gridDesc.available) { + missingGrid = true; + break; + } + } + if (!missingGrid) { + foundRegisteredTransformWithAllGridsAvailable = + true; + break; + } + } + } + if (!foundRegisteredTransformWithAllGridsAvailable && + srcGeogCRS && + !srcGeogCRS->_isEquivalentTo( + geogDst, util::IComparable::Criterion::EQUIVALENT) && + !srcGeogCRS->is2DPartOf3D(NN_NO_CHECK(geogDst))) { + auto verticalTransformsTmp = createOperations( + componentsSrc[1], NN_NO_CHECK(srcGeogCRS), context); + bool foundRegisteredTransform = false; + foundRegisteredTransformWithAllGridsAvailable = false; + for (const auto &op : verticalTransformsTmp) { + if (!op->identifiers().empty() && authFactory) { + bool missingGrid = false; + const auto gridsNeeded = + op->gridsNeeded(authFactory->databaseContext()); + for (const auto &gridDesc : gridsNeeded) { + if (!gridDesc.available) { + missingGrid = true; + break; + } + } + foundRegisteredTransform = true; + if (!missingGrid) { + foundRegisteredTransformWithAllGridsAvailable = + true; + break; + } + } + } + if (foundRegisteredTransformWithAllGridsAvailable) { + verticalTransforms = verticalTransformsTmp; + } else if (foundRegisteredTransform) { + verticalTransforms.insert(verticalTransforms.end(), + verticalTransformsTmp.begin(), + verticalTransformsTmp.end()); + } + } } if (!horizTransforms.empty() && !verticalTransforms.empty()) { for (const auto &horizTransform : horizTransforms) { diff --git a/src/iso19111/factory.cpp b/src/iso19111/factory.cpp index 4bf5427d..3ed69ae9 100644 --- a/src/iso19111/factory.cpp +++ b/src/iso19111/factory.cpp @@ -3312,9 +3312,9 @@ AuthorityFactory::createFromCoordinateReferenceSystemCodes( bool discardSuperseded) const { auto cacheKey(d->authority()); - cacheKey += sourceCRSAuthName; + cacheKey += sourceCRSAuthName.empty() ? "{empty}" : sourceCRSAuthName; cacheKey += sourceCRSCode; - cacheKey += targetCRSAuthName; + cacheKey += targetCRSAuthName.empty() ? "{empty}" : targetCRSAuthName; cacheKey += targetCRSCode; cacheKey += (usePROJAlternativeGridNames ? '1' : '0'); cacheKey += (discardIfMissingGrid ? '1' : '0'); @@ -3326,27 +3326,30 @@ AuthorityFactory::createFromCoordinateReferenceSystemCodes( return list; } - // Look-up first for conversion which is the most precise. - std::string sql("SELECT conversion_auth_name, " - "geodetic_crs_auth_name, geodetic_crs_code FROM " - "projected_crs WHERE auth_name = ? AND code = ?"); - auto params = ListOfParams{targetCRSAuthName, targetCRSCode}; - auto res = d->run(sql, params); - if (!res.empty()) { - const auto &row = res.front(); - bool ok = row[1] == sourceCRSAuthName && row[2] == sourceCRSCode; - if (ok && d->hasAuthorityRestriction()) { - ok = row[0] == d->authority(); - } - if (ok) { - auto targetCRS = d->createFactory(targetCRSAuthName) - ->createProjectedCRS(targetCRSCode); - auto conv = targetCRS->derivingConversion(); - list.emplace_back(conv); - d->context()->d->cache(cacheKey, list); - return list; + if (!targetCRSAuthName.empty()) { + // Look-up first for conversion which is the most precise. + std::string sql("SELECT conversion_auth_name, " + "geodetic_crs_auth_name, geodetic_crs_code FROM " + "projected_crs WHERE auth_name = ? AND code = ?"); + auto params = ListOfParams{targetCRSAuthName, targetCRSCode}; + auto res = d->run(sql, params); + if (!res.empty()) { + const auto &row = res.front(); + bool ok = row[1] == sourceCRSAuthName && row[2] == sourceCRSCode; + if (ok && d->hasAuthorityRestriction()) { + ok = row[0] == d->authority(); + } + if (ok) { + auto targetCRS = d->createFactory(targetCRSAuthName) + ->createProjectedCRS(targetCRSCode); + auto conv = targetCRS->derivingConversion(); + list.emplace_back(conv); + d->context()->d->cache(cacheKey, list); + return list; + } } } + std::string sql; if (discardSuperseded) { sql = "SELECT cov.auth_name, cov.code, cov.table_name, " "ss.replacement_auth_name, ss.replacement_code FROM " @@ -3358,20 +3361,26 @@ AuthorityFactory::createFromCoordinateReferenceSystemCodes( "ss.superseded_auth_name = cov.auth_name AND " "ss.superseded_code = cov.code AND " "ss.superseded_table_name = ss.replacement_table_name " - "WHERE source_crs_auth_name = ? AND source_crs_code = ? AND " - "target_crs_auth_name = ? AND target_crs_code = ? AND " - "cov.deprecated = 0"; + "WHERE "; } else { sql = "SELECT cov.auth_name, cov.code, cov.table_name FROM " "coordinate_operation_view cov JOIN area " "ON cov.area_of_use_auth_name = area.auth_name AND " "cov.area_of_use_code = area.code " - "WHERE source_crs_auth_name = ? AND source_crs_code = ? AND " - "target_crs_auth_name = ? AND target_crs_code = ? AND " - "cov.deprecated = 0"; + "WHERE "; + } + ListOfParams params; + if (!sourceCRSAuthName.empty()) { + sql += "source_crs_auth_name = ? AND source_crs_code = ? AND "; + params.emplace_back(sourceCRSAuthName); + params.emplace_back(sourceCRSCode); + } + if (!targetCRSAuthName.empty()) { + sql += "target_crs_auth_name = ? AND target_crs_code = ? AND "; + params.emplace_back(targetCRSAuthName); + params.emplace_back(targetCRSCode); } - params = {sourceCRSAuthName, sourceCRSCode, targetCRSAuthName, - targetCRSCode}; + sql += "cov.deprecated = 0"; if (d->hasAuthorityRestriction()) { sql += " AND cov.auth_name = ?"; params.emplace_back(d->authority()); @@ -3379,7 +3388,7 @@ AuthorityFactory::createFromCoordinateReferenceSystemCodes( sql += " ORDER BY pseudo_area_from_swne(south_lat, west_lon, north_lat, " "east_lon) DESC, " "(CASE WHEN accuracy is NULL THEN 1 ELSE 0 END), accuracy"; - res = d->run(sql, params); + auto res = d->run(sql, params); std::set> setTransf; if (discardSuperseded) { for (const auto &row : res) { diff --git a/test/unit/test_operation.cpp b/test/unit/test_operation.cpp index a3b49ba9..e3eb4b7c 100644 --- a/test/unit/test_operation.cpp +++ b/test/unit/test_operation.cpp @@ -6859,23 +6859,28 @@ TEST(operation, compoundCRS_to_geogCRS_3D_context) { { auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0); + ctxt->setGridAvailabilityUse( + CoordinateOperationContext::GridAvailabilityUse:: + IGNORE_GRID_AVAILABILITY); auto list = CoordinateOperationFactory::create()->createOperations( authFactory->createCoordinateReferenceSystem( "5500"), // NAD83(NSRS2007) + NAVD88 height authFactory->createCoordinateReferenceSystem("4979"), // WGS 84 ctxt); ASSERT_GE(list.size(), 1U); - EXPECT_TRUE(list[0]->hasBallparkTransformation()); EXPECT_EQ(list[0]->nameStr(), - "NAD83(NSRS2007) to WGS 84 (1) + Transformation from NAVD88 " - "height to WGS 84 (ballpark vertical transformation, without " - "ellipsoid height to vertical height correction)"); + "NAD83(NSRS2007) to WGS 84 (1) + " + "Inverse of NAD83(2011) to NAVD88 height (1)"); EXPECT_EQ(list[0]->exportToPROJString( PROJStringFormatter::create( PROJStringFormatter::Convention::PROJ_5, authFactory->databaseContext()) .get()), - "+proj=noop"); + "+proj=pipeline +step +proj=axisswap +order=2,1 " + "+step +proj=unitconvert +xy_in=deg +xy_out=rad " + "+step +proj=vgridshift +grids=g2012bu0.gtx +multiplier=1 " + "+step +proj=unitconvert +xy_in=rad +xy_out=deg " + "+step +proj=axisswap +order=2,1"); } } -- cgit v1.2.3 From e6eae43cf2310c77a466fee257d9974b14ee85fd Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Thu, 12 Sep 2019 22:31:07 +0200 Subject: createOperations(): when tranforming from a compoundCRS whose vertical component is a BoundCRS, do not apply the horizontal transformation twice --- src/iso19111/coordinateoperation.cpp | 58 +++++++++++++++++++++++++----------- test/unit/test_operation.cpp | 38 ++++++++++++++++++++++- 2 files changed, 77 insertions(+), 19 deletions(-) diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp index aad86410..aea8400c 100644 --- a/src/iso19111/coordinateoperation.cpp +++ b/src/iso19111/coordinateoperation.cpp @@ -10173,6 +10173,7 @@ struct CoordinateOperationFactory::Private { bool inCreateOperationsWithDatumPivotAntiRecursion = false; bool inCreateOperationsThroughPreferredHub = false; bool inCreateOperationsGeogToVertWithIntermediate = false; + bool skipHorizontalTransformation = false; Context(const crs::CRSNNPtr &sourceCRSIn, const crs::CRSNNPtr &targetCRSIn, @@ -12784,27 +12785,32 @@ CoordinateOperationFactory::Private::createOperations( hubSrcGeog->coordinateSystem()->axisList().size() == 3 && geogDst->coordinateSystem()->axisList().size() == 3) { auto opsFirst = createOperations(sourceCRS, hubSrc, context); - auto opsSecond = createOperations(hubSrc, targetCRS, context); - if (!opsFirst.empty() && !opsSecond.empty()) { - for (const auto &opFirst : opsFirst) { - for (const auto &opLast : opsSecond) { - // Exclude artificial transformations from the hub - // to the target CRS - if (!opLast->hasBallparkTransformation()) { - try { - res.emplace_back( - ConcatenatedOperation:: - createComputeMetadata( - {opFirst, opLast}, - !allowEmptyIntersection)); - } catch ( - const InvalidOperationEmptyIntersection &) { + if (context.skipHorizontalTransformation) { + if (!opsFirst.empty()) + return opsFirst; + } else { + auto opsSecond = createOperations(hubSrc, targetCRS, context); + if (!opsFirst.empty() && !opsSecond.empty()) { + for (const auto &opFirst : opsFirst) { + for (const auto &opLast : opsSecond) { + // Exclude artificial transformations from the hub + // to the target CRS + if (!opLast->hasBallparkTransformation()) { + try { + res.emplace_back( + ConcatenatedOperation:: + createComputeMetadata( + {opFirst, opLast}, + !allowEmptyIntersection)); + } catch ( + const InvalidOperationEmptyIntersection &) { + } } } } - } - if (!res.empty()) { - return res; + if (!res.empty()) { + return res; + } } } } @@ -13026,6 +13032,22 @@ CoordinateOperationFactory::Private::createOperations( std::vector verticalTransforms; if (componentsSrc.size() >= 2 && componentsSrc[1]->extractVerticalCRS()) { + + struct SetSkipHorizontalTransform { + Context &context; + + explicit SetSkipHorizontalTransform(Context &contextIn) + : context(contextIn) { + assert(!context.skipHorizontalTransformation); + context.skipHorizontalTransformation = true; + } + + ~SetSkipHorizontalTransform() { + context.skipHorizontalTransformation = false; + } + }; + SetSkipHorizontalTransform setSkipHorizontalTransform(context); + verticalTransforms = createOperations(componentsSrc[1], targetCRS, context); bool foundRegisteredTransformWithAllGridsAvailable = false; diff --git a/test/unit/test_operation.cpp b/test/unit/test_operation.cpp index e3eb4b7c..735b8b64 100644 --- a/test/unit/test_operation.cpp +++ b/test/unit/test_operation.cpp @@ -6081,8 +6081,9 @@ TEST(operation, ETRS89_3D_to_proj_string_with_geoidgrids_nadgrids) { src, NN_NO_CHECK(dst), ctxt); ASSERT_EQ(list.size(), 1U); EXPECT_EQ(list[0]->exportToPROJString(PROJStringFormatter::create().get()), - "+proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad " + "+proj=pipeline " "+step +proj=axisswap +order=2,1 " + "+step +proj=unitconvert +xy_in=deg +xy_out=rad " "+step +inv +proj=vgridshift +grids=naptrans2008.gtx " "+multiplier=1 " "+step +inv +proj=hgridshift +grids=rdtrans2008.gsb " @@ -6093,6 +6094,41 @@ TEST(operation, ETRS89_3D_to_proj_string_with_geoidgrids_nadgrids) { // --------------------------------------------------------------------------- +TEST(operation, WGS84_G1762_to_compoundCRS_with_bound_vertCRS) { + auto authFactoryEPSG = + AuthorityFactory::create(DatabaseContext::create(), "EPSG"); + // WGS 84 (G1762) 3D + auto src = authFactoryEPSG->createCoordinateReferenceSystem("7665"); + auto objDst = PROJStringParser().createFromPROJString( + "+proj=longlat +datum=NAD83 +geoidgrids=@foo.gtx +type=crs"); + auto dst = nn_dynamic_pointer_cast(objDst); + ASSERT_TRUE(dst != nullptr); + auto authFactory = + AuthorityFactory::create(DatabaseContext::create(), std::string()); + auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0); + ctxt->setSpatialCriterion( + CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION); + ctxt->setGridAvailabilityUse( + CoordinateOperationContext::GridAvailabilityUse:: + IGNORE_GRID_AVAILABILITY); + auto list = CoordinateOperationFactory::create()->createOperations( + src, NN_NO_CHECK(dst), ctxt); + ASSERT_GE(list.size(), 53U); + EXPECT_EQ(list[0]->nameStr(), + "Inverse of unknown to WGS84 ellipsoidal height + " + "Inverse of WGS 84 to WGS 84 (G1762) + " + "Inverse of NAD83 to WGS 84 (1) + " + "Inverse of axis order change (2D)"); + 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 " + "+step +inv +proj=vgridshift +grids=@foo.gtx +multiplier=1 " + "+step +proj=unitconvert +xy_in=rad +xy_out=deg"); +} + +// --------------------------------------------------------------------------- + static VerticalCRSNNPtr createVerticalCRS() { PropertyMap propertiesVDatum; propertiesVDatum.set(Identifier::CODESPACE_KEY, "EPSG") -- cgit v1.2.3