Skip to content

Commit

Permalink
Merge pull request #4379 from rouault/proj_geoidgrid_extent
Browse files Browse the repository at this point in the history 
createOperations(): use more appropriate operation when using a 'PROJ {grid_name}' geoid model, based on matching the vertical datum
  • Loading branch information
@rouault
rouault authored Jan 15, 2025
2 parents b047dd2 + 2564b59 commit 98dcc57
Show file tree
Hide file tree
Showing 2 changed files with 134 additions and 9 deletions.
38 changes: 29 additions & 9 deletions src/iso19111/operation/coordinateoperationfactory.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4253,14 +4253,39 @@ CoordinateOperationFactory::Private::createOperationsGeogToVertFromGeoid(
std::vector<metadata::PositionalAccuracyNNPtr> accuracies;
const auto &modelAccuracies =
model->coordinateOperationAccuracies();
std::vector<CoordinateOperationNNPtr> transformationsForGrid;
std::vector<CoordinateOperationNNPtr> 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<CoordinateOperationNNPtr>
transformationsForGridMatchingDatum;
for (const auto &op : transformationsForGrid) {
const auto opTargetCRS =
dynamic_cast<const crs::VerticalCRS *>(
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<size_t>(-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);
Expand All @@ -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<size_t>(-1)) {
const auto &transf = transformationsForGrid[idx];
extent = getExtent(transf, true, dummy);
Expand Down
105 changes: 105 additions & 0 deletions test/unit/test_operationfactory.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<CRS>(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");
Expand Down

0 comments on commit 98dcc57

Please sign in to comment.