From 2564b59ed24674aa8c8c22bff6cba840ab0e527d Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Tue, 14 Jan 2025 16:13:41 +0100 Subject: [PATCH] createOperations(): use more appropriate operation when using a 'PROJ {grid_name}' geoid model, based on matching the vertical datum Helps when transforming from 'NAD83(2011) / Puerto Rico and Virgin Is. + "VIVD09 height + PROJ us_noaa_g2012bp0.tif' to NAD83(2011). Without that fix, PROJ would use the EPSG record for the NAD83(2011) <--> PRVD02 height (Puerto Rico) to get the extent of the transformation, instead of the one for NAD83(2011) <--> VIVD09 height (Virgin Islands), causing the resulting global transformation to not be usable in practice. --- .../operation/coordinateoperationfactory.cpp | 38 +++++-- test/unit/test_operationfactory.cpp | 105 ++++++++++++++++++ 2 files changed, 134 insertions(+), 9 deletions(-) diff --git a/src/iso19111/operation/coordinateoperationfactory.cpp b/src/iso19111/operation/coordinateoperationfactory.cpp index 2059aa59a0..fa9a9be56c 100644 --- a/src/iso19111/operation/coordinateoperationfactory.cpp +++ b/src/iso19111/operation/coordinateoperationfactory.cpp @@ -4253,14 +4253,39 @@ CoordinateOperationFactory::Private::createOperationsGeogToVertFromGeoid( std::vector accuracies; const auto &modelAccuracies = model->coordinateOperationAccuracies(); - std::vector transformationsForGrid; + std::vector transformationsForGrid = + io::DatabaseContext::getTransformationsForGridName( + authFactory->databaseContext(), projFilename); + + // Only select transformations whose datum of the target vertical + // CRS match the one of the target vertical CRS of interest (when + // there's such match) Helps for example if specifying GEOID + // g2012bp0 whose has a record for Puerto Rico and another one for + // Virgin Islands. + { + std::vector + transformationsForGridMatchingDatum; + for (const auto &op : transformationsForGrid) { + const auto opTargetCRS = + dynamic_cast( + op->targetCRS().get()); + if (opTargetCRS && + opTargetCRS->datum()->_isEquivalentTo( + vertDst->datum().get(), + util::IComparable::Criterion::EQUIVALENT)) { + transformationsForGridMatchingDatum.push_back(op); + } + } + if (!transformationsForGridMatchingDatum.empty()) { + transformationsForGrid = + std::move(transformationsForGridMatchingDatum); + } + } + double accuracy = -1; size_t idx = static_cast(-1); if (modelAccuracies.empty()) { if (authFactory) { - transformationsForGrid = - io::DatabaseContext::getTransformationsForGridName( - authFactory->databaseContext(), projFilename); for (size_t i = 0; i < transformationsForGrid.size(); ++i) { const auto &transf = transformationsForGrid[i]; const double transfAcc = getAccuracy(transf); @@ -4284,11 +4309,6 @@ CoordinateOperationFactory::Private::createOperationsGeogToVertFromGeoid( // Otherwise fallback to the extent of a transformation using // the grid. if (extent == nullptr && authFactory != nullptr) { - if (transformationsForGrid.empty()) { - transformationsForGrid = - io::DatabaseContext::getTransformationsForGridName( - authFactory->databaseContext(), projFilename); - } if (idx != static_cast(-1)) { const auto &transf = transformationsForGrid[idx]; extent = getExtent(transf, true, dummy); diff --git a/test/unit/test_operationfactory.cpp b/test/unit/test_operationfactory.cpp index 755c749017..7d4da1a879 100644 --- a/test/unit/test_operationfactory.cpp +++ b/test/unit/test_operationfactory.cpp @@ -8284,6 +8284,111 @@ TEST(operation, compoundCRS_of_vertCRS_with_geoid_model_by_id_to_geogCRS) { // --------------------------------------------------------------------------- +TEST( + operation, + compoundCRS_of_vertCRS_with_geoid_model_by_name_and_several_records_to_geogCRS) { + auto authFactory = + AuthorityFactory::create(DatabaseContext::create(), "EPSG"); + auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0); + ctxt->setSpatialCriterion( + CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION); + ctxt->setGridAvailabilityUse( + CoordinateOperationContext::GridAvailabilityUse:: + IGNORE_GRID_AVAILABILITY); + auto wkt = + "COMPOUNDCRS[\"Compound CRS NAD83(2011) / Puerto Rico and Virgin Is. + " + "VIVD09 height + PROJ us_noaa_g2012bp0.tif\",\n" + " PROJCRS[\"NAD83(2011) / Puerto Rico and Virgin Is.\",\n" + " BASEGEOGCRS[\"NAD83(2011)\",\n" + " DATUM[\"NAD83 (National Spatial Reference System " + "2011)\",\n" + " ELLIPSOID[\"GRS 1980\",6378137,298.257222101,\n" + " LENGTHUNIT[\"metre\",1]]],\n" + " PRIMEM[\"Greenwich\",0,\n" + " ANGLEUNIT[\"degree\",0.0174532925199433]],\n" + " ID[\"EPSG\",6318]],\n" + " CONVERSION[\"SPCS83Puerto Rico & Virgin Islands zone " + "(meter)\",\n" + " METHOD[\"Lambert Conic Conformal (2SP)\",\n" + " ID[\"EPSG\",9802]],\n" + " PARAMETER[\"Latitude of false origin\",17.8333333333333,\n" + " ANGLEUNIT[\"degree\",0.0174532925199433],\n" + " ID[\"EPSG\",8821]],\n" + " PARAMETER[\"Longitude of false " + "origin\",-66.4333333333333,\n" + " ANGLEUNIT[\"degree\",0.0174532925199433],\n" + " ID[\"EPSG\",8822]],\n" + " PARAMETER[\"Latitude of 1st standard " + "parallel\",18.4333333333333,\n" + " ANGLEUNIT[\"degree\",0.0174532925199433],\n" + " ID[\"EPSG\",8823]],\n" + " PARAMETER[\"Latitude of 2nd standard " + "parallel\",18.0333333333333,\n" + " ANGLEUNIT[\"degree\",0.0174532925199433],\n" + " ID[\"EPSG\",8824]],\n" + " PARAMETER[\"Easting at false origin\",200000,\n" + " LENGTHUNIT[\"metre\",1],\n" + " ID[\"EPSG\",8826]],\n" + " PARAMETER[\"Northing at false origin\",200000,\n" + " LENGTHUNIT[\"metre\",1],\n" + " ID[\"EPSG\",8827]],\n" + " ID[\"EPSG\",15230]],\n" + " CS[Cartesian,2],\n" + " AXIS[\"easting (X)\",east,\n" + " ORDER[1],\n" + " LENGTHUNIT[\"metre\",1,\n" + " ID[\"EPSG\",9001]]],\n" + " AXIS[\"northing (Y)\",north,\n" + " ORDER[2],\n" + " LENGTHUNIT[\"metre\",1,\n" + " ID[\"EPSG\",9001]]]],\n" + " VERTCRS[\"VIVD09 height + PROJ us_noaa_g2012bp0.tif\",\n" + " VDATUM[\"Virgin Islands Vertical Datum of 2009\",\n" + " ID[\"EPSG\",1124]],\n" + " CS[vertical,1],\n" + " AXIS[\"gravity-related height (H)\",up,\n" + " LENGTHUNIT[\"metre\",1,\n" + " ID[\"EPSG\",9001]]],\n" + " GEOIDMODEL[\"PROJ us_noaa_g2012bp0.tif\"]]]"; + auto srcObj = + createFromUserInput(wkt, authFactory->databaseContext(), false); + auto src = nn_dynamic_pointer_cast(srcObj); + ASSERT_TRUE(src != nullptr); + auto dst = + authFactory->createCoordinateReferenceSystem("6319")->promoteTo3D( + std::string(), authFactory->databaseContext()); // NAD83(2011) 3d + + auto list = CoordinateOperationFactory::create()->createOperations( + NN_NO_CHECK(src), dst, ctxt); + ASSERT_TRUE(!list.empty()); + EXPECT_STREQ(list[0]->nameStr().c_str(), + "Inverse of SPCS83Puerto Rico & Virgin Islands zone (meter) + " + "Transformation from VIVD09 height + " + "PROJ us_noaa_g2012bp0.tif to NAD83(2011)"); + auto op_proj = + list[0]->exportToPROJString(PROJStringFormatter::create().get()); + EXPECT_STREQ(op_proj.c_str(), + "+proj=pipeline " + "+step +inv +proj=lcc +lat_0=17.8333333333333 " + "+lon_0=-66.4333333333333 +lat_1=18.4333333333333 " + "+lat_2=18.0333333333333 +x_0=200000 +y_0=200000 +ellps=GRS80 " + "+step +proj=vgridshift +grids=us_noaa_g2012bp0.tif " + "+multiplier=1 " + "+step +proj=unitconvert +xy_in=rad +xy_out=deg " + "+step +proj=axisswap +order=2,1"); + ASSERT_EQ(list[0]->domains().size(), 1U); + auto domain = list[0]->domains()[0]; + ASSERT_TRUE(domain->domainOfValidity() != nullptr); + EXPECT_TRUE(domain->domainOfValidity()->description().has_value()); + // This is the thing we actually want to check, that the area of use + // is the one of Virgin Islands, and not Puerto Rico + EXPECT_STREQ( + domain->domainOfValidity()->description()->c_str(), + "US Virgin Islands - onshore - St Croix, St John, and St Thomas."); +} + +// --------------------------------------------------------------------------- + TEST(operation, compoundCRS_of_bound_horizCRS_and_bound_vertCRS_to_geogCRS) { auto authFactory = AuthorityFactory::create(DatabaseContext::create(), "EPSG");