aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--data/sql/commit.sql16
-rw-r--r--data/sql/customizations.sql17
-rw-r--r--data/sql/proj_db_table_defs.sql17
-rw-r--r--include/proj/io.hpp4
-rw-r--r--scripts/reference_exported_symbols.txt1
-rw-r--r--src/iso19111/c_api.cpp76
-rw-r--r--src/iso19111/coordinateoperation.cpp143
-rw-r--r--src/iso19111/factory.cpp22
-rw-r--r--src/proj_experimental.h13
-rw-r--r--test/unit/test_c_api.cpp107
-rw-r--r--test/unit/test_operation.cpp82
11 files changed, 489 insertions, 9 deletions
diff --git a/data/sql/commit.sql b/data/sql/commit.sql
index eb49828a..a708df0f 100644
--- a/data/sql/commit.sql
+++ b/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;
diff --git a/data/sql/customizations.sql b/data/sql/customizations.sql
index 99891dc1..b0804ebd 100644
--- a/data/sql/customizations.sql
+++ b/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;
diff --git a/data/sql/proj_db_table_defs.sql b/data/sql/proj_db_table_defs.sql
index 8a02acf4..a97e75a8 100644
--- a/data/sql/proj_db_table_defs.sql
+++ b/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',
diff --git a/include/proj/io.hpp b/include/proj/io.hpp
index f8ee1c97..17f0fea5 100644
--- a/include/proj/io.hpp
+++ b/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:
diff --git a/scripts/reference_exported_symbols.txt b/scripts/reference_exported_symbols.txt
index cb568566..67d4837f 100644
--- a/scripts/reference_exported_symbols.txt
+++ b/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
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,
diff --git a/test/unit/test_c_api.cpp b/test/unit/test_c_api.cpp
index 5e6d1924..cb59563e 100644
--- a/test/unit/test_c_api.cpp
+++ b/test/unit/test_c_api.cpp
@@ -4213,4 +4213,111 @@ TEST_F(
EXPECT_NEAR(outcoord.xyzt.z, -32.5823, 1e-3);
}
+// ---------------------------------------------------------------------------
+
+TEST_F(CApi, proj_create_vertical_crs_ex) {
+
+ // NAD83(2011) / UTM zone 11N
+ auto horiz_crs = proj_create_from_database(m_ctxt, "EPSG", "6340",
+ PJ_CATEGORY_CRS, false, nullptr);
+ ObjectKeeper keeper_horiz_crs(horiz_crs);
+ ASSERT_NE(horiz_crs, 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);
+ ObjectKeeper keeper_vert_crs(vert_crs);
+ ASSERT_NE(vert_crs, nullptr);
+
+ auto compound =
+ proj_create_compound_crs(m_ctxt, "Compound", horiz_crs, vert_crs);
+ ObjectKeeper keeper_compound(compound);
+ ASSERT_NE(compound, nullptr);
+
+ // NAD83(2011) 3D
+ PJ *geog_crs = proj_create(m_ctxt, "EPSG:6319");
+ ObjectKeeper keeper_geog_crs(geog_crs);
+ ASSERT_NE(geog_crs, nullptr);
+
+ auto P = proj_create_crs_to_crs_from_pj(m_ctxt, compound, geog_crs, nullptr,
+ nullptr);
+ ObjectKeeper keeper_P(P);
+ ASSERT_NE(P, nullptr);
+
+ auto name = proj_get_name(P);
+ ASSERT_TRUE(name != nullptr);
+ EXPECT_EQ(name,
+ std::string("Inverse of UTM zone 11N + "
+ "Transformation from myVertCRS (ftUS) to myVertCRS + "
+ "Transformation from myVertCRS to NAD83(2011)"));
+
+ auto proj_5 = proj_as_proj_string(m_ctxt, P, PJ_PROJ_5, nullptr);
+ ASSERT_NE(proj_5, nullptr);
+ EXPECT_EQ(std::string(proj_5),
+ "+proj=pipeline "
+ "+step +inv +proj=utm +zone=11 +ellps=GRS80 "
+ "+step +proj=unitconvert +z_in=us-ft +z_out=m "
+ "+step +proj=vgridshift +grids=@foo.gtx +multiplier=1 "
+ "+step +proj=unitconvert +xy_in=rad +xy_out=deg "
+ "+step +proj=axisswap +order=2,1");
+}
+
+// ---------------------------------------------------------------------------
+
+TEST_F(CApi, proj_create_vertical_crs_ex_with_geog_crs) {
+
+ // NAD83(2011) / UTM zone 11N
+ auto horiz_crs = proj_create_from_database(m_ctxt, "EPSG", "6340",
+ PJ_CATEGORY_CRS, false, nullptr);
+ ObjectKeeper keeper_horiz_crs(horiz_crs);
+ ASSERT_NE(horiz_crs, nullptr);
+
+ // WGS84
+ PJ *wgs84 = proj_create(m_ctxt, "EPSG:4979");
+ ObjectKeeper keeper_wgs84(wgs84);
+ ASSERT_NE(wgs84, nullptr);
+
+ auto vert_crs = proj_create_vertical_crs_ex(
+ m_ctxt, "myVertCRS", "myVertDatum", nullptr, nullptr, "US survey foot",
+ 0.304800609601219, "PROJ @foo.gtx", nullptr, nullptr, wgs84, nullptr);
+ ObjectKeeper keeper_vert_crs(vert_crs);
+ ASSERT_NE(vert_crs, nullptr);
+
+ auto compound =
+ proj_create_compound_crs(m_ctxt, "Compound", horiz_crs, vert_crs);
+ ObjectKeeper keeper_compound(compound);
+ ASSERT_NE(compound, nullptr);
+
+ // NAD83(2011) 3D
+ PJ *geog_crs = proj_create(m_ctxt, "EPSG:6319");
+ ObjectKeeper keeper_geog_crs(geog_crs);
+ ASSERT_NE(geog_crs, nullptr);
+
+ auto P = proj_create_crs_to_crs_from_pj(m_ctxt, compound, geog_crs, nullptr,
+ nullptr);
+ ObjectKeeper keeper_P(P);
+ ASSERT_NE(P, nullptr);
+
+ auto name = proj_get_name(P);
+ ASSERT_TRUE(name != nullptr);
+ EXPECT_EQ(
+ name,
+ std::string("Inverse of UTM zone 11N + "
+ "Ballpark geographic offset from NAD83(2011) to WGS 84 + "
+ "Transformation from myVertCRS to myVertCRS (metre) + "
+ "Transformation from myVertCRS (metre) to WGS 84 + "
+ "Ballpark geographic offset from WGS 84 to NAD83(2011)"));
+
+ auto proj_5 = proj_as_proj_string(m_ctxt, P, PJ_PROJ_5, nullptr);
+ ASSERT_NE(proj_5, nullptr);
+ EXPECT_EQ(std::string(proj_5),
+ "+proj=pipeline "
+ "+step +inv +proj=utm +zone=11 +ellps=GRS80 "
+ "+step +proj=unitconvert +z_in=us-ft +z_out=m "
+ "+step +proj=vgridshift +grids=@foo.gtx +multiplier=1 "
+ "+step +proj=unitconvert +xy_in=rad +xy_out=deg "
+ "+step +proj=axisswap +order=2,1");
+}
+
} // namespace
diff --git a/test/unit/test_operation.cpp b/test/unit/test_operation.cpp
index 7d457dab..5e04bf23 100644
--- a/test/unit/test_operation.cpp
+++ b/test/unit/test_operation.cpp
@@ -7580,6 +7580,88 @@ TEST(
// ---------------------------------------------------------------------------
+TEST(operation, compoundCRS_of_vertCRS_with_geoid_model_to_geogCRS) {
+ auto authFactory =
+ AuthorityFactory::create(DatabaseContext::create(), "EPSG");
+ auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0);
+ ctxt->setSpatialCriterion(
+ CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION);
+ ctxt->setGridAvailabilityUse(
+ CoordinateOperationContext::GridAvailabilityUse::
+ IGNORE_GRID_AVAILABILITY);
+ auto wkt =
+ "COMPOUNDCRS[\"NAD83 / Pennsylvania South + NAVD88 height\",\n"
+ " PROJCRS[\"NAD83 / Pennsylvania South\",\n"
+ " BASEGEOGCRS[\"NAD83\",\n"
+ " DATUM[\"North American Datum 1983\",\n"
+ " ELLIPSOID[\"GRS 1980\",6378137,298.257222101,\n"
+ " LENGTHUNIT[\"metre\",1]]],\n"
+ " PRIMEM[\"Greenwich\",0,\n"
+ " ANGLEUNIT[\"degree\",0.0174532925199433]]],\n"
+ " CONVERSION[\"SPCS83 Pennsylvania South zone (meters)\",\n"
+ " METHOD[\"Lambert Conic Conformal (2SP)\",\n"
+ " ID[\"EPSG\",9802]],\n"
+ " PARAMETER[\"Latitude of false origin\",39.3333333333333,\n"
+ " ANGLEUNIT[\"degree\",0.0174532925199433],\n"
+ " ID[\"EPSG\",8821]],\n"
+ " PARAMETER[\"Longitude of false origin\",-77.75,\n"
+ " ANGLEUNIT[\"degree\",0.0174532925199433],\n"
+ " ID[\"EPSG\",8822]],\n"
+ " PARAMETER[\"Latitude of 1st standard "
+ "parallel\",40.9666666666667,\n"
+ " ANGLEUNIT[\"degree\",0.0174532925199433],\n"
+ " ID[\"EPSG\",8823]],\n"
+ " PARAMETER[\"Latitude of 2nd standard "
+ "parallel\",39.9333333333333,\n"
+ " ANGLEUNIT[\"degree\",0.0174532925199433],\n"
+ " ID[\"EPSG\",8824]],\n"
+ " PARAMETER[\"Easting at false origin\",600000,\n"
+ " LENGTHUNIT[\"metre\",1],\n"
+ " ID[\"EPSG\",8826]],\n"
+ " PARAMETER[\"Northing at false origin\",0,\n"
+ " LENGTHUNIT[\"metre\",1],\n"
+ " ID[\"EPSG\",8827]]],\n"
+ " CS[Cartesian,2],\n"
+ " AXIS[\"easting (X)\",east,\n"
+ " ORDER[1],\n"
+ " LENGTHUNIT[\"metre\",1]],\n"
+ " AXIS[\"northing (Y)\",north,\n"
+ " ORDER[2],\n"
+ " LENGTHUNIT[\"metre\",1]]],\n"
+ " VERTCRS[\"NAVD88 height\",\n"
+ " VDATUM[\"North American Vertical Datum 1988\"],\n"
+ " CS[vertical,1],\n"
+ " AXIS[\"gravity-related height (H)\",up,\n"
+ " LENGTHUNIT[\"metre\",1]],\n"
+ " GEOIDMODEL[\"GEOID12B\"]]]";
+ auto srcObj =
+ createFromUserInput(wkt, authFactory->databaseContext(), false);
+ auto src = nn_dynamic_pointer_cast<CRS>(srcObj);
+ ASSERT_TRUE(src != nullptr);
+ auto dst = authFactory->createCoordinateReferenceSystem("4269"); // NAD83
+
+ auto list = CoordinateOperationFactory::create()->createOperations(
+ NN_NO_CHECK(src), dst, ctxt);
+ ASSERT_TRUE(!list.empty());
+ EXPECT_EQ(list[0]->nameStr(),
+ "Inverse of SPCS83 Pennsylvania South zone (meters) + "
+ "Ballpark geographic offset from NAD83 to NAD83(2011) + "
+ "Inverse of NAD83(2011) to NAVD88 height (1) + "
+ "Ballpark geographic offset from NAD83(2011) to NAD83");
+ auto op_proj =
+ list[0]->exportToPROJString(PROJStringFormatter::create().get());
+ EXPECT_EQ(op_proj,
+ "+proj=pipeline "
+ "+step +inv +proj=lcc +lat_0=39.3333333333333 +lon_0=-77.75 "
+ "+lat_1=40.9666666666667 +lat_2=39.9333333333333 +x_0=600000 "
+ "+y_0=0 +ellps=GRS80 "
+ "+step +proj=vgridshift +grids=g2012bu0.gtx +multiplier=1 "
+ "+step +proj=unitconvert +xy_in=rad +xy_out=deg "
+ "+step +proj=axisswap +order=2,1");
+}
+
+// ---------------------------------------------------------------------------
+
TEST(operation, compoundCRS_from_WKT2_to_geogCRS_3D_context) {
auto authFactory =
AuthorityFactory::create(DatabaseContext::create(), "EPSG");
generated by cgit v1.2.3 (git 2.25.1) at 2025-12-15 21:46:57 +0000