diff options
| -rw-r--r-- | data/sql/customizations.sql | 17 | ||||
| -rw-r--r-- | data/sql/proj_db_table_defs.sql | 13 | ||||
| -rw-r--r-- | include/proj/io.hpp | 5 | ||||
| -rw-r--r-- | scripts/reference_exported_symbols.txt | 1 | ||||
| -rw-r--r-- | src/iso19111/c_api.cpp | 46 | ||||
| -rw-r--r-- | src/iso19111/coordinateoperation.cpp | 544 | ||||
| -rw-r--r-- | src/iso19111/factory.cpp | 92 | ||||
| -rw-r--r-- | src/proj_experimental.h | 4 | ||||
| -rw-r--r-- | test/unit/test_c_api.cpp | 31 | ||||
| -rw-r--r-- | test/unit/test_operation.cpp | 174 |
10 files changed, 797 insertions, 130 deletions
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<datum::GeodeticReferenceFrameNNPtr> + getPreferredHubGeodeticReferenceFrames( + const std::string &geodeticReferenceFrameCode) const; + //! @endcond protected: 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<VerticalCRS>(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<PositionalAccuracyNNPtr>()); + 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/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp index 768ef76f..aea8400c 100644 --- a/src/iso19111/coordinateoperation.cpp +++ b/src/iso19111/coordinateoperation.cpp @@ -56,7 +56,9 @@ #include <string> #include <vector> -#ifdef DEBUG +// #define DEBUG +// #define DEBUG_SORT +#if defined(DEBUG) || defined(DEBUG_SORT) #include <iostream> #endif @@ -10169,6 +10171,9 @@ struct CoordinateOperationFactory::Private { const crs::CRSNNPtr &targetCRS; const CoordinateOperationContextNNPtr &context; bool inCreateOperationsWithDatumPivotAntiRecursion = false; + bool inCreateOperationsThroughPreferredHub = false; + bool inCreateOperationsGeogToVertWithIntermediate = false; + bool skipHorizontalTransformation = false; Context(const crs::CRSNNPtr &sourceCRSIn, const crs::CRSNNPtr &targetCRSIn, @@ -10193,6 +10198,17 @@ struct CoordinateOperationFactory::Private { const crs::GeodeticCRS *geodSrc, const crs::GeodeticCRS *geodDst, Context &context); + static void createOperationsThroughPreferredHub( + std::vector<CoordinateOperationNNPtr> &res, + const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS, + const crs::GeodeticCRS *geodSrc, const crs::GeodeticCRS *geodDst, + Context &context); + + static std::vector<CoordinateOperationNNPtr> + createOperationsGeogToVertWithIntermediate(const crs::CRSNNPtr &sourceCRS, + const crs::CRSNNPtr &targetCRS, + Context &context); + static bool hasPerfectAccuracyResult(const std::vector<CoordinateOperationNNPtr> &res, const Context &context); @@ -10309,16 +10325,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 +10692,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 } // ---------------------------------------------------------------------- @@ -10917,18 +10957,21 @@ applyInverse(const std::vector<CoordinateOperationNNPtr> &list) { //! @cond Doxygen_Suppress -static void buildSourceAndTargetCRSIds( - const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS, - const CoordinateOperationContextNNPtr &context, - std::list<std::pair<std::string, std::string>> &sourceIds, - std::list<std::pair<std::string, std::string>> &targetIds) { +static void buildCRSIds(const crs::CRSNNPtr &crs, + const CoordinateOperationContextNNPtr &context, + std::list<std::pair<std::string, std::string>> &ids) { const auto &authFactory = context->getAuthorityFactory(); assert(authFactory); const auto &authFactoryName = authFactory->getAuthority(); - const auto findObjectInDB = [&authFactory, &authFactoryName]( - const crs::CRSNNPtr &crs, - std::list<std::pair<std::string, std::string>> &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(), @@ -10954,35 +10997,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<std::string> +getCandidateAuthorities(const io::AuthorityFactoryPtr &authFactory, + const std::string &srcAuthName, + const std::string &targetAuthName) { + const auto &authFactoryName = authFactory->getAuthority(); + std::vector<std::string> 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; } // --------------------------------------------------------------------------- @@ -10994,12 +11039,11 @@ findOpsInRegistryDirect(const crs::CRSNNPtr &sourceCRS, const CoordinateOperationContextNNPtr &context) { const auto &authFactory = context->getAuthorityFactory(); assert(authFactory); - const auto &authFactoryName = authFactory->getAuthority(); std::list<std::pair<std::string, std::string>> sourceIds; std::list<std::pair<std::string, std::string>> 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; @@ -11008,20 +11052,8 @@ findOpsInRegistryDirect(const crs::CRSNNPtr &sourceCRS, const auto &targetAuthName = idTarget.first; const auto &targetCode = idTarget.second; - std::vector<std::string> 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(), @@ -11042,6 +11074,81 @@ findOpsInRegistryDirect(const crs::CRSNNPtr &sourceCRS, } return std::vector<CoordinateOperationNNPtr>(); } + +// --------------------------------------------------------------------------- + +// Look in the authority registry for operations from sourceCRS +static std::vector<CoordinateOperationNNPtr> +findOpsInRegistryDirectFrom(const crs::CRSNNPtr &sourceCRS, + const CoordinateOperationContextNNPtr &context) { + const auto &authFactory = context->getAuthorityFactory(); + assert(authFactory); + + std::list<std::pair<std::string, std::string>> 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<CoordinateOperationNNPtr>(); +} + +// --------------------------------------------------------------------------- + +// Look in the authority registry for operations to targetCRS +static std::vector<CoordinateOperationNNPtr> +findOpsInRegistryDirectTo(const crs::CRSNNPtr &targetCRS, + const CoordinateOperationContextNNPtr &context) { + const auto &authFactory = context->getAuthorityFactory(); + assert(authFactory); + + std::list<std::pair<std::string, std::string>> 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<CoordinateOperationNNPtr>(); +} + //! @endcond // --------------------------------------------------------------------------- @@ -11056,12 +11163,11 @@ static std::vector<CoordinateOperationNNPtr> findsOpsInRegistryWithIntermediate( const auto &authFactory = context->getAuthorityFactory(); assert(authFactory); - const auto &authFactoryName = authFactory->getAuthority(); std::list<std::pair<std::string, std::string>> sourceIds; std::list<std::pair<std::string, std::string>> 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; @@ -11070,20 +11176,8 @@ static std::vector<CoordinateOperationNNPtr> findsOpsInRegistryWithIntermediate( const auto &targetAuthName = idTarget.first; const auto &targetCode = idTarget.second; - std::vector<std::string> 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(), @@ -11990,6 +12084,167 @@ void CoordinateOperationFactory::Private::createOperationsWithDatumPivot( // --------------------------------------------------------------------------- +void CoordinateOperationFactory::Private::createOperationsThroughPreferredHub( + std::vector<CoordinateOperationNNPtr> &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<CoordinateOperationNNPtr> 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<crs::CRSNNPtr> candidatesIntermCRS; + for (const auto &datumHub : dstPreferredHubs) { + auto candidates = + findCandidateGeodCRSForDatum(authFactory, datumHub.get()); + bool addedGeog2D = false; + for (const auto &intermCRS : candidates) { + auto geogCRS = dynamic_cast<crs::GeographicCRS *>(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 &) { + } + } + } + } + } +} + +// --------------------------------------------------------------------------- + +std::vector<CoordinateOperationNNPtr> +CoordinateOperationFactory::Private::createOperationsGeogToVertWithIntermediate( + const crs::CRSNNPtr & /*sourceCRS*/, // geographic CRS + const crs::CRSNNPtr &targetCRS, // vertical CRS + Private::Context &context) { + + std::vector<CoordinateOperationNNPtr> 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<const crs::GeographicCRS *>(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) { @@ -12201,6 +12456,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) { @@ -12236,7 +12506,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(); + } } } @@ -12505,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; + } } } } @@ -12747,8 +13032,79 @@ CoordinateOperationFactory::Private::createOperations( std::vector<CoordinateOperationNNPtr> 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; + 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 1f40f1f0..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<std::pair<std::string, std::string>> setTransf; if (discardSuperseded) { for (const auto &row : res) { @@ -4982,6 +4991,29 @@ AuthorityFactory::createCompoundCRSFromExisting( // --------------------------------------------------------------------------- //! @cond Doxygen_Suppress +std::list<datum::GeodeticReferenceFrameNNPtr> +AuthorityFactory::getPreferredHubGeodeticReferenceFrames( + const std::string &geodeticReferenceFrameCode) const { + std::list<datum::GeodeticReferenceFrameNNPtr> 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/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 diff --git a/test/unit/test_operation.cpp b/test/unit/test_operation.cpp index 21275284..735b8b64 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,100 @@ 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, 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"); @@ -4588,6 +4691,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 +4721,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 +5362,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( @@ -5969,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 " @@ -5981,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<CRS>(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") @@ -6535,6 +6683,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 +6864,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) @@ -6741,23 +6895,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"); } } @@ -6924,6 +7083,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( |