From 31fd3de9c2b2f823c01b3c2dbadddf4a7101fa16 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Sat, 2 Nov 2019 15:43:42 +0100 Subject: WKT and PROJJSON: add import/export of geoid model of VertCRS --- src/iso19111/io.cpp | 35 ++++++++++++++++++++++++++++++++--- 1 file changed, 32 insertions(+), 3 deletions(-) (limited to 'src/iso19111/io.cpp') diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp index fc66b6c9..116a8b78 100644 --- a/src/iso19111/io.cpp +++ b/src/iso19111/io.cpp @@ -89,7 +89,7 @@ static const std::string emptyString{}; // If changing that value, change it in data/projjson.schema.json as well #define PROJJSON_CURRENT_VERSION \ - "https://proj.org/schemas/v0.1/projjson.schema.json" + "https://proj.org/schemas/v0.2/projjson.schema.json" //! @endcond #if 0 @@ -3820,8 +3820,22 @@ CRSNNPtr WKTParser::Private::buildVerticalCRS(const WKTNodeNNPtr &node) { ThrowNotExpectedCSType("vertical"); } + auto &props = buildProperties(node); + auto &geoidModelNode = nodeP->lookForChild(WKTConstants::GEOIDMODEL); + if (!isNull(geoidModelNode)) { + auto &propsModel = buildProperties(geoidModelNode); + const auto dummyCRS = VerticalCRS::create( + PropertyMap(), datum, datumEnsemble, NN_NO_CHECK(verticalCS)); + const auto model(Transformation::create( + propsModel, dummyCRS, dummyCRS, nullptr, + OperationMethod::create(PropertyMap(), + std::vector()), + {}, {})); + props.set("GEOID_MODEL", model); + } + auto crs = nn_static_pointer_cast(VerticalCRS::create( - buildProperties(node), datum, datumEnsemble, NN_NO_CHECK(verticalCS))); + props, datum, datumEnsemble, NN_NO_CHECK(verticalCS))); if (!isNull(datumNode)) { auto &extensionNode = datumNode->lookForChild(WKTConstants::EXTENSION); @@ -4987,7 +5001,22 @@ VerticalCRSNNPtr JSONParser::buildVerticalCRS(const json &j) { if (!verticalCS) { throw ParsingException("expected a vertical CS"); } - return VerticalCRS::create(buildProperties(j), datum, datumEnsemble, + + auto props = buildProperties(j); + if (j.contains("geoid_model")) { + auto geoidModelJ = getObject(j, "geoid_model"); + auto propsModel = buildProperties(geoidModelJ); + const auto dummyCRS = VerticalCRS::create( + PropertyMap(), datum, datumEnsemble, NN_NO_CHECK(verticalCS)); + const auto model(Transformation::create( + propsModel, dummyCRS, dummyCRS, nullptr, + OperationMethod::create(PropertyMap(), + std::vector()), + {}, {})); + props.set("GEOID_MODEL", model); + } + + return VerticalCRS::create(props, datum, datumEnsemble, NN_NO_CHECK(verticalCS)); } -- cgit v1.2.3 From ff69c7e4535e869d0d9c6831c50d2f0fd2d288ca Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Sun, 3 Nov 2019 10:46:36 +0100 Subject: Import from WKT: add tweaks for Lidar WKT1 VERT_CS that embeds geoid model in CRS name --- src/iso19111/io.cpp | 58 +++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 58 insertions(+) (limited to 'src/iso19111/io.cpp') diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp index 116a8b78..d4e6132f 100644 --- a/src/iso19111/io.cpp +++ b/src/iso19111/io.cpp @@ -3821,6 +3821,64 @@ CRSNNPtr WKTParser::Private::buildVerticalCRS(const WKTNodeNNPtr &node) { } auto &props = buildProperties(node); + + // Deal with Lidar WKT1 VertCRS that embeds geoid model in CRS name, + // following conventions from + // https://pubs.usgs.gov/tm/11b4/pdf/tm11-B4.pdf + // page 9 + if (ci_equal(nodeValue, WKTConstants::VERT_CS)) { + std::string name; + if (props.getStringValue(IdentifiedObject::NAME_KEY, name)) { + std::string geoidName; + for (const char *prefix : + {"NAVD88 - ", "NAVD88 via ", "NAVD88 height - ", + "NAVD88 height (ftUS) - "}) { + if (starts_with(name, prefix)) { + geoidName = name.substr(strlen(prefix)); + auto pos = geoidName.find_first_of(" ("); + if (pos != std::string::npos) { + geoidName.resize(pos); + } + break; + } + } + if (!geoidName.empty()) { + const auto &axis = verticalCS->axisList()[0]; + const auto &dir = axis->direction(); + if (dir == cs::AxisDirection::UP) { + if (axis->unit() == common::UnitOfMeasure::METRE) { + props.set(IdentifiedObject::NAME_KEY, "NAVD88 height"); + props.set(Identifier::CODE_KEY, 5703); + props.set(Identifier::CODESPACE_KEY, Identifier::EPSG); + } else if (axis->unit().name() == "US survey foot") { + props.set(IdentifiedObject::NAME_KEY, + "NAVD88 height (ftUS)"); + props.set(Identifier::CODE_KEY, 6360); + props.set(Identifier::CODESPACE_KEY, Identifier::EPSG); + } + } + PropertyMap propsModel; + propsModel.set(IdentifiedObject::NAME_KEY, toupper(geoidName)); + PropertyMap propsDatum; + propsDatum.set(IdentifiedObject::NAME_KEY, + "North American Vertical Datum 1988"); + propsDatum.set(Identifier::CODE_KEY, 5103); + propsDatum.set(Identifier::CODESPACE_KEY, Identifier::EPSG); + datum = + VerticalReferenceFrame::create(propsDatum).as_nullable(); + const auto dummyCRS = + VerticalCRS::create(PropertyMap(), datum, datumEnsemble, + NN_NO_CHECK(verticalCS)); + const auto model(Transformation::create( + propsModel, dummyCRS, dummyCRS, nullptr, + OperationMethod::create( + PropertyMap(), std::vector()), + {}, {})); + props.set("GEOID_MODEL", model); + } + } + } + auto &geoidModelNode = nodeP->lookForChild(WKTConstants::GEOIDMODEL); if (!isNull(geoidModelNode)) { auto &propsModel = buildProperties(geoidModelNode); -- cgit v1.2.3 From cff22050faa67017eb4c5b44cb383b265c07d156 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Sun, 3 Nov 2019 11:09:14 +0100 Subject: Import from WKT of VERT_CS: remove outdated EXTENSION.PROJ4_GRIDS for NAVD88 --- src/iso19111/io.cpp | 36 +++++++++++++++++++++++------------- 1 file changed, 23 insertions(+), 13 deletions(-) (limited to 'src/iso19111/io.cpp') diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp index d4e6132f..aab1db4f 100644 --- a/src/iso19111/io.cpp +++ b/src/iso19111/io.cpp @@ -3900,20 +3900,30 @@ CRSNNPtr WKTParser::Private::buildVerticalCRS(const WKTNodeNNPtr &node) { const auto &extensionChildren = extensionNode->GP()->children(); if (extensionChildren.size() == 2) { if (ci_equal(stripQuotes(extensionChildren[0]), "PROJ4_GRIDS")) { - std::string transformationName(crs->nameStr()); - if (!ends_with(transformationName, " height")) { - transformationName += " height"; + const auto gridName(stripQuotes(extensionChildren[1])); + // This is the expansion of EPSG:5703 by old GDAL versions. + // See + // https://trac.osgeo.org/metacrs/changeset?reponame=&new=2281%40geotiff%2Ftrunk%2Flibgeotiff%2Fcsv%2Fvertcs.override.csv&old=1893%40geotiff%2Ftrunk%2Flibgeotiff%2Fcsv%2Fvertcs.override.csv + // It is unlikely that the user really explictly wants this. + if (gridName != "g2003conus.gtx,g2003alaska.gtx," + "g2003h01.gtx,g2003p01.gtx" && + gridName != "g2012a_conus.gtx,g2012a_alaska.gtx," + "g2012a_guam.gtx,g2012a_hawaii.gtx," + "g2012a_puertorico.gtx,g2012a_samoa.gtx") { + std::string transformationName(crs->nameStr()); + if (!ends_with(transformationName, " height")) { + transformationName += " height"; + } + transformationName += " to WGS84 ellipsoidal height"; + auto transformation = Transformation:: + createGravityRelatedHeightToGeographic3D( + PropertyMap().set(IdentifiedObject::NAME_KEY, + transformationName), + crs, GeographicCRS::EPSG_4979, nullptr, gridName, + std::vector()); + return nn_static_pointer_cast(BoundCRS::create( + crs, GeographicCRS::EPSG_4979, transformation)); } - transformationName += " to WGS84 ellipsoidal height"; - auto transformation = - Transformation::createGravityRelatedHeightToGeographic3D( - PropertyMap().set(IdentifiedObject::NAME_KEY, - transformationName), - crs, GeographicCRS::EPSG_4979, nullptr, - stripQuotes(extensionChildren[1]), - std::vector()); - return nn_static_pointer_cast(BoundCRS::create( - crs, GeographicCRS::EPSG_4979, transformation)); } } } -- cgit v1.2.3