diff options
| -rw-r--r-- | include/proj/crs.hpp | 10 | ||||
| -rw-r--r-- | src/apps/projinfo.cpp | 4 | ||||
| -rw-r--r-- | src/iso19111/crs.cpp | 142 | ||||
| -rw-r--r-- | src/iso19111/io.cpp | 17 | ||||
| -rw-r--r-- | test/unit/test_io.cpp | 142 |
5 files changed, 285 insertions, 30 deletions
diff --git a/include/proj/crs.hpp b/include/proj/crs.hpp index 30134cb4..61844f79 100644 --- a/include/proj/crs.hpp +++ b/include/proj/crs.hpp @@ -138,6 +138,8 @@ class PROJ_GCC_DLL CRS : public common::ObjectUsage, PROJ_INTERNAL CRSNNPtr normalizeForVisualization() const; + PROJ_INTERNAL CRSNNPtr allowNonConformantWKT1Export() const; + //! @endcond protected: @@ -892,6 +894,14 @@ class PROJ_GCC_DLL CompoundCRS final : public CRS, const std::vector<CRSNNPtr> &components); // throw InvalidCompoundCRSException + //! @cond Doxygen_Suppress + PROJ_INTERNAL static CRSNNPtr + createLax(const util::PropertyMap &properties, + const std::vector<CRSNNPtr> &components, + const io::DatabaseContextPtr + &dbContext); // throw InvalidCompoundCRSException + //! @endcond + protected: // relaxed: standard say SingleCRSNNPtr PROJ_INTERNAL explicit CompoundCRS(const std::vector<CRSNNPtr> &components); diff --git a/src/apps/projinfo.cpp b/src/apps/projinfo.cpp index ed732aaa..cd04433d 100644 --- a/src/apps/projinfo.cpp +++ b/src/apps/projinfo.cpp @@ -479,8 +479,8 @@ static void outputObject( std::cout << "WKT1:GDAL string:" << std::endl; } - auto formatter = - WKTFormatter::create(WKTFormatter::Convention::WKT1_GDAL); + auto formatter = WKTFormatter::create( + WKTFormatter::Convention::WKT1_GDAL, dbContext); if (outputOpt.singleLine) { formatter->setMultiLine(false); } diff --git a/src/iso19111/crs.cpp b/src/iso19111/crs.cpp index b0f5bcf9..77261211 100644 --- a/src/iso19111/crs.cpp +++ b/src/iso19111/crs.cpp @@ -93,6 +93,7 @@ struct CRS::Private { BoundCRSPtr canonicalBoundCRS_{}; std::string extensionProj4_{}; bool implicitCS_ = false; + bool allowNonConformantWKT1Export_ = false; void setImplicitCS(const util::PropertyMap &properties) { const auto pVal = properties.get("IMPLICIT_CS"); @@ -559,6 +560,18 @@ CRSNNPtr CRS::shallowClone() const { return _shallowClone(); } //! @cond Doxygen_Suppress +CRSNNPtr CRS::allowNonConformantWKT1Export() const { + auto crs = shallowClone(); + crs->d->allowNonConformantWKT1Export_ = true; + return crs; +} + +//! @endcond + +// --------------------------------------------------------------------------- + +//! @cond Doxygen_Suppress + CRSNNPtr CRS::alterName(const std::string &newName) const { auto crs = shallowClone(); auto newNameMod(newName); @@ -1334,30 +1347,49 @@ void GeodeticCRS::_exportToWKT(io::WKTFormatter *formatter) const { const bool isWKT2 = formatter->version() == io::WKTFormatter::Version::WKT2; const bool isGeographic = dynamic_cast<const GeographicCRS *>(this) != nullptr; - formatter->startNode(isWKT2 - ? ((formatter->use2019Keywords() && isGeographic) - ? io::WKTConstants::GEOGCRS - : io::WKTConstants::GEODCRS) - : isGeocentric() ? io::WKTConstants::GEOCCS - : io::WKTConstants::GEOGCS, - !identifiers().empty()); - auto l_name = nameStr(); + const auto &cs = coordinateSystem(); const auto &axisList = cs->axisList(); - const auto oldAxisOutputRule = formatter->outputAxis(); + auto l_name = nameStr(); + const auto &dbContext = formatter->databaseContext(); if (formatter->useESRIDialect()) { if (axisList.size() != 2) { io::FormattingException::Throw( "Only export of Geographic 2D CRS is supported in WKT1_ESRI"); } + } + + if (!isWKT2 && formatter->isStrict() && isGeographic && + axisList.size() != 2 && + oldAxisOutputRule != io::WKTFormatter::OutputAxisRule::NO) { + if (CRS::getPrivate()->allowNonConformantWKT1Export_) { + formatter->startNode(io::WKTConstants::COMPD_CS, false); + formatter->addQuotedString(l_name + " + " + l_name); + auto geogCRS = demoteTo2D(std::string(), dbContext); + geogCRS->_exportToWKT(formatter); + geogCRS->_exportToWKT(formatter); + formatter->endNode(); + return; + } + io::FormattingException::Throw( + "WKT1 does not support Geographic 3D CRS."); + } + formatter->startNode(isWKT2 + ? ((formatter->use2019Keywords() && isGeographic) + ? io::WKTConstants::GEOGCRS + : io::WKTConstants::GEODCRS) + : isGeocentric() ? io::WKTConstants::GEOCCS + : io::WKTConstants::GEOGCS, + !identifiers().empty()); + + if (formatter->useESRIDialect()) { if (l_name == "WGS 84") { l_name = "GCS_WGS_1984"; } else { bool aliasFound = false; - const auto &dbContext = formatter->databaseContext(); if (dbContext) { auto l_alias = dbContext->getAliasFromOfficialName( l_name, "geodetic_crs", "ESRI"); @@ -1373,11 +1405,6 @@ void GeodeticCRS::_exportToWKT(io::WKTFormatter *formatter) const { } } } - } else if (!isWKT2 && formatter->isStrict() && isGeographic && - axisList.size() != 2 && - oldAxisOutputRule != io::WKTFormatter::OutputAxisRule::NO) { - io::FormattingException::Throw( - "WKT1 does not support Geographic 3D CRS."); } if (!isWKT2 && !formatter->useESRIDialect() && isDeprecated()) { @@ -3163,6 +3190,36 @@ void ProjectedCRS::_exportToWKT(io::WKTFormatter *formatter) const { const auto &dbContext = formatter->databaseContext(); auto l_name = nameStr(); + const auto &l_coordinateSystem = d->coordinateSystem(); + const auto &axisList = l_coordinateSystem->axisList(); + if (axisList.size() == 3 && !(isWKT2 && formatter->use2019Keywords())) { + if (!formatter->useESRIDialect() && + CRS::getPrivate()->allowNonConformantWKT1Export_) { + formatter->startNode(io::WKTConstants::COMPD_CS, false); + formatter->addQuotedString(l_name + " + " + baseCRS()->nameStr()); + auto projCRS2D = demoteTo2D(std::string(), dbContext); + if (dbContext) { + const auto res = + projCRS2D->identify(io::AuthorityFactory::create( + NN_NO_CHECK(dbContext), "EPSG")); + if (res.size() == 1) { + const auto &front = res.front(); + if (front.second == 100) { + projCRS2D = front.first; + } + } + } + projCRS2D->_exportToWKT(formatter); + baseCRS() + ->demoteTo2D(std::string(), dbContext) + ->_exportToWKT(formatter); + formatter->endNode(); + return; + } + io::FormattingException::Throw( + "Projected 3D CRS can only be exported since WKT2:2019"); + } + std::string l_alias; if (formatter->useESRIDialect() && dbContext) { l_alias = dbContext->getAliasFromOfficialName(l_name, "projected_crs", @@ -3214,13 +3271,6 @@ void ProjectedCRS::_exportToWKT(io::WKTFormatter *formatter) const { } } - const auto &l_coordinateSystem = d->coordinateSystem(); - const auto &axisList = l_coordinateSystem->axisList(); - if (axisList.size() == 3 && !(isWKT2 && formatter->use2019Keywords())) { - io::FormattingException::Throw( - "Projected 3D CRS can only be exported since WKT2:2019"); - } - const auto exportAxis = [&l_coordinateSystem, &axisList, &formatter]() { const auto oldAxisOutputRule = formatter->outputAxis(); if (oldAxisOutputRule == @@ -4091,6 +4141,54 @@ CompoundCRSNNPtr CompoundCRS::create(const util::PropertyMap &properties, // --------------------------------------------------------------------------- //! @cond Doxygen_Suppress + +/** \brief Instantiate a CompoundCRS, a Geographic 3D CRS or a Projected CRS + * from a vector of CRS. + * + * Be a bit "lax", in allowing formulations like EPSG:4326+4326 or + * EPSG:32631+4326 to express Geographic 3D CRS / Projected3D CRS. + * + * @param properties See \ref general_properties. + * At minimum the name should be defined. + * @param components the component CRS of the CompoundCRS. + * @return new CRS. + * @throw InvalidCompoundCRSException + */ +CRSNNPtr CompoundCRS::createLax(const util::PropertyMap &properties, + const std::vector<CRSNNPtr> &components, + const io::DatabaseContextPtr &dbContext) { + + if (components.size() == 2) { + auto comp0 = components[0].get(); + auto comp1 = components[1].get(); + auto comp0Geog = dynamic_cast<const GeographicCRS *>(comp0); + auto comp0Proj = dynamic_cast<const ProjectedCRS *>(comp0); + auto comp1Geog = dynamic_cast<const GeographicCRS *>(comp1); + if ((comp0Geog != nullptr || comp0Proj != nullptr) && + comp1Geog != nullptr) { + const auto horizGeog = + (comp0Proj != nullptr) + ? comp0Proj->baseCRS().as_nullable().get() + : comp0Geog; + if (horizGeog->_isEquivalentTo( + comp1Geog->demoteTo2D(std::string(), dbContext).get())) { + return components[0] + ->promoteTo3D(std::string(), dbContext) + ->allowNonConformantWKT1Export(); + } + throw InvalidCompoundCRSException( + "The 'vertical' geographic CRS is not equivalent to the " + "geographic CRS of the horizontal part"); + } + } + + return create(properties, components); +} +//! @endcond + +// --------------------------------------------------------------------------- + +//! @cond Doxygen_Suppress void CompoundCRS::_exportToWKT(io::WKTFormatter *formatter) const { const bool isWKT2 = formatter->version() == io::WKTFormatter::Version::WKT2; formatter->startNode(isWKT2 ? io::WKTConstants::COMPOUNDCRS diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp index e516e4b5..59116ad6 100644 --- a/src/iso19111/io.cpp +++ b/src/iso19111/io.cpp @@ -1334,7 +1334,7 @@ struct WKTParser::Private { DerivedVerticalCRSNNPtr buildDerivedVerticalCRS(const WKTNodeNNPtr &node); - CompoundCRSNNPtr buildCompoundCRS(const WKTNodeNNPtr &node); + CRSNNPtr buildCompoundCRS(const WKTNodeNNPtr &node); BoundCRSNNPtr buildBoundCRS(const WKTNodeNNPtr &node); @@ -4157,8 +4157,7 @@ WKTParser::Private::buildDerivedVerticalCRS(const WKTNodeNNPtr &node) { // --------------------------------------------------------------------------- -CompoundCRSNNPtr -WKTParser::Private::buildCompoundCRS(const WKTNodeNNPtr &node) { +CRSNNPtr WKTParser::Private::buildCompoundCRS(const WKTNodeNNPtr &node) { std::vector<CRSNNPtr> components; for (const auto &child : node->GP()->children()) { auto crs = buildCRS(child); @@ -4166,7 +4165,13 @@ WKTParser::Private::buildCompoundCRS(const WKTNodeNNPtr &node) { components.push_back(NN_NO_CHECK(crs)); } } - return CompoundCRS::create(buildProperties(node), components); + + if (ci_equal(node->GP()->value(), WKTConstants::COMPD_CS)) { + return CompoundCRS::createLax(buildProperties(node), components, + dbContext_); + } else { + return CompoundCRS::create(buildProperties(node), components); + } } // --------------------------------------------------------------------------- @@ -5851,11 +5856,11 @@ static BaseObjectNNPtr createFromUserInput(const std::string &text, tokensCode[0], false)); auto crs2(factory->createCoordinateReferenceSystem( tokensCode[1], false)); - return CompoundCRS::create( + return CompoundCRS::createLax( util::PropertyMap().set( IdentifiedObject::NAME_KEY, crs1->nameStr() + " + " + crs2->nameStr()), - {crs1, crs2}); + {crs1, crs2}, dbContext); } throw; } diff --git a/test/unit/test_io.cpp b/test/unit/test_io.cpp index a227ff31..da26aa2d 100644 --- a/test/unit/test_io.cpp +++ b/test/unit/test_io.cpp @@ -2430,6 +2430,53 @@ TEST(wkt_parse, COMPD_CS) { // --------------------------------------------------------------------------- +TEST(wkt_parse, COMPD_CS_non_conformant_horizontal_plus_horizontal_as_in_LAS) { + auto obj = WKTParser().createFromWKT( + "COMPD_CS[\"horizontal + vertical\",\n" + " PROJCS[\"WGS 84 / UTM zone 31N\",\n" + " GEOGCS[\"WGS 84\",\n" + " DATUM[\"World Geodetic System 1984\",\n" + " SPHEROID[\"WGS 84\",6378137,298.257223563,\n" + " AUTHORITY[\"EPSG\",\"7030\"]],\n" + " AUTHORITY[\"EPSG\",\"6326\"]],\n" + " PRIMEM[\"Greenwich\",0,\n" + " AUTHORITY[\"EPSG\",\"8901\"]],\n" + " UNIT[\"degree\",0.0174532925199433,\n" + " AUTHORITY[\"EPSG\",\"9122\"]],\n" + " AXIS[\"Latitude\",NORTH],\n" + " AXIS[\"Longitude\",EAST],\n" + " AUTHORITY[\"EPSG\",\"4326\"]],\n" + " PROJECTION[\"Transverse_Mercator\"],\n" + " PARAMETER[\"latitude_of_origin\",0],\n" + " PARAMETER[\"central_meridian\",3],\n" + " PARAMETER[\"scale_factor\",0.9996],\n" + " PARAMETER[\"false_easting\",500000],\n" + " PARAMETER[\"false_northing\",0],\n" + " UNIT[\"metre\",1,\n" + " AUTHORITY[\"EPSG\",\"9001\"]],\n" + " AXIS[\"Easting\",EAST],\n" + " AXIS[\"Northing\",NORTH],\n" + " AUTHORITY[\"EPSG\",\"32631\"]],\n" + " GEOGCS[\"WGS 84\",\n" + " DATUM[\"World Geodetic System 1984\",\n" + " SPHEROID[\"WGS 84\",6378137,298.257223563,\n" + " AUTHORITY[\"EPSG\",\"7030\"]],\n" + " AUTHORITY[\"EPSG\",\"6326\"]],\n" + " PRIMEM[\"Greenwich\",0,\n" + " AUTHORITY[\"EPSG\",\"8901\"]],\n" + " UNIT[\"degree\",0.0174532925199433,\n" + " AUTHORITY[\"EPSG\",\"9122\"]],\n" + " AXIS[\"Latitude\",NORTH],\n" + " AXIS[\"Longitude\",EAST],\n" + " AUTHORITY[\"EPSG\",\"4326\"]]]"); + auto crs = nn_dynamic_pointer_cast<ProjectedCRS>(obj); + ASSERT_TRUE(crs != nullptr); + EXPECT_EQ(crs->nameStr(), "WGS 84 / UTM zone 31N"); + EXPECT_EQ(crs->coordinateSystem()->axisList().size(), 3U); +} + +// --------------------------------------------------------------------------- + TEST(wkt_parse, COORDINATEOPERATION) { std::string src_wkt; @@ -9490,6 +9537,101 @@ TEST(io, createFromUserInput) { EXPECT_EQ(crs->derivingConversion()->getEPSGCode(), 16031); } + // We accept non-conformant EPSG:4326+4326 + { + auto obj = createFromUserInput("EPSG:4326+4326", dbContext); + auto crs = nn_dynamic_pointer_cast<GeographicCRS>(obj); + ASSERT_TRUE(crs != nullptr); + EXPECT_EQ(crs->nameStr(), "WGS 84"); + EXPECT_EQ(crs->getEPSGCode(), 4979); + EXPECT_EQ(crs->coordinateSystem()->axisList().size(), 3U); + + const auto wkt = + "COMPD_CS[\"WGS 84 + WGS 84\",\n" + " GEOGCS[\"WGS 84\",\n" + " DATUM[\"WGS_1984\",\n" + " SPHEROID[\"WGS 84\",6378137,298.257223563,\n" + " AUTHORITY[\"EPSG\",\"7030\"]],\n" + " AUTHORITY[\"EPSG\",\"6326\"]],\n" + " PRIMEM[\"Greenwich\",0,\n" + " AUTHORITY[\"EPSG\",\"8901\"]],\n" + " UNIT[\"degree\",0.0174532925199433,\n" + " AUTHORITY[\"EPSG\",\"9122\"]],\n" + " AUTHORITY[\"EPSG\",\"4326\"]],\n" + " GEOGCS[\"WGS 84\",\n" + " DATUM[\"WGS_1984\",\n" + " SPHEROID[\"WGS 84\",6378137,298.257223563,\n" + " AUTHORITY[\"EPSG\",\"7030\"]],\n" + " AUTHORITY[\"EPSG\",\"6326\"]],\n" + " PRIMEM[\"Greenwich\",0,\n" + " AUTHORITY[\"EPSG\",\"8901\"]],\n" + " UNIT[\"degree\",0.0174532925199433,\n" + " AUTHORITY[\"EPSG\",\"9122\"]],\n" + " AUTHORITY[\"EPSG\",\"4326\"]]]"; + + EXPECT_EQ( + crs->exportToWKT(WKTFormatter::create( + WKTFormatter::Convention::WKT1_GDAL, dbContext) + .get()), + wkt); + } + + // Non consistent + EXPECT_THROW(createFromUserInput("EPSG:4326+4258", dbContext), + InvalidCompoundCRSException); + + // We accept non-conformant EPSG:32631+4326 + { + auto obj = createFromUserInput("EPSG:32631+4326", dbContext); + auto crs = nn_dynamic_pointer_cast<ProjectedCRS>(obj); + ASSERT_TRUE(crs != nullptr); + EXPECT_EQ(crs->nameStr(), "WGS 84 / UTM zone 31N"); + EXPECT_EQ(crs->baseCRS()->getEPSGCode(), 4979); + EXPECT_EQ(crs->baseCRS()->coordinateSystem()->axisList().size(), 3U); + EXPECT_EQ(crs->coordinateSystem()->axisList().size(), 3U); + + const auto wkt = + "COMPD_CS[\"WGS 84 / UTM zone 31N + WGS 84\",\n" + " PROJCS[\"WGS 84 / UTM zone 31N\",\n" + " GEOGCS[\"WGS 84\",\n" + " DATUM[\"WGS_1984\",\n" + " SPHEROID[\"WGS 84\",6378137,298.257223563,\n" + " AUTHORITY[\"EPSG\",\"7030\"]],\n" + " AUTHORITY[\"EPSG\",\"6326\"]],\n" + " PRIMEM[\"Greenwich\",0,\n" + " AUTHORITY[\"EPSG\",\"8901\"]],\n" + " UNIT[\"degree\",0.0174532925199433,\n" + " AUTHORITY[\"EPSG\",\"9122\"]],\n" + " AUTHORITY[\"EPSG\",\"4326\"]],\n" + " PROJECTION[\"Transverse_Mercator\"],\n" + " PARAMETER[\"latitude_of_origin\",0],\n" + " PARAMETER[\"central_meridian\",3],\n" + " PARAMETER[\"scale_factor\",0.9996],\n" + " PARAMETER[\"false_easting\",500000],\n" + " PARAMETER[\"false_northing\",0],\n" + " UNIT[\"metre\",1,\n" + " AUTHORITY[\"EPSG\",\"9001\"]],\n" + " AXIS[\"Easting\",EAST],\n" + " AXIS[\"Northing\",NORTH],\n" + " AUTHORITY[\"EPSG\",\"32631\"]],\n" + " GEOGCS[\"WGS 84\",\n" + " DATUM[\"WGS_1984\",\n" + " SPHEROID[\"WGS 84\",6378137,298.257223563,\n" + " AUTHORITY[\"EPSG\",\"7030\"]],\n" + " AUTHORITY[\"EPSG\",\"6326\"]],\n" + " PRIMEM[\"Greenwich\",0,\n" + " AUTHORITY[\"EPSG\",\"8901\"]],\n" + " UNIT[\"degree\",0.0174532925199433,\n" + " AUTHORITY[\"EPSG\",\"9122\"]],\n" + " AUTHORITY[\"EPSG\",\"4326\"]]]"; + + EXPECT_EQ( + crs->exportToWKT(WKTFormatter::create( + WKTFormatter::Convention::WKT1_GDAL, dbContext) + .get()), + wkt); + } + EXPECT_THROW(createFromUserInput( "urn:ogc:def:crs,crs:EPSG::4979," "cs:PROJ::ENh," |
