createOperations(): fix bug involving geoidmodel and non-metre vertical unit

2 files changed, 463 insertions, 5 deletions

diff --git a/src/iso19111/operation/coordinateoperationfactory.cpp b/src/iso19111/operation/coordinateoperationfactory.cpp
index 902efb4c..51262ff6 100644
--- a/src/iso19111/operation/coordinateoperationfactory.cpp
+++ b/src/iso19111/operation/coordinateoperationfactory.cpp
@@ -3285,21 +3285,62 @@ CoordinateOperationFactory::Private::createOperationsGeogToVertFromGeoid(
 
     ENTER_FUNCTION();
 
-    const auto useTransf = [&targetCRS, &context,
+    const auto useTransf = [&sourceCRS, &targetCRS, &context,
                             vertDst](const CoordinateOperationNNPtr &op) {
+
+        // If the source geographic CRS has a non-metre vertical unit, we need
+        // to create an intermediate and operation to do the vertical unit
+        // conversion from that vertical unit to the one of the geographic CRS
+        // of the source of the operation
+        const auto geogCRS =
+            dynamic_cast<const crs::GeographicCRS *>(sourceCRS.get());
+        assert(geogCRS);
+        const auto &srcAxisList = geogCRS->coordinateSystem()->axisList();
+        CoordinateOperationPtr opPtr;
+        const auto opSourceCRSGeog =
+            dynamic_cast<const crs::GeographicCRS *>(op->sourceCRS().get());
+        // I assume opSourceCRSGeog should always be null in practice...
+        if (opSourceCRSGeog && srcAxisList.size() == 3 &&
+            srcAxisList[2]->unit().conversionToSI() != 1) {
+            const auto &authFactory = context.context->getAuthorityFactory();
+            const auto dbContext =
+                authFactory ? authFactory->databaseContext().as_nullable()
+                            : nullptr;
+            auto tmpCRSWithSrcZ =
+                opSourceCRSGeog->demoteTo2D(std::string(), dbContext)
+                    ->promoteTo3D(std::string(), dbContext, srcAxisList[2]);
+
+            std::vector<CoordinateOperationNNPtr> opsUnitConvert;
+            createOperationsGeogToGeog(
+                opsUnitConvert, tmpCRSWithSrcZ, NN_NO_CHECK(op->sourceCRS()),
+                context,
+                dynamic_cast<const crs::GeographicCRS *>(tmpCRSWithSrcZ.get()),
+                opSourceCRSGeog);
+            assert(opsUnitConvert.size() == 1);
+            opPtr = opsUnitConvert.front().as_nullable();
+        }
+
+        std::vector<CoordinateOperationNNPtr> ops;
+        if (opPtr)
+            ops.emplace_back(NN_NO_CHECK(opPtr));
+        ops.emplace_back(op);
+
         const auto targetOp =
             dynamic_cast<const crs::VerticalCRS *>(op->targetCRS().get());
         assert(targetOp);
         if (targetOp->_isEquivalentTo(
                 vertDst, util::IComparable::Criterion::EQUIVALENT)) {
-            return op;
+            auto ret = ConcatenatedOperation::createComputeMetadata(
+                ops, disallowEmptyIntersection);
+            return ret;
         }
         std::vector<CoordinateOperationNNPtr> tmp;
         createOperationsVertToVert(NN_NO_CHECK(op->targetCRS()), targetCRS,
                                    context, targetOp, vertDst, tmp);
         assert(!tmp.empty());
+        ops.emplace_back(tmp.front());
         auto ret = ConcatenatedOperation::createComputeMetadata(
-            {op, tmp.front()}, disallowEmptyIntersection);
+            ops, disallowEmptyIntersection);
         return ret;
     };
 
@@ -3331,10 +3372,16 @@ CoordinateOperationFactory::Private::createOperationsGeogToVertFromGeoid(
         };
 
         const auto &axis = vertDst->coordinateSystem()->axisList()[0];
+        const auto &authFactory = context.context->getAuthorityFactory();
+        const auto dbContext =
+            authFactory ? authFactory->databaseContext().as_nullable()
+                        : nullptr;
+
         const auto geogSrcCRS =
             dynamic_cast<crs::GeographicCRS *>(model->interpolationCRS().get())
                 ? NN_NO_CHECK(model->interpolationCRS())
-                : sourceCRS;
+                : sourceCRS->demoteTo2D(std::string(), dbContext)
+                      ->promoteTo3D(std::string(), dbContext);
         const auto vertCRSMetre =
             axis->unit() == common::UnitOfMeasure::METRE &&
                     axis->direction() == cs::AxisDirection::UP
@@ -3356,7 +3403,6 @@ CoordinateOperationFactory::Private::createOperationsGeogToVertFromGeoid(
         std::vector<metadata::PositionalAccuracyNNPtr> accuracies;
         const auto &modelAccuracies = model->coordinateOperationAccuracies();
         if (modelAccuracies.empty()) {
-            const auto &authFactory = context.context->getAuthorityFactory();
             if (authFactory) {
                 const auto transformationsForGrid =
                     io::DatabaseContext::getTransformationsForGridName(
diff --git a/test/unit/test_operationfactory.cpp b/test/unit/test_operationfactory.cpp
index 2f0e97b5..82bdd2fd 100644
--- a/test/unit/test_operationfactory.cpp
+++ b/test/unit/test_operationfactory.cpp
@@ -5898,6 +5898,418 @@ TEST(operation, compoundCRS_to_proj_string_with_non_metre_height) {
 
 // ---------------------------------------------------------------------------
 
+TEST(operation, compoundCRS_to_PROJJSON_with_non_metre_height) {
+    auto srcPROJJSON =
+        "{\n"
+        "  \"$schema\": "
+        "\"https://proj.org/schemas/v0.2/projjson.schema.json\",\n"
+        "  \"type\": \"CompoundCRS\",\n"
+        "  \"name\": \"Compound CRS NAD83(2011) / Nebraska (ftUS) + North "
+        "American Vertical Datum 1988 + PROJ us_noaa_g2012bu0.tif\",\n"
+        "  \"components\": [\n"
+        "    {\n"
+        "      \"type\": \"ProjectedCRS\",\n"
+        "      \"name\": \"NAD83(2011) / Nebraska (ftUS)\",\n"
+        "      \"base_crs\": {\n"
+        "        \"name\": \"NAD83(2011)\",\n"
+        "        \"datum\": {\n"
+        "          \"type\": \"GeodeticReferenceFrame\",\n"
+        "          \"name\": \"NAD83 (National Spatial Reference System "
+        "2011)\",\n"
+        "          \"ellipsoid\": {\n"
+        "            \"name\": \"GRS 1980\",\n"
+        "            \"semi_major_axis\": 6378137,\n"
+        "            \"inverse_flattening\": 298.257222101\n"
+        "          }\n"
+        "        },\n"
+        "        \"coordinate_system\": {\n"
+        "          \"subtype\": \"ellipsoidal\",\n"
+        "          \"axis\": [\n"
+        "            {\n"
+        "              \"name\": \"Geodetic latitude\",\n"
+        "              \"abbreviation\": \"Lat\",\n"
+        "              \"direction\": \"north\",\n"
+        "              \"unit\": \"degree\"\n"
+        "            },\n"
+        "            {\n"
+        "              \"name\": \"Geodetic longitude\",\n"
+        "              \"abbreviation\": \"Lon\",\n"
+        "              \"direction\": \"east\",\n"
+        "              \"unit\": \"degree\"\n"
+        "            }\n"
+        "          ]\n"
+        "        },\n"
+        "        \"id\": {\n"
+        "          \"authority\": \"EPSG\",\n"
+        "          \"code\": 6318\n"
+        "        }\n"
+        "      },\n"
+        "      \"conversion\": {\n"
+        "        \"name\": \"SPCS83 Nebraska zone (US Survey feet)\",\n"
+        "        \"method\": {\n"
+        "          \"name\": \"Lambert Conic Conformal (2SP)\",\n"
+        "          \"id\": {\n"
+        "            \"authority\": \"EPSG\",\n"
+        "            \"code\": 9802\n"
+        "          }\n"
+        "        },\n"
+        "        \"parameters\": [\n"
+        "          {\n"
+        "            \"name\": \"Latitude of false origin\",\n"
+        "            \"value\": 39.8333333333333,\n"
+        "            \"unit\": \"degree\",\n"
+        "            \"id\": {\n"
+        "              \"authority\": \"EPSG\",\n"
+        "              \"code\": 8821\n"
+        "            }\n"
+        "          },\n"
+        "          {\n"
+        "            \"name\": \"Longitude of false origin\",\n"
+        "            \"value\": -100,\n"
+        "            \"unit\": \"degree\",\n"
+        "            \"id\": {\n"
+        "              \"authority\": \"EPSG\",\n"
+        "              \"code\": 8822\n"
+        "            }\n"
+        "          },\n"
+        "          {\n"
+        "            \"name\": \"Latitude of 1st standard parallel\",\n"
+        "            \"value\": 43,\n"
+        "            \"unit\": \"degree\",\n"
+        "            \"id\": {\n"
+        "              \"authority\": \"EPSG\",\n"
+        "              \"code\": 8823\n"
+        "            }\n"
+        "          },\n"
+        "          {\n"
+        "            \"name\": \"Latitude of 2nd standard parallel\",\n"
+        "            \"value\": 40,\n"
+        "            \"unit\": \"degree\",\n"
+        "            \"id\": {\n"
+        "              \"authority\": \"EPSG\",\n"
+        "              \"code\": 8824\n"
+        "            }\n"
+        "          },\n"
+        "          {\n"
+        "            \"name\": \"Easting at false origin\",\n"
+        "            \"value\": 1640416.6667,\n"
+        "            \"unit\": {\n"
+        "              \"type\": \"LinearUnit\",\n"
+        "              \"name\": \"US survey foot\",\n"
+        "              \"conversion_factor\": 0.304800609601219\n"
+        "            },\n"
+        "            \"id\": {\n"
+        "              \"authority\": \"EPSG\",\n"
+        "              \"code\": 8826\n"
+        "            }\n"
+        "          },\n"
+        "          {\n"
+        "            \"name\": \"Northing at false origin\",\n"
+        "            \"value\": 0,\n"
+        "            \"unit\": {\n"
+        "              \"type\": \"LinearUnit\",\n"
+        "              \"name\": \"US survey foot\",\n"
+        "              \"conversion_factor\": 0.304800609601219\n"
+        "            },\n"
+        "            \"id\": {\n"
+        "              \"authority\": \"EPSG\",\n"
+        "              \"code\": 8827\n"
+        "            }\n"
+        "          }\n"
+        "        ],\n"
+        "        \"id\": {\n"
+        "          \"authority\": \"EPSG\",\n"
+        "          \"code\": 15396\n"
+        "        }\n"
+        "      },\n"
+        "      \"coordinate_system\": {\n"
+        "        \"subtype\": \"Cartesian\",\n"
+        "        \"axis\": [\n"
+        "          {\n"
+        "            \"name\": \"Easting\",\n"
+        "            \"abbreviation\": \"X\",\n"
+        "            \"direction\": \"east\",\n"
+        "            \"unit\": {\n"
+        "              \"type\": \"LinearUnit\",\n"
+        "              \"name\": \"US survey foot\",\n"
+        "              \"conversion_factor\": 0.304800609601219\n"
+        "            }\n"
+        "          },\n"
+        "          {\n"
+        "            \"name\": \"Northing\",\n"
+        "            \"abbreviation\": \"Y\",\n"
+        "            \"direction\": \"north\",\n"
+        "            \"unit\": {\n"
+        "              \"type\": \"LinearUnit\",\n"
+        "              \"name\": \"US survey foot\",\n"
+        "              \"conversion_factor\": 0.304800609601219\n"
+        "            }\n"
+        "          }\n"
+        "        ]\n"
+        "      }\n"
+        "    },\n"
+        "    {\n"
+        "      \"type\": \"VerticalCRS\",\n"
+        "      \"name\": \"North American Vertical Datum 1988 + PROJ "
+        "us_noaa_g2012bu0.tif\",\n"
+        "      \"datum\": {\n"
+        "        \"type\": \"VerticalReferenceFrame\",\n"
+        "        \"name\": \"North American Vertical Datum 1988\",\n"
+        "        \"id\": {\n"
+        "          \"authority\": \"EPSG\",\n"
+        "          \"code\": 5103\n"
+        "        }\n"
+        "      },\n"
+        "      \"coordinate_system\": {\n"
+        "        \"subtype\": \"vertical\",\n"
+        "        \"axis\": [\n"
+        "          {\n"
+        "            \"name\": \"Gravity-related height\",\n"
+        "            \"abbreviation\": \"H\",\n"
+        "            \"direction\": \"up\",\n"
+        "            \"unit\": {\n"
+        "              \"type\": \"LinearUnit\",\n"
+        "              \"name\": \"US survey foot\",\n"
+        "              \"conversion_factor\": 0.304800609601219\n"
+        "            }\n"
+        "          }\n"
+        "        ]\n"
+        "      },\n"
+        "      \"geoid_model\": {\n"
+        "        \"name\": \"PROJ us_noaa_g2012bu0.tif\",\n"
+        "        \"interpolation_crs\": {\n"
+        "          \"type\": \"GeographicCRS\",\n"
+        "          \"name\": \"NAD83(2011)\",\n"
+        "          \"datum\": {\n"
+        "            \"type\": \"GeodeticReferenceFrame\",\n"
+        "            \"name\": \"NAD83 (National Spatial Reference System "
+        "2011)\",\n"
+        "            \"ellipsoid\": {\n"
+        "              \"name\": \"GRS 1980\",\n"
+        "              \"semi_major_axis\": 6378137,\n"
+        "              \"inverse_flattening\": 298.257222101\n"
+        "            }\n"
+        "          },\n"
+        "          \"coordinate_system\": {\n"
+        "            \"subtype\": \"ellipsoidal\",\n"
+        "            \"axis\": [\n"
+        "              {\n"
+        "                \"name\": \"Geodetic latitude\",\n"
+        "                \"abbreviation\": \"Lat\",\n"
+        "                \"direction\": \"north\",\n"
+        "                \"unit\": \"degree\"\n"
+        "              },\n"
+        "              {\n"
+        "                \"name\": \"Geodetic longitude\",\n"
+        "                \"abbreviation\": \"Lon\",\n"
+        "                \"direction\": \"east\",\n"
+        "                \"unit\": \"degree\"\n"
+        "              },\n"
+        "              {\n"
+        "                \"name\": \"Ellipsoidal height\",\n"
+        "                \"abbreviation\": \"h\",\n"
+        "                \"direction\": \"up\",\n"
+        "                \"unit\": \"metre\"\n"
+        "              }\n"
+        "            ]\n"
+        "          },\n"
+        "          \"id\": {\n"
+        "            \"authority\": \"EPSG\",\n"
+        "            \"code\": 6319\n"
+        "          }\n"
+        "        }\n"
+        "      }\n"
+        "    }\n"
+        "  ]\n"
+        "}";
+    auto objSrc =
+        createFromUserInput(srcPROJJSON, DatabaseContext::create(), false);
+    auto src = nn_dynamic_pointer_cast<CRS>(objSrc);
+    ASSERT_TRUE(src != nullptr);
+
+    // The untypical potentially a bit buggy thing (and what caused a bug)
+    // is the US-ft unit for the vertical axis of the base CRS ...
+    // When outputing that to WKT, and
+    // re-exporting to PROJJSON, one gets metre, which conforms more to the
+    // official definition of NAD83(2011) 3D.
+    // The vertical unit of the base CRS shouldn't matter much anyway, so this
+    // is valid.
+    auto dstPROJJSON =
+        "{\n"
+        "  \"$schema\": "
+        "\"https://proj.org/schemas/v0.2/projjson.schema.json\",\n"
+        "  \"type\": \"ProjectedCRS\",\n"
+        "  \"name\": \"Projected CRS NAD83(2011) / UTM zone 14N with "
+        "ellipsoidal NAD83(2011) height\",\n"
+        "  \"base_crs\": {\n"
+        "    \"name\": \"NAD83(2011)\",\n"
+        "    \"datum\": {\n"
+        "      \"type\": \"GeodeticReferenceFrame\",\n"
+        "      \"name\": \"NAD83 (National Spatial Reference System 2011)\",\n"
+        "      \"ellipsoid\": {\n"
+        "        \"name\": \"GRS 1980\",\n"
+        "        \"semi_major_axis\": 6378137,\n"
+        "        \"inverse_flattening\": 298.257222101\n"
+        "      },\n"
+        "      \"id\": {\n"
+        "        \"authority\": \"EPSG\",\n"
+        "        \"code\": 1116\n"
+        "      }\n"
+        "    },\n"
+        "    \"coordinate_system\": {\n"
+        "      \"subtype\": \"ellipsoidal\",\n"
+        "      \"axis\": [\n"
+        "        {\n"
+        "          \"name\": \"Geodetic latitude\",\n"
+        "          \"abbreviation\": \"Lat\",\n"
+        "          \"direction\": \"north\",\n"
+        "          \"unit\": \"degree\"\n"
+        "        },\n"
+        "        {\n"
+        "          \"name\": \"Geodetic longitude\",\n"
+        "          \"abbreviation\": \"Lon\",\n"
+        "          \"direction\": \"east\",\n"
+        "          \"unit\": \"degree\"\n"
+        "        },\n"
+        "        {\n"
+        "          \"name\": \"Ellipsoidal height\",\n"
+        "          \"abbreviation\": \"h\",\n"
+        "          \"direction\": \"up\",\n"
+        "          \"unit\": {\n"
+        "            \"type\": \"LinearUnit\",\n"
+        "            \"name\": \"US survey foot\",\n"
+        "            \"conversion_factor\": 0.304800609601219\n"
+        "          }\n"
+        "        }\n"
+        "      ]\n"
+        "    }\n"
+        "  },\n"
+        "  \"conversion\": {\n"
+        "    \"name\": \"UTM zone 14N\",\n"
+        "    \"method\": {\n"
+        "      \"name\": \"Transverse Mercator\",\n"
+        "      \"id\": {\n"
+        "        \"authority\": \"EPSG\",\n"
+        "        \"code\": 9807\n"
+        "      }\n"
+        "    },\n"
+        "    \"parameters\": [\n"
+        "      {\n"
+        "        \"name\": \"Latitude of natural origin\",\n"
+        "        \"value\": 0,\n"
+        "        \"unit\": \"degree\",\n"
+        "        \"id\": {\n"
+        "          \"authority\": \"EPSG\",\n"
+        "          \"code\": 8801\n"
+        "        }\n"
+        "      },\n"
+        "      {\n"
+        "        \"name\": \"Longitude of natural origin\",\n"
+        "        \"value\": -99,\n"
+        "        \"unit\": \"degree\",\n"
+        "        \"id\": {\n"
+        "          \"authority\": \"EPSG\",\n"
+        "          \"code\": 8802\n"
+        "        }\n"
+        "      },\n"
+        "      {\n"
+        "        \"name\": \"Scale factor at natural origin\",\n"
+        "        \"value\": 0.9996,\n"
+        "        \"unit\": \"unity\",\n"
+        "        \"id\": {\n"
+        "          \"authority\": \"EPSG\",\n"
+        "          \"code\": 8805\n"
+        "        }\n"
+        "      },\n"
+        "      {\n"
+        "        \"name\": \"False easting\",\n"
+        "        \"value\": 500000,\n"
+        "        \"unit\": \"metre\",\n"
+        "        \"id\": {\n"
+        "          \"authority\": \"EPSG\",\n"
+        "          \"code\": 8806\n"
+        "        }\n"
+        "      },\n"
+        "      {\n"
+        "        \"name\": \"False northing\",\n"
+        "        \"value\": 0,\n"
+        "        \"unit\": \"metre\",\n"
+        "        \"id\": {\n"
+        "          \"authority\": \"EPSG\",\n"
+        "          \"code\": 8807\n"
+        "        }\n"
+        "      }\n"
+        "    ],\n"
+        "    \"id\": {\n"
+        "      \"authority\": \"EPSG\",\n"
+        "      \"code\": 16014\n"
+        "    }\n"
+        "  },\n"
+        "  \"coordinate_system\": {\n"
+        "    \"subtype\": \"Cartesian\",\n"
+        "    \"axis\": [\n"
+        "      {\n"
+        "        \"name\": \"Easting\",\n"
+        "        \"abbreviation\": \"E\",\n"
+        "        \"direction\": \"east\",\n"
+        "        \"unit\": {\n"
+        "          \"type\": \"LinearUnit\",\n"
+        "          \"name\": \"US survey foot\",\n"
+        "          \"conversion_factor\": 0.304800609601219\n"
+        "        }\n"
+        "      },\n"
+        "      {\n"
+        "        \"name\": \"Northing\",\n"
+        "        \"abbreviation\": \"N\",\n"
+        "        \"direction\": \"north\",\n"
+        "        \"unit\": {\n"
+        "          \"type\": \"LinearUnit\",\n"
+        "          \"name\": \"US survey foot\",\n"
+        "          \"conversion_factor\": 0.304800609601219\n"
+        "        }\n"
+        "      },\n"
+        "      {\n"
+        "        \"name\": \"Ellipsoidal height\",\n"
+        "        \"abbreviation\": \"h\",\n"
+        "        \"direction\": \"up\",\n"
+        "        \"unit\": {\n"
+        "          \"type\": \"LinearUnit\",\n"
+        "          \"name\": \"US survey foot\",\n"
+        "          \"conversion_factor\": 0.304800609601219\n"
+        "        }\n"
+        "      }\n"
+        "    ]\n"
+        "  }\n"
+        "}";
+
+    auto objDst =
+        createFromUserInput(dstPROJJSON, DatabaseContext::create(), false);
+    auto dst = nn_dynamic_pointer_cast<CRS>(objDst);
+    ASSERT_TRUE(dst != nullptr);
+
+    auto authFactory =
+        AuthorityFactory::create(DatabaseContext::create(), "EPSG");
+    auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0);
+    ctxt->setSpatialCriterion(
+        CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION);
+    auto list = CoordinateOperationFactory::create()->createOperations(
+        NN_NO_CHECK(src), NN_NO_CHECK(dst), ctxt);
+    ASSERT_GT(list.size(), 1U);
+    // What is important to check here is the vertical unit conversion
+    EXPECT_EQ(
+        list[0]->exportToPROJString(PROJStringFormatter::create().get()),
+        "+proj=pipeline "
+        "+step +proj=unitconvert +xy_in=us-ft +xy_out=m "
+        "+step +inv +proj=lcc +lat_0=39.8333333333333 +lon_0=-100 +lat_1=43 "
+        "+lat_2=40 +x_0=500000.00001016 +y_0=0 +ellps=GRS80 "
+        "+step +proj=unitconvert +z_in=us-ft +z_out=m "
+        "+step +proj=vgridshift +grids=us_noaa_g2012bu0.tif +multiplier=1 "
+        "+step +proj=utm +zone=14 +ellps=GRS80 "
+        "+step +proj=unitconvert +xy_in=m +z_in=m +xy_out=us-ft +z_out=us-ft");
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(operation, createOperation_ossfuzz_18587) {
     auto objSrc =
         createFromUserInput("EPSG:4326", DatabaseContext::create(), false);