Add API and WKT mapping for 'nsper' to EPSG Vertical Persepective method

Relates to https://github.com/OSGeo/gdal/issues/1856

10 files changed, 253 insertions, 1 deletions

diff --git a/include/proj/coordinateoperation.hpp b/include/proj/coordinateoperation.hpp
index 2f6095f7..fb57238d 100644
--- a/include/proj/coordinateoperation.hpp
+++ b/include/proj/coordinateoperation.hpp
@@ -1311,6 +1311,15 @@ class PROJ_GCC_DLL Conversion : public SingleOperation {
                      const common::Length &falseNorthing);
 
     PROJ_DLL static ConversionNNPtr
+    createVerticalPerspective(const util::PropertyMap &properties,
+                              const common::Angle &topoOriginLat,
+                              const common::Angle &topoOriginLong,
+                              const common::Length &topoOriginHeight,
+                              const common::Length &viewPointHeight,
+                              const common::Length &falseEasting,
+                              const common::Length &falseNorthing);
+
+    PROJ_DLL static ConversionNNPtr
     createChangeVerticalUnit(const util::PropertyMap &properties,
                              const common::Scale &factor);
 
diff --git a/include/proj/internal/coordinateoperation_constants.hpp b/include/proj/internal/coordinateoperation_constants.hpp
index 9d65dab8..f1925c9b 100644
--- a/include/proj/internal/coordinateoperation_constants.hpp
+++ b/include/proj/internal/coordinateoperation_constants.hpp
@@ -490,6 +490,35 @@ static const ParamMapping *const paramsLabordeObliqueMercator[] = {
     &paramFalseNorthing,
     nullptr};
 
+static const ParamMapping paramLatTopoOrigin = {
+    EPSG_NAME_PARAMETER_LATITUDE_TOPOGRAPHIC_ORIGIN,
+    EPSG_CODE_PARAMETER_LATITUDE_TOPOGRAPHIC_ORIGIN, nullptr,
+    common::UnitOfMeasure::Type::ANGULAR, lat_0};
+
+static const ParamMapping paramLonTopoOrigin = {
+    EPSG_NAME_PARAMETER_LONGITUDE_TOPOGRAPHIC_ORIGIN,
+    EPSG_CODE_PARAMETER_LONGITUDE_TOPOGRAPHIC_ORIGIN, nullptr,
+    common::UnitOfMeasure::Type::ANGULAR, lon_0};
+
+static const ParamMapping paramHeightTopoOrigin = {
+    EPSG_NAME_PARAMETER_ELLIPSOIDAL_HEIGHT_TOPOCENTRIC_ORIGIN,
+    EPSG_CODE_PARAMETER_ELLIPSOIDAL_HEIGHT_TOPOCENTRIC_ORIGIN, nullptr,
+    common::UnitOfMeasure::Type::LINEAR,
+    nullptr}; // unsupported by PROJ right now
+
+static const ParamMapping paramViewpointHeight = {
+    EPSG_NAME_PARAMETER_VIEWPOINT_HEIGHT, EPSG_CODE_PARAMETER_VIEWPOINT_HEIGHT,
+    nullptr, common::UnitOfMeasure::Type::LINEAR, "h"};
+
+static const ParamMapping *const paramsVerticalPerspective[] = {
+    &paramLatTopoOrigin,
+    &paramLonTopoOrigin,
+    &paramHeightTopoOrigin,
+    &paramViewpointHeight,
+    &paramFalseEasting,  // PROJ addition
+    &paramFalseNorthing, // PROJ addition
+    nullptr};
+
 static const MethodMapping projectionMethodMappings[] = {
     {EPSG_NAME_METHOD_TRANSVERSE_MERCATOR, EPSG_CODE_METHOD_TRANSVERSE_MERCATOR,
      "Transverse_Mercator", "tmerc", nullptr, paramsNatOriginScaleK},
@@ -756,6 +785,9 @@ static const MethodMapping projectionMethodMappings[] = {
      EPSG_CODE_METHOD_LABORDE_OBLIQUE_MERCATOR, "Laborde_Oblique_Mercator",
      "labrd", nullptr, paramsLabordeObliqueMercator},
 
+    {EPSG_NAME_METHOD_VERTICAL_PERSPECTIVE,
+     EPSG_CODE_METHOD_VERTICAL_PERSPECTIVE, nullptr, "nsper", nullptr,
+     paramsVerticalPerspective},
 };
 
 #define METHOD_NAME_CODE(method)                                               \
