From 82695cce869e1bcf2eb2a8ff078b679b6a21c663 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Thu, 8 Oct 2020 17:08:27 +0200 Subject: [PATCH 1/4] When reading from database, possibly return Geographic/GeodeticCRS with a DatumEnsemble, typically for WGS 84 and ETRS89 ('breaking change') --- include/proj/io.hpp | 6 ++ src/iso19111/factory.cpp | 158 ++++++++++++++++++++++++--------- src/iso19111/io.cpp | 4 + test/cli/testprojinfo_out.dist | 116 +++++++++++++++++++++--- test/unit/test_c_api.cpp | 11 ++- test/unit/test_factory.cpp | 18 ++-- 6 files changed, 248 insertions(+), 65 deletions(-) diff --git a/include/proj/io.hpp b/include/proj/io.hpp index 20dcedfe06..35e924ecdd 100644 --- a/include/proj/io.hpp +++ b/include/proj/io.hpp @@ -1225,6 +1225,12 @@ class PROJ_GCC_DLL AuthorityFactory { private: PROJ_OPAQUE_PRIVATE_DATA + + PROJ_INTERNAL void + createGeodeticDatumOrEnsemble(const std::string &code, + datum::GeodeticReferenceFramePtr &outDatum, + datum::DatumEnsemblePtr &outDatumEnsemble, + bool turnEnsembleAsDatum) const; }; // --------------------------------------------------------------------------- diff --git a/src/iso19111/factory.cpp b/src/iso19111/factory.cpp index 5d02aeead4..cce80b8a89 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 f4ec7740d6..fde697ad05 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"; } diff --git a/test/cli/testprojinfo_out.dist b/test/cli/testprojinfo_out.dist index 7da5add7f6..2d4990ee85 100644 --- a/test/cli/testprojinfo_out.dist +++ b/test/cli/testprojinfo_out.dist @@ -4,9 +4,16 @@ PROJ.4 string: WKT2:2019 string: GEOGCRS["WGS 84", - DATUM["World Geodetic System 1984", + ENSEMBLE["World Geodetic System 1984 ensemble", + MEMBER["World Geodetic System 1984 (Transit)"], + MEMBER["World Geodetic System 1984 (G730)"], + MEMBER["World Geodetic System 1984 (G873)"], + MEMBER["World Geodetic System 1984 (G1150)"], + MEMBER["World Geodetic System 1984 (G1674)"], + MEMBER["World Geodetic System 1984 (G1762)"], ELLIPSOID["WGS 84",6378137,298.257223563, - LENGTHUNIT["metre",1]]], + LENGTHUNIT["metre",1]], + ENSEMBLEACCURACY[2.0]], PRIMEM["Greenwich",0, ANGLEUNIT["degree",0.0174532925199433]], CS[ellipsoidal,2], @@ -59,9 +66,16 @@ GEODCRS["WGS 84", Testing projinfo -o WKT2_2019 EPSG:4326 WKT2:2019 string: GEOGCRS["WGS 84", - DATUM["World Geodetic System 1984", + ENSEMBLE["World Geodetic System 1984 ensemble", + MEMBER["World Geodetic System 1984 (Transit)"], + MEMBER["World Geodetic System 1984 (G730)"], + MEMBER["World Geodetic System 1984 (G873)"], + MEMBER["World Geodetic System 1984 (G1150)"], + MEMBER["World Geodetic System 1984 (G1674)"], + MEMBER["World Geodetic System 1984 (G1762)"], ELLIPSOID["WGS 84",6378137,298.257223563, - LENGTHUNIT["metre",1]]], + LENGTHUNIT["metre",1]], + ENSEMBLEACCURACY[2.0]], PRIMEM["Greenwich",0, ANGLEUNIT["degree",0.0174532925199433]], CS[ellipsoidal,2], @@ -102,9 +116,16 @@ GEODCRS["WGS 84", WKT2:2019 string: GEOGCRS["WGS 84", - DATUM["World Geodetic System 1984", + ENSEMBLE["World Geodetic System 1984 ensemble", + MEMBER["World Geodetic System 1984 (Transit)"], + MEMBER["World Geodetic System 1984 (G730)"], + MEMBER["World Geodetic System 1984 (G873)"], + MEMBER["World Geodetic System 1984 (G1150)"], + MEMBER["World Geodetic System 1984 (G1674)"], + MEMBER["World Geodetic System 1984 (G1762)"], ELLIPSOID["WGS 84",6378137,298.257223563, - LENGTHUNIT["metre",1]]], + LENGTHUNIT["metre",1]], + ENSEMBLEACCURACY[2.0]], PRIMEM["Greenwich",0, ANGLEUNIT["degree",0.0174532925199433]], CS[ellipsoidal,2], @@ -142,13 +163,61 @@ PROJJSON: "$schema": "https://proj.org/schemas/v0.2/projjson.schema.json", "type": "GeographicCRS", "name": "WGS 84", - "datum": { - "type": "GeodeticReferenceFrame", - "name": "World Geodetic System 1984", + "datum_ensemble": { + "name": "World Geodetic System 1984 ensemble", + "members": [ + { + "name": "World Geodetic System 1984 (Transit)", + "id": { + "authority": "EPSG", + "code": 1166 + } + }, + { + "name": "World Geodetic System 1984 (G730)", + "id": { + "authority": "EPSG", + "code": 1152 + } + }, + { + "name": "World Geodetic System 1984 (G873)", + "id": { + "authority": "EPSG", + "code": 1153 + } + }, + { + "name": "World Geodetic System 1984 (G1150)", + "id": { + "authority": "EPSG", + "code": 1154 + } + }, + { + "name": "World Geodetic System 1984 (G1674)", + "id": { + "authority": "EPSG", + "code": 1155 + } + }, + { + "name": "World Geodetic System 1984 (G1762)", + "id": { + "authority": "EPSG", + "code": 1156 + } + } + ], "ellipsoid": { "name": "WGS 84", "semi_major_axis": 6378137, "inverse_flattening": 298.257223563 + }, + "accuracy": "2.0", + "id": { + "authority": "EPSG", + "code": 6326 } }, "coordinate_system": { @@ -1023,9 +1092,16 @@ PROJ.4 string: WKT2:2019 string: PROJCRS["WGS 84 / UTM zone 31N", BASEGEOGCRS["WGS 84", - DATUM["World Geodetic System 1984", + ENSEMBLE["World Geodetic System 1984 ensemble", + MEMBER["World Geodetic System 1984 (Transit)"], + MEMBER["World Geodetic System 1984 (G730)"], + MEMBER["World Geodetic System 1984 (G873)"], + MEMBER["World Geodetic System 1984 (G1150)"], + MEMBER["World Geodetic System 1984 (G1674)"], + MEMBER["World Geodetic System 1984 (G1762)"], ELLIPSOID["WGS 84",6378137,298.257223563, - LENGTHUNIT["metre",1]]], + LENGTHUNIT["metre",1]], + ENSEMBLEACCURACY[2.0]], PRIMEM["Greenwich",0, ANGLEUNIT["degree",0.0174532925199433]], ID["EPSG",4979]], @@ -1198,9 +1274,23 @@ Testing -k operation EPSG:8457 -o PROJ -q Testing D_WGS_1984 WKT2:2019 string: -DATUM["World Geodetic System 1984", +ENSEMBLE["World Geodetic System 1984 ensemble", + MEMBER["World Geodetic System 1984 (Transit)", + ID["EPSG",1166]], + MEMBER["World Geodetic System 1984 (G730)", + ID["EPSG",1152]], + MEMBER["World Geodetic System 1984 (G873)", + ID["EPSG",1153]], + MEMBER["World Geodetic System 1984 (G1150)", + ID["EPSG",1154]], + MEMBER["World Geodetic System 1984 (G1674)", + ID["EPSG",1155]], + MEMBER["World Geodetic System 1984 (G1762)", + ID["EPSG",1156]], ELLIPSOID["WGS 84",6378137,298.257223563, - LENGTHUNIT["metre",1]], + LENGTHUNIT["metre",1], + ID["EPSG",7030]], + ENSEMBLEACCURACY[2.0], ID["EPSG",6326]] Testing -k datum D_WGS_1984 diff --git a/test/unit/test_c_api.cpp b/test/unit/test_c_api.cpp index 5da6b369f7..a72aa72efd 100644 --- a/test/unit/test_c_api.cpp +++ b/test/unit/test_c_api.cpp @@ -4841,9 +4841,16 @@ TEST_F(CApi, proj_create_derived_geographic_crs) { const char *expected_wkt = "GEOGCRS[\"my rotated CRS\",\n" " BASEGEOGCRS[\"WGS 84\",\n" - " DATUM[\"World Geodetic System 1984\",\n" + " ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n" + " MEMBER[\"World Geodetic System 1984 (Transit)\"],\n" + " MEMBER[\"World Geodetic System 1984 (G730)\"],\n" + " MEMBER[\"World Geodetic System 1984 (G873)\"],\n" + " MEMBER[\"World Geodetic System 1984 (G1150)\"],\n" + " MEMBER[\"World Geodetic System 1984 (G1674)\"],\n" + " MEMBER[\"World Geodetic System 1984 (G1762)\"],\n" " ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n" - " LENGTHUNIT[\"metre\",1]]],\n" + " LENGTHUNIT[\"metre\",1]],\n" + " ENSEMBLEACCURACY[2.0]],\n" " PRIMEM[\"Greenwich\",0,\n" " ANGLEUNIT[\"degree\",0.0174532925199433]]],\n" " DERIVINGCONVERSION[\"Pole rotation (GRIB convention)\",\n" diff --git a/test/unit/test_factory.cpp b/test/unit/test_factory.cpp index 1005d49bf3..c869fa50d5 100644 --- a/test/unit/test_factory.cpp +++ b/test/unit/test_factory.cpp @@ -531,8 +531,10 @@ TEST(factory, AuthorityFactory_createGeodeticCRS_geographic2D) { EXPECT_EQ(gcrs->identifiers()[0]->code(), "4326"); EXPECT_EQ(*(gcrs->identifiers()[0]->codeSpace()), "EPSG"); EXPECT_EQ(*(gcrs->name()->description()), "WGS 84"); - EXPECT_TRUE( - gcrs->datum()->isEquivalentTo(factory->createDatum("6326").get())); + ASSERT_TRUE(gcrs->datum() == nullptr); + ASSERT_TRUE(gcrs->datumEnsemble() != nullptr); + EXPECT_TRUE(gcrs->datumEnsemble()->isEquivalentTo( + factory->createDatumEnsemble("6326").get())); EXPECT_TRUE(gcrs->coordinateSystem()->isEquivalentTo( factory->createCoordinateSystem("6422").get())); auto domain = crs->domains()[0]; @@ -566,8 +568,10 @@ TEST(factory, AuthorityFactory_createGeodeticCRS_geographic3D) { EXPECT_EQ(gcrs->identifiers()[0]->code(), "4979"); EXPECT_EQ(*(gcrs->identifiers()[0]->codeSpace()), "EPSG"); EXPECT_EQ(*(gcrs->name()->description()), "WGS 84"); - EXPECT_TRUE( - gcrs->datum()->isEquivalentTo(factory->createDatum("6326").get())); + ASSERT_TRUE(gcrs->datum() == nullptr); + ASSERT_TRUE(gcrs->datumEnsemble() != nullptr); + EXPECT_TRUE(gcrs->datumEnsemble()->isEquivalentTo( + factory->createDatumEnsemble("6326").get())); EXPECT_TRUE(gcrs->coordinateSystem()->isEquivalentTo( factory->createCoordinateSystem("6423").get())); } @@ -582,8 +586,10 @@ TEST(factory, AuthorityFactory_createGeodeticCRS_geocentric) { EXPECT_EQ(crs->identifiers()[0]->code(), "4978"); EXPECT_EQ(*(crs->identifiers()[0]->codeSpace()), "EPSG"); EXPECT_EQ(*(crs->name()->description()), "WGS 84"); - EXPECT_TRUE( - crs->datum()->isEquivalentTo(factory->createDatum("6326").get())); + ASSERT_TRUE(crs->datum() == nullptr); + ASSERT_TRUE(crs->datumEnsemble() != nullptr); + EXPECT_TRUE(crs->datumEnsemble()->isEquivalentTo( + factory->createDatumEnsemble("6326").get())); EXPECT_TRUE(crs->coordinateSystem()->isEquivalentTo( factory->createCoordinateSystem("6500").get())); } From c2b0dcc468b4e722e46fe10fca93fe70a95fcb8e Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Thu, 8 Oct 2020 18:41:57 +0200 Subject: [PATCH 2/4] 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' --- include/proj/io.hpp | 6 +++ src/iso19111/factory.cpp | 84 ++++++++++++++++++++++++++++---------- test/unit/test_factory.cpp | 14 +++++++ 3 files changed, 83 insertions(+), 21 deletions(-) diff --git a/include/proj/io.hpp b/include/proj/io.hpp index 35e924ecdd..a30b04dcd9 100644 --- a/include/proj/io.hpp +++ b/include/proj/io.hpp @@ -1231,6 +1231,12 @@ class PROJ_GCC_DLL AuthorityFactory { datum::GeodeticReferenceFramePtr &outDatum, datum::DatumEnsemblePtr &outDatumEnsemble, bool turnEnsembleAsDatum) const; + + PROJ_INTERNAL void + createVerticalDatumOrEnsemble(const std::string &code, + datum::VerticalReferenceFramePtr &outDatum, + datum::DatumEnsemblePtr &outDatumEnsemble, + bool turnEnsembleAsDatum) const; }; // --------------------------------------------------------------------------- diff --git a/src/iso19111/factory.cpp b/src/iso19111/factory.cpp index cce80b8a89..301b88be0a 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; } diff --git a/test/unit/test_factory.cpp b/test/unit/test_factory.cpp index c869fa50d5..681ac81014 100644 --- a/test/unit/test_factory.cpp +++ b/test/unit/test_factory.cpp @@ -619,6 +619,20 @@ TEST(factory, AuthorityFactory_createVerticalCRS) { // --------------------------------------------------------------------------- +TEST(factory, AuthorityFactory_createVerticalCRS_with_datum_ensemble) { + auto factory = AuthorityFactory::create(DatabaseContext::create(), "EPSG"); + EXPECT_THROW(factory->createVerticalCRS("-1"), + NoSuchAuthorityCodeException); + + auto crs = factory->createVerticalCRS("9451"); // BI height + ASSERT_TRUE(crs->datum() == nullptr); + ASSERT_TRUE(crs->datumEnsemble() != nullptr); + EXPECT_TRUE(crs->datumEnsemble()->isEquivalentTo( + factory->createDatumEnsemble("1288").get())); +} + +// --------------------------------------------------------------------------- + TEST(factory, AuthorityFactory_createConversion) { auto factory = AuthorityFactory::create(DatabaseContext::create(), "EPSG"); EXPECT_THROW(factory->createConversion("-1"), NoSuchAuthorityCodeException); From 1e5acb00a0c0fc2533b9bce2e5803da10ed1d8d6 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Thu, 8 Oct 2020 20:59:19 +0200 Subject: [PATCH 3/4] projinfo / createObjectsFromName(): support returning a datum ensemble --- docs/source/apps/projinfo.rst | 4 +-- include/proj/io.hpp | 2 ++ src/apps/projinfo.cpp | 18 ++++++++---- src/iso19111/c_api.cpp | 2 +- src/iso19111/datum.cpp | 8 ++++++ src/iso19111/factory.cpp | 51 +++++++++++++++++++++++++++++++++- src/iso19111/io.cpp | 5 ++++ test/cli/testprojinfo | 4 +++ test/cli/testprojinfo_out.dist | 21 ++++++++++++++ test/unit/test_factory.cpp | 24 ++++++++++++++++ test/unit/test_io.cpp | 8 ++++++ 11 files changed, 137 insertions(+), 10 deletions(-) diff --git a/docs/source/apps/projinfo.rst b/docs/source/apps/projinfo.rst index 05184ef99a..a07cf70951 100644 --- a/docs/source/apps/projinfo.rst +++ b/docs/source/apps/projinfo.rst @@ -16,7 +16,7 @@ Synopsis ******** | **projinfo** - | [-o formats] [-k crs|operation|datum|ellipsoid] [--summary] [-q] + | [-o formats] [-k crs|operation|datum|ensemble|ellipsoid] [--summary] [-q] | [[--area name_or_code] | [--bbox west_long,south_lat,east_long,north_lat]] | [--spatial-test contains|intersects] | [--crs-extent-use none|both|intersection|smallest] @@ -86,7 +86,7 @@ The following control parameters can appear in any order: .. note:: Before PROJ 6.3.0, WKT1:GDAL was implicitly calling --boundcrs-to-wgs84. This is no longer the case. -.. option:: -k crs|operation|datum|ellipsoid +.. option:: -k crs|operation|datum|ensemble|ellipsoid When used to query a single object with a AUTHORITY:CODE, determines the (k)ind of the object in case there are CRS, coordinate operations or ellipsoids with the same CODE. diff --git a/include/proj/io.hpp b/include/proj/io.hpp index a30b04dcd9..de8a90fca7 100644 --- a/include/proj/io.hpp +++ b/include/proj/io.hpp @@ -1032,6 +1032,8 @@ class PROJ_GCC_DLL AuthorityFactory { DYNAMIC_GEODETIC_REFERENCE_FRAME, /** Object of type datum::DynamicVerticalReferenceFrame */ DYNAMIC_VERTICAL_REFERENCE_FRAME, + /** Object of type datum::DatumEnsemble */ + DATUM_ENSEMBLE, }; PROJ_DLL std::set diff --git a/src/apps/projinfo.cpp b/src/apps/projinfo.cpp index 37b76346dd..6ee5925a69 100644 --- a/src/apps/projinfo.cpp +++ b/src/apps/projinfo.cpp @@ -79,7 +79,8 @@ struct OutputOptions { static void usage() { std::cerr - << "usage: projinfo [-o formats] [-k crs|operation|datum|ellipsoid] " + << "usage: projinfo [-o formats] " + "[-k crs|operation|datum|ensemble|ellipsoid] " "[--summary] [-q]" << std::endl << " ([--area name_or_code] | " @@ -184,11 +185,11 @@ static BaseObjectNNPtr buildObject( auto urn = "urn:ogc:def:coordinateOperation:" + tokens[0] + "::" + tokens[1]; obj = createFromUserInput(urn, dbContext).as_nullable(); - } else if (kind == "ellipsoid" && tokens.size() == 2) { - auto urn = "urn:ogc:def:ellipsoid:" + tokens[0] + "::" + tokens[1]; - obj = createFromUserInput(urn, dbContext).as_nullable(); - } else if (kind == "datum" && tokens.size() == 2) { - auto urn = "urn:ogc:def:datum:" + tokens[0] + "::" + tokens[1]; + } else if ((kind == "ellipsoid" || kind == "datum" || + kind == "ensemble") && + tokens.size() == 2) { + auto urn = + "urn:ogc:def:" + kind + ":" + tokens[0] + "::" + tokens[1]; obj = createFromUserInput(urn, dbContext).as_nullable(); } else { // Convenience to be able to use C escaped strings... @@ -222,6 +223,9 @@ static BaseObjectNNPtr buildObject( AuthorityFactory::ObjectType::ELLIPSOID); else if (kind == "datum") allowedTypes.push_back(AuthorityFactory::ObjectType::DATUM); + else if (kind == "ensemble") + allowedTypes.push_back( + AuthorityFactory::ObjectType::DATUM_ENSEMBLE); constexpr size_t limitResultCount = 10; auto factory = AuthorityFactory::create(NN_NO_CHECK(dbContext), std::string()); @@ -929,6 +933,8 @@ int main(int argc, char **argv) { objectKind = "ellipsoid"; } else if (ci_equal(kind, "datum")) { objectKind = "datum"; + } else if (ci_equal(kind, "ensemble")) { + objectKind = "ensemble"; } else { std::cerr << "Unrecognized value for option -k: " << kind << std::endl; diff --git a/src/iso19111/c_api.cpp b/src/iso19111/c_api.cpp index 90a414c68e..4fb73d3199 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 5bc8074cbd..e29f63194f 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 301b88be0a..0ccfefc129 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 fde697ad05..b2e5822b65 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); diff --git a/test/cli/testprojinfo b/test/cli/testprojinfo index 62713af581..c31cfef078 100755 --- a/test/cli/testprojinfo +++ b/test/cli/testprojinfo @@ -171,6 +171,10 @@ echo 'Testing -k datum EPSG:6326' >> ${OUT} $EXE -k datum EPSG:6326 >>${OUT} 2>&1 echo "" >>${OUT} +echo 'Testing -k ensemble WGS84' >> ${OUT} +$EXE -k ensemble WGS84 >>${OUT} 2>&1 +echo "" >>${OUT} + echo 'Testing -k operation EPSG:8457 -o PROJ -q' >> ${OUT} $EXE -k operation EPSG:8457 -o PROJ -q >>${OUT} 2>&1 echo "" >>${OUT} diff --git a/test/cli/testprojinfo_out.dist b/test/cli/testprojinfo_out.dist index 2d4990ee85..174f11d331 100644 --- a/test/cli/testprojinfo_out.dist +++ b/test/cli/testprojinfo_out.dist @@ -1260,6 +1260,27 @@ DATUM["World Geodetic System 1984", LENGTHUNIT["metre",1]], ID["EPSG",6326]] +Testing -k ensemble WGS84 +WKT2:2019 string: +ENSEMBLE["World Geodetic System 1984 ensemble", + MEMBER["World Geodetic System 1984 (Transit)", + ID["EPSG",1166]], + MEMBER["World Geodetic System 1984 (G730)", + ID["EPSG",1152]], + MEMBER["World Geodetic System 1984 (G873)", + ID["EPSG",1153]], + MEMBER["World Geodetic System 1984 (G1150)", + ID["EPSG",1154]], + MEMBER["World Geodetic System 1984 (G1674)", + ID["EPSG",1155]], + MEMBER["World Geodetic System 1984 (G1762)", + ID["EPSG",1156]], + ELLIPSOID["WGS 84",6378137,298.257223563, + LENGTHUNIT["metre",1], + ID["EPSG",7030]], + ENSEMBLEACCURACY[2.0], + ID["EPSG",6326]] + Testing -k operation EPSG:8457 -o PROJ -q +proj=pipeline +step +proj=axisswap +order=2,1 diff --git a/test/unit/test_factory.cpp b/test/unit/test_factory.cpp index 681ac81014..8e9b7ab69f 100644 --- a/test/unit/test_factory.cpp +++ b/test/unit/test_factory.cpp @@ -3186,6 +3186,29 @@ TEST(factory, createObjectsFromName) { .size(), 1U); + { + auto res = factory->createObjectsFromName( + "World Geodetic System 1984 ensemble", + {AuthorityFactory::ObjectType::DATUM_ENSEMBLE}, false); + EXPECT_EQ(res.size(), 1U); + if (!res.empty()) { + EXPECT_EQ(res.front()->getEPSGCode(), 6326); + EXPECT_TRUE(dynamic_cast(res.front().get()) != + nullptr); + } + } + + { + auto res = factory->createObjectsFromName( + "World Geodetic System 1984 ensemble", {}, false); + EXPECT_EQ(res.size(), 1U); + if (!res.empty()) { + EXPECT_EQ(res.front()->getEPSGCode(), 6326); + EXPECT_TRUE(dynamic_cast(res.front().get()) != + nullptr); + } + } + const auto types = std::vector{ AuthorityFactory::ObjectType::PRIME_MERIDIAN, AuthorityFactory::ObjectType::ELLIPSOID, @@ -3207,6 +3230,7 @@ TEST(factory, createObjectsFromName) { AuthorityFactory::ObjectType::CONVERSION, AuthorityFactory::ObjectType::TRANSFORMATION, AuthorityFactory::ObjectType::CONCATENATED_OPERATION, + AuthorityFactory::ObjectType::DATUM_ENSEMBLE, }; for (const auto type : types) { factory->createObjectsFromName("i_dont_exist", {type}, false, 1); diff --git a/test/unit/test_io.cpp b/test/unit/test_io.cpp index 555d115969..3d8df9982f 100644 --- a/test/unit/test_io.cpp +++ b/test/unit/test_io.cpp @@ -10444,6 +10444,14 @@ TEST(io, createFromUserInput) { ParsingException); EXPECT_THROW(createFromUserInput("foobar + EGM96 height", dbContext), ParsingException); + + { + auto obj = createFromUserInput("World Geodetic System 1984 ensemble", + dbContext); + auto ensemble = nn_dynamic_pointer_cast(obj); + ASSERT_TRUE(ensemble != nullptr); + EXPECT_EQ(ensemble->identifiers().size(), 1U); + } } // --------------------------------------------------------------------------- From 119888b041258267768d632b89395e7074323326 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Sat, 10 Oct 2020 15:35:29 +0200 Subject: [PATCH 4/4] proj.h: add PJ_CATEGORY_DATUM_ENSEMBLE for proj_create_from_database() --- src/iso19111/c_api.cpp | 3 +++ src/proj.h | 3 ++- test/unit/test_c_api.cpp | 7 +++++++ 3 files changed, 12 insertions(+), 1 deletion(-) diff --git a/src/iso19111/c_api.cpp b/src/iso19111/c_api.cpp index 4fb73d3199..cbbdcaa8d3 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) { diff --git a/src/proj.h b/src/proj.h index e77a59841c..4185ddbf73 100644 --- a/src/proj.h +++ b/src/proj.h @@ -694,7 +694,8 @@ typedef enum PJ_CATEGORY_PRIME_MERIDIAN, PJ_CATEGORY_DATUM, PJ_CATEGORY_CRS, - PJ_CATEGORY_COORDINATE_OPERATION + PJ_CATEGORY_COORDINATE_OPERATION, + PJ_CATEGORY_DATUM_ENSEMBLE } PJ_CATEGORY; /** \brief Object type. */ diff --git a/test/unit/test_c_api.cpp b/test/unit/test_c_api.cpp index a72aa72efd..a7c1eb535e 100644 --- a/test/unit/test_c_api.cpp +++ b/test/unit/test_c_api.cpp @@ -844,6 +844,13 @@ TEST_F(CApi, proj_create_from_database) { ObjectKeeper keeper(datum); EXPECT_EQ(proj_get_type(datum), PJ_TYPE_GEODETIC_REFERENCE_FRAME); } + { + auto ensemble = proj_create_from_database( + m_ctxt, "EPSG", "6326", PJ_CATEGORY_DATUM_ENSEMBLE, false, nullptr); + ASSERT_NE(ensemble, nullptr); + ObjectKeeper keeper(ensemble); + EXPECT_EQ(proj_get_type(ensemble), PJ_TYPE_DATUM_ENSEMBLE); + } { // International Terrestrial Reference Frame 2008 auto datum = proj_create_from_database(