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

5 files changed, 116 insertions, 6 deletions

diff --git a/include/proj/io.hpp b/include/proj/io.hpp
index 014d9987..19ed292d 100644
--- a/include/proj/io.hpp
+++ b/include/proj/io.hpp
@@ -873,6 +873,10 @@ class PROJ_GCC_DLL DatabaseContext {
     getNonDeprecated(const std::string &tableName, const std::string &authName,
                      const std::string &code) const;
 
+    PROJ_INTERNAL static std::vector<operation::CoordinateOperationNNPtr>
+    getTransformationsForGridName(const DatabaseContextNNPtr &databaseContext,
+                                  const std::string &gridName);
+
     //! @endcond
 
   protected:
@@ -1184,6 +1188,7 @@ class PROJ_GCC_DLL AuthorityFactory {
                                 std::vector<ObjectType>(),
                             bool approximateMatch = true,
                             size_t limitResultCount = 0) const;
+
     //! @endcond
 
   protected:
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
 
 // ---------------------------------------------------------------------------
diff --git a/test/unit/test_c_api.cpp b/test/unit/test_c_api.cpp
index 947b993c..f6ba708f 100644
--- a/test/unit/test_c_api.cpp
+++ b/test/unit/test_c_api.cpp
@@ -4571,10 +4571,11 @@ TEST_F(CApi, proj_create_vertical_crs_ex) {
     ObjectKeeper keeper_horiz_crs(horiz_crs);
     ASSERT_NE(horiz_crs, nullptr);
 
+    const char *options[] = {"ACCURACY=123", nullptr};
     auto vert_crs = proj_create_vertical_crs_ex(
         m_ctxt, "myVertCRS (ftUS)", "myVertDatum", nullptr, nullptr,
         "US survey foot", 0.304800609601219, "PROJ @foo.gtx", nullptr, nullptr,
-        nullptr, nullptr);
+        nullptr, options);
     ObjectKeeper keeper_vert_crs(vert_crs);
     ASSERT_NE(vert_crs, nullptr);
 
@@ -4609,6 +4610,8 @@ TEST_F(CApi, proj_create_vertical_crs_ex) {
               "+step +proj=vgridshift +grids=@foo.gtx +multiplier=1 "
               "+step +proj=unitconvert +xy_in=rad +xy_out=deg "
               "+step +proj=axisswap +order=2,1");
+
+    ASSERT_EQ(proj_coordoperation_get_accuracy(m_ctxt, P), 123.0);
 }
 
 // ---------------------------------------------------------------------------
@@ -4691,6 +4694,54 @@ TEST_F(CApi, proj_create_vertical_crs_ex_with_geog_crs) {
 
 // ---------------------------------------------------------------------------
 
+TEST_F(CApi, proj_create_vertical_crs_ex_implied_accuracy) {
+
+    PJ *crsH = proj_create(m_ctxt, "EPSG:4283"); // GDA94
+    ASSERT_NE(crsH, nullptr);
+    ObjectKeeper keeper_crsH(crsH);
+    PJ *crsV = proj_create(m_ctxt, "EPSG:5711"); // AHD height
+    ASSERT_NE(crsV, nullptr);
+    ObjectKeeper keeper_crsV(crsV);
+    PJ *crsGeoid = proj_create(m_ctxt, "EPSG:4939"); // GDA94 3D
+    ASSERT_NE(crsGeoid, nullptr);
+    ObjectKeeper keeper_crsGeoid(crsGeoid);
+
+    PJ *vertDatum = proj_crs_get_datum(m_ctxt, crsV);
+    ObjectKeeper keeper_vertDatum(vertDatum);
+    const char *vertDatumName = proj_get_name(vertDatum);
+    const char *vertDatumAuthority = proj_get_id_auth_name(vertDatum, 0);
+    const char *vertDatumCode = proj_get_id_code(vertDatum, 0);
+    PJ *crsVGeoid = proj_create_vertical_crs_ex(
+        m_ctxt, "Vertical", vertDatumName, vertDatumAuthority, vertDatumCode,
+        "metre", 1.0, "PROJ au_ga_AUSGeoid09_V1.01.tif", nullptr, nullptr,
+        crsGeoid, nullptr);
+    ObjectKeeper keeper_crsVGeoid(crsVGeoid);
+    PJ *crsCompoundGeoid = proj_create_compound_crs(
+        m_ctxt, "Compound with geoid", crsH, crsVGeoid);
+    ObjectKeeper keeper_crsCompoundGeoid(crsCompoundGeoid);
+
+    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, crsCompoundGeoid, crsGeoid, ctxt);
+    ASSERT_NE(operations, nullptr);
+    ObjListKeeper keeper_operations(operations);
+    EXPECT_GE(proj_list_get_count(operations), 1);
+    PJ *transform = proj_list_get(m_ctxt, operations, 0);
+    ObjectKeeper keeper_transform(transform);
+
+    // This is the accuracy of operations EPSG:5656 / 5657
+    ASSERT_EQ(proj_coordoperation_get_accuracy(m_ctxt, transform), 0.03);
+}
+
+// ---------------------------------------------------------------------------
+
 TEST_F(CApi, proj_create_derived_geographic_crs) {
 
     PJ *crs_4326 = proj_create(m_ctxt, "EPSG:4326");