@@ -791,6 +823,7 @@ static const struct MethodNameCode {
     METHOD_NAME_CODE(POLAR_STEREOGRAPHIC_VARIANT_A),
     METHOD_NAME_CODE(POLAR_STEREOGRAPHIC_VARIANT_B),
     METHOD_NAME_CODE(EQUAL_EARTH), METHOD_NAME_CODE(LABORDE_OBLIQUE_MERCATOR),
+    METHOD_NAME_CODE(VERTICAL_PERSPECTIVE),
     // Other conversions
     METHOD_NAME_CODE(CHANGE_VERTICAL_UNIT),
     METHOD_NAME_CODE(AXIS_ORDER_REVERSAL_2D),
diff --git a/scripts/create_c_api_projections.py b/scripts/create_c_api_projections.py
index e99c36a8..1d6c638d 100755
--- a/scripts/create_c_api_projections.py
+++ b/scripts/create_c_api_projections.py
@@ -157,7 +157,7 @@ for sectiondef in compounddef.iter('sectiondef'):
                     cppfile.write(", Scale(" + param[1] + ")")
 
             cppfile.write(");\n")
-            cppfile.write("    return proj_create_conversion(conv);\n")
+            cppfile.write("    return proj_create_conversion(ctx, conv);\n")
             cppfile.write("  } catch (const std::exception &e) {\n");
             cppfile.write("    proj_log_error(ctx, __FUNCTION__, e.what());\n")
             cppfile.write("  }\n")
diff --git a/scripts/reference_exported_symbols.txt b/scripts/reference_exported_symbols.txt
index c9f44ead..040ef88e 100644
--- a/scripts/reference_exported_symbols.txt
+++ b/scripts/reference_exported_symbols.txt
@@ -531,6 +531,7 @@ osgeo::proj::operation::Conversion::createTunisiaMappingGrid(osgeo::proj::util::
 osgeo::proj::operation::Conversion::createTwoPointEquidistant(osgeo::proj::util::PropertyMap const&, osgeo::proj::common::Angle const&, osgeo::proj::common::Angle const&, osgeo::proj::common::Angle const&, osgeo::proj::common::Angle const&, osgeo::proj::common::Length const&, osgeo::proj::common::Length const&)
 osgeo::proj::operation::Conversion::createUTM(osgeo::proj::util::PropertyMap const&, int, bool)
 osgeo::proj::operation::Conversion::createVanDerGrinten(osgeo::proj::util::PropertyMap const&, osgeo::proj::common::Angle const&, osgeo::proj::common::Length const&, osgeo::proj::common::Length const&)
+osgeo::proj::operation::Conversion::createVerticalPerspective(osgeo::proj::util::PropertyMap const&, osgeo::proj::common::Angle const&, osgeo::proj::common::Angle const&, osgeo::proj::common::Length const&, osgeo::proj::common::Length const&, osgeo::proj::common::Length const&, osgeo::proj::common::Length const&)
 osgeo::proj::operation::Conversion::createWagnerIII(osgeo::proj::util::PropertyMap const&, osgeo::proj::common::Angle const&, osgeo::proj::common::Angle const&, osgeo::proj::common::Length const&, osgeo::proj::common::Length const&)
 osgeo::proj::operation::Conversion::createWagnerII(osgeo::proj::util::PropertyMap const&, osgeo::proj::common::Angle const&, osgeo::proj::common::Length const&, osgeo::proj::common::Length const&)
 osgeo::proj::operation::Conversion::createWagnerI(osgeo::proj::util::PropertyMap const&, osgeo::proj::common::Angle const&, osgeo::proj::common::Length const&, osgeo::proj::common::Length const&)
@@ -869,6 +870,7 @@ proj_create_conversion_tunisia_mapping_grid
 proj_create_conversion_two_point_equidistant
 proj_create_conversion_utm
 proj_create_conversion_van_der_grinten
+proj_create_conversion_vertical_perspective
 proj_create_conversion_wagner_i
 proj_create_conversion_wagner_ii
 proj_create_conversion_wagner_iii
diff --git a/src/iso19111/c_api.cpp b/src/iso19111/c_api.cpp
index 090d59a1..27a727dc 100644
--- a/src/iso19111/c_api.cpp
+++ b/src/iso19111/c_api.cpp
@@ -6266,6 +6266,46 @@ PJ *proj_create_conversion_equal_earth(PJ_CONTEXT *ctx, double center_long,
     }
     return nullptr;
 }
+
+// ---------------------------------------------------------------------------
+
+/** \brief Instantiate a conversion based on the Vertical Perspective projection
+ * method.
+ *
+ * See osgeo::proj::operation::Conversion::createVerticalPerspective().
+ *
+ * Linear parameters are expressed in (linear_unit_name,
+ * linear_unit_conv_factor).
+ * Angular parameters are expressed in (ang_unit_name, ang_unit_conv_factor).
+ *
+ * @since 7.0
+ */
+PJ *proj_create_conversion_vertical_perspective(
+    PJ_CONTEXT *ctx, double topo_origin_lat, double topo_origin_long,
+    double topo_origin_height, double view_point_height, double false_easting,
+    double false_northing, const char *ang_unit_name,
+    double ang_unit_conv_factor, const char *linear_unit_name,
+    double linear_unit_conv_factor) {
+    SANITIZE_CTX(ctx);
+    try {
+        UnitOfMeasure linearUnit(
+            createLinearUnit(linear_unit_name, linear_unit_conv_factor));
+        UnitOfMeasure angUnit(
+            createAngularUnit(ang_unit_name, ang_unit_conv_factor));
+        auto conv = Conversion::createVerticalPerspective(
+            PropertyMap(), Angle(topo_origin_lat, angUnit),
+            Angle(topo_origin_long, angUnit),
+            Length(topo_origin_height, linearUnit),
+            Length(view_point_height, linearUnit),
+            Length(false_easting, linearUnit),
+            Length(false_northing, linearUnit));
+        return proj_create_conversion(ctx, conv);
+    } catch (const std::exception &e) {
+        proj_log_error(ctx, __FUNCTION__, e.what());
+    }
+    return nullptr;
+}
+
 /* END: Generated by scripts/create_c_api_projections.py*/
 
 // ---------------------------------------------------------------------------
diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp
index d4f2903f..6e6d4cec 100644
--- a/src/iso19111/coordinateoperation.cpp
+++ b/src/iso19111/coordinateoperation.cpp
@@ -4577,6 +4577,52 @@ ConversionNNPtr Conversion::createEqualEarth(
 
 // ---------------------------------------------------------------------------
 
+/** \brief Instantiate a conversion based on the [Vertical Perspective]
+ * (https://proj.org/operations/projections/nsper.html) projection method.
+ *
+ * This method is defined as [EPSG:9838]
+ * (https://www.epsg-registry.org/export.htm?gml=urn:ogc:def:method:EPSG::9838)
+ *
+ * The PROJ implementation of the EPSG Vertical Perspective has the current
+ * limitations with respect to the method described in EPSG:
+ * <ul>
+ * <li> it is a 2D-only method, ignoring the ellipsoidal height of the point to
+ *      project.</li>
+ * <li> it has only a spherical development.</li>
+ * <li> the height of the topocentric origin is ignored, and thus assumed to be
+ * 0.</li>
+ * </ul>
+ *
+ * For completness, PROJ adds the falseEasting and falseNorthing parameter,
+ * which are not described in EPSG. They should usually be set to 0.
+ *
+ * @param properties See \ref general_properties of the conversion. If the name
+ * is not provided, it is automatically set.
+ * @param topoOriginLat Latitude of topocentric origin
+ * @param topoOriginLong Longitude of topocentric origin
+ * @param topoOriginHeight Ellipsoidal height of topocentric origin. Ignored by
+ * PROJ (that is assumed to be 0)
+ * @param viewPointHeight Viewpoint height with respect to the
+ * topocentric/mapping plane. In the case where topoOriginHeight = 0, this is
+ * the height above the ellipsoid surface at topoOriginLat, topoOriginLong.
+ * @param falseEasting See \ref false_easting . (not in EPSG)
+ * @param falseNorthing See \ref false_northing . (not in EPSG)
+ * @return a new Conversion.
+ *
+ * @since 7.0
+ */
+ConversionNNPtr Conversion::createVerticalPerspective(
+    const util::PropertyMap &properties, const common::Angle &topoOriginLat,
+    const common::Angle &topoOriginLong, const common::Length &topoOriginHeight,
+    const common::Length &viewPointHeight, const common::Length &falseEasting,
+    const common::Length &falseNorthing) {
+    return create(properties, EPSG_CODE_METHOD_VERTICAL_PERSPECTIVE,
+                  createParams(topoOriginLat, topoOriginLong, topoOriginHeight,
+                               viewPointHeight, falseEasting, falseNorthing));
+}
+
+// ---------------------------------------------------------------------------
+
 //! @cond Doxygen_Suppress
 
 static OperationParameterNNPtr createOpParamNameEPSGCode(int code) {
@@ -5450,6 +5496,27 @@ void Conversion::_exportToWKT(io::WKTFormatter *formatter) const {
                     }
                 }
             }
+            // Same for false easting / false northing for Vertical Perspective
+            else if (methodEPSGCode == EPSG_CODE_METHOD_VERTICAL_PERSPECTIVE) {
+                auto opParamvalue =
+                    dynamic_cast<const OperationParameterValue *>(
+                        genOpParamvalue.get());
+                if (opParamvalue) {
+                    const auto paramEPSGCode =
+                        opParamvalue->parameter()->getEPSGCode();
+                    if (paramEPSGCode == EPSG_CODE_PARAMETER_FALSE_EASTING ||
+                        paramEPSGCode == EPSG_CODE_PARAMETER_FALSE_NORTHING) {
+                        const auto &paramValue = opParamvalue->parameterValue();
+                        if (paramValue->type() ==
+                            ParameterValue::Type::MEASURE) {
+                            const auto &measure = paramValue->value();
+                            if (measure.getSIValue() == 0) {
+                                continue;
+                            }
+                        }
+                    }
+                }
+            }
             genOpParamvalue->_exportToWKT(formatter, mapping);
         }
     }
