From c08da1a6604f9cb422bb71e34b6006ede6146642 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Sat, 31 Oct 2020 10:52:44 +0100 Subject: proj_create_ellipsoidal_2D/3D_cs(): doc improvement --- src/iso19111/c_api.cpp | 21 +++++++++++++-------- 1 file changed, 13 insertions(+), 8 deletions(-) (limited to 'src/iso19111') diff --git a/src/iso19111/c_api.cpp b/src/iso19111/c_api.cpp index 7582f37f..f4f6e6a2 100644 --- a/src/iso19111/c_api.cpp +++ b/src/iso19111/c_api.cpp @@ -4438,8 +4438,9 @@ PJ *proj_create_cartesian_2D_cs(PJ_CONTEXT *ctx, PJ_CARTESIAN_CS_2D_TYPE type, * * @param ctx PROJ context, or NULL for default context * @param type Coordinate system type. - * @param unit_name Unit name. - * @param unit_conv_factor Unit conversion factor to SI. + * @param unit_name Name of the angular units. Or NULL for Degree + * @param unit_conv_factor Conversion factor from the angular unit to radian. + * Or 0 for Degree if unit_name == NULL. Otherwise should be not NULL * * @return Object that must be unreferenced with * proj_destroy(), or NULL in case of error. @@ -4478,13 +4479,17 @@ PJ *proj_create_ellipsoidal_2D_cs(PJ_CONTEXT *ctx, * * @param ctx PROJ context, or NULL for default context * @param type Coordinate system type. - * @param horizontal_angular_unit_name Horizontal angular unit name. - * @param horizontal_angular_unit_conv_factor Horizontal angular unit conversion - * factor to SI. - * @param vertical_linear_unit_name Vertical linear unit name. + * @param horizontal_angular_unit_name Name of the angular units. Or NULL for + * Degree. + * @param horizontal_angular_unit_conv_factor Conversion factor from the angular + * unit to radian. Or 0 for Degree if horizontal_angular_unit_name == NULL. + * Otherwise should be not NULL + * @param vertical_linear_unit_name Vertical linear unit name. Or NULL for + * Metre. * @param vertical_linear_unit_conv_factor Vertical linear unit conversion - * factor to SI. - * + * factor to metre. Or 0 for Metre if vertical_linear_unit_name == NULL. + * Otherwise should be not NULL + * @return Object that must be unreferenced with * proj_destroy(), or NULL in case of error. * @since 6.3 -- cgit v1.2.3 From baed8dae4980698873500dec16b158638df71318 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Sat, 31 Oct 2020 13:59:25 +0100 Subject: Make sure that remarks is preserved when substituting (EPSG) grid name by PROJ one --- src/iso19111/coordinateoperation.cpp | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) (limited to 'src/iso19111') diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp index be15b3e0..599e04bd 100644 --- a/src/iso19111/coordinateoperation.cpp +++ b/src/iso19111/coordinateoperation.cpp @@ -8836,11 +8836,16 @@ createSimilarPropertiesTransformation(TransformationNNPtr obj) { // The domain(s) are unchanged addDomains(map, obj.get()); - std::string forwardName = obj->nameStr(); + const std::string &forwardName = obj->nameStr(); if (!forwardName.empty()) { map.set(common::IdentifiedObject::NAME_KEY, forwardName); } + const std::string &remarks = obj->remarks(); + if (!remarks.empty()) { + map.set(common::IdentifiedObject::REMARKS_KEY, remarks); + } + addModifiedIdentifier(map, obj.get(), false, true); return map; -- cgit v1.2.3 From 67718853e1248656cc35caa5a605ac102ec7e985 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Sat, 31 Oct 2020 15:47:52 +0100 Subject: Fill remarks for PROJ-based operation mixing horizontal and vertical transformations --- src/iso19111/coordinateoperation.cpp | 63 ++++++++++++++++++++++++++++++++++++ 1 file changed, 63 insertions(+) (limited to 'src/iso19111') diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp index 599e04bd..e9d73f52 100644 --- a/src/iso19111/coordinateoperation.cpp +++ b/src/iso19111/coordinateoperation.cpp @@ -12620,6 +12620,55 @@ createGeodToGeodPROJBased(const crs::CRSNNPtr &geodSrc, // --------------------------------------------------------------------------- +static std::string +getRemarks(const std::vector &ops) { + std::string remarks; + for (const auto &op : ops) { + const auto &opRemarks = op->remarks(); + if (!opRemarks.empty()) { + if (!remarks.empty()) { + remarks += '\n'; + } + + std::string opName(op->nameStr()); + if (starts_with(opName, INVERSE_OF)) { + opName = opName.substr(INVERSE_OF.size()); + } + + remarks += "For "; + remarks += opName; + + const auto &ids = op->identifiers(); + if (!ids.empty()) { + std::string authority(*ids.front()->codeSpace()); + if (starts_with(authority, "INVERSE(") && + authority.back() == ')') { + authority = authority.substr(strlen("INVERSE("), + authority.size() - 1 - + strlen("INVERSE(")); + } + if (starts_with(authority, "DERIVED_FROM(") && + authority.back() == ')') { + authority = authority.substr(strlen("DERIVED_FROM("), + authority.size() - 1 - + strlen("DERIVED_FROM(")); + } + + remarks += " ("; + remarks += authority; + remarks += ':'; + remarks += ids.front()->code(); + remarks += ')'; + } + remarks += ": "; + remarks += opRemarks; + } + } + return remarks; +} + +// --------------------------------------------------------------------------- + static CoordinateOperationNNPtr createHorizVerticalPROJBased( const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS, const operation::CoordinateOperationNNPtr &horizTransform, @@ -12645,6 +12694,10 @@ static CoordinateOperationNNPtr createHorizVerticalPROJBased( properties.set(common::ObjectUsage::DOMAIN_OF_VALIDITY_KEY, NN_NO_CHECK(extent)); } + const auto &remarks = verticalTransform->remarks(); + if (!remarks.empty()) { + properties.set(common::IdentifiedObject::REMARKS_KEY, remarks); + } return createPROJBased( properties, exportable, sourceCRS, targetCRS, nullptr, verticalTransform->coordinateOperationAccuracies(), @@ -12669,6 +12722,11 @@ static CoordinateOperationNNPtr createHorizVerticalPROJBased( NN_NO_CHECK(extent)); } + const auto remarks = getRemarks(ops); + if (!remarks.empty()) { + properties.set(common::IdentifiedObject::REMARKS_KEY, remarks); + } + std::vector accuracies; const double accuracy = getAccuracy(ops); if (accuracy >= 0.0) { @@ -12729,6 +12787,11 @@ static CoordinateOperationNNPtr createHorizVerticalHorizPROJBased( NN_NO_CHECK(extent)); } + const auto remarks = getRemarks(ops); + if (!remarks.empty()) { + properties.set(common::IdentifiedObject::REMARKS_KEY, remarks); + } + std::vector accuracies; const double accuracy = getAccuracy(ops); if (accuracy >= 0.0) { -- cgit v1.2.3 From 6e7124c8b16dcc00e4a3f4336d1535f9655f1aa3 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Sat, 31 Oct 2020 15:53:15 +0100 Subject: Doxygen: improve doc of functions related to Datum/DatumEnsemble --- src/iso19111/c_api.cpp | 8 ++++++++ 1 file changed, 8 insertions(+) (limited to 'src/iso19111') diff --git a/src/iso19111/c_api.cpp b/src/iso19111/c_api.cpp index 7582f37f..f463dc2a 100644 --- a/src/iso19111/c_api.cpp +++ b/src/iso19111/c_api.cpp @@ -2049,6 +2049,8 @@ PJ *proj_get_ellipsoid(PJ_CONTEXT *ctx, const PJ *obj) { // --------------------------------------------------------------------------- /** \brief Get the horizontal datum from a CRS + * + * This function may return a Datum or DatumEnsemble object. * * The returned object must be unreferenced with proj_destroy() after * use. @@ -7986,6 +7988,9 @@ double proj_coordoperation_get_accuracy(PJ_CONTEXT *ctx, // --------------------------------------------------------------------------- /** \brief Returns the datum of a SingleCRS. + * + * If that function returns NULL, @see proj_crs_get_datum_ensemble() to + * potentially get a DatumEnsemble instead. * * The returned object must be unreferenced with proj_destroy() after * use. @@ -8017,6 +8022,9 @@ PJ *proj_crs_get_datum(PJ_CONTEXT *ctx, const PJ *crs) { // --------------------------------------------------------------------------- /** \brief Returns the datum ensemble of a SingleCRS. + * + * If that function returns NULL, @see proj_crs_get_datum() to + * potentially get a Datum instead. * * The returned object must be unreferenced with proj_destroy() after * use. -- cgit v1.2.3 From 82695cce869e1bcf2eb2a8ff078b679b6a21c663 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Thu, 8 Oct 2020 17:08:27 +0200 Subject: When reading from database, possibly return Geographic/GeodeticCRS with a DatumEnsemble, typically for WGS 84 and ETRS89 ('breaking change') --- src/iso19111/factory.cpp | 158 ++++++++++++++++++++++++++++++++++------------- src/iso19111/io.cpp | 4 ++ 2 files changed, 118 insertions(+), 44 deletions(-) (limited to 'src/iso19111') diff --git a/src/iso19111/factory.cpp b/src/iso19111/factory.cpp index 5d02aeea..cce80b8a 100644 --- a/src/iso19111/factory.cpp +++ b/src/iso19111/factory.cpp @@ -192,6 +192,13 @@ struct DatabaseContext::Private { void cache(const std::string &code, const datum::GeodeticReferenceFrameNNPtr &datum); + datum::DatumEnsemblePtr + // cppcheck-suppress functionStatic + getDatumEnsembleFromCache(const std::string &code); + // cppcheck-suppress functionStatic + void cache(const std::string &code, + const datum::DatumEnsembleNNPtr &datumEnsemble); + datum::EllipsoidPtr // cppcheck-suppress functionStatic getEllipsoidFromCache(const std::string &code); @@ -258,6 +265,7 @@ struct DatabaseContext::Private { LRUCacheOfObjects cacheCRS_{CACHE_SIZE}; LRUCacheOfObjects cacheEllipsoid_{CACHE_SIZE}; LRUCacheOfObjects cacheGeodeticDatum_{CACHE_SIZE}; + LRUCacheOfObjects cacheDatumEnsemble_{CACHE_SIZE}; LRUCacheOfObjects cachePrimeMeridian_{CACHE_SIZE}; LRUCacheOfObjects cacheCS_{CACHE_SIZE}; LRUCacheOfObjects cacheExtent_{CACHE_SIZE}; @@ -417,6 +425,22 @@ void DatabaseContext::Private::cache( // --------------------------------------------------------------------------- +datum::DatumEnsemblePtr +DatabaseContext::Private::getDatumEnsembleFromCache(const std::string &code) { + util::BaseObjectPtr obj; + getFromCache(cacheDatumEnsemble_, code, obj); + return std::static_pointer_cast(obj); +} + +// --------------------------------------------------------------------------- + +void DatabaseContext::Private::cache( + const std::string &code, const datum::DatumEnsembleNNPtr &datumEnsemble) { + insertIntoCache(cacheDatumEnsemble_, code, datumEnsemble.as_nullable()); +} + +// --------------------------------------------------------------------------- + datum::EllipsoidPtr DatabaseContext::Private::getEllipsoidFromCache(const std::string &code) { util::BaseObjectPtr obj; @@ -1064,18 +1088,6 @@ std::string DatabaseContext::getOldProjGridName(const std::string &gridName) { // --------------------------------------------------------------------------- -// FIXME: as we don't support datum ensemble yet, add it from name -static std::string removeEnsembleSuffix(const std::string &name) { - if (name == "World Geodetic System 1984 ensemble") { - return "World Geodetic System 1984"; - } else if (name == "European Terrestrial Reference System 1989 ensemble") { - return "European Terrestrial Reference System 1989"; - } - return name; -} - -// --------------------------------------------------------------------------- - /** \brief Gets the alias name from an official name. * * @param officialName Official name. Mandatory @@ -2012,18 +2024,38 @@ AuthorityFactory::createEllipsoid(const std::string &code) const { datum::GeodeticReferenceFrameNNPtr AuthorityFactory::createGeodeticDatum(const std::string &code) const { + + datum::GeodeticReferenceFramePtr datum; + datum::DatumEnsemblePtr datumEnsemble; + constexpr bool turnEnsembleAsDatum = true; + createGeodeticDatumOrEnsemble(code, datum, datumEnsemble, + turnEnsembleAsDatum); + return NN_NO_CHECK(datum); +} + +// --------------------------------------------------------------------------- + +void AuthorityFactory::createGeodeticDatumOrEnsemble( + const std::string &code, datum::GeodeticReferenceFramePtr &outDatum, + datum::DatumEnsemblePtr &outDatumEnsemble, bool turnEnsembleAsDatum) const { const auto cacheKey(d->authority() + code); { - auto datum = d->context()->d->getGeodeticDatumFromCache(cacheKey); - if (datum) { - return NN_NO_CHECK(datum); + outDatumEnsemble = d->context()->d->getDatumEnsembleFromCache(cacheKey); + if (outDatumEnsemble) { + if (!turnEnsembleAsDatum) + return; + outDatumEnsemble = nullptr; + } + outDatum = d->context()->d->getGeodeticDatumFromCache(cacheKey); + if (outDatum) { + return; } } auto res = d->runWithCodeParam("SELECT name, ellipsoid_auth_name, ellipsoid_code, " "prime_meridian_auth_name, prime_meridian_code, " "publication_date, frame_reference_epoch, " - "deprecated FROM geodetic_datum " + "ensemble_accuracy, deprecated FROM geodetic_datum " "WHERE " "auth_name = ? AND code = ?", code); @@ -2040,29 +2072,63 @@ AuthorityFactory::createGeodeticDatum(const std::string &code) const { const auto &prime_meridian_code = row[4]; const auto &publication_date = row[5]; const auto &frame_reference_epoch = row[6]; - const bool deprecated = row[7] == "1"; - auto ellipsoid = d->createFactory(ellipsoid_auth_name) - ->createEllipsoid(ellipsoid_code); - auto pm = d->createFactory(prime_meridian_auth_name) - ->createPrimeMeridian(prime_meridian_code); - auto props = d->createPropertiesSearchUsages( - "geodetic_datum", code, removeEnsembleSuffix(name), deprecated); - auto anchor = util::optional(); - if (!publication_date.empty()) { - props.set("PUBLICATION_DATE", publication_date); + const auto &ensemble_accuracy = row[7]; + const bool deprecated = row[8] == "1"; + + std::string massagedName = name; + if (turnEnsembleAsDatum) { + if (name == "World Geodetic System 1984 ensemble") { + massagedName = "World Geodetic System 1984"; + } else if (name == + "European Terrestrial Reference System 1989 ensemble") { + massagedName = "European Terrestrial Reference System 1989"; + } + } + auto props = d->createPropertiesSearchUsages("geodetic_datum", code, + massagedName, deprecated); + + if (!turnEnsembleAsDatum && !ensemble_accuracy.empty()) { + auto resMembers = + d->run("SELECT member_auth_name, member_code FROM " + "geodetic_datum_ensemble_member WHERE " + "ensemble_auth_name = ? AND ensemble_code = ? " + "ORDER BY sequence", + {d->authority(), code}); + + std::vector members; + for (const auto &memberRow : resMembers) { + members.push_back( + d->createFactory(memberRow[0])->createDatum(memberRow[1])); + } + auto datumEnsemble = datum::DatumEnsemble::create( + props, std::move(members), + metadata::PositionalAccuracy::create(ensemble_accuracy)); + d->context()->d->cache(cacheKey, datumEnsemble); + outDatumEnsemble = datumEnsemble.as_nullable(); + } else { + auto ellipsoid = d->createFactory(ellipsoid_auth_name) + ->createEllipsoid(ellipsoid_code); + auto pm = d->createFactory(prime_meridian_auth_name) + ->createPrimeMeridian(prime_meridian_code); + + auto anchor = util::optional(); + if (!publication_date.empty()) { + props.set("PUBLICATION_DATE", publication_date); + } + auto datum = frame_reference_epoch.empty() + ? datum::GeodeticReferenceFrame::create( + props, ellipsoid, anchor, pm) + : util::nn_static_pointer_cast< + datum::GeodeticReferenceFrame>( + datum::DynamicGeodeticReferenceFrame::create( + props, ellipsoid, anchor, pm, + common::Measure( + c_locale_stod(frame_reference_epoch), + common::UnitOfMeasure::YEAR), + util::optional())); + d->context()->d->cache(cacheKey, datum); + outDatum = datum.as_nullable(); } - auto datum = - frame_reference_epoch.empty() - ? datum::GeodeticReferenceFrame::create(props, ellipsoid, - anchor, pm) - : util::nn_static_pointer_cast( - datum::DynamicGeodeticReferenceFrame::create( - props, ellipsoid, anchor, pm, - common::Measure(c_locale_stod(frame_reference_epoch), - common::UnitOfMeasure::YEAR), - util::optional())); - d->context()->d->cache(cacheKey, datum); - return datum; } catch (const std::exception &ex) { throw buildFactoryException("geodetic reference frame", code, ex); } @@ -2472,20 +2538,24 @@ AuthorityFactory::createGeodeticCRS(const std::string &code, auto cs = d->createFactory(cs_auth_name)->createCoordinateSystem(cs_code); - auto datum = - d->createFactory(datum_auth_name)->createGeodeticDatum(datum_code); + datum::GeodeticReferenceFramePtr datum; + datum::DatumEnsemblePtr datumEnsemble; + constexpr bool turnEnsembleAsDatum = false; + d->createFactory(datum_auth_name) + ->createGeodeticDatumOrEnsemble(datum_code, datum, datumEnsemble, + turnEnsembleAsDatum); auto ellipsoidalCS = util::nn_dynamic_pointer_cast(cs); if ((type == GEOG_2D || type == GEOG_3D) && ellipsoidalCS) { auto crsRet = crs::GeographicCRS::create( - props, datum, NN_NO_CHECK(ellipsoidalCS)); + props, datum, datumEnsemble, NN_NO_CHECK(ellipsoidalCS)); d->context()->d->cache(cacheKey, crsRet); return crsRet; } auto geocentricCS = util::nn_dynamic_pointer_cast(cs); if (type == GEOCENTRIC && geocentricCS) { - auto crsRet = crs::GeodeticCRS::create(props, datum, + auto crsRet = crs::GeodeticCRS::create(props, datum, datumEnsemble, NN_NO_CHECK(geocentricCS)); d->context()->d->cache(cacheKey, crsRet); return crsRet; @@ -5617,7 +5687,7 @@ std::string AuthorityFactory::getOfficialNameFromAlias( if (res.empty()) { // shouldn't happen normally return std::string(); } - return removeEnsembleSuffix(res.front()[0]); + return res.front()[0]; } } return std::string(); @@ -5667,7 +5737,7 @@ std::string AuthorityFactory::getOfficialNameFromAlias( outTableName = row[1]; outAuthName = row[2]; outCode = row[3]; - return removeEnsembleSuffix(row[0]); + return row[0]; } } diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp index f4ec7740..fde697ad 100644 --- a/src/iso19111/io.cpp +++ b/src/iso19111/io.cpp @@ -2110,6 +2110,10 @@ GeodeticReferenceFrameNNPtr WKTParser::Private::buildGeodeticReferenceFrame( name = "World Geodetic System 1984"; authNameFromAlias = Identifier::EPSG; codeFromAlias = "6326"; + } else if (name == "D_ETRS_1989") { + name = "European Terrestrial Reference System 1989"; + authNameFromAlias = Identifier::EPSG; + codeFromAlias = "6258"; } else { tableNameForAlias = "geodetic_datum"; } -- cgit v1.2.3 From c2b0dcc468b4e722e46fe10fca93fe70a95fcb8e Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Thu, 8 Oct 2020 18:41:57 +0200 Subject: When reading from database, possibly return VerticalCRS with a DatumEnsemble Only occurence for now is EPSG:9451 'BI height' using the 'British Isles height ensemble' --- src/iso19111/factory.cpp | 84 ++++++++++++++++++++++++++++++++++++------------ 1 file changed, 63 insertions(+), 21 deletions(-) (limited to 'src/iso19111') diff --git a/src/iso19111/factory.cpp b/src/iso19111/factory.cpp index cce80b8a..301b88be 100644 --- a/src/iso19111/factory.cpp +++ b/src/iso19111/factory.cpp @@ -2146,9 +2146,23 @@ void AuthorityFactory::createGeodeticDatumOrEnsemble( datum::VerticalReferenceFrameNNPtr AuthorityFactory::createVerticalDatum(const std::string &code) const { + datum::VerticalReferenceFramePtr datum; + datum::DatumEnsemblePtr datumEnsemble; + constexpr bool turnEnsembleAsDatum = true; + createVerticalDatumOrEnsemble(code, datum, datumEnsemble, + turnEnsembleAsDatum); + return NN_NO_CHECK(datum); +} + +// --------------------------------------------------------------------------- + +void AuthorityFactory::createVerticalDatumOrEnsemble( + const std::string &code, datum::VerticalReferenceFramePtr &outDatum, + datum::DatumEnsemblePtr &outDatumEnsemble, bool turnEnsembleAsDatum) const { auto res = d->runWithCodeParam("SELECT name, publication_date, " - "frame_reference_epoch, deprecated FROM " + "frame_reference_epoch, ensemble_accuracy, " + "deprecated FROM " "vertical_datum WHERE auth_name = ? AND code = ?", code); if (res.empty()) { @@ -2160,24 +2174,49 @@ AuthorityFactory::createVerticalDatum(const std::string &code) const { const auto &name = row[0]; const auto &publication_date = row[1]; const auto &frame_reference_epoch = row[2]; - const bool deprecated = row[3] == "1"; + const auto &ensemble_accuracy = row[3]; + const bool deprecated = row[4] == "1"; auto props = d->createPropertiesSearchUsages("vertical_datum", code, name, deprecated); - if (!publication_date.empty()) { - props.set("PUBLICATION_DATE", publication_date); - } - if (d->authority() == "ESRI" && starts_with(code, "from_geogdatum_")) { - props.set("VERT_DATUM_TYPE", "2002"); - } - auto anchor = util::optional(); - if (frame_reference_epoch.empty()) { - return datum::VerticalReferenceFrame::create(props, anchor); + if (!turnEnsembleAsDatum && !ensemble_accuracy.empty()) { + auto resMembers = + d->run("SELECT member_auth_name, member_code FROM " + "vertical_datum_ensemble_member WHERE " + "ensemble_auth_name = ? AND ensemble_code = ? " + "ORDER BY sequence", + {d->authority(), code}); + + std::vector members; + for (const auto &memberRow : resMembers) { + members.push_back( + d->createFactory(memberRow[0])->createDatum(memberRow[1])); + } + auto datumEnsemble = datum::DatumEnsemble::create( + props, std::move(members), + metadata::PositionalAccuracy::create(ensemble_accuracy)); + outDatumEnsemble = datumEnsemble.as_nullable(); } else { - return datum::DynamicVerticalReferenceFrame::create( - props, anchor, util::optional(), - common::Measure(c_locale_stod(frame_reference_epoch), - common::UnitOfMeasure::YEAR), - util::optional()); + if (!publication_date.empty()) { + props.set("PUBLICATION_DATE", publication_date); + } + if (d->authority() == "ESRI" && + starts_with(code, "from_geogdatum_")) { + props.set("VERT_DATUM_TYPE", "2002"); + } + auto anchor = util::optional(); + if (frame_reference_epoch.empty()) { + outDatum = datum::VerticalReferenceFrame::create(props, anchor) + .as_nullable(); + } else { + outDatum = + datum::DynamicVerticalReferenceFrame::create( + props, anchor, + util::optional(), + common::Measure(c_locale_stod(frame_reference_epoch), + common::UnitOfMeasure::YEAR), + util::optional()) + .as_nullable(); + } } } catch (const std::exception &ex) { throw buildFactoryException("vertical reference frame", code, ex); @@ -2609,16 +2648,19 @@ AuthorityFactory::createVerticalCRS(const std::string &code) const { const bool deprecated = row[5] == "1"; auto cs = d->createFactory(cs_auth_name)->createCoordinateSystem(cs_code); - auto datum = - d->createFactory(datum_auth_name)->createVerticalDatum(datum_code); - + datum::VerticalReferenceFramePtr datum; + datum::DatumEnsemblePtr datumEnsemble; + constexpr bool turnEnsembleAsDatum = false; + d->createFactory(datum_auth_name) + ->createVerticalDatumOrEnsemble(datum_code, datum, datumEnsemble, + turnEnsembleAsDatum); auto props = d->createPropertiesSearchUsages("vertical_crs", code, name, deprecated); auto verticalCS = util::nn_dynamic_pointer_cast(cs); if (verticalCS) { - auto crsRet = - crs::VerticalCRS::create(props, datum, NN_NO_CHECK(verticalCS)); + auto crsRet = crs::VerticalCRS::create(props, datum, datumEnsemble, + NN_NO_CHECK(verticalCS)); d->context()->d->cache(cacheKey, crsRet); return crsRet; } -- cgit v1.2.3 From 1e5acb00a0c0fc2533b9bce2e5803da10ed1d8d6 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Thu, 8 Oct 2020 20:59:19 +0200 Subject: projinfo / createObjectsFromName(): support returning a datum ensemble --- src/iso19111/c_api.cpp | 2 +- src/iso19111/datum.cpp | 8 ++++++++ src/iso19111/factory.cpp | 51 +++++++++++++++++++++++++++++++++++++++++++++++- src/iso19111/io.cpp | 5 +++++ 4 files changed, 64 insertions(+), 2 deletions(-) (limited to 'src/iso19111') diff --git a/src/iso19111/c_api.cpp b/src/iso19111/c_api.cpp index 90a414c6..4fb73d31 100644 --- a/src/iso19111/c_api.cpp +++ b/src/iso19111/c_api.cpp @@ -916,7 +916,7 @@ convertPJObjectTypeToObjectType(PJ_TYPE type, bool &valid) { break; case PJ_TYPE_DATUM_ENSEMBLE: - cppType = AuthorityFactory::ObjectType::DATUM; + cppType = AuthorityFactory::ObjectType::DATUM_ENSEMBLE; break; case PJ_TYPE_TEMPORAL_DATUM: diff --git a/src/iso19111/datum.cpp b/src/iso19111/datum.cpp index 5bc8074c..e29f6319 100644 --- a/src/iso19111/datum.cpp +++ b/src/iso19111/datum.cpp @@ -1731,6 +1731,14 @@ void DatumEnsemble::_exportToWKT( formatter->startNode(io::WKTConstants::ENSEMBLEACCURACY, false); formatter->add(positionalAccuracy()->value()); formatter->endNode(); + + // In theory, we should do the following, but currently the WKT grammar + // doesn't allow this + // ObjectUsage::baseExportToWKT(formatter); + if (formatter->outputId()) { + formatID(formatter); + } + formatter->endNode(); } //! @endcond diff --git a/src/iso19111/factory.cpp b/src/iso19111/factory.cpp index 301b88be..0ccfefc1 100644 --- a/src/iso19111/factory.cpp +++ b/src/iso19111/factory.cpp @@ -5452,6 +5452,11 @@ AuthorityFactory::getAuthorityCodes(const ObjectType &type, case ObjectType::CONCATENATED_OPERATION: sql = "SELECT code FROM concatenated_operation WHERE "; break; + case ObjectType::DATUM_ENSEMBLE: + sql = "SELECT code FROM object_view WHERE table_name IN " + "('geodetic_datum', 'vertical_datum') AND " + "type = 'ensemble' AND "; + break; } sql += "auth_name = ?"; @@ -5961,12 +5966,28 @@ AuthorityFactory::createObjectsFromNameEx( res.emplace_back( TableType("concatenated_operation", std::string())); break; + case ObjectType::DATUM_ENSEMBLE: + res.emplace_back(TableType("geodetic_datum", "ensemble")); + res.emplace_back(TableType("vertical_datum", "ensemble")); + break; } } } return res; }; + bool datumEnsembleAllowed = false; + if (allowedObjectTypes.empty()) { + datumEnsembleAllowed = true; + } else { + for (const auto type : allowedObjectTypes) { + if (type == ObjectType::DATUM_ENSEMBLE) { + datumEnsembleAllowed = true; + break; + } + } + } + const auto listTableNameType = getTableAndTypeConstraints(); bool first = true; ListOfParams params; @@ -5984,6 +6005,8 @@ AuthorityFactory::createObjectsFromNameEx( if (!tableNameTypePair.second.empty()) { if (tableNameTypePair.second == "frame_reference_epoch") { sql += "AND frame_reference_epoch IS NOT NULL "; + } else if (tableNameTypePair.second == "ensemble") { + sql += "AND ensemble_accuracy IS NOT NULL "; } else { sql += "AND type = '"; sql += tableNameTypePair.second; @@ -6018,6 +6041,8 @@ AuthorityFactory::createObjectsFromNameEx( if (!tableNameTypePair.second.empty()) { if (tableNameTypePair.second == "frame_reference_epoch") { sql += "AND ov.frame_reference_epoch IS NOT NULL "; + } else if (tableNameTypePair.second == "ensemble") { + sql += "AND ov.ensemble_accuracy IS NOT NULL "; } else { sql += "AND ov.type = '"; sql += tableNameTypePair.second; @@ -6165,7 +6190,7 @@ AuthorityFactory::createObjectsFromNameEx( break; } auto factory = d->createFactory(auth_name); - auto getObject = [&factory]( + auto getObject = [&factory, datumEnsembleAllowed]( const std::string &l_table_name, const std::string &l_code) -> common::IdentifiedObjectNNPtr { if (l_table_name == "prime_meridian") { @@ -6173,8 +6198,32 @@ AuthorityFactory::createObjectsFromNameEx( } else if (l_table_name == "ellipsoid") { return factory->createEllipsoid(l_code); } else if (l_table_name == "geodetic_datum") { + if (datumEnsembleAllowed) { + datum::GeodeticReferenceFramePtr datum; + datum::DatumEnsemblePtr datumEnsemble; + constexpr bool turnEnsembleAsDatum = false; + factory->createGeodeticDatumOrEnsemble( + l_code, datum, datumEnsemble, turnEnsembleAsDatum); + if (datum) { + return NN_NO_CHECK(datum); + } + assert(datumEnsemble); + return NN_NO_CHECK(datumEnsemble); + } return factory->createGeodeticDatum(l_code); } else if (l_table_name == "vertical_datum") { + if (datumEnsembleAllowed) { + datum::VerticalReferenceFramePtr datum; + datum::DatumEnsemblePtr datumEnsemble; + constexpr bool turnEnsembleAsDatum = false; + factory->createVerticalDatumOrEnsemble( + l_code, datum, datumEnsemble, turnEnsembleAsDatum); + if (datum) { + return NN_NO_CHECK(datum); + } + assert(datumEnsemble); + return NN_NO_CHECK(datumEnsemble); + } return factory->createVerticalDatum(l_code); } else if (l_table_name == "geodetic_crs") { return factory->createGeodeticCRS(l_code); diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp index fde697ad..b2e5822b 100644 --- a/src/iso19111/io.cpp +++ b/src/iso19111/io.cpp @@ -6267,6 +6267,9 @@ static BaseObjectNNPtr createFromUserInput(const std::string &text, if (type == "datum") { return factory->createDatum(code); } + if (type == "ensemble") { + return factory->createDatumEnsemble(code); + } if (type == "ellipsoid") { return factory->createEllipsoid(code); } @@ -6384,6 +6387,8 @@ static BaseObjectNNPtr createFromUserInput(const std::string &text, ELLIPSOID, AuthorityFactory::ObjectType:: DATUM, + AuthorityFactory::ObjectType:: + DATUM_ENSEMBLE, AuthorityFactory::ObjectType:: COORDINATE_OPERATION}, goOn); -- cgit v1.2.3 From 119888b041258267768d632b89395e7074323326 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Sat, 10 Oct 2020 15:35:29 +0200 Subject: proj.h: add PJ_CATEGORY_DATUM_ENSEMBLE for proj_create_from_database() --- src/iso19111/c_api.cpp | 3 +++ 1 file changed, 3 insertions(+) (limited to 'src/iso19111') diff --git a/src/iso19111/c_api.cpp b/src/iso19111/c_api.cpp index 4fb73d31..cbbdcaa8 100644 --- a/src/iso19111/c_api.cpp +++ b/src/iso19111/c_api.cpp @@ -675,6 +675,9 @@ PJ *proj_create_from_database(PJ_CONTEXT *ctx, const char *auth_name, codeStr, usePROJAlternativeGridNames != 0) .as_nullable(); break; + case PJ_CATEGORY_DATUM_ENSEMBLE: + obj = factory->createDatumEnsemble(codeStr).as_nullable(); + break; } return pj_obj_create(ctx, NN_NO_CHECK(obj)); } catch (const std::exception &e) { -- cgit v1.2.3 From 848ee7b6a8bf56a9d967f104ba2ef2ad92dc48d5 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Wed, 4 Nov 2020 15:24:44 +0100 Subject: createBoundCRSToWGS84IfPossible(): make it return same result with a CRS built from EPSG code or WKT1 Related to https://github.com/OSGeo/gdal/issues/3144 --- src/iso19111/coordinateoperation.cpp | 296 ++++++++++++++++++----------------- src/iso19111/crs.cpp | 6 +- 2 files changed, 156 insertions(+), 146 deletions(-) (limited to 'src/iso19111') diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp index e9d73f52..4cf1ae89 100644 --- a/src/iso19111/coordinateoperation.cpp +++ b/src/iso19111/coordinateoperation.cpp @@ -15545,149 +15545,6 @@ void CoordinateOperationFactory::Private::createOperationsBoundToCompound( // --------------------------------------------------------------------------- -static crs::CRSNNPtr -getResolvedCRS(const crs::CRSNNPtr &crs, - const CoordinateOperationContextNNPtr &context, - metadata::ExtentPtr &extentOut) { - const auto &authFactory = context->getAuthorityFactory(); - const auto &ids = crs->identifiers(); - const auto &name = crs->nameStr(); - - bool approxExtent; - extentOut = getExtentPossiblySynthetized(crs, approxExtent); - - // We try to "identify" the provided CRS with the ones of the database, - // but in a more restricted way that what identify() does. - // If we get a match from id in priority, and from name as a fallback, and - // that they are equivalent to the input CRS, then use the identified CRS. - // Even if they aren't equivalent, we update extentOut with the one of the - // identified CRS if our input one is absent/not reliable. - - const auto tryToIdentifyByName = [&crs, &name, &authFactory, approxExtent, - &extentOut]( - io::AuthorityFactory::ObjectType objectType) { - if (name != "unknown" && name != "unnamed") { - auto matches = authFactory->createObjectsFromName( - name, {objectType}, false, 2); - if (matches.size() == 1) { - const auto match = - util::nn_static_pointer_cast(matches.front()); - if (approxExtent || !extentOut) { - extentOut = getExtent(match); - } - if (match->isEquivalentTo( - crs.get(), util::IComparable::Criterion::EQUIVALENT)) { - return match; - } - } - } - return crs; - }; - - auto geogCRS = dynamic_cast(crs.get()); - if (geogCRS && authFactory) { - if (!ids.empty()) { - const auto tmpAuthFactory = io::AuthorityFactory::create( - authFactory->databaseContext(), *ids.front()->codeSpace()); - try { - auto resolvedCrs( - tmpAuthFactory->createGeographicCRS(ids.front()->code())); - if (approxExtent || !extentOut) { - extentOut = getExtent(resolvedCrs); - } - if (resolvedCrs->isEquivalentTo( - crs.get(), util::IComparable::Criterion::EQUIVALENT)) { - return util::nn_static_pointer_cast(resolvedCrs); - } - } catch (const std::exception &) { - } - } else { - return tryToIdentifyByName( - geogCRS->coordinateSystem()->axisList().size() == 2 - ? io::AuthorityFactory::ObjectType::GEOGRAPHIC_2D_CRS - : io::AuthorityFactory::ObjectType::GEOGRAPHIC_3D_CRS); - } - } - - auto projectedCrs = dynamic_cast(crs.get()); - if (projectedCrs && authFactory) { - if (!ids.empty()) { - const auto tmpAuthFactory = io::AuthorityFactory::create( - authFactory->databaseContext(), *ids.front()->codeSpace()); - try { - auto resolvedCrs( - tmpAuthFactory->createProjectedCRS(ids.front()->code())); - if (approxExtent || !extentOut) { - extentOut = getExtent(resolvedCrs); - } - if (resolvedCrs->isEquivalentTo( - crs.get(), util::IComparable::Criterion::EQUIVALENT)) { - return util::nn_static_pointer_cast(resolvedCrs); - } - } catch (const std::exception &) { - } - } else { - return tryToIdentifyByName( - io::AuthorityFactory::ObjectType::PROJECTED_CRS); - } - } - - auto compoundCrs = dynamic_cast(crs.get()); - if (compoundCrs && authFactory) { - if (!ids.empty()) { - const auto tmpAuthFactory = io::AuthorityFactory::create( - authFactory->databaseContext(), *ids.front()->codeSpace()); - try { - auto resolvedCrs( - tmpAuthFactory->createCompoundCRS(ids.front()->code())); - if (approxExtent || !extentOut) { - extentOut = getExtent(resolvedCrs); - } - if (resolvedCrs->isEquivalentTo( - crs.get(), util::IComparable::Criterion::EQUIVALENT)) { - return util::nn_static_pointer_cast(resolvedCrs); - } - } catch (const std::exception &) { - } - } else { - auto outCrs = tryToIdentifyByName( - io::AuthorityFactory::ObjectType::COMPOUND_CRS); - const auto &components = compoundCrs->componentReferenceSystems(); - if (outCrs.get() != crs.get()) { - bool hasGeoid = false; - if (components.size() == 2) { - auto vertCRS = - dynamic_cast(components[1].get()); - if (vertCRS && !vertCRS->geoidModel().empty()) { - hasGeoid = true; - } - } - if (!hasGeoid) { - return outCrs; - } - } - if (approxExtent || !extentOut) { - // If we still did not get a reliable extent, then try to - // resolve the components of the compoundCRS, and take the - // intersection of their extent. - extentOut = metadata::ExtentPtr(); - for (const auto &component : components) { - metadata::ExtentPtr componentExtent; - getResolvedCRS(component, context, componentExtent); - if (extentOut && componentExtent) - extentOut = extentOut->intersection( - NN_NO_CHECK(componentExtent)); - else if (componentExtent) - extentOut = componentExtent; - } - } - } - } - return crs; -} - -// --------------------------------------------------------------------------- - /** \brief Find a list of CoordinateOperation from sourceCRS to targetCRS. * * The operations are sorted with the most relevant ones first: by @@ -15723,13 +15580,14 @@ CoordinateOperationFactory::createOperations( const auto &targetBoundCRS = targetCRS->canonicalBoundCRS(); auto l_sourceCRS = srcBoundCRS ? NN_NO_CHECK(srcBoundCRS) : sourceCRS; auto l_targetCRS = targetBoundCRS ? NN_NO_CHECK(targetBoundCRS) : targetCRS; + const auto &authFactory = context->getAuthorityFactory(); metadata::ExtentPtr sourceCRSExtent; auto l_resolvedSourceCRS = - getResolvedCRS(l_sourceCRS, context, sourceCRSExtent); + crs::CRS::getResolvedCRS(l_sourceCRS, authFactory, sourceCRSExtent); metadata::ExtentPtr targetCRSExtent; auto l_resolvedTargetCRS = - getResolvedCRS(l_targetCRS, context, targetCRSExtent); + crs::CRS::getResolvedCRS(l_targetCRS, authFactory, targetCRSExtent); Private::Context contextPrivate(sourceCRSExtent, targetCRSExtent, context); if (context->getSourceAndTargetCRSExtentUse() == @@ -16054,4 +15912,152 @@ PROJBasedOperation::gridsNeeded(const io::DatabaseContextPtr &databaseContext, // --------------------------------------------------------------------------- } // namespace operation + +namespace crs { +// --------------------------------------------------------------------------- + +//! @cond Doxygen_Suppress + +crs::CRSNNPtr CRS::getResolvedCRS(const crs::CRSNNPtr &crs, + const io::AuthorityFactoryPtr &authFactory, + metadata::ExtentPtr &extentOut) { + const auto &ids = crs->identifiers(); + const auto &name = crs->nameStr(); + + bool approxExtent; + extentOut = getExtentPossiblySynthetized(crs, approxExtent); + + // We try to "identify" the provided CRS with the ones of the database, + // but in a more restricted way that what identify() does. + // If we get a match from id in priority, and from name as a fallback, and + // that they are equivalent to the input CRS, then use the identified CRS. + // Even if they aren't equivalent, we update extentOut with the one of the + // identified CRS if our input one is absent/not reliable. + + const auto tryToIdentifyByName = [&crs, &name, &authFactory, approxExtent, + &extentOut]( + io::AuthorityFactory::ObjectType objectType) { + if (name != "unknown" && name != "unnamed") { + auto matches = authFactory->createObjectsFromName( + name, {objectType}, false, 2); + if (matches.size() == 1) { + const auto match = + util::nn_static_pointer_cast(matches.front()); + if (approxExtent || !extentOut) { + extentOut = getExtent(match); + } + if (match->isEquivalentTo( + crs.get(), util::IComparable::Criterion::EQUIVALENT)) { + return match; + } + } + } + return crs; + }; + + auto geogCRS = dynamic_cast(crs.get()); + if (geogCRS && authFactory) { + if (!ids.empty()) { + const auto tmpAuthFactory = io::AuthorityFactory::create( + authFactory->databaseContext(), *ids.front()->codeSpace()); + try { + auto resolvedCrs( + tmpAuthFactory->createGeographicCRS(ids.front()->code())); + if (approxExtent || !extentOut) { + extentOut = getExtent(resolvedCrs); + } + if (resolvedCrs->isEquivalentTo( + crs.get(), util::IComparable::Criterion::EQUIVALENT)) { + return util::nn_static_pointer_cast(resolvedCrs); + } + } catch (const std::exception &) { + } + } else { + return tryToIdentifyByName( + geogCRS->coordinateSystem()->axisList().size() == 2 + ? io::AuthorityFactory::ObjectType::GEOGRAPHIC_2D_CRS + : io::AuthorityFactory::ObjectType::GEOGRAPHIC_3D_CRS); + } + } + + auto projectedCrs = dynamic_cast(crs.get()); + if (projectedCrs && authFactory) { + if (!ids.empty()) { + const auto tmpAuthFactory = io::AuthorityFactory::create( + authFactory->databaseContext(), *ids.front()->codeSpace()); + try { + auto resolvedCrs( + tmpAuthFactory->createProjectedCRS(ids.front()->code())); + if (approxExtent || !extentOut) { + extentOut = getExtent(resolvedCrs); + } + if (resolvedCrs->isEquivalentTo( + crs.get(), util::IComparable::Criterion::EQUIVALENT)) { + return util::nn_static_pointer_cast(resolvedCrs); + } + } catch (const std::exception &) { + } + } else { + return tryToIdentifyByName( + io::AuthorityFactory::ObjectType::PROJECTED_CRS); + } + } + + auto compoundCrs = dynamic_cast(crs.get()); + if (compoundCrs && authFactory) { + if (!ids.empty()) { + const auto tmpAuthFactory = io::AuthorityFactory::create( + authFactory->databaseContext(), *ids.front()->codeSpace()); + try { + auto resolvedCrs( + tmpAuthFactory->createCompoundCRS(ids.front()->code())); + if (approxExtent || !extentOut) { + extentOut = getExtent(resolvedCrs); + } + if (resolvedCrs->isEquivalentTo( + crs.get(), util::IComparable::Criterion::EQUIVALENT)) { + return util::nn_static_pointer_cast(resolvedCrs); + } + } catch (const std::exception &) { + } + } else { + auto outCrs = tryToIdentifyByName( + io::AuthorityFactory::ObjectType::COMPOUND_CRS); + const auto &components = compoundCrs->componentReferenceSystems(); + if (outCrs.get() != crs.get()) { + bool hasGeoid = false; + if (components.size() == 2) { + auto vertCRS = + dynamic_cast(components[1].get()); + if (vertCRS && !vertCRS->geoidModel().empty()) { + hasGeoid = true; + } + } + if (!hasGeoid) { + return outCrs; + } + } + if (approxExtent || !extentOut) { + // If we still did not get a reliable extent, then try to + // resolve the components of the compoundCRS, and take the + // intersection of their extent. + extentOut = metadata::ExtentPtr(); + for (const auto &component : components) { + metadata::ExtentPtr componentExtent; + getResolvedCRS(component, authFactory, componentExtent); + if (extentOut && componentExtent) + extentOut = extentOut->intersection( + NN_NO_CHECK(componentExtent)); + else if (componentExtent) + extentOut = componentExtent; + } + } + } + } + return crs; +} + +//! @endcond + +} // namespace crs NS_PROJ_END diff --git a/src/iso19111/crs.cpp b/src/iso19111/crs.cpp index edc8a71f..84b98984 100644 --- a/src/iso19111/crs.cpp +++ b/src/iso19111/crs.cpp @@ -533,8 +533,12 @@ CRSNNPtr CRS::createBoundCRSToWGS84IfPossible( auto authFactory = io::AuthorityFactory::create( NN_NO_CHECK(dbContext), authority == "any" ? std::string() : authority); + metadata::ExtentPtr extentResolved(extent); + if (!extent) { + getResolvedCRS(thisAsCRS, authFactory, extentResolved); + } auto ctxt = operation::CoordinateOperationContext::create( - authFactory, extent, 0.0); + authFactory, extentResolved, 0.0); ctxt->setAllowUseIntermediateCRS(allowIntermediateCRSUse); // ctxt->setSpatialCriterion( // operation::CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION); -- cgit v1.2.3 From 63f9f4e09ab26f1a52751cec15f394e88660b8bf Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Mon, 16 Nov 2020 13:57:09 +0100 Subject: WKT2 parsing: several fixes related to map projection parameter units - WKT2 grammar: accept PARAMETER[name,value,id] without unit - Recognize "Ellipsoid scaling factor" as a parameter with a scaling unit, and defaults to Unity when not specified - WKT2 parsing: implement the requirement of 18-010r7.html#80, that is when a map projection parameter has no explicit unit, use metre/degree/unity as the default unit (contrary to WKT1 parsing rules where they are deduced from the GeogCRS angular unit and coordinate system) --- src/iso19111/io.cpp | 15 ++++++++++++--- 1 file changed, 12 insertions(+), 3 deletions(-) (limited to 'src/iso19111') diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp index b2e5822b..623ad6f9 100644 --- a/src/iso19111/io.cpp +++ b/src/iso19111/io.cpp @@ -2944,7 +2944,8 @@ UnitOfMeasure WKTParser::Private::guessUnitForParameter( const UnitOfMeasure &defaultAngularUnit) { UnitOfMeasure unit; // scale must be first because of 'Scale factor on pseudo standard parallel' - if (ci_find(paramName, "scale") != std::string::npos) { + if (ci_find(paramName, "scale") != std::string::npos || + ci_find(paramName, "scaling factor") != std::string::npos) { unit = UnitOfMeasure::SCALE_UNITY; } else if (ci_find(paramName, "latitude") != std::string::npos || ci_find(paramName, "longitude") != std::string::npos || @@ -3876,8 +3877,16 @@ WKTParser::Private::buildProjectedCRS(const WKTNodeNNPtr &node) { return createPseudoMercator(props, NN_NO_CHECK(cartesianCS)); } - auto linearUnit = buildUnitInSubNode(node, UnitOfMeasure::Type::LINEAR); - auto angularUnit = baseGeodCRS->coordinateSystem()->axisList()[0]->unit(); + // For WKT2, if there is no explicit parameter unit, use metre for linear + // units and degree for angular units + auto linearUnit = + !isNull(conversionNode) + ? UnitOfMeasure::METRE + : buildUnitInSubNode(node, UnitOfMeasure::Type::LINEAR); + auto angularUnit = + !isNull(conversionNode) + ? UnitOfMeasure::DEGREE + : baseGeodCRS->coordinateSystem()->axisList()[0]->unit(); auto conversion = !isNull(conversionNode) -- cgit v1.2.3 From 9419b42ade1f1d35a61c05f0f7ce9aab2eeff3c7 Mon Sep 17 00:00:00 2001 From: Kristian Evers Date: Sun, 1 Nov 2020 13:34:32 +0100 Subject: Remove proj_api.h Removes proj_api.h from the public API. The contents of the header file has been moved to proj_internal.h verbatim and any references to proj_api.h has been changed to proj_internal.h. The documentation of proj_api.h has been removed. The only exception to this is the API migration guides which still mention the old API. Fixes #837 --- src/iso19111/datum.cpp | 1 - src/iso19111/factory.cpp | 1 - src/iso19111/io.cpp | 1 - 3 files changed, 3 deletions(-) (limited to 'src/iso19111') diff --git a/src/iso19111/datum.cpp b/src/iso19111/datum.cpp index e29f6319..28dbc462 100644 --- a/src/iso19111/datum.cpp +++ b/src/iso19111/datum.cpp @@ -43,7 +43,6 @@ // clang-format off #include "proj.h" #include "proj_internal.h" -#include "proj_api.h" // clang-format on #include "proj_json_streaming_writer.hpp" diff --git a/src/iso19111/factory.cpp b/src/iso19111/factory.cpp index 0ccfefc1..03fd5267 100644 --- a/src/iso19111/factory.cpp +++ b/src/iso19111/factory.cpp @@ -65,7 +65,6 @@ // clang-format off #include "proj.h" #include "proj_internal.h" -#include "proj_api.h" // clang-format on #include diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp index 623ad6f9..8c815f36 100644 --- a/src/iso19111/io.cpp +++ b/src/iso19111/io.cpp @@ -70,7 +70,6 @@ // clang-format off #include "proj.h" #include "proj_internal.h" -#include "proj_api.h" // clang-format on using namespace NS_PROJ::common; -- cgit v1.2.3 From aa46197d66ce70ece382bf955326c46b13f35864 Mon Sep 17 00:00:00 2001 From: Kristian Evers Date: Sat, 14 Nov 2020 22:55:31 +0100 Subject: Weed out proj_api.h datatypes and replace them with their proj.h counterparts --- src/iso19111/datum.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'src/iso19111') diff --git a/src/iso19111/datum.cpp b/src/iso19111/datum.cpp index 28dbc462..4d7e78d3 100644 --- a/src/iso19111/datum.cpp +++ b/src/iso19111/datum.cpp @@ -417,7 +417,7 @@ std::string PrimeMeridian::getPROJStringWellKnownName(const common::Angle &angle) { const double valRad = angle.getSIValue(); std::string projPMName; - projCtx ctxt = pj_ctx_alloc(); + PJ_CONTEXT *ctxt = pj_ctx_alloc(); auto proj_pm = proj_list_prime_meridians(); for (int i = 0; proj_pm[i].id != nullptr; ++i) { double valRefRad = dmstor_ctx(ctxt, proj_pm[i].defn, nullptr); -- cgit v1.2.3 From a74b985b5006c2d279353a245cfcb850cf7fcc94 Mon Sep 17 00:00:00 2001 From: Kristian Evers Date: Tue, 17 Nov 2020 12:54:24 +0100 Subject: Remove pj_ctx_* functions and use their proj_context counterparts --- src/iso19111/c_api.cpp | 4 ++-- src/iso19111/datum.cpp | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) (limited to 'src/iso19111') diff --git a/src/iso19111/c_api.cpp b/src/iso19111/c_api.cpp index cbbdcaa8..d6308c49 100644 --- a/src/iso19111/c_api.cpp +++ b/src/iso19111/c_api.cpp @@ -79,10 +79,10 @@ static void PROJ_NO_INLINE proj_log_error(PJ_CONTEXT *ctx, const char *function, msg += ": "; msg += text; ctx->logger(ctx->logger_app_data, PJ_LOG_ERROR, msg.c_str()); - auto previous_errno = pj_ctx_get_errno(ctx); + auto previous_errno = proj_context_errno(ctx); if (previous_errno == 0) { // only set errno if it wasn't set deeper down the call stack - pj_ctx_set_errno(ctx, PJD_ERR_GENERIC_ERROR); + proj_context_errno_set(ctx, PJD_ERR_GENERIC_ERROR); } } diff --git a/src/iso19111/datum.cpp b/src/iso19111/datum.cpp index 4d7e78d3..60244779 100644 --- a/src/iso19111/datum.cpp +++ b/src/iso19111/datum.cpp @@ -417,7 +417,7 @@ std::string PrimeMeridian::getPROJStringWellKnownName(const common::Angle &angle) { const double valRad = angle.getSIValue(); std::string projPMName; - PJ_CONTEXT *ctxt = pj_ctx_alloc(); + PJ_CONTEXT *ctxt = proj_context_create(); auto proj_pm = proj_list_prime_meridians(); for (int i = 0; proj_pm[i].id != nullptr; ++i) { double valRefRad = dmstor_ctx(ctxt, proj_pm[i].defn, nullptr); @@ -426,7 +426,7 @@ PrimeMeridian::getPROJStringWellKnownName(const common::Angle &angle) { break; } } - pj_ctx_free(ctxt); + proj_context_destroy(ctxt); return projPMName; } //! @endcond -- cgit v1.2.3 From 046270a85faf20f38d01e02942d197db74bb8542 Mon Sep 17 00:00:00 2001 From: Kristian Evers Date: Fri, 20 Nov 2020 16:37:12 +0100 Subject: Remove old pj_ memory (de)allocation functions Gone are pj_malloc, pj_calloc, pj_dalloc and pj_dealloc. Their primary function as API memory functions in proj_api.h is no longer there and the other use as a workaround for old errno problems is no longer valid either. Replaced with malloc and free across the codebase. --- src/iso19111/c_api.cpp | 20 ++++++++++---------- src/iso19111/io.cpp | 4 ++-- 2 files changed, 12 insertions(+), 12 deletions(-) (limited to 'src/iso19111') diff --git a/src/iso19111/c_api.cpp b/src/iso19111/c_api.cpp index d6308c49..8fdebcbd 100644 --- a/src/iso19111/c_api.cpp +++ b/src/iso19111/c_api.cpp @@ -2690,11 +2690,11 @@ proj_get_crs_info_list_from_database(PJ_CONTEXT *ctx, const char *auth_name, void proj_crs_info_list_destroy(PROJ_CRS_INFO **list) { if (list) { for (int i = 0; list[i] != nullptr; i++) { - pj_dalloc(list[i]->auth_name); - pj_dalloc(list[i]->code); - pj_dalloc(list[i]->name); - pj_dalloc(list[i]->area_name); - pj_dalloc(list[i]->projection_method_name); + free(list[i]->auth_name); + free(list[i]->code); + free(list[i]->name); + free(list[i]->area_name); + free(list[i]->projection_method_name); delete list[i]; } delete[] list; @@ -2784,11 +2784,11 @@ PROJ_UNIT_INFO **proj_get_units_from_database(PJ_CONTEXT *ctx, void proj_unit_list_destroy(PROJ_UNIT_INFO **list) { if (list) { for (int i = 0; list[i] != nullptr; i++) { - pj_dalloc(list[i]->auth_name); - pj_dalloc(list[i]->code); - pj_dalloc(list[i]->name); - pj_dalloc(list[i]->category); - pj_dalloc(list[i]->proj_short_name); + free(list[i]->auth_name); + free(list[i]->code); + free(list[i]->name); + free(list[i]->category); + free(list[i]->proj_short_name); delete list[i]; } delete[] list; diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp index 8c815f36..6c710697 100644 --- a/src/iso19111/io.cpp +++ b/src/iso19111/io.cpp @@ -9778,7 +9778,7 @@ PROJStringParser::createFromPROJString(const std::string &projString) { paralist *list = pj_expand_init(ctx, init); ctx->projStringParserCreateFromPROJStringRecursionCounter--; if (!list) { - pj_dealloc(init); + free(init); throw ParsingException("cannot expand " + projString); } std::string expanded; @@ -9801,7 +9801,7 @@ PROJStringParser::createFromPROJString(const std::string &projString) { } auto n = t->next; - pj_dealloc(t); + free(t); t = n; } for (const auto &pair : d->steps_[0].paramValues) { -- cgit v1.2.3 From 3dea81d9be4712aa90ba79f333338c0b3ecb21e2 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Sat, 21 Nov 2020 13:55:11 +0100 Subject: Make GeographicCRS/GeodeticCRS::isEquivalentTo() work properly when comparing to a DerivedGeographicCRS/DerivedGeodeticCRS --- src/iso19111/crs.cpp | 27 ++++++++++++++++++--------- 1 file changed, 18 insertions(+), 9 deletions(-) (limited to 'src/iso19111') diff --git a/src/iso19111/crs.cpp b/src/iso19111/crs.cpp index 84b98984..de882105 100644 --- a/src/iso19111/crs.cpp +++ b/src/iso19111/crs.cpp @@ -1924,13 +1924,21 @@ getStandardCriterion(util::IComparable::Criterion criterion) { //! @cond Doxygen_Suppress bool GeodeticCRS::_isEquivalentTo( + const util::IComparable *other, util::IComparable::Criterion criterion, + const io::DatabaseContextPtr &dbContext) const { + if (other == nullptr || !util::isOfExactType(*other)) { + return false; + } + return _isEquivalentToNoTypeCheck(other, criterion, dbContext); +} + +bool GeodeticCRS::_isEquivalentToNoTypeCheck( const util::IComparable *other, util::IComparable::Criterion criterion, const io::DatabaseContextPtr &dbContext) const { const auto standardCriterion = getStandardCriterion(criterion); - auto otherGeodCRS = dynamic_cast(other); + // TODO test velocityModel - return otherGeodCRS != nullptr && - SingleCRS::baseIsEquivalentTo(other, standardCriterion, dbContext); + return SingleCRS::baseIsEquivalentTo(other, standardCriterion, dbContext); } //! @endcond @@ -2486,12 +2494,13 @@ bool GeographicCRS::is2DPartOf3D(util::nn other, bool GeographicCRS::_isEquivalentTo( const util::IComparable *other, util::IComparable::Criterion criterion, const io::DatabaseContextPtr &dbContext) const { - auto otherGeogCRS = dynamic_cast(other); - if (otherGeogCRS == nullptr) { + if (other == nullptr || !util::isOfExactType(*other)) { return false; } + const auto standardCriterion = getStandardCriterion(criterion); - if (GeodeticCRS::_isEquivalentTo(other, standardCriterion, dbContext)) { + if (GeodeticCRS::_isEquivalentToNoTypeCheck(other, standardCriterion, + dbContext)) { return true; } if (criterion != @@ -2510,7 +2519,8 @@ bool GeographicCRS::_isEquivalentTo( cs::EllipsoidalCS::AxisOrder::LONG_EAST_LAT_NORTH ? cs::EllipsoidalCS::createLatitudeLongitude(unit) : cs::EllipsoidalCS::createLongitudeLatitude(unit)) - ->GeodeticCRS::_isEquivalentTo(other, standardCriterion, dbContext); + ->GeodeticCRS::_isEquivalentToNoTypeCheck(other, standardCriterion, + dbContext); } return false; } @@ -3889,8 +3899,7 @@ ProjectedCRS::create(const util::PropertyMap &properties, bool ProjectedCRS::_isEquivalentTo( const util::IComparable *other, util::IComparable::Criterion criterion, const io::DatabaseContextPtr &dbContext) const { - auto otherProjCRS = dynamic_cast(other); - return otherProjCRS != nullptr && + return other != nullptr && util::isOfExactType(*other) && DerivedCRS::_isEquivalentTo(other, criterion, dbContext); } -- cgit v1.2.3 From c2edd467c591f936f4eb81dfadd62a20f6bff371 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Sat, 21 Nov 2020 13:56:18 +0100 Subject: createOperation(): make it work properly when one of the CRS is a BoundCRS of a DerivedGeographicCRS (+proj=ob_tran +o_proj=lonlat +towgs84=....) --- src/iso19111/coordinateoperation.cpp | 53 +++++++++++++++++++++++++++++++++++- src/iso19111/io.cpp | 4 +-- 2 files changed, 54 insertions(+), 3 deletions(-) (limited to 'src/iso19111') diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp index 4cf1ae89..1d3b7a90 100644 --- a/src/iso19111/coordinateoperation.cpp +++ b/src/iso19111/coordinateoperation.cpp @@ -64,6 +64,30 @@ // #define DEBUG_CONCATENATED_OPERATION #if defined(DEBUG_SORT) || defined(DEBUG_CONCATENATED_OPERATION) #include + +void dumpWKT(const NS_PROJ::crs::CRS *crs); +void dumpWKT(const NS_PROJ::crs::CRS *crs) { + auto f(NS_PROJ::io::WKTFormatter::create( + NS_PROJ::io::WKTFormatter::Convention::WKT2_2019)); + std::cerr << crs->exportToWKT(f.get()) << std::endl; +} + +void dumpWKT(const NS_PROJ::crs::CRSPtr &crs); +void dumpWKT(const NS_PROJ::crs::CRSPtr &crs) { dumpWKT(crs.get()); } + +void dumpWKT(const NS_PROJ::crs::CRSNNPtr &crs); +void dumpWKT(const NS_PROJ::crs::CRSNNPtr &crs) { + dumpWKT(crs.as_nullable().get()); +} + +void dumpWKT(const NS_PROJ::crs::GeographicCRSPtr &crs); +void dumpWKT(const NS_PROJ::crs::GeographicCRSPtr &crs) { dumpWKT(crs.get()); } + +void dumpWKT(const NS_PROJ::crs::GeographicCRSNNPtr &crs); +void dumpWKT(const NS_PROJ::crs::GeographicCRSNNPtr &crs) { + dumpWKT(crs.as_nullable().get()); +} + #endif using namespace NS_PROJ::internal; @@ -9184,6 +9208,14 @@ static void setupPROJGeodeticSourceCRS(io::PROJStringFormatter *formatter, formatter->startInversion(); sourceCRSGeog->_exportToPROJString(formatter); formatter->stopInversion(); + if (util::isOfExactType( + *(sourceCRSGeog.get()))) { + // The export of a DerivedGeographicCRS in non-CRS mode adds + // unit conversion and axis swapping. We must compensate for that + formatter->startInversion(); + sourceCRSGeog->addAngularUnitConvertAndAxisSwap(formatter); + formatter->stopInversion(); + } if (addPushV3) { formatter->addStep("push"); @@ -9217,7 +9249,12 @@ static void setupPROJGeodeticTargetCRS(io::PROJStringFormatter *formatter, formatter->addStep("pop"); formatter->addParam("v_3"); } - + if (util::isOfExactType( + *(targetCRSGeog.get()))) { + // The export of a DerivedGeographicCRS in non-CRS mode adds + // unit conversion and axis swapping. We must compensate for that + targetCRSGeog->addAngularUnitConvertAndAxisSwap(formatter); + } targetCRSGeog->_exportToPROJString(formatter); } else { auto targetCRSGeod = dynamic_cast(crs.get()); @@ -14290,6 +14327,20 @@ void CoordinateOperationFactory::Private::createOperationsBoundToGeog( const auto &hubSrc = boundSrc->hubCRS(); auto hubSrcGeog = dynamic_cast(hubSrc.get()); auto geogCRSOfBaseOfBoundSrc = boundSrc->baseCRS()->extractGeographicCRS(); + { + // If geogCRSOfBaseOfBoundSrc is a DerivedGeographicCRS, use its base + // instead (if it is a GeographicCRS) + auto derivedGeogCRS = + std::dynamic_pointer_cast( + geogCRSOfBaseOfBoundSrc); + if (derivedGeogCRS) { + auto baseCRS = std::dynamic_pointer_cast( + derivedGeogCRS->baseCRS().as_nullable()); + if (baseCRS) { + geogCRSOfBaseOfBoundSrc = baseCRS; + } + } + } const auto &authFactory = context.context->getAuthorityFactory(); const auto dbContext = diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp index 623ad6f9..fa359d39 100644 --- a/src/iso19111/io.cpp +++ b/src/iso19111/io.cpp @@ -9478,8 +9478,8 @@ CRSNNPtr PROJStringParser::Private::buildProjectedCRS( std::string methodName = "PROJ " + step.name; for (const auto ¶m : step.paramValues) { if (is_in_stringlist(param.key, - "wktext,no_defs,datum,ellps,a,b,R,towgs84," - "nadgrids,geoidgrids," + "wktext,no_defs,datum,ellps,a,b,R,f,rf," + "towgs84,nadgrids,geoidgrids," "units,to_meter,vunits,vto_meter,type")) { continue; } -- cgit v1.2.3 From eb66679cde834096ff18f40b5b1d4bc10e3f4c1d Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Sat, 21 Nov 2020 22:58:05 +0100 Subject: WKT parsing: fix ingestion of WKT with a Geocentric CRS as the base of the projected CRS --- src/iso19111/io.cpp | 26 ++++++++++++++++++++++++++ 1 file changed, 26 insertions(+) (limited to 'src/iso19111') diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp index fa359d39..284729cd 100644 --- a/src/iso19111/io.cpp +++ b/src/iso19111/io.cpp @@ -2794,6 +2794,32 @@ WKTParser::Private::buildGeodeticCRS(const WKTNodeNNPtr &node) { .as_nullable() : nullptr; auto cs = buildCS(csNode, node, angularUnit); + + // If there's no CS[] node, typically for a BASEGEODCRS of a projected CRS, + // in a few rare cases, this might be a Geocentric CRS, and thus a + // Cartesian CS, and not the ellipsoidalCS we assumed above. The only way + // to figure that is to resolve the CRS from its code... + if (isNull(csNode) && dbContext_ && + ci_equal(nodeName, WKTConstants::BASEGEODCRS)) { + const auto &nodeChildren = nodeP->children(); + for (const auto &subNode : nodeChildren) { + const auto &subNodeName(subNode->GP()->value()); + if (ci_equal(subNodeName, WKTConstants::ID) || + ci_equal(subNodeName, WKTConstants::AUTHORITY)) { + auto id = buildId(subNode, true, false); + if (id) { + try { + auto authFactory = AuthorityFactory::create( + NN_NO_CHECK(dbContext_), *id->codeSpace()); + auto dbCRS = authFactory->createGeodeticCRS(id->code()); + cs = dbCRS->coordinateSystem(); + } catch (const util::Exception &) { + } + } + } + } + } + auto ellipsoidalCS = nn_dynamic_pointer_cast(cs); if (ellipsoidalCS) { if (ci_equal(nodeName, WKTConstants::GEOCCS)) { -- cgit v1.2.3 From 650e37f159b96d18f5ff03784466784a769688e4 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Sat, 21 Nov 2020 23:59:33 +0100 Subject: PROJJSON parsing: fix parsing of a Geodetic CRS with a DatumEnsemble, and fix parsing of a ProjectedCRS whose base is a Geocentric CRS --- src/iso19111/io.cpp | 33 ++++++++++++++++++++++++--------- 1 file changed, 24 insertions(+), 9 deletions(-) (limited to 'src/iso19111') diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp index 284729cd..713a471d 100644 --- a/src/iso19111/io.cpp +++ b/src/iso19111/io.cpp @@ -4980,6 +4980,10 @@ class JSONParser { TransformationNNPtr buildTransformation(const json &j); ConcatenatedOperationNNPtr buildConcatenatedOperation(const json &j); + void buildGeodeticDatumOrDatumEnsemble(const json &j, + GeodeticReferenceFramePtr &datum, + DatumEnsemblePtr &datumEnsemble); + static util::optional getAnchor(const json &j) { util::optional anchor; if (j.contains("anchor")) { @@ -5439,9 +5443,9 @@ BaseObjectNNPtr JSONParser::create(const json &j) // --------------------------------------------------------------------------- -GeographicCRSNNPtr JSONParser::buildGeographicCRS(const json &j) { - GeodeticReferenceFramePtr datum; - DatumEnsemblePtr datumEnsemble; +void JSONParser::buildGeodeticDatumOrDatumEnsemble( + const json &j, GeodeticReferenceFramePtr &datum, + DatumEnsemblePtr &datumEnsemble) { if (j.contains("datum")) { auto datumJ = getObject(j, "datum"); datum = util::nn_dynamic_pointer_cast( @@ -5454,6 +5458,14 @@ GeographicCRSNNPtr JSONParser::buildGeographicCRS(const json &j) { datumEnsemble = buildDatumEnsemble(getObject(j, "datum_ensemble")).as_nullable(); } +} + +// --------------------------------------------------------------------------- + +GeographicCRSNNPtr JSONParser::buildGeographicCRS(const json &j) { + GeodeticReferenceFramePtr datum; + DatumEnsemblePtr datumEnsemble; + buildGeodeticDatumOrDatumEnsemble(j, datum, datumEnsemble); auto csJ = getObject(j, "coordinate_system"); auto ellipsoidalCS = util::nn_dynamic_pointer_cast(buildCS(csJ)); @@ -5467,12 +5479,9 @@ GeographicCRSNNPtr JSONParser::buildGeographicCRS(const json &j) { // --------------------------------------------------------------------------- GeodeticCRSNNPtr JSONParser::buildGeodeticCRS(const json &j) { - auto datumJ = getObject(j, "datum"); - if (getType(datumJ) != "GeodeticReferenceFrame") { - throw ParsingException("Unsupported type for datum."); - } - auto datum = buildGeodeticReferenceFrame(datumJ); + GeodeticReferenceFramePtr datum; DatumEnsemblePtr datumEnsemble; + buildGeodeticDatumOrDatumEnsemble(j, datum, datumEnsemble); auto csJ = getObject(j, "coordinate_system"); auto cs = buildCS(csJ); auto props = buildProperties(j); @@ -5507,7 +5516,13 @@ GeodeticCRSNNPtr JSONParser::buildGeodeticCRS(const json &j) { // --------------------------------------------------------------------------- ProjectedCRSNNPtr JSONParser::buildProjectedCRS(const json &j) { - auto baseCRS = buildGeographicCRS(getObject(j, "base_crs")); + auto jBaseCRS = getObject(j, "base_crs"); + auto jBaseCS = getObject(jBaseCRS, "coordinate_system"); + auto baseCS = buildCS(jBaseCS); + auto baseCRS = dynamic_cast(baseCS.get()) != nullptr + ? util::nn_static_pointer_cast( + buildGeographicCRS(jBaseCRS)) + : buildGeodeticCRS(jBaseCRS); auto csJ = getObject(j, "coordinate_system"); auto cartesianCS = util::nn_dynamic_pointer_cast(buildCS(csJ)); if (!cartesianCS) { -- cgit v1.2.3 From b503c2e5ae4b90da5e390be929ab5a37ab534232 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Mon, 23 Nov 2020 16:54:48 +0100 Subject: GeographicCRS::_isEquivalentTo(EQUIVALENT_EXCEPT_AXIS_ORDER_GEOGCRS): make it work when comparing easting,northing,up and northing,easting,up --- src/iso19111/crs.cpp | 21 +++++++++++++++++++++ 1 file changed, 21 insertions(+) (limited to 'src/iso19111') diff --git a/src/iso19111/crs.cpp b/src/iso19111/crs.cpp index de882105..5d2c6d7b 100644 --- a/src/iso19111/crs.cpp +++ b/src/iso19111/crs.cpp @@ -2522,6 +2522,27 @@ bool GeographicCRS::_isEquivalentTo( ->GeodeticCRS::_isEquivalentToNoTypeCheck(other, standardCriterion, dbContext); } + if (axisOrder == + cs::EllipsoidalCS::AxisOrder::LONG_EAST_LAT_NORTH_HEIGHT_UP || + axisOrder == + cs::EllipsoidalCS::AxisOrder::LAT_NORTH_LONG_EAST_HEIGHT_UP) { + const auto &angularUnit = coordinateSystem()->axisList()[0]->unit(); + const auto &linearUnit = coordinateSystem()->axisList()[2]->unit(); + return GeographicCRS::create( + util::PropertyMap().set(common::IdentifiedObject::NAME_KEY, + nameStr()), + datum(), datumEnsemble(), + axisOrder == cs::EllipsoidalCS::AxisOrder:: + LONG_EAST_LAT_NORTH_HEIGHT_UP + ? cs::EllipsoidalCS:: + createLatitudeLongitudeEllipsoidalHeight( + angularUnit, linearUnit) + : cs::EllipsoidalCS:: + createLongitudeLatitudeEllipsoidalHeight( + angularUnit, linearUnit)) + ->GeodeticCRS::_isEquivalentToNoTypeCheck(other, standardCriterion, + dbContext); + } return false; } //! @endcond -- cgit v1.2.3 From d8fe9964bcbc6d1eeaddf6f5d47ca6d444dc8744 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Sun, 22 Nov 2020 15:08:41 +0100 Subject: Add +proj=topocentric geocentric->topocentric conversion (fixes #500) --- src/iso19111/coordinateoperation.cpp | 55 ++++++++++++++++++++++++++++++------ 1 file changed, 47 insertions(+), 8 deletions(-) (limited to 'src/iso19111') diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp index 1d3b7a90..dfd349b5 100644 --- a/src/iso19111/coordinateoperation.cpp +++ b/src/iso19111/coordinateoperation.cpp @@ -6034,28 +6034,39 @@ void Conversion::_exportToPROJString( !isHeightDepthReversal; bool applyTargetCRSModifiers = applySourceCRSModifiers; + if (formatter->getCRSExport()) { + if (methodEPSGCode == EPSG_CODE_METHOD_GEOCENTRIC_TOPOCENTRIC || + methodEPSGCode == EPSG_CODE_METHOD_GEOGRAPHIC_TOPOCENTRIC) { + throw io::FormattingException("Transformation cannot be exported " + "as a PROJ.4 string (but can be part " + "of a PROJ pipeline)"); + } + } + auto l_sourceCRS = sourceCRS(); + crs::GeographicCRSPtr srcGeogCRS; if (!formatter->getCRSExport() && l_sourceCRS && applySourceCRSModifiers) { - crs::CRS *horiz = l_sourceCRS.get(); - const auto compound = dynamic_cast(horiz); + crs::CRSPtr horiz = l_sourceCRS; + const auto compound = + dynamic_cast(l_sourceCRS.get()); if (compound) { const auto &components = compound->componentReferenceSystems(); if (!components.empty()) { - horiz = components.front().get(); + horiz = components.front().as_nullable(); } } - auto geogCRS = dynamic_cast(horiz); - if (geogCRS) { + srcGeogCRS = std::dynamic_pointer_cast(horiz); + if (srcGeogCRS) { formatter->setOmitProjLongLatIfPossible(true); formatter->startInversion(); - geogCRS->_exportToPROJString(formatter); + srcGeogCRS->_exportToPROJString(formatter); formatter->stopInversion(); formatter->setOmitProjLongLatIfPossible(false); } - auto projCRS = dynamic_cast(horiz); + auto projCRS = dynamic_cast(horiz.get()); if (projCRS) { formatter->startInversion(); formatter->pushOmitZUnitConversion(); @@ -6325,6 +6336,30 @@ void Conversion::_exportToPROJString( } bConversionDone = true; bEllipsoidParametersDone = true; + } else if (methodEPSGCode == EPSG_CODE_METHOD_GEOGRAPHIC_TOPOCENTRIC) { + if (!srcGeogCRS) { + throw io::FormattingException( + "Export of Geographic/Topocentric conversion to a PROJ string " + "requires an input geographic CRS"); + } + + formatter->addStep("cart"); + srcGeogCRS->ellipsoid()->_exportToPROJString(formatter); + + formatter->addStep("topocentric"); + const auto latOrigin = parameterValueNumeric( + EPSG_CODE_PARAMETER_LATITUDE_TOPOGRAPHIC_ORIGIN, + common::UnitOfMeasure::DEGREE); + const auto lonOrigin = parameterValueNumeric( + EPSG_CODE_PARAMETER_LONGITUDE_TOPOGRAPHIC_ORIGIN, + common::UnitOfMeasure::DEGREE); + const auto heightOrigin = parameterValueNumeric( + EPSG_CODE_PARAMETER_ELLIPSOIDAL_HEIGHT_TOPOCENTRIC_ORIGIN, + common::UnitOfMeasure::METRE); + formatter->addParam("lat_0", latOrigin); + formatter->addParam("lon_0", lonOrigin); + formatter->addParam("h_0", heightOrigin); + bConversionDone = true; } auto l_targetCRS = targetCRS(); @@ -6449,7 +6484,9 @@ void Conversion::_exportToPROJString( } if (!bEllipsoidParametersDone) { - auto targetGeogCRS = horiz->extractGeographicCRS(); + auto targetGeodCRS = horiz->extractGeodeticCRS(); + auto targetGeogCRS = + std::dynamic_pointer_cast(targetGeodCRS); if (targetGeogCRS) { if (formatter->getCRSExport()) { targetGeogCRS->addDatumInfoToPROJString(formatter); @@ -6458,6 +6495,8 @@ void Conversion::_exportToPROJString( targetGeogCRS->primeMeridian()->_exportToPROJString( formatter); } + } else if (targetGeodCRS) { + targetGeodCRS->ellipsoid()->_exportToPROJString(formatter); } } -- cgit v1.2.3 From d6f6862352d0caeb827c0c3ec98ccf9aa2f9a1d2 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Tue, 24 Nov 2020 16:41:06 +0100 Subject: createOperation(): add a ballpark vertical transformation when dealing with GEOIDMODEL[] --- src/iso19111/coordinateoperation.cpp | 99 +++++++++++++++++++++++++++++++----- 1 file changed, 85 insertions(+), 14 deletions(-) (limited to 'src/iso19111') diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp index 1d3b7a90..0f1c680b 100644 --- a/src/iso19111/coordinateoperation.cpp +++ b/src/iso19111/coordinateoperation.cpp @@ -11233,6 +11233,12 @@ struct CoordinateOperationFactory::Private { const crs::GeographicCRS *geogDst, std::vector &res); + static void createOperationsVertToGeogBallpark( + const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS, + Private::Context &context, const crs::VerticalCRS *vertSrc, + const crs::GeographicCRS *geogDst, + std::vector &res); + static void createOperationsBoundToBound( const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS, Private::Context &context, const crs::BoundCRS *boundSrc, @@ -11334,6 +11340,7 @@ struct PrecomputedOpCharacteristics { bool gridsKnown_ = false; size_t stepCount_ = 0; bool isApprox_ = false; + bool hasBallparkVertical_ = false; bool isNullTransformation_ = false; PrecomputedOpCharacteristics() = default; @@ -11341,10 +11348,12 @@ struct PrecomputedOpCharacteristics { bool isPROJExportable, bool hasGrids, bool gridsAvailable, bool gridsKnown, size_t stepCount, bool isApprox, + bool hasBallparkVertical, bool isNullTransformation) : area_(area), accuracy_(accuracy), isPROJExportable_(isPROJExportable), hasGrids_(hasGrids), gridsAvailable_(gridsAvailable), gridsKnown_(gridsKnown), stepCount_(stepCount), isApprox_(isApprox), + hasBallparkVertical_(hasBallparkVertical), isNullTransformation_(isNullTransformation) {} }; @@ -11388,6 +11397,15 @@ struct SortFunction { return false; } + if (!iterA->second.hasBallparkVertical_ && + iterB->second.hasBallparkVertical_) { + return true; + } + if (iterA->second.hasBallparkVertical_ && + !iterB->second.hasBallparkVertical_) { + return false; + } + if (!iterA->second.isNullTransformation_ && iterB->second.isNullTransformation_) { return true; @@ -11654,7 +11672,9 @@ struct FilterResults { ? CoordinateOperationContext::SpatialCriterion:: STRICT_CONTAINMENT : context->getSpatialCriterion(); - bool hasFoundOpWithExtent = false; + bool hasOnlyBallpark = true; + bool hasNonBallparkWithoutExtent = false; + bool hasNonBallparkOpWithExtent = false; const bool allowBallpark = context->getAllowBallparkTransformations(); for (const auto &op : sourceList) { if (desiredAccuracy != 0) { @@ -11669,9 +11689,15 @@ struct FilterResults { if (areaOfInterest) { bool emptyIntersection = false; auto extent = getExtent(op, true, emptyIntersection); - if (!extent) + if (!extent) { + if (!op->hasBallparkTransformation()) { + hasNonBallparkWithoutExtent = true; + } continue; - hasFoundOpWithExtent = true; + } + if (!op->hasBallparkTransformation()) { + hasNonBallparkOpWithExtent = true; + } bool extentContains = extent->contains(NN_NO_CHECK(areaOfInterest)); if (!hasOpThatContainsAreaOfInterestAndNoGrid && @@ -11698,9 +11724,15 @@ struct FilterResults { BOTH) { bool emptyIntersection = false; auto extent = getExtent(op, true, emptyIntersection); - if (!extent) + if (!extent) { + if (!op->hasBallparkTransformation()) { + hasNonBallparkWithoutExtent = true; + } continue; - hasFoundOpWithExtent = true; + } + if (!op->hasBallparkTransformation()) { + hasNonBallparkOpWithExtent = true; + } bool extentContainsExtent1 = !extent1 || extent->contains(NN_NO_CHECK(extent1)); bool extentContainsExtent2 = @@ -11730,12 +11762,16 @@ struct FilterResults { } } } + if (!op->hasBallparkTransformation()) { + hasOnlyBallpark = false; + } res.emplace_back(op); } // In case no operation has an extent and no result is found, // retain all initial operations that match accuracy criterion. - if (res.empty() && !hasFoundOpWithExtent) { + if ((res.empty() && !hasNonBallparkOpWithExtent) || + (hasOnlyBallpark && hasNonBallparkWithoutExtent)) { for (const auto &op : sourceList) { if (desiredAccuracy != 0) { const double accuracy = getAccuracy(op); @@ -11839,6 +11875,8 @@ struct FilterResults { area, getAccuracy(op), isPROJExportable, hasGrids, gridsAvailable, gridsKnown, stepCount, op->hasBallparkTransformation(), + op->nameStr().find("ballpark vertical transformation") != + std::string::npos, isNullTransformation(op->nameStr())); } @@ -13647,11 +13685,17 @@ bool CoordinateOperationFactory::Private::createOperationsFromDatabase( ENTER_FUNCTION(); if (geogSrc && vertDst) { - res = createOperationsGeogToVertFromGeoid(sourceCRS, targetCRS, vertDst, - context); + createOperationsFromDatabase(targetCRS, sourceCRS, context, geodDst, + geodSrc, geogDst, geogSrc, vertDst, + vertSrc, res); + res = applyInverse(res); } else if (geogDst && vertSrc) { res = applyInverse(createOperationsGeogToVertFromGeoid( targetCRS, sourceCRS, vertSrc, context)); + if (!res.empty()) { + createOperationsVertToGeogBallpark(sourceCRS, targetCRS, context, + vertSrc, geogDst, res); + } } if (!res.empty()) { @@ -14762,17 +14806,32 @@ void CoordinateOperationFactory::Private::createOperationsVertToGeog( match.get(), util::IComparable::Criterion::EQUIVALENT) && !match->identifiers().empty()) { - res = createOperations( + auto resTmp = createOperations( NN_NO_CHECK( util::nn_dynamic_pointer_cast( match)), targetCRS, context); + res.insert(res.end(), resTmp.begin(), resTmp.end()); return; } } } } + createOperationsVertToGeogBallpark(sourceCRS, targetCRS, context, vertSrc, + geogDst, res); +} + +// --------------------------------------------------------------------------- + +void CoordinateOperationFactory::Private::createOperationsVertToGeogBallpark( + const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS, + Private::Context &, const crs::VerticalCRS *vertSrc, + const crs::GeographicCRS *geogDst, + std::vector &res) { + + ENTER_FUNCTION(); + const auto &srcAxis = vertSrc->coordinateSystem()->axisList()[0]; const double convSrc = srcAxis->unit().conversionToSI(); double convDst = 1.0; @@ -14791,12 +14850,24 @@ void CoordinateOperationFactory::Private::createOperationsVertToGeog( ((srcIsUp && dstIsDown) || (srcIsDown && dstIsUp)); const double factor = convSrc / convDst; - auto conv = Transformation::createChangeVerticalUnit( - util::PropertyMap().set( - common::IdentifiedObject::NAME_KEY, + + const auto &sourceCRSExtent = getExtent(sourceCRS); + const auto &targetCRSExtent = getExtent(targetCRS); + const bool sameExtent = + sourceCRSExtent && targetCRSExtent && + sourceCRSExtent->_isEquivalentTo( + targetCRSExtent.get(), util::IComparable::Criterion::EQUIVALENT); + + util::PropertyMap map; + map.set(common::IdentifiedObject::NAME_KEY, buildTransfName(sourceCRS->nameStr(), targetCRS->nameStr()) + - BALLPARK_VERTICAL_TRANSFORMATION_NO_ELLIPSOID_VERT_HEIGHT), - sourceCRS, targetCRS, + BALLPARK_VERTICAL_TRANSFORMATION_NO_ELLIPSOID_VERT_HEIGHT) + .set(common::ObjectUsage::DOMAIN_OF_VALIDITY_KEY, + sameExtent ? NN_NO_CHECK(sourceCRSExtent) + : metadata::Extent::WORLD); + + auto conv = Transformation::createChangeVerticalUnit( + map, sourceCRS, targetCRS, common::Scale(heightDepthReversal ? -factor : factor), {}); conv->setHasBallparkTransformation(true); res.push_back(conv); -- cgit v1.2.3 From 4caf32aedd4da6b1fd1b1ce0e04a1a08dc1e3f33 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Tue, 24 Nov 2020 19:16:12 +0100 Subject: Add option to allow export of Geographic/Projected 3D CRS in WKT1_GDAL as CompoundCRS with a VerticalCRS being an ellipsoidal height, which is not conformant. But needed for LAS 1.4 that only supports WKT1 --- src/iso19111/c_api.cpp | 12 ++++++++++++ src/iso19111/crs.cpp | 43 +++++++++++++++++++++++++++++++++++++++++++ src/iso19111/io.cpp | 25 +++++++++++++++++++++++++ 3 files changed, 80 insertions(+) (limited to 'src/iso19111') diff --git a/src/iso19111/c_api.cpp b/src/iso19111/c_api.cpp index cbbdcaa8..d3c01e8f 100644 --- a/src/iso19111/c_api.cpp +++ b/src/iso19111/c_api.cpp @@ -1421,6 +1421,13 @@ const char *proj_get_id_code(const PJ *obj, int index) { * variants, for WKT1_GDAL for ProjectedCRS with easting/northing ordering * (otherwise stripped), but not for WKT1_ESRI. Setting to YES will output * them unconditionally, and to NO will omit them unconditionally. + *
  • STRICT=YES/NO. Default is YES. If NO, a Geographic 3D CRS can be for + * example exported as WKT1_GDAL with 3 axes, whereas this is normally not + * allowed.
    • + *
  • ALLOW_ELLIPSOIDAL_HEIGHT_AS_VERTICAL_CRS=YES/NO. Default is NO. If set + * to YES and type == PJ_WKT1_GDAL, a Geographic 3D CRS or a Projected 3D CRS + * will be exported as a compound CRS whose vertical part represents an + * ellipsoidal height (for example for use with LAS 1.4 WKT1).
    • * * @return a string, or NULL in case of error. */ @@ -1471,6 +1478,11 @@ const char *proj_as_wkt(PJ_CONTEXT *ctx, const PJ *obj, PJ_WKT_TYPE type, } } else if ((value = getOptionValue(*iter, "STRICT="))) { formatter->setStrict(ci_equal(value, "YES")); + } else if ((value = getOptionValue( + *iter, + "ALLOW_ELLIPSOIDAL_HEIGHT_AS_VERTICAL_CRS="))) { + formatter->setAllowEllipsoidalHeightAsVerticalCRS( + ci_equal(value, "YES")); } else { std::string msg("Unknown option :"); msg += *iter; diff --git a/src/iso19111/crs.cpp b/src/iso19111/crs.cpp index 5d2c6d7b..573dd6db 100644 --- a/src/iso19111/crs.cpp +++ b/src/iso19111/crs.cpp @@ -1620,6 +1620,34 @@ static bool exportAsESRIWktCompoundCRSWithEllipsoidalHeight( vertCRSList.front()->_exportToWKT(formatter); return true; } + +// --------------------------------------------------------------------------- + +// Try to format a Geographic/ProjectedCRS 3D CRS as a +// GEOGCS[]/PROJCS[],VERTCS["Ellipsoid (metre)",DATUM["Ellipsoid",2002],...] +static bool exportAsWKT1CompoundCRSWithEllipsoidalHeight( + const CRSNNPtr &base2DCRS, + const cs::CoordinateSystemAxisNNPtr &verticalAxis, + io::WKTFormatter *formatter) { + std::string verticalCRSName = "Ellipsoid ("; + verticalCRSName += verticalAxis->unit().name(); + verticalCRSName += ')'; + auto vertDatum = datum::VerticalReferenceFrame::create( + util::PropertyMap() + .set(common::IdentifiedObject::NAME_KEY, "Ellipsoid") + .set("VERT_DATUM_TYPE", "2002")); + auto vertCRS = VerticalCRS::create( + util::PropertyMap().set(common::IdentifiedObject::NAME_KEY, + verticalCRSName), + vertDatum.as_nullable(), nullptr, + cs::VerticalCS::create(util::PropertyMap(), verticalAxis)); + formatter->startNode(io::WKTConstants::COMPD_CS, false); + formatter->addQuotedString(base2DCRS->nameStr() + " + " + verticalCRSName); + base2DCRS->_exportToWKT(formatter); + vertCRS->_exportToWKT(formatter); + formatter->endNode(); + return true; +} //! @endcond // --------------------------------------------------------------------------- @@ -1687,6 +1715,13 @@ void GeodeticCRS::_exportToWKT(io::WKTFormatter *formatter) const { return; } + if (formatter->isAllowedEllipsoidalHeightAsVerticalCRS()) { + if (exportAsWKT1CompoundCRSWithEllipsoidalHeight( + geogCRS2D, axisList[2], formatter)) { + return; + } + } + io::FormattingException::Throw( "WKT1 does not support Geographic 3D CRS."); } @@ -3640,6 +3675,14 @@ void ProjectedCRS::_exportToWKT(io::WKTFormatter *formatter) const { return; } + if (!formatter->useESRIDialect() && + formatter->isAllowedEllipsoidalHeightAsVerticalCRS()) { + if (exportAsWKT1CompoundCRSWithEllipsoidalHeight( + projCRS2D, axisList[2], formatter)) { + return; + } + } + io::FormattingException::Throw( "Projected 3D CRS can only be exported since WKT2:2019"); } diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp index 713a471d..b6a09bb6 100644 --- a/src/iso19111/io.cpp +++ b/src/iso19111/io.cpp @@ -140,6 +140,7 @@ struct WKTFormatter::Private { bool primeMeridianInDegree_ = false; bool use2019Keywords_ = false; bool useESRIDialect_ = false; + bool allowEllipsoidalHeightAsVerticalCRS_ = false; OutputAxisRule outputAxis_ = WKTFormatter::OutputAxisRule::YES; }; Params params_{}; @@ -251,6 +252,8 @@ WKTFormatter::setOutputAxis(OutputAxisRule outputAxisIn) noexcept { * * The default is strict mode, in which case a FormattingException can be * thrown. + * In non-strict mode, a Geographic 3D CRS can be for example exported as + * WKT1_GDAL with 3 axes, whereas this is normally not allowed. */ WKTFormatter &WKTFormatter::setStrict(bool strictIn) noexcept { d->params_.strict_ = strictIn; @@ -264,6 +267,28 @@ bool WKTFormatter::isStrict() const noexcept { return d->params_.strict_; } // --------------------------------------------------------------------------- +/** \brief Set whether the formatter should export, in WKT1, a Geographic or + * Projected 3D CRS as a compound CRS whose vertical part represents an + * ellipsoidal height. + */ +WKTFormatter & +WKTFormatter::setAllowEllipsoidalHeightAsVerticalCRS(bool allow) noexcept { + d->params_.allowEllipsoidalHeightAsVerticalCRS_ = allow; + return *this; +} + +// --------------------------------------------------------------------------- + +/** \brief Return whether the formatter should export, in WKT1, a Geographic or + * Projected 3D CRS as a compound CRS whose vertical part represents an + * ellipsoidal height. + */ +bool WKTFormatter::isAllowedEllipsoidalHeightAsVerticalCRS() const noexcept { + return d->params_.allowEllipsoidalHeightAsVerticalCRS_; +} + +// --------------------------------------------------------------------------- + /** Returns the WKT string from the formatter. */ const std::string &WKTFormatter::toString() const { if (d->indentLevel_ > 0 || d->level_ > 0) { -- cgit v1.2.3 From b27efedef0b3deb488d322d42b3d5b652ab55662 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Thu, 26 Nov 2020 18:49:25 +0100 Subject: PRIMEM WKT handling: fixes on import for 'sexagesimal DMS' or from WKT1:GDAL/ESRI when GEOGCS UNIT != Degree; morph to ESRI the PRIMEM name on export --- src/iso19111/datum.cpp | 17 +++++++++++ src/iso19111/io.cpp | 76 +++++++++++++++++++++++++++++++++++++++++++++++--- 2 files changed, 89 insertions(+), 4 deletions(-) (limited to 'src/iso19111') diff --git a/src/iso19111/datum.cpp b/src/iso19111/datum.cpp index 60244779..8e1b8257 100644 --- a/src/iso19111/datum.cpp +++ b/src/iso19111/datum.cpp @@ -351,6 +351,23 @@ void PrimeMeridian::_exportToWKT( if (!(isWKT2 && formatter->primeMeridianOmittedIfGreenwich() && l_name == "Greenwich")) { formatter->startNode(io::WKTConstants::PRIMEM, !identifiers().empty()); + + if (formatter->useESRIDialect()) { + bool aliasFound = false; + const auto &dbContext = formatter->databaseContext(); + if (dbContext) { + auto l_alias = dbContext->getAliasFromOfficialName( + l_name, "prime_meridian", "ESRI"); + if (!l_alias.empty()) { + l_name = l_alias; + aliasFound = true; + } + } + if (!aliasFound) { + l_name = io::WKTFormatter::morphNameToESRI(l_name); + } + } + formatter->addQuotedString(l_name); const auto &l_long = longitude(); if (formatter->primeMeridianInDegree()) { diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp index 867a0d53..52d46d27 100644 --- a/src/iso19111/io.cpp +++ b/src/iso19111/io.cpp @@ -1978,14 +1978,80 @@ PrimeMeridianNNPtr WKTParser::Private::buildPrimeMeridian( try { double angleValue = asDouble(children[1]); - // Correct for GDAL WKT1 departure + // Correct for GDAL WKT1 and WKT1-ESRI departure if (name == "Paris" && std::fabs(angleValue - 2.33722917) < 1e-8 && - unit == UnitOfMeasure::GRAD) { + unit._isEquivalentTo(UnitOfMeasure::GRAD, + util::IComparable::Criterion::EQUIVALENT)) { angleValue = 2.5969213; + } else { + static const struct { + const char *name; + int deg; + int min; + double sec; + } primeMeridiansDMS[] = { + {"Lisbon", -9, 7, 54.862}, {"Bogota", -74, 4, 51.3}, + {"Madrid", -3, 41, 14.55}, {"Rome", 12, 27, 8.4}, + {"Bern", 7, 26, 22.5}, {"Jakarta", 106, 48, 27.79}, + {"Ferro", -17, 40, 0}, {"Brussels", 4, 22, 4.71}, + {"Stockholm", 18, 3, 29.8}, {"Athens", 23, 42, 58.815}, + {"Oslo", 10, 43, 22.5}, {"Paris RGS", 2, 20, 13.95}, + {"Paris_RGS", 2, 20, 13.95}}; + + // Current epsg.org output may use the EPSG:9110 "sexagesimal DMS" + // unit and a DD.MMSSsss value, but this will likely be changed to + // use decimal degree. + // Or WKT1 may for example use the Paris RGS decimal degree value + // but with a GEOGCS with UNIT["Grad"] + for (const auto &pmDef : primeMeridiansDMS) { + if (name == pmDef.name) { + double dmsAsDecimalValue = + (pmDef.deg >= 0 ? 1 : -1) * + (std::abs(pmDef.deg) + pmDef.min / 100. + + pmDef.sec / 10000.); + double dmsAsDecimalDegreeValue = + (pmDef.deg >= 0 ? 1 : -1) * + (std::abs(pmDef.deg) + pmDef.min / 60. + + pmDef.sec / 3600.); + if (std::fabs(angleValue - dmsAsDecimalValue) < 1e-8 || + std::fabs(angleValue - dmsAsDecimalDegreeValue) < + 1e-8) { + angleValue = dmsAsDecimalDegreeValue; + unit = UnitOfMeasure::DEGREE; + } + break; + } + } + } + + auto &properties = buildProperties(node); + if (dbContext_ && esriStyle_) { + std::string outTableName; + std::string codeFromAlias; + std::string authNameFromAlias; + auto authFactory = AuthorityFactory::create(NN_NO_CHECK(dbContext_), + std::string()); + auto officialName = authFactory->getOfficialNameFromAlias( + name, "prime_meridian", "ESRI", false, outTableName, + authNameFromAlias, codeFromAlias); + if (!officialName.empty()) { + properties.set(IdentifiedObject::NAME_KEY, officialName); + if (!authNameFromAlias.empty()) { + auto identifiers = ArrayOfBaseObject::create(); + identifiers->add(Identifier::create( + codeFromAlias, + PropertyMap() + .set(Identifier::CODESPACE_KEY, authNameFromAlias) + .set(Identifier::AUTHORITY_KEY, + authNameFromAlias))); + properties.set(IdentifiedObject::IDENTIFIERS_KEY, + identifiers); + } + } } Angle angle(angleValue, unit); - return PrimeMeridian::create(buildProperties(node), angle); + return PrimeMeridian::create(properties, angle); } catch (const std::exception &e) { throw buildRethrow(__FUNCTION__, e); } @@ -2737,6 +2803,9 @@ WKTParser::Private::buildGeodeticCRS(const WKTNodeNNPtr &node) { throw ParsingException("Missing DATUM or ENSEMBLE node"); } + // Do that now so that esriStyle_ can be set before buildPrimeMeridian() + auto props = buildProperties(node); + auto &dynamicNode = nodeP->lookForChild(WKTConstants::DYNAMIC); auto &csNode = nodeP->lookForChild(WKTConstants::CS_); @@ -2774,7 +2843,6 @@ WKTParser::Private::buildGeodeticCRS(const WKTNodeNNPtr &node) { angularUnit = primeMeridian->longitude().unit(); } - auto props = buildProperties(node); addExtensionProj4ToProp(nodeP, props); // No explicit AXIS node ? (WKT1) -- cgit v1.2.3 From 4de6f0bd24051ed60d9c0a98dbba1da376457f12 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Sat, 28 Nov 2020 00:07:18 +0100 Subject: createObjectsFromName(): in exact match, make looking for 'ETRS89 / UTM zone 32N' return only the exact match --- src/iso19111/factory.cpp | 31 ++++++++++++++++++++++++++++++- 1 file changed, 30 insertions(+), 1 deletion(-) (limited to 'src/iso19111') diff --git a/src/iso19111/factory.cpp b/src/iso19111/factory.cpp index 03fd5267..7c9eac00 100644 --- a/src/iso19111/factory.cpp +++ b/src/iso19111/factory.cpp @@ -6155,6 +6155,8 @@ AuthorityFactory::createObjectsFromNameEx( auto sqlRes = d->run(sql, params); bool isFirst = true; bool firstIsDeprecated = false; + bool foundExactMatch = false; + std::size_t hashCodeFirstMatch = 0; for (const auto &row : sqlRes) { const auto &name = row[3]; if (approximateMatch) { @@ -6242,11 +6244,38 @@ AuthorityFactory::createObjectsFromNameEx( } throw std::runtime_error("Unsupported table_name"); }; - res.emplace_back(PairObjectName(getObject(table_name, code), name)); + const auto obj = getObject(table_name, code); + if (metadata::Identifier::canonicalizeName(obj->nameStr()) == + canonicalizedSearchedName) { + foundExactMatch = true; + } + + const auto objPtr = obj.get(); + if (res.empty()) { + hashCodeFirstMatch = typeid(*objPtr).hash_code(); + } else if (hashCodeFirstMatch != typeid(*objPtr).hash_code()) { + hashCodeFirstMatch = 0; + } + + res.emplace_back(PairObjectName(obj, name)); if (limitResultCount > 0 && res.size() == limitResultCount) { break; } } + + // If we found a name that is an exact match, and all objects have the + // same type, and we are not in approximate mode, only keep the objet(s) + // with the exact name match. + if (foundExactMatch && hashCodeFirstMatch != 0 && !approximateMatch) { + std::list resTmp; + for (const auto &pair : res) { + if (metadata::Identifier::canonicalizeName( + pair.first->nameStr()) == canonicalizedSearchedName) { + resTmp.emplace_back(pair); + } + } + res = std::move(resTmp); + } } auto sortLambda = [](const PairObjectName &a, const PairObjectName &b) { -- cgit v1.2.3 From c6dda7637c0c49ece9216395caffa2bbd6184a92 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Sat, 28 Nov 2020 20:02:15 +0100 Subject: WKT1 import: better deal with apps.epsg.org output w.r.t datum ensemble names --- src/iso19111/io.cpp | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) (limited to 'src/iso19111') diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp index 52d46d27..125105ec 100644 --- a/src/iso19111/io.cpp +++ b/src/iso19111/io.cpp @@ -2163,9 +2163,15 @@ GeodeticReferenceFrameNNPtr WKTParser::Private::buildGeodeticReferenceFrame( return false; }; - if (name == "WGS_1984") { + // Remap GDAL WGS_1984 to EPSG v9 "World Geodetic System 1984" official + // name. + // Also remap EPSG v10 datum ensemble names to non-ensemble EPSG v9 + if (name == "WGS_1984" || name == "World Geodetic System 1984 ensemble") { properties.set(IdentifiedObject::NAME_KEY, GeodeticReferenceFrame::EPSG_6326->nameStr()); + } else if (name == "European Terrestrial Reference System 1989 ensemble") { + properties.set(IdentifiedObject::NAME_KEY, + "European Terrestrial Reference System 1989"); } else if (starts_with(name, "D_")) { esriStyle_ = true; const char *tableNameForAlias = nullptr; -- cgit v1.2.3 From 87270ae6f1840872d5f49c70258941197b0ce809 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Sat, 28 Nov 2020 23:39:33 +0100 Subject: createOperations(): get a '+proj=noop' instead of '+proj=affine +s33=-1' when attempting (non-sensical) EPSG:3855 to EPSG:4326 --- src/iso19111/coordinateoperation.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'src/iso19111') diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp index 0f1c680b..5ef32d31 100644 --- a/src/iso19111/coordinateoperation.cpp +++ b/src/iso19111/coordinateoperation.cpp @@ -14837,7 +14837,7 @@ void CoordinateOperationFactory::Private::createOperationsVertToGeogBallpark( double convDst = 1.0; const auto &geogAxis = geogDst->coordinateSystem()->axisList(); bool dstIsUp = true; - bool dstIsDown = true; + bool dstIsDown = false; if (geogAxis.size() == 3) { const auto &dstAxis = geogAxis[2]; convDst = dstAxis->unit().conversionToSI(); -- cgit v1.2.3 From eef73e3eca2959768563d2dbb95a591c68b1d807 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Mon, 30 Nov 2020 14:00:46 +0100 Subject: Database: add metadata with the version number of the database layout and check it in the code This is aimed at detecting running a PROJ version against a proj.db that is not meant to work with it. This happens sometimes in complex setups mixing PROJ versions. Hopefully this will help spotting the issue earlier. --- src/iso19111/factory.cpp | 64 ++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 64 insertions(+) (limited to 'src/iso19111') diff --git a/src/iso19111/factory.cpp b/src/iso19111/factory.cpp index 7c9eac00..2a03fd4e 100644 --- a/src/iso19111/factory.cpp +++ b/src/iso19111/factory.cpp @@ -93,6 +93,12 @@ namespace io { #define GEOG_3D_SINGLE_QUOTED "'geographic 3D'" #define GEOCENTRIC_SINGLE_QUOTED "'geocentric'" +// See data/sql/metadata.sql for the semantics of those constants +constexpr int DATABASE_LAYOUT_VERSION_MAJOR = 1; +// If the code depends on the new additions, then DATABASE_LAYOUT_VERSION_MINOR +// must be incremented. +constexpr int DATABASE_LAYOUT_VERSION_MINOR = 0; + // --------------------------------------------------------------------------- struct SQLValues { @@ -277,6 +283,8 @@ struct DatabaseContext::Private { lru11::Cache> cacheAliasNames_{ CACHE_SIZE}; + void checkDatabaseLayout(); + static void insertIntoCache(LRUCacheOfObjects &cache, const std::string &code, const util::BaseObjectPtr &obj); @@ -571,6 +579,61 @@ void DatabaseContext::Private::open(const std::string &databasePath, // --------------------------------------------------------------------------- +void DatabaseContext::Private::checkDatabaseLayout() { + auto res = run("SELECT key, value FROM metadata WHERE key IN " + "('DATABASE.LAYOUT.VERSION.MAJOR', " + "'DATABASE.LAYOUT.VERSION.MINOR')"); + if (res.size() != 2) { + // The database layout of PROJ 7.2 that shipped with EPSG v10.003 is + // at the time of writing still compatible of the one we support. + static_assert( + // cppcheck-suppress knownConditionTrueFalse + DATABASE_LAYOUT_VERSION_MAJOR == 1 && + // cppcheck-suppress knownConditionTrueFalse + DATABASE_LAYOUT_VERSION_MINOR == 0, + "remove that assertion and below lines next time we upgrade " + "database structure"); + res = run("SELECT 1 FROM metadata WHERE key = 'EPSG.VERSION' AND " + "value = 'v10.003'"); + if (!res.empty()) { + return; + } + + throw FactoryException( + databasePath_ + + " lacks DATABASE.LAYOUT.VERSION.MAJOR / " + "DATABASE.LAYOUT.VERSION.MINOR " + "metadata. It comes from another PROJ installation."); + } + int nMajor = 0; + int nMinor = 0; + for (const auto &row : res) { + if (row[0] == "DATABASE.LAYOUT.VERSION.MAJOR") { + nMajor = atoi(row[1].c_str()); + } else if (row[0] == "DATABASE.LAYOUT.VERSION.MINOR") { + nMinor = atoi(row[1].c_str()); + } + } + if (nMajor != DATABASE_LAYOUT_VERSION_MAJOR) { + throw FactoryException(databasePath_ + + " contains DATABASE.LAYOUT.VERSION.MAJOR = " + + toString(nMajor) + " whereas " + + toString(DATABASE_LAYOUT_VERSION_MAJOR) + + " is expected. " + "It comes from another PROJ installation."); + } + if (nMinor < DATABASE_LAYOUT_VERSION_MINOR) { + throw FactoryException(databasePath_ + + " contains DATABASE.LAYOUT.VERSION.MINOR = " + + toString(nMinor) + " whereas a number >= " + + toString(DATABASE_LAYOUT_VERSION_MINOR) + + " is expected. " + "It comes from another PROJ installation."); + } +} + +// --------------------------------------------------------------------------- + void DatabaseContext::Private::setHandle(sqlite3 *sqlite_handle) { assert(sqlite_handle); @@ -888,6 +951,7 @@ DatabaseContext::create(const std::string &databasePath, if (!auxiliaryDatabasePaths.empty()) { dbCtx->getPrivate()->attachExtraDatabases(auxiliaryDatabasePaths); } + dbCtx->getPrivate()->checkDatabaseLayout(); return dbCtx; } -- cgit v1.2.3