Merge pull request #1710 from rouault/geoid_model

Add support for GEOIDMODEL

7 files changed, 391 insertions, 27 deletions

diff --git a/src/iso19111/c_api.cpp b/src/iso19111/c_api.cpp
index 337de313..f3badb51 100644
--- a/src/iso19111/c_api.cpp
+++ b/src/iso19111/c_api.cpp
@@ -2615,13 +2615,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);
 }
 
@@ -2946,15 +2952,76 @@ 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.
+ *
+ * This is an extented (_ex) version of proj_create_vertical_crs() that adds
+ * the capability of defining a geoid model.
+ *
+ * @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 0983c471..c51b4d3d 100644
--- a/src/iso19111/coordinateoperation.cpp
+++ b/src/iso19111/coordinateoperation.cpp
@@ -10480,6 +10480,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);
@@ -12813,6 +12819,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
@@ -12936,6 +12954,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,
@@ -14261,16 +14389,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/crs.cpp b/src/iso19111/crs.cpp
index 26725e24..b9b38c80 100644
--- a/src/iso19111/crs.cpp
+++ b/src/iso19111/crs.cpp
@@ -2478,6 +2478,14 @@ void VerticalCRS::_exportToWKT(io::WKTFormatter *formatter) const {
     cs->_exportToWKT(formatter);
     formatter->setOutputAxis(oldAxisOutputRule);
 
+    if (isWKT2 && formatter->use2019Keywords() && !d->geoidModel.empty()) {
+        const auto &model = d->geoidModel[0];
+        formatter->startNode(io::WKTConstants::GEOIDMODEL, false);
+        formatter->addQuotedString(model->nameStr());
+        model->formatID(formatter);
+        formatter->endNode();
+    }
+
     ObjectUsage::baseExportToWKT(formatter);
     formatter->endNode();
 }
@@ -2539,6 +2547,15 @@ void VerticalCRS::_exportToJSON(
     formatter->setOmitTypeInImmediateChild();
     coordinateSystem()->_exportToJSON(formatter);
 
+    if (!d->geoidModel.empty()) {
+        const auto &model = d->geoidModel[0];
+        writer.AddObjKey("geoid_model");
+        auto objectContext2(formatter->MakeObjectContext(nullptr, false));
+        writer.AddObjKey("name");
+        writer.Add(model->nameStr());
+        model->formatID(formatter);
+    }
+
     ObjectUsage::baseExportToJSON(formatter);
 }
 //! @endcond
@@ -2572,7 +2589,8 @@ void VerticalCRS::addLinearUnitConvert(
  * cs::VerticalCS.
  *
  * @param properties See \ref general_properties.
- * At minimum the name should be defined.
+ * At minimum the name should be defined. The GEOID_MODEL property can be set
+ * to a TransformationNNPtr object.
  * @param datumIn The datum of the CRS.
  * @param csIn a VerticalCS.
  * @return new VerticalCRS.
@@ -2592,7 +2610,8 @@ VerticalCRS::create(const util::PropertyMap &properties,
  * One and only one of datum or datumEnsemble should be set to a non-null value.
  *
  * @param properties See \ref general_properties.
- * At minimum the name should be defined.
+ * At minimum the name should be defined. The GEOID_MODEL property can be set
+ * to a TransformationNNPtr object.
  * @param datumIn The datum of the CRS, or nullptr
  * @param datumEnsembleIn The datum ensemble of the CRS, or nullptr.
  * @param csIn a VerticalCS.
@@ -2607,6 +2626,14 @@ VerticalCRS::create(const util::PropertyMap &properties,
                                                       csIn));
     crs->assignSelf(crs);
     crs->setProperties(properties);
+    const auto geoidModelPtr = properties.get("GEOID_MODEL");
+    if (geoidModelPtr) {
+        auto transf = util::nn_dynamic_pointer_cast<operation::Transformation>(
+            *geoidModelPtr);
+        if (transf) {
+            crs->d->geoidModel.emplace_back(NN_NO_CHECK(transf));
+        }
+    }
     return crs;
 }
 
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/iso19111/io.cpp b/src/iso19111/io.cpp
index 50ad5a87..645bec0b 100644
--- a/src/iso19111/io.cpp
+++ b/src/iso19111/io.cpp
@@ -89,7 +89,7 @@ static const std::string emptyString{};
 
 // If changing that value, change it in data/projjson.schema.json as well
 #define PROJJSON_CURRENT_VERSION                                               \