diff --git a/src/proj_constants.h b/src/proj_constants.h
index d3f33c73..619d9d91 100644
--- a/src/proj_constants.h
+++ b/src/proj_constants.h
@@ -212,6 +212,9 @@
 #define EPSG_NAME_METHOD_LABORDE_OBLIQUE_MERCATOR "Laborde Oblique Mercator"
 #define EPSG_CODE_METHOD_LABORDE_OBLIQUE_MERCATOR 9813
 
+#define EPSG_NAME_METHOD_VERTICAL_PERSPECTIVE "Vertical Perspective"
+#define EPSG_CODE_METHOD_VERTICAL_PERSPECTIVE 9838
+
 /* ------------------------------------------------------------------------ */
 
 /* Projection parameters */
@@ -301,6 +304,18 @@
 #define EPSG_NAME_PARAMETER_ELLIPSOID_SCALE_FACTOR "Ellipsoid scaling factor"
 #define EPSG_CODE_PARAMETER_ELLIPSOID_SCALE_FACTOR 1038
 
+#define EPSG_NAME_PARAMETER_LATITUDE_TOPOGRAPHIC_ORIGIN "Latitude of topocentric origin"
+#define EPSG_CODE_PARAMETER_LATITUDE_TOPOGRAPHIC_ORIGIN 8834
+
+#define EPSG_NAME_PARAMETER_LONGITUDE_TOPOGRAPHIC_ORIGIN "Longitude of topocentric origin"
+#define EPSG_CODE_PARAMETER_LONGITUDE_TOPOGRAPHIC_ORIGIN 8835
+
+#define EPSG_NAME_PARAMETER_ELLIPSOIDAL_HEIGHT_TOPOCENTRIC_ORIGIN "Ellipsoidal height of topocentric origin"
+#define EPSG_CODE_PARAMETER_ELLIPSOIDAL_HEIGHT_TOPOCENTRIC_ORIGIN 8836
+
+#define EPSG_NAME_PARAMETER_VIEWPOINT_HEIGHT "Viewpoint height"
+#define EPSG_CODE_PARAMETER_VIEWPOINT_HEIGHT 8840
+
 /* ------------------------------------------------------------------------ */
 
 /* Other conversions and transformations */
