Add support for GEOIDMODEL

data/sql/commit.sql

@@ -29,6 +29,22 @@ FOR EACH ROW BEGIN
         WHERE (SELECT 1 FROM object_view LIMIT 1) = 0;
     SELECT RAISE(ABORT, 'corrupt definition of authority_list')
         WHERE (SELECT 1 FROM authority_list LIMIT 1) = 0;
+
+    -- check geoid_model table
+    SELECT RAISE(ABORT, 'missing GEOID99 in geoid_model')
+        WHERE NOT EXISTS(SELECT 1 FROM geoid_model WHERE name = 'GEOID99');
+    SELECT RAISE(ABORT, 'missing GEOID03 in geoid_model')
+        WHERE NOT EXISTS(SELECT 1 FROM geoid_model WHERE name = 'GEOID03');
+    SELECT RAISE(ABORT, 'missing GEOID06 in geoid_model')
+        WHERE NOT EXISTS(SELECT 1 FROM geoid_model WHERE name = 'GEOID06');
+    SELECT RAISE(ABORT, 'missing GEOID09 in geoid_model')
+        WHERE NOT EXISTS(SELECT 1 FROM geoid_model WHERE name = 'GEOID09');
+    SELECT RAISE(ABORT, 'missing GEOID12A in geoid_model')
+        WHERE NOT EXISTS(SELECT 1 FROM geoid_model WHERE name = 'GEOID12A');
+    SELECT RAISE(ABORT, 'missing GEOID12B in geoid_model')
+        WHERE NOT EXISTS(SELECT 1 FROM geoid_model WHERE name = 'GEOID12B');
+    SELECT RAISE(ABORT, 'missing GEOID18 in geoid_model')
+        WHERE NOT EXISTS(SELECT 1 FROM geoid_model WHERE name = 'GEOID18');
 END;
 INSERT INTO dummy DEFAULT VALUES;
 DROP TRIGGER final_checks;

data/sql/customizations.sql

@@ -95,3 +95,20 @@ INSERT INTO "helmert_transformation" VALUES('PROJ','WGS84_TO_WGS84_TRANSIT','WGS
 DELETE FROM supersession WHERE superseded_table_name = 'grid_transformation' AND superseded_auth_name = 'EPSG' AND superseded_code = '6326';
 DELETE FROM supersession WHERE superseded_table_name = 'grid_transformation' AND superseded_auth_name = 'EPSG' AND superseded_code = '7646';
 DELETE FROM supersession WHERE superseded_table_name = 'grid_transformation' AND superseded_auth_name = 'EPSG' AND superseded_code = '7647';
+
+---- Geoid models -----
+
+INSERT INTO "geoid_model" SELECT 'GEOID99', auth_name, code FROM grid_transformation WHERE auth_name = 'EPSG' AND grid_name LIKE 'g1999%' AND deprecated = 0;
+
+INSERT INTO "geoid_model" SELECT 'GEOID03', auth_name, code FROM grid_transformation WHERE auth_name = 'EPSG' AND grid_name LIKE 'geoid03%' AND deprecated = 0;
+
+INSERT INTO "geoid_model" SELECT 'GEOID06', auth_name, code FROM grid_transformation WHERE auth_name = 'EPSG' AND grid_name LIKE 'geoid06%' AND deprecated = 0;
+
+INSERT INTO "geoid_model" SELECT 'GEOID09', auth_name, code FROM grid_transformation WHERE auth_name = 'EPSG' AND grid_name LIKE 'geoid09%' AND deprecated = 0;
+
+-- Geoid12A and Geoid12B are identical
+INSERT INTO "geoid_model" SELECT 'GEOID12A', auth_name, code FROM grid_transformation WHERE auth_name = 'EPSG' AND grid_name LIKE 'g2012b%' AND deprecated = 0;
+
+INSERT INTO "geoid_model" SELECT 'GEOID12B', auth_name, code FROM grid_transformation WHERE auth_name = 'EPSG' AND grid_name LIKE 'g2012b%' AND deprecated = 0;
+
+INSERT INTO "geoid_model" SELECT 'GEOID18', auth_name, code FROM grid_transformation WHERE auth_name = 'EPSG' AND grid_name LIKE 'g2018%' AND deprecated = 0;

data/sql/proj_db_table_defs.sql

@@ -1365,6 +1365,23 @@ FOR EACH ROW BEGIN
 
 END;
 
+
+CREATE TABLE geoid_model(
+    name TEXT NOT NULL,
+    operation_auth_name TEXT NOT NULL,
+    operation_code TEXT NOT NULL,
+    CONSTRAINT pk_geoid_model PRIMARY KEY (name, operation_auth_name, operation_code)
+    -- CONSTRATINT fk_geoid_model_operation FOREIGN KEY (operation_auth_name, operation_code) REFERENCES coordinate_operation(auth_name, code)
+);
+
+CREATE TRIGGER geoid_model_insert_trigger
+BEFORE INSERT ON geoid_model
+FOR EACH ROW BEGIN
+    SELECT RAISE(ABORT, 'insert on geoid_model violates constraint: (operation_auth_name, operation_code) must already exist in coordinate_operation_with_conversion_view')
+        WHERE NOT EXISTS (SELECT 1 FROM coordinate_operation_with_conversion_view covwv WHERE covwv.auth_name = NEW.operation_auth_name AND covwv.code = NEW.operation_code);
+END;
+
+
 CREATE TABLE alias_name(
     table_name TEXT NOT NULL CHECK (table_name IN (
         'unit_of_measure', 'celestial_body', 'ellipsoid', 

include/proj/io.hpp

@@ -1118,6 +1118,10 @@ class PROJ_GCC_DLL AuthorityFactory {
     getPreferredHubGeodeticReferenceFrames(
         const std::string &geodeticReferenceFrameCode) const;
 
+    PROJ_INTERNAL std::vector<operation::CoordinateOperationNNPtr>
+    getTransformationsForGeoid(const std::string &geoidName,
+                               bool usePROJAlternativeGridNames) const;
+
     //! @endcond
 
   protected:

scripts/reference_exported_symbols.txt

@@ -899,6 +899,7 @@ proj_create_operations
 proj_create_projected_crs
 proj_create_transformation
 proj_create_vertical_crs
+proj_create_vertical_crs_ex
 proj_crs_alter_cs_angular_unit
 proj_crs_alter_cs_linear_unit
 proj_crs_alter_geodetic_crs

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());

src/iso19111/coordinateoperation.cpp

@@ -10441,6 +10441,12 @@ struct CoordinateOperationFactory::Private {
         const crs::VerticalCRS *vertDst,
         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,
@@ -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);

src/iso19111/factory.cpp

@@ -5379,6 +5379,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) {}