-    "https://proj.org/schemas/v0.1/projjson.schema.json"
+    "https://proj.org/schemas/v0.2/projjson.schema.json"
 //! @endcond
 
 #if 0
@@ -3873,28 +3873,110 @@ CRSNNPtr WKTParser::Private::buildVerticalCRS(const WKTNodeNNPtr &node) {
         ThrowNotExpectedCSType("vertical");
     }
 
+    auto &props = buildProperties(node);
+
+    // Deal with Lidar WKT1 VertCRS that embeds geoid model in CRS name,
+    // following conventions from
+    // https://pubs.usgs.gov/tm/11b4/pdf/tm11-B4.pdf
+    // page 9
+    if (ci_equal(nodeValue, WKTConstants::VERT_CS)) {
+        std::string name;
+        if (props.getStringValue(IdentifiedObject::NAME_KEY, name)) {
+            std::string geoidName;
+            for (const char *prefix :
+                 {"NAVD88 - ", "NAVD88 via ", "NAVD88 height - ",
+                  "NAVD88 height (ftUS) - "}) {
+                if (starts_with(name, prefix)) {
+                    geoidName = name.substr(strlen(prefix));
+                    auto pos = geoidName.find_first_of(" (");
+                    if (pos != std::string::npos) {
+                        geoidName.resize(pos);
+                    }
+                    break;
+                }
+            }
+            if (!geoidName.empty()) {
+                const auto &axis = verticalCS->axisList()[0];
+                const auto &dir = axis->direction();
+                if (dir == cs::AxisDirection::UP) {
+                    if (axis->unit() == common::UnitOfMeasure::METRE) {
+                        props.set(IdentifiedObject::NAME_KEY, "NAVD88 height");
+                        props.set(Identifier::CODE_KEY, 5703);
+                        props.set(Identifier::CODESPACE_KEY, Identifier::EPSG);
+                    } else if (axis->unit().name() == "US survey foot") {
+                        props.set(IdentifiedObject::NAME_KEY,
+                                  "NAVD88 height (ftUS)");
+                        props.set(Identifier::CODE_KEY, 6360);
+                        props.set(Identifier::CODESPACE_KEY, Identifier::EPSG);
+                    }
+                }
+                PropertyMap propsModel;
+                propsModel.set(IdentifiedObject::NAME_KEY, toupper(geoidName));
+                PropertyMap propsDatum;
+                propsDatum.set(IdentifiedObject::NAME_KEY,
+                               "North American Vertical Datum 1988");
+                propsDatum.set(Identifier::CODE_KEY, 5103);
+                propsDatum.set(Identifier::CODESPACE_KEY, Identifier::EPSG);
+                datum =
+                    VerticalReferenceFrame::create(propsDatum).as_nullable();
+                const auto dummyCRS =
+                    VerticalCRS::create(PropertyMap(), datum, datumEnsemble,
+                                        NN_NO_CHECK(verticalCS));
+                const auto model(Transformation::create(
+                    propsModel, dummyCRS, dummyCRS, nullptr,
+                    OperationMethod::create(
+                        PropertyMap(), std::vector<OperationParameterNNPtr>()),
+                    {}, {}));
+                props.set("GEOID_MODEL", model);
+            }
+        }
+    }
+
+    auto &geoidModelNode = nodeP->lookForChild(WKTConstants::GEOIDMODEL);
+    if (!isNull(geoidModelNode)) {
+        auto &propsModel = buildProperties(geoidModelNode);
+        const auto dummyCRS = VerticalCRS::create(
+            PropertyMap(), datum, datumEnsemble, NN_NO_CHECK(verticalCS));
+        const auto model(Transformation::create(
+            propsModel, dummyCRS, dummyCRS, nullptr,
+            OperationMethod::create(PropertyMap(),
+                                    std::vector<OperationParameterNNPtr>()),
+            {}, {}));
+        props.set("GEOID_MODEL", model);
+    }
+
     auto crs = nn_static_pointer_cast<CRS>(VerticalCRS::create(
-        buildProperties(node), datum, datumEnsemble, NN_NO_CHECK(verticalCS)));
+        props, datum, datumEnsemble, NN_NO_CHECK(verticalCS)));
 
     if (!isNull(datumNode)) {
         auto &extensionNode = datumNode->lookForChild(WKTConstants::EXTENSION);
         const auto &extensionChildren = extensionNode->GP()->children();
         if (extensionChildren.size() == 2) {
             if (ci_equal(stripQuotes(extensionChildren[0]), "PROJ4_GRIDS")) {
-                std::string transformationName(crs->nameStr());
-                if (!ends_with(transformationName, " height")) {
-                    transformationName += " height";
+                const auto gridName(stripQuotes(extensionChildren[1]));
+                // This is the expansion of EPSG:5703 by old GDAL versions.
+                // See
+                // https://trac.osgeo.org/metacrs/changeset?reponame=&new=2281%40geotiff%2Ftrunk%2Flibgeotiff%2Fcsv%2Fvertcs.override.csv&old=1893%40geotiff%2Ftrunk%2Flibgeotiff%2Fcsv%2Fvertcs.override.csv
+                // It is unlikely that the user really explictly wants this.
+                if (gridName != "g2003conus.gtx,g2003alaska.gtx,"
+                                "g2003h01.gtx,g2003p01.gtx" &&
+                    gridName != "g2012a_conus.gtx,g2012a_alaska.gtx,"
+                                "g2012a_guam.gtx,g2012a_hawaii.gtx,"
+                                "g2012a_puertorico.gtx,g2012a_samoa.gtx") {
+                    std::string transformationName(crs->nameStr());
+                    if (!ends_with(transformationName, " height")) {
+                        transformationName += " height";
+                    }
+                    transformationName += " to WGS84 ellipsoidal height";
+                    auto transformation = Transformation::
+                        createGravityRelatedHeightToGeographic3D(
+                            PropertyMap().set(IdentifiedObject::NAME_KEY,
+                                              transformationName),
+                            crs, GeographicCRS::EPSG_4979, nullptr, gridName,
+                            std::vector<PositionalAccuracyNNPtr>());
+                    return nn_static_pointer_cast<CRS>(BoundCRS::create(
+                        crs, GeographicCRS::EPSG_4979, transformation));
                 }
-                transformationName += " to WGS84 ellipsoidal height";
-                auto transformation =
-                    Transformation::createGravityRelatedHeightToGeographic3D(
-                        PropertyMap().set(IdentifiedObject::NAME_KEY,
-                                          transformationName),
-                        crs, GeographicCRS::EPSG_4979, nullptr,
-                        stripQuotes(extensionChildren[1]),
-                        std::vector<PositionalAccuracyNNPtr>());
-                return nn_static_pointer_cast<CRS>(BoundCRS::create(
-                    crs, GeographicCRS::EPSG_4979, transformation));
             }
         }
     }
@@ -5040,7 +5122,22 @@ VerticalCRSNNPtr JSONParser::buildVerticalCRS(const json &j) {
     if (!verticalCS) {
         throw ParsingException("expected a vertical CS");
     }
-    return VerticalCRS::create(buildProperties(j), datum, datumEnsemble,
+
+    auto props = buildProperties(j);
+    if (j.contains("geoid_model")) {
+        auto geoidModelJ = getObject(j, "geoid_model");
+        auto propsModel = buildProperties(geoidModelJ);
+        const auto dummyCRS = VerticalCRS::create(
+            PropertyMap(), datum, datumEnsemble, NN_NO_CHECK(verticalCS));
+        const auto model(Transformation::create(
+            propsModel, dummyCRS, dummyCRS, nullptr,
+            OperationMethod::create(PropertyMap(),
+                                    std::vector<OperationParameterNNPtr>()),
+            {}, {}));
+        props.set("GEOID_MODEL", model);
+    }
+
+    return VerticalCRS::create(props, datum, datumEnsemble,
                                NN_NO_CHECK(verticalCS));
 }
 
diff --git a/src/iso19111/static.cpp b/src/iso19111/static.cpp
index 824047f0..0ed08c95 100644
--- a/src/iso19111/static.cpp
+++ b/src/iso19111/static.cpp
@@ -274,6 +274,7 @@ DEFINE_WKT_CONSTANT(BASEENGCRS);
 DEFINE_WKT_CONSTANT(BASEPARAMCRS);
 DEFINE_WKT_CONSTANT(BASETIMECRS);
 DEFINE_WKT_CONSTANT(VERSION);
+DEFINE_WKT_CONSTANT(GEOIDMODEL);
 
 DEFINE_WKT_CONSTANT(GEODETICCRS);
 DEFINE_WKT_CONSTANT(GEODETICDATUM);
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,