diff --git a/src/proj_experimental.h b/src/proj_experimental.h
index dcb2e888..9a1a6045 100644
--- a/src/proj_experimental.h
+++ b/src/proj_experimental.h
@@ -918,6 +918,17 @@ PJ PROJ_DLL *proj_create_conversion_equal_earth(
     const char* ang_unit_name, double ang_unit_conv_factor,
     const char* linear_unit_name, double linear_unit_conv_factor);
 
+PJ PROJ_DLL *proj_create_conversion_vertical_perspective(
+    PJ_CONTEXT *ctx,
+    double topo_origin_lat,
+    double topo_origin_long,
+    double topo_origin_height,
+    double view_point_height,
+    double false_easting,
+    double false_northing,
+    const char* ang_unit_name, double ang_unit_conv_factor,
+    const char* linear_unit_name, double linear_unit_conv_factor);
+
 /* END: Generated by scripts/create_c_api_projections.py*/
 
 /**@}*/
diff --git a/test/unit/test_c_api.cpp b/test/unit/test_c_api.cpp
index 483ecdbd..b4c620ce 100644
--- a/test/unit/test_c_api.cpp
+++ b/test/unit/test_c_api.cpp
@@ -2402,6 +2402,13 @@ TEST_F(CApi, proj_create_projections) {
         ObjectKeeper keeper_projCRS(projCRS);
         ASSERT_NE(projCRS, nullptr);
     }
