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

1 files changed, 412 insertions, 0 deletions

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