Add proj_get_geoid_models_from_database() (#2681)

to list all geoid model names that apply to a vertical CRS

7 files changed, 212 insertions, 0 deletions

diff --git a/include/proj/io.hpp b/include/proj/io.hpp
index e51dd19a..e56e0994 100644
--- a/include/proj/io.hpp
+++ b/include/proj/io.hpp
@@ -1000,6 +1000,9 @@ class PROJ_GCC_DLL AuthorityFactory {
         const std::string &sourceCRSCode,
         const std::string &targetCRSCode) const;
 
+    PROJ_DLL std::list<std::string>
+    getGeoidModels(const std::string &code) const;
+
     PROJ_DLL const std::string &getAuthority() PROJ_PURE_DECL;
 
     /** Object type. */
diff --git a/scripts/reference_exported_symbols.txt b/scripts/reference_exported_symbols.txt
index 6364ed67..4f853090 100644
--- a/scripts/reference_exported_symbols.txt
+++ b/scripts/reference_exported_symbols.txt
@@ -350,6 +350,7 @@ osgeo::proj::io::AuthorityFactory::getAuthority() const
 osgeo::proj::io::AuthorityFactory::getCelestialBodyList() const
 osgeo::proj::io::AuthorityFactory::getCRSInfoList() const
 osgeo::proj::io::AuthorityFactory::getDescriptionText(std::string const&) const
+osgeo::proj::io::AuthorityFactory::getGeoidModels(std::string const&) const
 osgeo::proj::io::AuthorityFactory::getOfficialNameFromAlias(std::string const&, std::string const&, std::string const&, bool, std::string&, std::string&, std::string&) const
 osgeo::proj::io::AuthorityFactory::getUnitList() const
 osgeo::proj::io::AuthorityFactory::identifyBodyFromSemiMajorAxis(double, double) const
@@ -970,6 +971,7 @@ proj_get_crs_info_list_from_database
 proj_get_crs_list_parameters_create
 proj_get_crs_list_parameters_destroy
 proj_get_ellipsoid
+proj_get_geoid_models_from_database
 proj_get_id_auth_name
 proj_get_id_code
 proj_get_insert_statements
diff --git a/src/iso19111/c_api.cpp b/src/iso19111/c_api.cpp
index ce7f2ae2..f8d9e685 100644
--- a/src/iso19111/c_api.cpp
+++ b/src/iso19111/c_api.cpp
@@ -9162,3 +9162,46 @@ PROJ_STRING_LIST proj_get_insert_statements(
     }
     return nullptr;
 }