+    {
+        auto projCRS = proj_create_conversion_vertical_perspective(
+            m_ctxt, 0, 0, 0, 0, 0, 0, "Degree", 0.0174532925199433, "Metre",
+            1.0);
+        ObjectKeeper keeper_projCRS(projCRS);
+        ASSERT_NE(projCRS, nullptr);
+    }
     /* END: Generated by scripts/create_c_api_projections.py*/
 }
 
diff --git a/test/unit/test_operation.cpp b/test/unit/test_operation.cpp
index ddbb1cc9..fcbbb566 100644
--- a/test/unit/test_operation.cpp
+++ b/test/unit/test_operation.cpp
@@ -4091,6 +4091,74 @@ TEST(operation, eqearth_export) {
 
 // ---------------------------------------------------------------------------
 
+TEST(operation, vertical_perspective_export) {
+
+    auto conv = Conversion::createVerticalPerspective(
+        PropertyMap(), Angle(1), Angle(2), Length(3), Length(4), Length(5),
+        Length(6));
+
+    EXPECT_EQ(conv->exportToPROJString(PROJStringFormatter::create().get()),
+              "+proj=nsper +lat_0=1 +lon_0=2 +h=4 +x_0=5 +y_0=6");
+
+    EXPECT_EQ(conv->exportToWKT(WKTFormatter::create().get()),
+              "CONVERSION[\"Vertical Perspective\",\n"
+              "    METHOD[\"Vertical Perspective\",\n"
+              "        ID[\"EPSG\",9838]],\n"
+              "    PARAMETER[\"Latitude of topocentric origin\",1,\n"
+              "        ANGLEUNIT[\"degree\",0.0174532925199433],\n"
+              "        ID[\"EPSG\",8834]],\n"
+              "    PARAMETER[\"Longitude of topocentric origin\",2,\n"
+              "        ANGLEUNIT[\"degree\",0.0174532925199433],\n"
+              "        ID[\"EPSG\",8835]],\n"
+              "    PARAMETER[\"Ellipsoidal height of topocentric origin\",3,\n"
+              "        LENGTHUNIT[\"metre\",1],\n"
+              "        ID[\"EPSG\",8836]],\n"
+              "    PARAMETER[\"Viewpoint height\",4,\n"
+              "        LENGTHUNIT[\"metre\",1],\n"
+              "        ID[\"EPSG\",8840]],\n"
+              "    PARAMETER[\"False easting\",5,\n"
+              "        LENGTHUNIT[\"metre\",1],\n"
+              "        ID[\"EPSG\",8806]],\n"
+              "    PARAMETER[\"False northing\",6,\n"
+              "        LENGTHUNIT[\"metre\",1],\n"
+              "        ID[\"EPSG\",8807]]]");
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(
+    operation,
+    vertical_perspective_export_no_topocentric_height_and_false_easting_northing) {
+
+    auto conv = Conversion::createVerticalPerspective(
+        PropertyMap(), Angle(1), Angle(2), Length(0), Length(4), Length(0),
+        Length(0));
+
+    EXPECT_EQ(conv->exportToPROJString(PROJStringFormatter::create().get()),
+              "+proj=nsper +lat_0=1 +lon_0=2 +h=4 +x_0=0 +y_0=0");
+
+    // Check that False esting and False northing are not exported, when they
+    // are 0.
+    EXPECT_EQ(conv->exportToWKT(WKTFormatter::create().get()),
+              "CONVERSION[\"Vertical Perspective\",\n"
+              "    METHOD[\"Vertical Perspective\",\n"
+              "        ID[\"EPSG\",9838]],\n"
+              "    PARAMETER[\"Latitude of topocentric origin\",1,\n"
+              "        ANGLEUNIT[\"degree\",0.0174532925199433],\n"
+              "        ID[\"EPSG\",8834]],\n"
+              "    PARAMETER[\"Longitude of topocentric origin\",2,\n"
+              "        ANGLEUNIT[\"degree\",0.0174532925199433],\n"
+              "        ID[\"EPSG\",8835]],\n"
+              "    PARAMETER[\"Ellipsoidal height of topocentric origin\",0,\n"
+              "        LENGTHUNIT[\"metre\",1],\n"
+              "        ID[\"EPSG\",8836]],\n"
+              "    PARAMETER[\"Viewpoint height\",4,\n"
+              "        LENGTHUNIT[\"metre\",1],\n"
+              "        ID[\"EPSG\",8840]]]");
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(operation, laborde_oblique_mercator) {
 
     // Content of EPSG:29701 "Tananarive (Paris) / Laborde Grid"