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/io.cpp | 4 ++++ 1 file changed, 4 insertions(+) (limited to 'src/iso19111/io.cpp') 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 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/io.cpp | 5 +++++ 1 file changed, 5 insertions(+) (limited to 'src/iso19111/io.cpp') 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 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/io.cpp') 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/io.cpp | 1 - 1 file changed, 1 deletion(-) (limited to 'src/iso19111/io.cpp') 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 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/io.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) (limited to 'src/iso19111/io.cpp') 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 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/io.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) (limited to 'src/iso19111/io.cpp') 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/io.cpp') 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/io.cpp') 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 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/io.cpp | 25 +++++++++++++++++++++++++ 1 file changed, 25 insertions(+) (limited to 'src/iso19111/io.cpp') 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/io.cpp | 76 ++++++++++++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 72 insertions(+), 4 deletions(-) (limited to 'src/iso19111/io.cpp') 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 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/io.cpp') 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