+
+// ---------------------------------------------------------------------------
+
+/** \brief Returns a list of geoid models available for that crs
+ *
+ * The list includes the geoid models connected directly with the crs,
+ * or via "Height Depth Reversal" or "Change of Vertical Unit" transformations.
+ * The returned list is NULL terminated and must be freed with
+ * proj_string_list_destroy().
+ *
+ * @param ctx Context, or NULL for default context.
+ * @param auth_name Authority name (must not be NULL)
+ * @param code Object code (must not be NULL)
+ * @param options should be set to NULL for now
+ * @return list of geoid models names (to be freed with
+ * proj_string_list_destroy()), or NULL in case of error.
+ * @since 8.1
+ */
+PROJ_STRING_LIST
+proj_get_geoid_models_from_database(PJ_CONTEXT *ctx, const char *auth_name,
+                                    const char *code,
+                                    const char *const *options) {
+    SANITIZE_CTX(ctx);
+    if (!auth_name || !code) {
+        proj_context_errno_set(ctx, PROJ_ERR_OTHER_API_MISUSE);
+        proj_log_error(ctx, __FUNCTION__, "missing required input");
+        return nullptr;
+    }
+    (void)options;
+    try {
+        const std::string codeStr(code);
+        auto factory = AuthorityFactory::create(getDBcontext(ctx), auth_name);
+        auto geoidModels = factory->getGeoidModels(codeStr);
+        ctx->safeAutoCloseDbIfNeeded();
+        return to_string_list(std::move(geoidModels));
+    } catch (const std::exception &e) {
+        proj_log_error(ctx, __FUNCTION__, e.what());
+    }
+    ctx->safeAutoCloseDbIfNeeded();
+    return nullptr;
+}
+
+// ---------------------------------------------------------------------------
diff --git a/src/iso19111/factory.cpp b/src/iso19111/factory.cpp
index fabf57e1..d95ccbf4 100644
--- a/src/iso19111/factory.cpp
+++ b/src/iso19111/factory.cpp
@@ -5765,6 +5765,88 @@ AuthorityFactory::createFromCoordinateReferenceSystemCodes(
 
 // ---------------------------------------------------------------------------
 
+/** \brief Returns a list of geoid models available for that crs
+ *
+ * The list includes the geoid models connected directly with the crs,
+ * or via "Height Depth Reversal" or "Change of Vertical Unit" transformations
+ *
+ * @param code crs code allocated by authority.
+ * @return list of geoid model names
+ * @throw FactoryException
+ */
+
+std::list<std::string>
+AuthorityFactory::getGeoidModels(const std::string &code) const {
+
+    ListOfParams params;
+    std::string sql;
+    sql += "SELECT DISTINCT GM0.name "
+           " FROM geoid_model GM0 "
+           "INNER JOIN grid_transformation GT0 "
+           " ON  GT0.code = GM0.operation_code "
+           " AND GT0.auth_name = GM0.operation_auth_name "
+           " AND GT0.target_crs_code = ? ";
+    params.emplace_back(code);
+    if (d->hasAuthorityRestriction()) {
+        sql += " AND GT0.target_crs_auth_name = ? ";
+        params.emplace_back(d->authority());
+    }
+
+    /// The second part of the query is for CRSs that use that geoid model via
+    /// Height Depth Reversal (EPSG:1068) or Change of Vertical Unit (EPSG:1069)
+    sql += "UNION "
+           "SELECT DISTINCT GM0.name "
+           " FROM geoid_model GM0 "
+           "INNER JOIN grid_transformation GT1 "
+           " ON  GT1.code = GM0.operation_code "
+           " AND GT1.auth_name = GM0.operation_auth_name "
+           "INNER JOIN other_transformation OT1 "
+           " ON  OT1.source_crs_code = GT1.target_crs_code "
+           " AND OT1.source_crs_auth_name = GT1.target_crs_auth_name "
+           " AND OT1.method_auth_name = 'EPSG' "
+           " AND OT1.method_code IN (1068, 1069, 1104) "
+           " AND OT1.target_crs_code = ? ";
+    params.emplace_back(code);
+    if (d->hasAuthorityRestriction()) {
+        sql += " AND OT1.target_crs_auth_name = ? ";
+        params.emplace_back(d->authority());
+    }
+
+    /// The third part of the query is for CRSs that use that geoid model via
+    /// other_transformation table twice, like transforming depth and feet
+    sql += "UNION "
+           "SELECT DISTINCT GM0.name "
+           " FROM geoid_model GM0 "
+           "INNER JOIN grid_transformation GT1 "
+           " ON  GT1.code = GM0.operation_code "
+           " AND GT1.auth_name = GM0.operation_auth_name "
+           "INNER JOIN other_transformation OT1 "
+           " ON  OT1.source_crs_code = GT1.target_crs_code "
+           " AND OT1.source_crs_auth_name = GT1.target_crs_auth_name "
+           " AND OT1.method_auth_name = 'EPSG' "
+           " AND OT1.method_code IN (1068, 1069, 1104) "
+           "INNER JOIN other_transformation OT2 "
+           " ON  OT2.source_crs_code = OT1.target_crs_code "
+           " AND OT2.source_crs_auth_name = OT1.target_crs_auth_name "
+           " AND OT2.method_code IN (1068, 1069, 1104) "
+           " AND OT2.target_crs_code = ? ";
+    params.emplace_back(code);
+    if (d->hasAuthorityRestriction()) {
+        sql += " AND OT2.target_crs_auth_name = ? ";
+        params.emplace_back(d->authority());
+    }
+    sql += " ORDER BY 1 ";
+
+    auto sqlRes = d->run(sql, params);
+    std::list<std::string> res;
+    for (const auto &row : sqlRes) {
+        res.push_back(row[0]);
+    }
+    return res;
+}
+
+// ---------------------------------------------------------------------------
+
 /** \brief Returns a list operation::CoordinateOperation between two CRS.
  *
  * The list is ordered with preferred operations first. No attempt is made
diff --git a/src/proj.h b/src/proj.h
index 49bf991c..0f78271a 100644
--- a/src/proj.h
+++ b/src/proj.h
@@ -1188,6 +1188,12 @@ PJ_OBJ_LIST PROJ_DLL *proj_identify(PJ_CONTEXT *ctx,
                                         const char* const *options,
                                         int **out_confidence);
 
+PROJ_STRING_LIST PROJ_DLL proj_get_geoid_models_from_database(
+                               PJ_CONTEXT *ctx,
+                               const char *auth_name,
+                               const char *code,
+                               const char *const *options);
+
 void PROJ_DLL proj_int_list_destroy(int* list);
 
 /* ------------------------------------------------------------------------- */
diff --git a/test/unit/test_c_api.cpp b/test/unit/test_c_api.cpp
index 8fa1b228..dfc92139 100644
--- a/test/unit/test_c_api.cpp
+++ b/test/unit/test_c_api.cpp
@@ -5619,5 +5619,25 @@ TEST_F(CApi, proj_get_insert_statements) {
         proj_insert_object_session_destroy(m_ctxt, session);
     }
 }
