Merge pull request #1737 from rouault/proj_create_derived_geographic_crs

Add proj_create_derived_geographic_crs() and proj_create_conversion_pole_rotation_grib_convention() to address GRIB datasets using a pole rotation method

5 files changed, 208 insertions, 1 deletions

diff --git a/src/iso19111/c_api.cpp b/src/iso19111/c_api.cpp
index edd4fbf3..9db9e5b3 100644
--- a/src/iso19111/c_api.cpp
+++ b/src/iso19111/c_api.cpp
@@ -2947,6 +2947,66 @@ PJ *proj_create_geocentric_crs_from_datum(PJ_CONTEXT *ctx, const char *crs_name,
 
 // ---------------------------------------------------------------------------
 
+/** \brief Create a DerivedGeograhicCRS.
+ *
+ * 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 base_geographic_crs Base Geographic CRS. Must not be NULL.
+ * @param conversion Conversion from the base Geographic to the
+ * DerivedGeograhicCRS. Must not be NULL.
+ * @param ellipsoidal_cs Coordinate system. Must not be NULL.
+ *
+ * @return Object of type GeodeticCRS that must be unreferenced with
+ * proj_destroy(), or NULL in case of error.
+ *
+ * @since 7.0
+ */
+PJ *proj_create_derived_geographic_crs(PJ_CONTEXT *ctx, const char *crs_name,
+                                       const PJ *base_geographic_crs,
+                                       const PJ *conversion,
+                                       const PJ *ellipsoidal_cs) {
+    SANITIZE_CTX(ctx);
+    auto base_crs =
+        std::dynamic_pointer_cast<GeographicCRS>(base_geographic_crs->iso_obj);
+    auto conversion_cpp =
+        std::dynamic_pointer_cast<Conversion>(conversion->iso_obj);
+    auto cs = std::dynamic_pointer_cast<EllipsoidalCS>(ellipsoidal_cs->iso_obj);
+    if (!base_crs || !conversion_cpp || !cs) {
+        return nullptr;
+    }
+    try {
+        auto derivedCRS = DerivedGeographicCRS::create(
+            createPropertyMapName(crs_name), NN_NO_CHECK(base_crs),
+            NN_NO_CHECK(conversion_cpp), NN_NO_CHECK(cs));
+        return pj_obj_create(ctx, derivedCRS);
+    } catch (const std::exception &e) {
+        proj_log_error(ctx, __FUNCTION__, e.what());
+    }
+    return nullptr;
+}
+
+// ---------------------------------------------------------------------------
+
+/** \brief Return whether a CRS is a Derived CRS.
+ *
+ * @param ctx PROJ context, or NULL for default context
+ * @param crs CRS. Must not be NULL.
+ *
+ * @return whether a CRS is a Derived CRS.
+ *
+ * @since 7.0
+ */
+int proj_is_derived_crs(PJ_CONTEXT *ctx, const PJ *crs) {
+    SANITIZE_CTX(ctx);
+    return dynamic_cast<DerivedCRS *>(crs->iso_obj.get()) != nullptr;
+}
+
+// ---------------------------------------------------------------------------
+
 /** \brief Create a VerticalCRS
  *
  * The returned object must be unreferenced with proj_destroy() after
@@ -6454,6 +6514,36 @@ PJ *proj_create_conversion_vertical_perspective(
     return nullptr;
 }
 
+// ---------------------------------------------------------------------------
+
+/** \brief Instantiate a conversion based on the Pole Rotation method, using the
+ * conventions of the GRIB 1 and GRIB 2 data formats.
+ *
+ * See osgeo::proj::operation::Conversion::createPoleRotationGRIBConvention().
+ *
+ * Linear parameters are expressed in (linear_unit_name,
+ * linear_unit_conv_factor).
+ * Angular parameters are expressed in (ang_unit_name, ang_unit_conv_factor).
+ */
+PJ *proj_create_conversion_pole_rotation_grib_convention(
+    PJ_CONTEXT *ctx, double south_pole_lat_in_unrotated_crs,
+    double south_pole_long_in_unrotated_crs, double axis_rotation,
+    const char *ang_unit_name, double ang_unit_conv_factor) {
+    SANITIZE_CTX(ctx);
+    try {
+        UnitOfMeasure angUnit(
+            createAngularUnit(ang_unit_name, ang_unit_conv_factor));
+        auto conv = Conversion::createPoleRotationGRIBConvention(
+            PropertyMap(), Angle(south_pole_lat_in_unrotated_crs, angUnit),
+            Angle(south_pole_long_in_unrotated_crs, angUnit),
+            Angle(axis_rotation, angUnit));
+        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 556f18b6..1df16e3d 100644
--- a/src/iso19111/coordinateoperation.cpp
+++ b/src/iso19111/coordinateoperation.cpp
@@ -296,6 +296,12 @@ const MethodMapping *getMapping(const char *wkt2_name) noexcept {
             return &mapping;
         }
     }
+    for (const auto &mapping : otherMethodMappings) {
+        if (metadata::Identifier::isEquivalentName(mapping.wkt2_name,
+                                                   wkt2_name)) {
+            return &mapping;
+        }
+    }
     return nullptr;
 }
 
@@ -1697,6 +1703,16 @@ double SingleOperation::parameterValueNumeric(
     return 0.0;
 }
 
+double SingleOperation::parameterValueNumeric(
+    const char *param_name, const common::UnitOfMeasure &targetUnit) const
+    noexcept {
+    const auto &val = parameterValue(param_name, 0);
+    if (val && val->type() == ParameterValue::Type::MEASURE) {
+        return val->value().convertToUnit(targetUnit);
+    }
+    return 0.0;
+}
+
 //! @endcond
 // ---------------------------------------------------------------------------
 
@@ -4669,6 +4685,58 @@ ConversionNNPtr Conversion::createVerticalPerspective(
 
 // ---------------------------------------------------------------------------
 
+/** \brief Instantiate a conversion based on the Pole Rotation method, using
+ * the conventions of the GRIB 1 and GRIB 2 data formats.
+ *
+ * Those are mentionned in the Note 2 of
+ * https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_temp3-1.shtml
+ *
+ * Several conventions for the pole rotation method exists.
+ * The parameters provided in this method are remapped to the PROJ ob_tran
+ * operation with:
+ * <pre>
+ * +proj=ob_tran +o_proj=longlat +o_lon_p=-rotationAngle
+ *                               +o_lat_p=-southPoleLatInUnrotatedCRS
+ *                               +lon_0=southPoleLongInUnrotatedCRS
+ * </pre>
+ *
+ * Another implementation of that convention is also in the netcdf-java library:
+ * https://github.com/Unidata/netcdf-java/blob/3ce72c0cd167609ed8c69152bb4a004d1daa9273/cdm/core/src/main/java/ucar/unidata/geoloc/projection/RotatedLatLon.java
+ *
+ * The PROJ implementation of this method assumes a spherical ellipsoid.
+ *
+ * @param properties See \ref general_properties of the conversion. If the name
+ * is not provided, it is automatically set.
+ * @param southPoleLatInUnrotatedCRS Latitude of the point from the unrotated
+ * CRS, expressed in the unrotated CRS, that will become the south pole of the
+ * rotated CRS.
+ * @param southPoleLongInUnrotatedCRS Longitude of the point from the unrotated
+ * CRS, expressed in the unrotated CRS, that will become the south pole of the
+ * rotated CRS.
+ * @param axisRotation The angle of rotation about the new polar
+ * axis (measured clockwise when looking from the southern to the northern pole)
+ * of the coordinate system, assuming the new axis to have been obtained by
+ * first rotating the sphere through southPoleLongInUnrotatedCRS degrees about
+ * the geographic polar axis and then rotating through
+ * (90 + southPoleLatInUnrotatedCRS) degrees so that the southern pole moved
+ * along the (previously rotated) Greenwich meridian.
+ * @return a new Conversion.
+ *
+ * @since 7.0
+ */
+ConversionNNPtr Conversion::createPoleRotationGRIBConvention(
+    const util::PropertyMap &properties,
+    const common::Angle &southPoleLatInUnrotatedCRS,
+    const common::Angle &southPoleLongInUnrotatedCRS,
+    const common::Angle &axisRotation) {
+    return create(properties,
+                  PROJ_WKT2_NAME_METHOD_POLE_ROTATION_GRIB_CONVENTION,
+                  createParams(southPoleLatInUnrotatedCRS,
+                               southPoleLongInUnrotatedCRS, axisRotation));
+}
+
+// ---------------------------------------------------------------------------
+
 //! @cond Doxygen_Suppress
 
 static OperationParameterNNPtr createOpParamNameEPSGCode(int code) {
@@ -6069,6 +6137,23 @@ void Conversion::_exportToPROJString(
     } else if (starts_with(methodName, "PROJ ")) {
         bConversionDone = true;
         createPROJExtensionFromCustomProj(this, formatter, false);
+    } else if (ci_equal(methodName,
+                        PROJ_WKT2_NAME_METHOD_POLE_ROTATION_GRIB_CONVENTION)) {
+        double southPoleLat = parameterValueNumeric(
+            PROJ_WKT2_NAME_PARAMETER_SOUTH_POLE_LATITUDE_GRIB_CONVENTION,
+            common::UnitOfMeasure::DEGREE);
+        double southPoleLon = parameterValueNumeric(
+            PROJ_WKT2_NAME_PARAMETER_SOUTH_POLE_LONGITUDE_GRIB_CONVENTION,
+            common::UnitOfMeasure::DEGREE);
+        double rotation = parameterValueNumeric(
+            PROJ_WKT2_NAME_PARAMETER_AXIS_ROTATION_GRIB_CONVENTION,
+            common::UnitOfMeasure::DEGREE);
+        formatter->addStep("ob_tran");
+        formatter->addParam("o_proj", "longlat");
+        formatter->addParam("o_lon_p", -rotation);
+        formatter->addParam("o_lat_p", -southPoleLat);
+        formatter->addParam("lon_0", southPoleLon);
+        bConversionDone = true;
     } else if (formatter->convention() ==
                    io::PROJStringFormatter::Convention::PROJ_5 &&
                isZUnitConversion) {
diff --git a/src/iso19111/crs.cpp b/src/iso19111/crs.cpp
index a87ba82d..49cc050f 100644
--- a/src/iso19111/crs.cpp
+++ b/src/iso19111/crs.cpp
@@ -45,6 +45,8 @@
 #include "proj/internal/internal.hpp"
 #include "proj/internal/io_internal.hpp"
 
+#include "proj_constants.h"
+
 #include <algorithm>
 #include <cassert>
 #include <cstring>
@@ -4892,7 +4894,9 @@ void DerivedGeographicCRS::_exportToPROJString(
     if (methodName == "PROJ ob_tran o_proj=longlat" ||
         methodName == "PROJ ob_tran o_proj=lonlat" ||
         methodName == "PROJ ob_tran o_proj=latlong" ||
-        methodName == "PROJ ob_tran o_proj=latlon") {
+        methodName == "PROJ ob_tran o_proj=latlon" ||
+        ci_equal(methodName,
+                 PROJ_WKT2_NAME_METHOD_POLE_ROTATION_GRIB_CONVENTION)) {
         l_conv->_exportToPROJString(formatter);
         return;
     }
diff --git a/src/proj_constants.h b/src/proj_constants.h
index 62cf94be..5c642862 100644
--- a/src/proj_constants.h
+++ b/src/proj_constants.h
@@ -215,6 +215,8 @@
 #define EPSG_NAME_METHOD_VERTICAL_PERSPECTIVE "Vertical Perspective"
 #define EPSG_CODE_METHOD_VERTICAL_PERSPECTIVE 9838
 
+#define PROJ_WKT2_NAME_METHOD_POLE_ROTATION_GRIB_CONVENTION "Pole rotation (GRIB convention)"
+
 /* ------------------------------------------------------------------------ */
 
 /* Projection parameters */
@@ -474,6 +476,15 @@
 #define \
     EPSG_NAME_PARAMETER_FLATTENING_DIFFERENCE "Flattening difference"
 
+#define PROJ_WKT2_NAME_PARAMETER_SOUTH_POLE_LATITUDE_GRIB_CONVENTION \
+    "Latitude of the southern pole (GRIB convention)"
+
+#define PROJ_WKT2_NAME_PARAMETER_SOUTH_POLE_LONGITUDE_GRIB_CONVENTION \
+    "Longitude of the southern pole (GRIB convention)"
+
+#define PROJ_WKT2_NAME_PARAMETER_AXIS_ROTATION_GRIB_CONVENTION \
+    "Axis rotation (GRIB convention)"
+
 /* ------------------------------------------------------------------------ */
 
 #define EPSG_CODE_METHOD_NTV1 9614
diff --git a/src/proj_experimental.h b/src/proj_experimental.h
index 0886ba28..50f1fb9f 100644
--- a/src/proj_experimental.h
+++ b/src/proj_experimental.h
@@ -191,6 +191,16 @@ PJ PROJ_DLL *proj_create_geocentric_crs_from_datum(
                             const char *linear_units,
                             double linear_units_conv);
 
+PJ PROJ_DLL *proj_create_derived_geographic_crs(
+                            PJ_CONTEXT *ctx,
+                            const char *crs_name,
+                            const PJ* base_geographic_crs,
+                            const PJ* conversion,
+                            const PJ* ellipsoidal_cs);
+
+int PROJ_DLL proj_is_derived_crs(PJ_CONTEXT *ctx,
+                                 const PJ* crs);
+
 PJ PROJ_DLL *proj_alter_name(PJ_CONTEXT *ctx,
                                      const PJ* obj, const char* name);
 
@@ -947,6 +957,13 @@ PJ PROJ_DLL *proj_create_conversion_vertical_perspective(
     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_pole_rotation_grib_convention(
+    PJ_CONTEXT *ctx,
+    double south_pole_lat_in_unrotated_crs,
+    double south_pole_long_in_unrotated_crs,
+    double axis_rotation,
+    const char* ang_unit_name, double ang_unit_conv_factor);
+
 /* END: Generated by scripts/create_c_api_projections.py*/
 
 /**@}*/