proj_create_vertical_crs_ex(): add a ACCURACY option to provide an explicit accuracy, or derive it from the grid name if it is known

3 files changed, 59 insertions, 5 deletions

diff --git a/src/iso19111/c_api.cpp b/src/iso19111/c_api.cpp
index 2e968985..55219004 100644
--- a/src/iso19111/c_api.cpp
+++ b/src/iso19111/c_api.cpp
@@ -3287,7 +3287,9 @@ PJ *proj_create_vertical_crs(PJ_CONTEXT *ctx, const char *crs_name,
  * @param geoid_model_code Code of the transformation for
  * the geoid model. or NULL
  * @param geoid_geog_crs Geographic CRS for the geoid transformation, or NULL.
- * @param options should be set to NULL for now
+ * @param options NULL-terminated list of strings with "KEY=VALUE" format. or
+ * NULL.
+ * The currently recognized option is ACCURACY=value, where value is in metre.
  * @return Object of type VerticalCRS that must be unreferenced with
  * proj_destroy(), or NULL in case of error.
  */
@@ -3299,7 +3301,6 @@ PJ *proj_create_vertical_crs_ex(
     const char *geoid_model_code, const PJ *geoid_geog_crs,
     const char *const *options) {
     SANITIZE_CTX(ctx);
-    (void)options;
     try {
         const UnitOfMeasure linearUnit(
             createLinearUnit(linear_units, linear_units_conv));
@@ -3317,13 +3318,22 @@ PJ *proj_create_vertical_crs_ex(
                                       geoid_geog_crs->iso_obj)
                     ? std::dynamic_pointer_cast<CRS>(geoid_geog_crs->iso_obj)
                     : nullptr;
+
+            std::vector<metadata::PositionalAccuracyNNPtr> accuracies;
+            for (auto iter = options; iter && iter[0]; ++iter) {
+                const char *value;
+                if ((value = getOptionValue(*iter, "ACCURACY="))) {
+                    accuracies.emplace_back(
+                        metadata::PositionalAccuracy::create(value));
+                }
+            }
             const auto model(Transformation::create(
                 propsModel, vertCRSWithoutGeoid,
                 GeographicCRS::EPSG_4979, // arbitrarily chosen. Ignored
                 interpCRS,
                 OperationMethod::create(PropertyMap(),
                                         std::vector<OperationParameterNNPtr>()),
-                {}, {}));
+                {}, accuracies));
             props.set("GEOID_MODEL", model);
         }
         auto vertCRS = VerticalCRS::create(props, datum, cs);
diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp
index 3d5287a3..327af7a1 100644
--- a/src/iso19111/coordinateoperation.cpp
+++ b/src/iso19111/coordinateoperation.cpp
@@ -13738,7 +13738,8 @@ CoordinateOperationFactory::Private::createOperationsGeogToVertFromGeoid(
         return ret;
     };
 
-    const auto getProjGeoidTransformation = [&sourceCRS, &targetCRS, &vertDst](
+    const auto getProjGeoidTransformation = [&sourceCRS, &targetCRS, &vertDst,
+                                             &context](
         const CoordinateOperationNNPtr &model,
         const std::string &projFilename) {
 
@@ -13784,8 +13785,32 @@ CoordinateOperationFactory::Private::createOperationsGeogToVertFromGeoid(
         const auto properties = util::PropertyMap().set(
             common::IdentifiedObject::NAME_KEY,
             buildOpName("Transformation", vertCRSMetre, geogSrcCRS));
+
+        // Try to find a representative value for the accuracy of this grid
+        // from the registered transformations.
+        std::vector<metadata::PositionalAccuracyNNPtr> accuracies;
+        const auto &modelAccuracies = model->coordinateOperationAccuracies();
+        if (modelAccuracies.empty()) {
+            const auto &authFactory = context.context->getAuthorityFactory();
+            if (authFactory) {
+                const auto transformationsForGrid =
+                    io::DatabaseContext::getTransformationsForGridName(
+                        authFactory->databaseContext(), projFilename);
+                double accuracy = -1;
+                for (const auto &transf : transformationsForGrid) {
+                    accuracy = std::max(accuracy, getAccuracy(transf));
+                }
+                if (accuracy >= 0) {
+                    accuracies.emplace_back(
+                        metadata::PositionalAccuracy::create(
+                            toString(accuracy)));
+                }
+            }
+        }
+
         return Transformation::createGravityRelatedHeightToGeographic3D(
-            properties, vertCRSMetre, geogSrcCRS, nullptr, projFilename, {});
+            properties, vertCRSMetre, geogSrcCRS, nullptr, projFilename,
+            !modelAccuracies.empty() ? modelAccuracies : accuracies);
     };
 
     std::vector<CoordinateOperationNNPtr> res;
diff --git a/src/iso19111/factory.cpp b/src/iso19111/factory.cpp
index a9d82268..a011f397 100644
--- a/src/iso19111/factory.cpp
+++ b/src/iso19111/factory.cpp
@@ -1256,6 +1256,25 @@ DatabaseContext::getNonDeprecated(const std::string &tableName,
     return res;
 }
 
+// ---------------------------------------------------------------------------
+
+std::vector<operation::CoordinateOperationNNPtr>
+DatabaseContext::getTransformationsForGridName(
+    const DatabaseContextNNPtr &databaseContext, const std::string &gridName) {
+    auto sqlRes = databaseContext->d->run(
+        "SELECT auth_name, code FROM grid_transformation "
+        "WHERE grid_name = ? OR grid_name = "
+        "(SELECT original_grid_name FROM grid_alternatives "
+        "WHERE proj_grid_name = ?)",
+        {gridName, gridName});
+    std::vector<operation::CoordinateOperationNNPtr> res;
+    for (const auto &row : sqlRes) {
+        res.emplace_back(AuthorityFactory::create(databaseContext, row[0])
+                             ->createCoordinateOperation(row[1], true));
+    }
+    return res;
+}
+
 //! @endcond
 
 // ---------------------------------------------------------------------------