+// ---------------------------------------------------------------------------
+
+TEST_F(CApi, proj_get_geoid_models_from_database) {
+    auto findInList = [](PROJ_STRING_LIST list, const std::string &ref) {
+        while (list && *list) {
+            if (std::string(*list) == ref) {
+                return true;
+            }
+            list++;
+        }
+        return false;
+    };
+
+    auto list =
+        proj_get_geoid_models_from_database(m_ctxt, "EPSG", "5703", nullptr);
+    ListFreer freer(list);
+    EXPECT_TRUE(findInList(list, "GEOID12B"));
+    EXPECT_TRUE(findInList(list, "GEOID18"));
+    EXPECT_FALSE(findInList(list, "OSGM15"));
+}
 
 } // namespace
diff --git a/test/unit/test_factory.cpp b/test/unit/test_factory.cpp
index c10c4697..37013775 100644
--- a/test/unit/test_factory.cpp
+++ b/test/unit/test_factory.cpp
@@ -2139,6 +2139,62 @@ TEST(
     }
 }
 
+TEST(factory, AuthorityFactory_getAvailableGeoidmodels) {
+
+    const std::string OSGM15{"OSGM15"};
+    const std::string GEOID12B{"GEOID12B"};
+    const std::string GEOID18{"GEOID18"};
+
+    auto checkNavd88 = [&](const std::list<std::string> &res) {
+        EXPECT_TRUE(res.end() != std::find(res.begin(), res.end(), GEOID12B));
+        EXPECT_TRUE(res.end() != std::find(res.begin(), res.end(), GEOID18));
+        EXPECT_FALSE(res.end() != std::find(res.begin(), res.end(), OSGM15));
+    };
+
+    auto checkOdn = [&](const std::list<std::string> &res) {
+        EXPECT_FALSE(res.end() != std::find(res.begin(), res.end(), GEOID12B));
+        EXPECT_FALSE(res.end() != std::find(res.begin(), res.end(), GEOID18));
+        EXPECT_TRUE(res.end() != std::find(res.begin(), res.end(), OSGM15));
+    };
+
+    auto factory = AuthorityFactory::create(DatabaseContext::create(), "EPSG");
+
+    {
+        auto res = factory->getGeoidModels("4326");
+        ASSERT_TRUE(res.empty());
+    }
+
+    {
+        auto res = factory->getGeoidModels("5703"); // "NAVD88 height"
+        checkNavd88(res);
+    }
+    {
+        auto res = factory->getGeoidModels("6360"); // "NAVD88 height (ftUS)"
+        checkNavd88(res);
+    }
+    {
+        auto res = factory->getGeoidModels("8228"); // "NAVD88 height (ft)"
+        checkNavd88(res);
+    }
+    {
+        auto res = factory->getGeoidModels("6357"); // "NAVD88 depth"
+        checkNavd88(res);
+    }
+    {
+        auto res = factory->getGeoidModels("6358"); // "NAVD88 depth (ftUS)"
+        checkNavd88(res);
+    }
+
+    {
+        auto res = factory->getGeoidModels("5701"); // "ODN height"
+        checkOdn(res);
+    }
+    {
+        auto res = factory->getGeoidModels("5732"); // "Belfast height"
+        checkOdn(res);
+    }
+}
+
 // ---------------------------------------------------------------------------
 
 TEST_F(FactoryWithTmpDatabase,