createOperations(): make sure to associate an extent to the transform…

… of a CRS with a GEOIDMODEL using a PROJ grid, so that it is later used instead of a ballpark operation (fixes #2768)
OSGeo · Jul 6, 2021 · aa44d91 · aa44d91
1 parent fc983b6
commit aa44d91
Show file tree

Hide file tree

Showing 2 changed files with 40 additions and 3 deletions.
diff --git a/src/iso19111/operation/coordinateoperationfactory.cpp b/src/iso19111/operation/coordinateoperationfactory.cpp
@@ -3476,7 +3476,7 @@ CoordinateOperationFactory::Private::createOperationsGeogToVertFromGeoid(
                               vertDst->datum(), vertDst->datumEnsemble(),
                               cs::VerticalCS::createGravityRelatedHeight(
                                   common::UnitOfMeasure::METRE)));
-            const auto properties = util::PropertyMap().set(
+            auto properties = util::PropertyMap().set(
                 common::IdentifiedObject::NAME_KEY,
                 buildOpName("Transformation", vertCRSMetre, geogSrcCRS));
 
@@ -3485,12 +3485,13 @@ CoordinateOperationFactory::Private::createOperationsGeogToVertFromGeoid(
             std::vector<metadata::PositionalAccuracyNNPtr> accuracies;
             const auto &modelAccuracies =
                 model->coordinateOperationAccuracies();
+            std::vector<CoordinateOperationNNPtr> transformationsForGrid;
+            double accuracy = -1;
             if (modelAccuracies.empty()) {
                 if (authFactory) {
-                    const auto transformationsForGrid =
+                    transformationsForGrid =
                         io::DatabaseContext::getTransformationsForGridName(
                             authFactory->databaseContext(), projFilename);
-                    double accuracy = -1;
                     for (const auto &transf : transformationsForGrid) {
                         accuracy = std::max(accuracy, getAccuracy(transf));
                     }
@@ -3502,6 +3503,30 @@ CoordinateOperationFactory::Private::createOperationsGeogToVertFromGeoid(
                 }
             }
 
+            // Set extent
+            bool dummy = false;
+            // Use in priority the one of the geoid model transformation
+            auto extent = getExtent(model, true, dummy);
+            // 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);
+                }
+                for (const auto &transf : transformationsForGrid) {
+                    if (accuracy == -1 || accuracy == getAccuracy(transf)) {
+                        extent = getExtent(transf, true, dummy);
+                        break;
+                    }
+                }
+            }
+            if (extent) {
+                properties.set(common::ObjectUsage::DOMAIN_OF_VALIDITY_KEY,
+                               NN_NO_CHECK(extent));
+            }
+
             return Transformation::createGravityRelatedHeightToGeographic3D(
                 properties, vertCRSMetre, geogSrcCRS, nullptr, projFilename,
                 !modelAccuracies.empty() ? modelAccuracies : accuracies);

diff --git a/test/unit/test_c_api.cpp b/test/unit/test_c_api.cpp
@@ -5083,6 +5083,18 @@ TEST_F(CApi, proj_create_vertical_crs_ex_implied_accuracy) {
 
     // This is the accuracy of operations EPSG:5656 / 5657
     ASSERT_EQ(proj_coordoperation_get_accuracy(m_ctxt, transform), 0.15);
+
+    // Check there's an asssociated area of use
+    double west_lon_degree = 0;
+    double south_lat_degree = 0;
+    double east_lon_degree = 0;
+    double north_lat_degree = 0;
+    ASSERT_EQ(proj_get_area_of_use(m_ctxt, transform, &west_lon_degree,
+                                   &south_lat_degree, &east_lon_degree,
+                                   &north_lat_degree, nullptr),
+              true);
+    EXPECT_LE(north_lat_degree, -10);
+    EXPECT_GE(west_lon_degree, 110);
 }
 
 // ---------------------------------------------------------------------------