Add a geoid_model name in database, use GEOIDMODEL for transformations, add a proj_create_vertical_crs_ex()

4 files changed, 245 insertions, 9 deletions

diff --git a/src/iso19111/c_api.cpp b/src/iso19111/c_api.cpp
index f5f7ba55..32c8df45 100644
--- a/src/iso19111/c_api.cpp
+++ b/src/iso19111/c_api.cpp
@@ -2603,13 +2603,19 @@ int proj_coordoperation_get_method_info(PJ_CONTEXT *ctx,
 // ---------------------------------------------------------------------------
 
 //! @cond Doxygen_Suppress
-static PropertyMap createPropertyMapName(const char *c_name) {
+static PropertyMap createPropertyMapName(const char *c_name,
+                                         const char *auth_name = nullptr,
+                                         const char *code = nullptr) {
     std::string name(c_name ? c_name : "unnamed");
     PropertyMap properties;
     if (ends_with(name, " (deprecated)")) {
         name.resize(name.size() - strlen(" (deprecated)"));
         properties.set(common::IdentifiedObject::DEPRECATED_KEY, true);
     }
+    if (auth_name && code) {
+        properties.set(metadata::Identifier::CODESPACE_KEY, auth_name);
+        properties.set(metadata::Identifier::CODE_KEY, code);
+    }
     return properties.set(common::IdentifiedObject::NAME_KEY, name);
 }
 
@@ -2934,15 +2940,73 @@ PJ *proj_create_vertical_crs(PJ_CONTEXT *ctx, const char *crs_name,
                              const char *datum_name, const char *linear_units,
                              double linear_units_conv) {
 
+    return proj_create_vertical_crs_ex(
+        ctx, crs_name, datum_name, nullptr, nullptr, linear_units,
+        linear_units_conv, nullptr, nullptr, nullptr, nullptr, nullptr);
+}
+
+// ---------------------------------------------------------------------------
+
+/** \brief Create a VerticalCRS
+ *
+ * The returned object must be unreferenced with proj_destroy() after
+ * use.
+ * It should be used by at most one thread at a time.
+ *
+ * @param ctx PROJ context, or NULL for default context
+ * @param crs_name Name of the GeographicCRS. Or NULL
+ * @param datum_name Name of the VerticalReferenceFrame. Or NULL
+ * @param datum_auth_name Authority name of the VerticalReferenceFrame. Or NULL
+ * @param datum_code Code of the VerticalReferenceFrame. Or NULL
+ * @param linear_units Name of the linear units. Or NULL for Metre
+ * @param linear_units_conv Conversion factor from the linear unit to metre. Or
+ * 0 for Metre if linear_units == NULL. Otherwise should be not NULL
+ * @param geoid_model_name Geoid model name, or NULL. Can be a name from the
+ * geoid_model name or a string "PROJ foo.gtx"
+ * @param geoid_model_auth_name Authority name of the transformation for
+ * the geoid model. or NULL
+ * @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
+ * @return Object of type VerticalCRS that must be unreferenced with
+ * proj_destroy(), or NULL in case of error.
+ */
+PJ *proj_create_vertical_crs_ex(
+    PJ_CONTEXT *ctx, const char *crs_name, const char *datum_name,
+    const char *datum_auth_name, const char *datum_code,
+    const char *linear_units, double linear_units_conv,
+    const char *geoid_model_name, const char *geoid_model_auth_name,
+    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));
-        auto datum =
-            VerticalReferenceFrame::create(createPropertyMapName(datum_name));
-        auto vertCRS = VerticalCRS::create(
-            createPropertyMapName(crs_name), datum,
-            cs::VerticalCS::createGravityRelatedHeight(linearUnit));
+        auto datum = VerticalReferenceFrame::create(
+            createPropertyMapName(datum_name, datum_auth_name, datum_code));
+        auto props = createPropertyMapName(crs_name);
+        auto cs = cs::VerticalCS::createGravityRelatedHeight(linearUnit);
+        if (geoid_model_name) {
+            auto propsModel = createPropertyMapName(
+                geoid_model_name, geoid_model_auth_name, geoid_model_code);
+            const auto vertCRSWithoutGeoid =
+                VerticalCRS::create(props, datum, cs);
+            const auto interpCRS =
+                geoid_geog_crs && std::dynamic_pointer_cast<GeographicCRS>(
+                                      geoid_geog_crs->iso_obj)
+                    ? std::dynamic_pointer_cast<CRS>(geoid_geog_crs->iso_obj)
+                    : nullptr;
+            const auto model(Transformation::create(
+                propsModel, vertCRSWithoutGeoid, GeographicCRS::EPSG_4979,
+                interpCRS,
+                OperationMethod::create(PropertyMap(),
+                                        std::vector<OperationParameterNNPtr>()),
+                {}, {}));
+            props.set("GEOID_MODEL", model);
+        }
+        auto vertCRS = VerticalCRS::create(props, datum, cs);
         return pj_obj_create(ctx, vertCRS);
     } catch (const std::exception &e) {
         proj_log_error(ctx, __FUNCTION__, e.what());
diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp
index d0a1761c..e8f195e3 100644
--- a/src/iso19111/coordinateoperation.cpp
+++ b/src/iso19111/coordinateoperation.cpp
@@ -10442,6 +10442,12 @@ struct CoordinateOperationFactory::Private {
         std::vector<CoordinateOperationNNPtr> &res);
 
     static std::vector<CoordinateOperationNNPtr>
+    createOperationsGeogToVertFromGeoid(const crs::CRSNNPtr &sourceCRS,
+                                        const crs::CRSNNPtr &targetCRS,
+                                        const crs::VerticalCRS *vertDst,
+                                        Context &context);
+
+    static std::vector<CoordinateOperationNNPtr>
     createOperationsGeogToVertWithIntermediateVert(
         const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
         const crs::VerticalCRS *vertDst, Context &context);
@@ -12775,6 +12781,18 @@ bool CoordinateOperationFactory::Private::createOperationsFromDatabase(
 
     ENTER_FUNCTION();
 
+    if (geogSrc && vertDst) {
+        res = createOperationsGeogToVertFromGeoid(sourceCRS, targetCRS, vertDst,
+                                                  context);
+    } else if (geogDst && vertSrc) {
+        res = applyInverse(createOperationsGeogToVertFromGeoid(
+            targetCRS, sourceCRS, vertSrc, context));
+    }
+
+    if (!res.empty()) {
+        return true;
+    }
+
     res = findOpsInRegistryDirect(sourceCRS, targetCRS, context);
 
     // If we get at least a result with perfect accuracy, do not
@@ -12898,6 +12916,116 @@ findCandidateVertCRSForDatum(const io::AuthorityFactoryPtr &authFactory,
 
 // ---------------------------------------------------------------------------
 
+std::vector<CoordinateOperationNNPtr>
+CoordinateOperationFactory::Private::createOperationsGeogToVertFromGeoid(
+    const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
+    const crs::VerticalCRS *vertDst, Private::Context &context) {
+
+    ENTER_FUNCTION();
+
+    const auto useTransf = [&targetCRS, &context,
+                            vertDst](const CoordinateOperationNNPtr &op) {
+        const auto targetOp =
+            dynamic_cast<const crs::VerticalCRS *>(op->targetCRS().get());
+        assert(targetOp);
+        if (targetOp->_isEquivalentTo(
+                vertDst, util::IComparable::Criterion::EQUIVALENT)) {
+            return op;
+        }
+        std::vector<CoordinateOperationNNPtr> tmp;
+        createOperationsVertToVert(NN_NO_CHECK(op->targetCRS()), targetCRS,
+                                   context, targetOp, vertDst, tmp);
+        assert(!tmp.empty());
+        auto ret = ConcatenatedOperation::createComputeMetadata(
+            {op, tmp.front()}, !allowEmptyIntersection);
+        return ret;
+    };
+
+    const auto getProjGeoidTransformation = [&sourceCRS, &targetCRS, &vertDst](
+        const CoordinateOperationNNPtr &model,
+        const std::string &projFilename) {
+
+        const auto getNameVertCRSMetre = [](const std::string &name) {
+            if (name.empty())
+                return std::string("unnamed");
+            auto ret(name);
+            bool haveOriginalUnit = false;
+            if (name.back() == ')') {
+                const auto pos = ret.rfind(" (");
+                if (pos != std::string::npos) {
+                    haveOriginalUnit = true;
+                    ret = ret.substr(0, pos);
+                }
+            }
+            const auto pos = ret.rfind(" depth");
+            if (pos != std::string::npos) {
+                ret = ret.substr(0, pos) + " height";
+            }
+            if (!haveOriginalUnit) {
+                ret += " (metre)";
+            }
+            return ret;
+        };
+
+        const auto &axis = vertDst->coordinateSystem()->axisList()[0];
+        const auto geogSrcCRS =
+            dynamic_cast<crs::GeographicCRS *>(model->interpolationCRS().get())
+                ? NN_NO_CHECK(model->interpolationCRS())
+                : sourceCRS;
+        const auto vertCRSMetre =
+            axis->unit() == common::UnitOfMeasure::METRE &&
+                    axis->direction() == cs::AxisDirection::UP
+                ? targetCRS
+                : util::nn_static_pointer_cast<crs::CRS>(
+                      crs::VerticalCRS::create(
+                          util::PropertyMap().set(
+                              common::IdentifiedObject::NAME_KEY,
+                              getNameVertCRSMetre(targetCRS->nameStr())),
+                          vertDst->datum(), vertDst->datumEnsemble(),
+                          cs::VerticalCS::createGravityRelatedHeight(
+                              common::UnitOfMeasure::METRE)));
+        const auto properties = util::PropertyMap().set(
+            common::IdentifiedObject::NAME_KEY,
+            buildOpName("Transformation", vertCRSMetre, geogSrcCRS));
+        return Transformation::createGravityRelatedHeightToGeographic3D(
+            properties, vertCRSMetre, geogSrcCRS, nullptr, projFilename, {});
+    };
+
+    std::vector<CoordinateOperationNNPtr> res;
+    const auto &authFactory = context.context->getAuthorityFactory();
+    if (authFactory) {
+        const auto &models = vertDst->geoidModel();
+        for (const auto &model : models) {
+            const auto &modelName = model->nameStr();
+            const auto transformations =
+                starts_with(modelName, "PROJ ")
+                    ? std::vector<
+                          CoordinateOperationNNPtr>{getProjGeoidTransformation(
+                          model, modelName.substr(strlen("PROJ ")))}
+                    : authFactory->getTransformationsForGeoid(
+                          modelName,
+                          context.context->getUsePROJAlternativeGridNames());
+            for (const auto &transf : transformations) {
+                if (dynamic_cast<crs::GeographicCRS *>(
+                        transf->sourceCRS().get()) &&
+                    dynamic_cast<crs::VerticalCRS *>(
+                        transf->targetCRS().get())) {
+                    res.push_back(useTransf(transf));
+                } else if (dynamic_cast<crs::GeographicCRS *>(
+                               transf->targetCRS().get()) &&
+                           dynamic_cast<crs::VerticalCRS *>(
+                               transf->sourceCRS().get())) {
+                    res.push_back(useTransf(transf->inverse()));
+                }
+            }
+        }
+    }
+
+    return res;
+}
+
+// ---------------------------------------------------------------------------
+
 std::vector<CoordinateOperationNNPtr> CoordinateOperationFactory::Private::
     createOperationsGeogToVertWithIntermediateVert(
         const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
@@ -14223,16 +14351,25 @@ getResolvedCRS(const crs::CRSNNPtr &crs,
         } else {
             auto outCrs = tryToIdentifyByName(
                 io::AuthorityFactory::ObjectType::COMPOUND_CRS);
+            const auto &components = compoundCrs->componentReferenceSystems();
             if (outCrs.get() != crs.get()) {
-                return outCrs;
+                bool hasGeoid = false;
+                if (components.size() == 2) {
+                    auto vertCRS =
+                        dynamic_cast<crs::VerticalCRS *>(components[1].get());
+                    if (vertCRS && !vertCRS->geoidModel().empty()) {
+                        hasGeoid = true;
+                    }
+                }
+                if (!hasGeoid) {
+                    return outCrs;
+                }
             }
             if (approxExtent || !extentOut) {
                 // If we still did not get a reliable extent, then try to
                 // resolve the components of the compoundCRS, and take the
                 // intersection of their extent.
                 extentOut = metadata::ExtentPtr();
-                const auto &components =
-                    compoundCrs->componentReferenceSystems();
                 for (const auto &component : components) {
                     metadata::ExtentPtr componentExtent;
                     getResolvedCRS(component, context, componentExtent);
diff --git a/src/iso19111/factory.cpp b/src/iso19111/factory.cpp
index 307c3a07..d9917996 100644
--- a/src/iso19111/factory.cpp
+++ b/src/iso19111/factory.cpp
@@ -5380,6 +5380,28 @@ AuthorityFactory::getPreferredHubGeodeticReferenceFrames(
 // ---------------------------------------------------------------------------
 
 //! @cond Doxygen_Suppress
+std::vector<operation::CoordinateOperationNNPtr>
+AuthorityFactory::getTransformationsForGeoid(
+    const std::string &geoidName, bool usePROJAlternativeGridNames) const {
+    std::vector<operation::CoordinateOperationNNPtr> res;
+
+    const std::string sql("SELECT operation_auth_name, operation_code FROM "
+                          "geoid_model WHERE name = ?");
+    auto sqlRes = d->run(sql, {geoidName});
+    for (const auto &row : sqlRes) {
+        const auto &auth_name = row[0];
+        const auto &code = row[1];
+        res.emplace_back(d->createFactory(auth_name)->createCoordinateOperation(
+            code, usePROJAlternativeGridNames));
+    }
+
+    return res;
+}
+//! @endcond
+
+// ---------------------------------------------------------------------------
+
+//! @cond Doxygen_Suppress
 FactoryException::FactoryException(const char *message) : Exception(message) {}
 
 // ---------------------------------------------------------------------------
diff --git a/src/proj_experimental.h b/src/proj_experimental.h
index c6d5bc45..0886ba28 100644
--- a/src/proj_experimental.h
+++ b/src/proj_experimental.h
@@ -248,6 +248,19 @@ PJ PROJ_DLL *proj_create_vertical_crs(PJ_CONTEXT *ctx,
                                               const char *linear_units,
                                               double linear_units_conv);
 
+PJ PROJ_DLL *proj_create_vertical_crs_ex(PJ_CONTEXT *ctx,
+                                              const char *crs_name,
+                                              const char *datum_name,
+                                              const char *datum_auth_name,
+                                              const char* datum_code,
+                                              const char *linear_units,
+                                              double linear_units_conv,
+                                              const char* geoid_model_name,
+                                              const char* geoid_model_auth_name,
+                                              const char* geoid_model_code,
+                                              const PJ* geoid_geog_crs,
+                                              const char *const *options);
+
 PJ PROJ_DLL *proj_create_compound_crs(PJ_CONTEXT *ctx,
                                               const char *crs_name,
                                               PJ* horiz_crs,