From d6f6862352d0caeb827c0c3ec98ccf9aa2f9a1d2 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Tue, 24 Nov 2020 16:41:06 +0100 Subject: [PATCH] createOperation(): add a ballpark vertical transformation when dealing with GEOIDMODEL[] --- src/iso19111/coordinateoperation.cpp | 99 ++++++++++++++++++++++++---- test/cli/testprojinfo_out.dist | 69 ++++++++++++++++++- test/unit/test_c_api.cpp | 49 ++++++++++---- test/unit/test_operation.cpp | 12 ++-- 4 files changed, 194 insertions(+), 35 deletions(-) diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp index 1d3b7a90d9..0f1c680beb 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); diff --git a/test/cli/testprojinfo_out.dist b/test/cli/testprojinfo_out.dist index 174f11d331..8e8ef294d4 100644 --- a/test/cli/testprojinfo_out.dist +++ b/test/cli/testprojinfo_out.dist @@ -990,7 +990,7 @@ PROJ.4 string: +proj=utm +zone=32 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs +type=crs Testing RH2000 height to SWEREF99: projinfo -s EPSG:5613 -t EPSG:4977 -Candidate operations found: 1 +Candidate operations found: 2 ------------------------------------- Operation No. 1: @@ -1042,8 +1042,55 @@ COORDINATEOPERATION["Inverse of SWEREF99 to RH2000 height", BBOX[55.28,10.93,69.07,24.17]], ID["INVERSE(DERIVED_FROM(PROJ))","EPSG_4977_TO_EPSG_5613"]] +------------------------------------- +Operation No. 2: + +unknown id, Transformation from RH2000 height to SWEREF99 (ballpark vertical transformation, without ellipsoid height to vertical height correction), unknown accuracy, World, has ballpark transformation + +PROJ string: ++proj=noop + +WKT2:2019 string: +COORDINATEOPERATION["Transformation from RH2000 height to SWEREF99 (ballpark vertical transformation, without ellipsoid height to vertical height correction)", + SOURCECRS[ + VERTCRS["RH2000 height", + DYNAMIC[ + FRAMEEPOCH[2000]], + VDATUM["Rikets hojdsystem 2000"], + CS[vertical,1], + AXIS["gravity-related height (H)",up, + LENGTHUNIT["metre",1]], + ID["EPSG",5613]]], + TARGETCRS[ + GEOGCRS["SWEREF99", + DATUM["SWEREF99", + ELLIPSOID["GRS 1980",6378137,298.257222101, + LENGTHUNIT["metre",1]]], + PRIMEM["Greenwich",0, + ANGLEUNIT["degree",0.0174532925199433]], + CS[ellipsoidal,3], + AXIS["geodetic latitude (Lat)",north, + ORDER[1], + ANGLEUNIT["degree",0.0174532925199433]], + AXIS["geodetic longitude (Lon)",east, + ORDER[2], + ANGLEUNIT["degree",0.0174532925199433]], + AXIS["ellipsoidal height (h)",up, + ORDER[3], + LENGTHUNIT["metre",1]], + ID["EPSG",4977]]], + METHOD["Change of Vertical Unit", + ID["EPSG",1069]], + PARAMETER["Unit conversion scalar",1, + SCALEUNIT["unity",1], + ID["EPSG",1051]], + USAGE[ + SCOPE["unknown"], + AREA["World"], + BBOX[-90,-180,90,180]]] + Testing NAD83(2011) + NAVD88 height -> NAD83(2011) : projinfo -s EPSG:6349 -t EPSG:6319 --spatial-test intersects -o PROJ -Candidate operations found: 2 +Candidate operations found: 3 ------------------------------------- Operation No. 1: @@ -1070,8 +1117,16 @@ PROJ string: +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1 +------------------------------------- +Operation No. 3: + +unknown id, Transformation from NAVD88 height to NAD83(2011) (ballpark vertical transformation, without ellipsoid height to vertical height correction), unknown accuracy, World, has ballpark transformation + +PROJ string: ++proj=noop + Testing NGF IGN69 height to RGF93: projinfo -s EPSG:5720 -t EPSG:4965 -o PROJ -Candidate operations found: 1 +Candidate operations found: 2 ------------------------------------- Operation No. 1: @@ -1085,6 +1140,14 @@ PROJ string: +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1 +------------------------------------- +Operation No. 2: + +unknown id, Transformation from NGF-IGN69 height to RGF93 (ballpark vertical transformation, without ellipsoid height to vertical height correction), unknown accuracy, World, has ballpark transformation + +PROJ string: ++proj=noop + Testing EPSG:32631 --3d PROJ.4 string: +proj=utm +zone=31 +datum=WGS84 +units=m +no_defs +type=crs diff --git a/test/unit/test_c_api.cpp b/test/unit/test_c_api.cpp index a7c1eb535e..9885d221ff 100644 --- a/test/unit/test_c_api.cpp +++ b/test/unit/test_c_api.cpp @@ -4650,10 +4650,21 @@ TEST_F(CApi, proj_create_vertical_crs_ex) { ObjectKeeper keeper_geog_crs(geog_crs); ASSERT_NE(geog_crs, nullptr); - auto P = proj_create_crs_to_crs_from_pj(m_ctxt, compound, geog_crs, nullptr, - nullptr); - ObjectKeeper keeper_P(P); - ASSERT_NE(P, nullptr); + PJ_OPERATION_FACTORY_CONTEXT *ctxt = + proj_create_operation_factory_context(m_ctxt, nullptr); + ASSERT_NE(ctxt, nullptr); + ContextKeeper keeper_ctxt(ctxt); + proj_operation_factory_context_set_grid_availability_use( + m_ctxt, ctxt, PROJ_GRID_AVAILABILITY_IGNORED); + proj_operation_factory_context_set_spatial_criterion( + m_ctxt, ctxt, PROJ_SPATIAL_CRITERION_PARTIAL_INTERSECTION); + PJ_OBJ_LIST *operations = + proj_create_operations(m_ctxt, compound, geog_crs, ctxt); + ASSERT_NE(operations, nullptr); + ObjListKeeper keeper_operations(operations); + EXPECT_GE(proj_list_get_count(operations), 1); + auto P = proj_list_get(m_ctxt, operations, 0); + ObjectKeeper keeper_transform(P); auto name = proj_get_name(P); ASSERT_TRUE(name != nullptr); @@ -4706,10 +4717,21 @@ TEST_F(CApi, proj_create_vertical_crs_ex_with_geog_crs) { ObjectKeeper keeper_geog_crs(geog_crs); ASSERT_NE(geog_crs, nullptr); - auto P = proj_create_crs_to_crs_from_pj(m_ctxt, compound, geog_crs, nullptr, - nullptr); - ObjectKeeper keeper_P(P); - ASSERT_NE(P, nullptr); + PJ_OPERATION_FACTORY_CONTEXT *ctxt = + proj_create_operation_factory_context(m_ctxt, nullptr); + ASSERT_NE(ctxt, nullptr); + ContextKeeper keeper_ctxt(ctxt); + proj_operation_factory_context_set_grid_availability_use( + m_ctxt, ctxt, PROJ_GRID_AVAILABILITY_IGNORED); + proj_operation_factory_context_set_spatial_criterion( + m_ctxt, ctxt, PROJ_SPATIAL_CRITERION_PARTIAL_INTERSECTION); + PJ_OBJ_LIST *operations = + proj_create_operations(m_ctxt, compound, geog_crs, ctxt); + ASSERT_NE(operations, nullptr); + ObjListKeeper keeper_operations(operations); + EXPECT_GE(proj_list_get_count(operations), 1); + auto P = proj_list_get(m_ctxt, operations, 0); + ObjectKeeper keeper_transform(P); auto name = proj_get_name(P); ASSERT_TRUE(name != nullptr); @@ -4739,10 +4761,13 @@ TEST_F(CApi, proj_create_vertical_crs_ex_with_geog_crs) { ObjectKeeper keeper_compound_from_projjson(compound_from_projjson); ASSERT_NE(compound_from_projjson, nullptr); - auto P2 = proj_create_crs_to_crs_from_pj(m_ctxt, compound_from_projjson, - geog_crs, nullptr, nullptr); - ObjectKeeper keeper_P2(P2); - ASSERT_NE(P2, nullptr); + PJ_OBJ_LIST *operations2 = + proj_create_operations(m_ctxt, compound_from_projjson, geog_crs, ctxt); + ASSERT_NE(operations2, nullptr); + ObjListKeeper keeper_operations2(operations2); + EXPECT_GE(proj_list_get_count(operations2), 1); + auto P2 = proj_list_get(m_ctxt, operations2, 0); + ObjectKeeper keeper_transform2(P2); auto name_bis = proj_get_name(P2); ASSERT_TRUE(name_bis != nullptr); diff --git a/test/unit/test_operation.cpp b/test/unit/test_operation.cpp index b6e555b1a4..9bf883d990 100644 --- a/test/unit/test_operation.cpp +++ b/test/unit/test_operation.cpp @@ -4867,7 +4867,7 @@ TEST(operation, vertCRS_to_geogCRS_context) { "3855"), // EGM2008 height authFactory->createCoordinateReferenceSystem("4979"), // WGS 84 ctxt); - ASSERT_EQ(list.size(), 2U); + ASSERT_EQ(list.size(), 3U); EXPECT_EQ( list[1]->exportToPROJString( PROJStringFormatter::create( @@ -4889,7 +4889,7 @@ TEST(operation, vertCRS_to_geogCRS_context) { "3855"), // EGM2008 height authFactory->createCoordinateReferenceSystem("4979"), // WGS 84 ctxt); - ASSERT_EQ(list.size(), 2U); + ASSERT_EQ(list.size(), 3U); EXPECT_EQ( list[0]->exportToPROJString(PROJStringFormatter::create().get()), "+proj=pipeline " @@ -4938,7 +4938,7 @@ TEST(operation, vertCRS_to_geogCRS_context) { authFactory->createCoordinateReferenceSystem("7839"), // NZGD2000 authFactory->createCoordinateReferenceSystem("4959"), ctxt); - ASSERT_EQ(list.size(), 1U); + ASSERT_EQ(list.size(), 2U); EXPECT_EQ( list[0]->exportToPROJString(PROJStringFormatter::create().get()), "+proj=pipeline " @@ -4964,7 +4964,7 @@ TEST(operation, vertCRS_to_geogCRS_context) { authFactory->createCoordinateReferenceSystem("8357"), // ETRS89 authFactory->createCoordinateReferenceSystem("4937"), ctxt); - ASSERT_EQ(list.size(), 1U); + ASSERT_EQ(list.size(), 2U); EXPECT_EQ( list[0]->exportToPROJString(PROJStringFormatter::create().get()), "+proj=pipeline " @@ -8912,7 +8912,7 @@ TEST(operation, compoundCRS_to_geogCRS_2D_promote_to_3D_context) { ctxt); // The checked value is not that important, but in case this changes, // likely due to a EPSG upgrade, worth checking - EXPECT_EQ(listCompoundToGeog2D.size(), 141U); + EXPECT_EQ(listCompoundToGeog2D.size(), 142U); auto listGeog2DToCompound = CoordinateOperationFactory::create()->createOperations(dst, nnSrc, @@ -11169,7 +11169,7 @@ TEST(operation, normalizeForVisualization) { src, authFactory->createCoordinateReferenceSystem("4979"), // WGS 84 3D ctxt); - ASSERT_EQ(list.size(), 2U); + ASSERT_EQ(list.size(), 3U); auto op = list[1]; auto opNormalized = op->normalizeForVisualization(); EXPECT_FALSE(opNormalized->_isEquivalentTo(op.get()));