Merge pull request #4228 from jjimenezshaw/local-geographic

Add new Conversion "Local Orthographic"

docs/source/operations/projections/ortho.rst

@@ -37,7 +37,10 @@ around a given latitude and longitude.
           ellipsoid must be forced to a sphere, for example by adding a ``+f=0``
           parameter.
 
-This projection method corresponds to ``EPSG:9840``.
+.. note:: Parameters ``k_0`` and ``alpha`` are added in PROJ 9.5.0
+
+This projection method corresponds to ``EPSG:9840``
+(or ``EPSG:1130`` with ``k_0`` or ``alpha``).
 
 Parameters
 ################################################################################
@@ -48,6 +51,22 @@ Parameters
 
 .. include:: ../options/lat_0.rst
 
+.. option:: +alpha=<value>
+
+    .. versionadded:: 9.5.0
+
+    Azimuth clockwise from north at the center of projection.
+
+    *Defaults to 0.0.*
+
+.. option:: +k_0=<value>
+
+   .. versionadded:: 9.5.0
+
+    Scale factor. Determines scale factor used in the projection.
+
+    *Defaults to 1.0.*
+
 .. include:: ../options/ellps.rst
 
 .. include:: ../options/R.rst

include/proj/coordinateoperation.hpp

@@ -1328,6 +1328,13 @@ class PROJ_GCC_DLL Conversion : public SingleOperation {
         const common::Angle &centerLong, const common::Length &falseEasting,
         const common::Length &falseNorthing);
 
+    PROJ_DLL static ConversionNNPtr createLocalOrthographic(
+        const util::PropertyMap &properties, const common::Angle &centerLat,
+        const common::Angle &centerLong,
+        const common::Angle &azimuthInitialLine, const common::Scale &scale,
+        const common::Length &falseEasting,
+        const common::Length &falseNorthing);
+
     PROJ_DLL static ConversionNNPtr createAmericanPolyconic(
         const util::PropertyMap &properties, const common::Angle &centerLat,
         const common::Angle &centerLong, const common::Length &falseEasting,

scripts/build_esri_projection_mapping.py

@@ -215,7 +215,7 @@
             - Standard_Parallel_2: EPSG_NAME_PARAMETER_LATITUDE_2ND_STD_PARALLEL
             - Latitude_Of_Origin: EPSG_NAME_PARAMETER_LATITUDE_FALSE_ORIGIN
 
-    # From GDAL autotest 
+    # From GDAL autotest
     -   WKT2_name: EPSG_NAME_METHOD_LAMBERT_CONIC_CONFORMAL_2SP
         Params:
             - False_Easting: EPSG_NAME_PARAMETER_EASTING_FALSE_ORIGIN
@@ -456,12 +456,12 @@
         - Latitude_Of_Center: EPSG_NAME_PARAMETER_LATITUDE_OF_NATURAL_ORIGIN
 
 - Local:
-    WKT2_name: EPSG_NAME_METHOD_ORTHOGRAPHIC
+    WKT2_name: EPSG_NAME_METHOD_LOCAL_ORTHOGRAPHIC
     Params:
         - False_Easting: EPSG_NAME_PARAMETER_FALSE_EASTING
         - False_Northing: EPSG_NAME_PARAMETER_FALSE_NORTHING
-        - Scale_Factor: 1.0
-        - Azimuth: 0.0
+        - Scale_Factor: EPSG_NAME_PARAMETER_SCALE_FACTOR_PROJECTION_CENTRE
+        - Azimuth: EPSG_NAME_PARAMETER_AZIMUTH_PROJECTION_CENTRE
         - Longitude_Of_Center: EPSG_NAME_PARAMETER_LONGITUDE_OF_NATURAL_ORIGIN
         - Latitude_Of_Center: EPSG_NAME_PARAMETER_LATITUDE_OF_NATURAL_ORIGIN
 

scripts/reference_exported_symbols.txt

@@ -595,6 +595,7 @@ osgeo::proj::operation::Conversion::createLambertConicConformal_2SP_Michigan(osg
 osgeo::proj::operation::Conversion::createLambertConicConformal_2SP(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::createLambertCylindricalEqualArea(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::createLambertCylindricalEqualAreaSpherical(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::createLocalOrthographic(osgeo::proj::util::PropertyMap const&, osgeo::proj::common::Angle const&, osgeo::proj::common::Angle const&, osgeo::proj::common::Angle const&, osgeo::proj::common::Scale const&, osgeo::proj::common::Length const&, osgeo::proj::common::Length const&)
 osgeo::proj::operation::Conversion::createMercatorSpherical(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::createMercatorVariantA(osgeo::proj::util::PropertyMap const&, osgeo::proj::common::Angle const&, osgeo::proj::common::Angle const&, osgeo::proj::common::Scale const&, osgeo::proj::common::Length const&, osgeo::proj::common::Length const&)
 osgeo::proj::operation::Conversion::createMercatorVariantB(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&)
@@ -938,6 +939,7 @@ proj_create_conversion_lambert_conic_conformal_2sp_belgium
 proj_create_conversion_lambert_conic_conformal_2sp_michigan
 proj_create_conversion_lambert_cylindrical_equal_area
 proj_create_conversion_lambert_cylindrical_equal_area_spherical
+proj_create_conversion_local_orthographic
 proj_create_conversion_mercator_variant_a
 proj_create_conversion_mercator_variant_b
 proj_create_conversion_miller_cylindrical

src/iso19111/c_api.cpp

@@ -6775,6 +6775,39 @@ PJ *proj_create_conversion_orthographic(
 }
 // ---------------------------------------------------------------------------
 
+/** \brief Instantiate a ProjectedCRS with a conversion based on the Local
+ * Orthographic projection method.
+ *
+ * See osgeo::proj::operation::Conversion::createLocalOrthographic().
+ *
+ * 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_local_orthographic(
+    PJ_CONTEXT *ctx, double center_lat, double center_long, double azimuth,
+    double scale, 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::createLocalOrthographic(
+            PropertyMap(), Angle(center_lat, angUnit),
+            Angle(center_long, angUnit), Angle(azimuth, angUnit), Scale(scale),
+            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;
+}
+// ---------------------------------------------------------------------------
+
 /** \brief Instantiate a ProjectedCRS with a conversion based on the American
  * Polyconic projection method.
  *

src/iso19111/io.cpp

@@ -11866,6 +11866,13 @@ PROJStringParser::Private::buildProjectedCRS(int iStep,
                 axisType = AxisType::SOUTH_POLE;
             }
         }
+    } else if (step.name == "ortho") {
+        const std::string &k = getParamValueK(step);
+        if ((!k.empty() && getNumericValue(k) != 1.0) ||
+            (hasParamValue(step, "alpha") &&
+             getNumericValue(getParamValue(step, "alpha")) != 0.0)) {
+            mapping = getMapping(EPSG_CODE_METHOD_LOCAL_ORTHOGRAPHIC);
+        }
     }
 
     UnitOfMeasure unit = buildUnit(step, "units", "to_meter");

src/iso19111/operation/conversion.cpp

@@ -1893,6 +1893,36 @@ ConversionNNPtr Conversion::createOrthographic(
 
 // ---------------------------------------------------------------------------
 
+/** \brief Instantiate a conversion based on the
+ * <a href="../../../operations/projections/ortho.html">
+ * Orthographic</a> projection method.
+ *
+ * This method is defined as
+ * <a href="https://epsg.org/coord-operation-method_1130/index.html">
+ * EPSG:1130</a>.
+ *
+ * @param properties See \ref general_properties of the conversion. If the name
+ * is not provided, it is automatically set.
+ * @param centerLat See \ref center_latitude
+ * @param centerLong See \ref center_longitude
+ * @param azimuthInitialLine See \ref azimuth_initial_line
+ * @param scale See \ref scale_factor_initial_line
+ * @param falseEasting See \ref false_easting
+ * @param falseNorthing See \ref false_northing
+ * @return a new Conversion.
+ */
+ConversionNNPtr Conversion::createLocalOrthographic(
+    const util::PropertyMap &properties, const common::Angle &centerLat,
+    const common::Angle &centerLong, const common::Angle &azimuthInitialLine,
+    const common::Scale &scale, const common::Length &falseEasting,
+    const common::Length &falseNorthing) {
+    return create(properties, EPSG_CODE_METHOD_LOCAL_ORTHOGRAPHIC,
+                  createParams(centerLat, centerLong, azimuthInitialLine, scale,
+                               falseEasting, falseNorthing));
+}
+
+// ---------------------------------------------------------------------------
+
 /** \brief Instantiate a conversion based on the
  * <a href="../../../operations/projections/poly.html">
  * American Polyconic</a> projection method.

src/iso19111/operation/esriparammappings.cpp

@@ -615,8 +615,10 @@ static const ESRIParamMapping paramsESRI_Local[] = {
      EPSG_CODE_PARAMETER_FALSE_EASTING, "0.0", false},
     {"False_Northing", EPSG_NAME_PARAMETER_FALSE_NORTHING,
      EPSG_CODE_PARAMETER_FALSE_NORTHING, "0.0", false},
-    {"Scale_Factor", nullptr, 0, "1.0", false},
-    {"Azimuth", nullptr, 0, "0.0", false},
+    {"Scale_Factor", EPSG_NAME_PARAMETER_SCALE_FACTOR_PROJECTION_CENTRE,
+     EPSG_CODE_PARAMETER_SCALE_FACTOR_PROJECTION_CENTRE, "1.0", false},
+    {"Azimuth", EPSG_NAME_PARAMETER_AZIMUTH_PROJECTION_CENTRE,
+     EPSG_CODE_PARAMETER_AZIMUTH_PROJECTION_CENTRE, "0.0", false},
     {"Longitude_Of_Center", EPSG_NAME_PARAMETER_LONGITUDE_OF_NATURAL_ORIGIN,
      EPSG_CODE_PARAMETER_LONGITUDE_OF_NATURAL_ORIGIN, "0.0", false},
     {"Latitude_Of_Center", EPSG_NAME_PARAMETER_LATITUDE_OF_NATURAL_ORIGIN,
@@ -1089,8 +1091,8 @@ static const ESRIMethodMapping esriMappings[] = {
      paramsESRI_New_Zealand_Map_Grid},
     {"Orthographic", PROJ_WKT2_NAME_ORTHOGRAPHIC_SPHERICAL, 0,
      paramsESRI_Orthographic},
-    {"Local", EPSG_NAME_METHOD_ORTHOGRAPHIC, EPSG_CODE_METHOD_ORTHOGRAPHIC,
-     paramsESRI_Local},
+    {"Local", EPSG_NAME_METHOD_LOCAL_ORTHOGRAPHIC,
+     EPSG_CODE_METHOD_LOCAL_ORTHOGRAPHIC, paramsESRI_Local},
     {"Winkel_Tripel", "Winkel Tripel", 0, paramsESRI_Winkel_Tripel},
     {"Aitoff", "Aitoff", 0, paramsESRI_Aitoff},
     {"Flat_Polar_Quartic", PROJ_WKT2_NAME_METHOD_FLAT_POLAR_QUARTIC, 0,

src/iso19111/operation/parammappings.cpp

@@ -367,6 +367,15 @@ static const ParamMapping *const paramsHomTwoPoint[] = {
     &paramFalseNorthingProjectionCentre,
     nullptr};
 
+static const ParamMapping *const paramsNatOriginAzimuthScale[] = {
+    &paramLatitudeNatOrigin,
+    &paramLongitudeNatOrigin,
+    &paramAzimuth,
+    &paramScaleFactorProjectionCentre,
+    &paramFalseEasting,
+    &paramFalseNorthing,
+    nullptr};
+
 static const ParamMapping *const paramsIMWP[] = {
     &paramLongitudeNatOrigin, &paramLatFirstPoint, &paramLatSecondPoint,
     &paramFalseEasting,       &paramFalseNorthing, nullptr};
@@ -822,6 +831,9 @@ static const MethodMapping projectionMethodMappings[] = {
     {EPSG_NAME_METHOD_ORTHOGRAPHIC, EPSG_CODE_METHOD_ORTHOGRAPHIC,
      "Orthographic", "ortho", nullptr, paramsNatOrigin},
 
+    {EPSG_NAME_METHOD_LOCAL_ORTHOGRAPHIC, EPSG_CODE_METHOD_LOCAL_ORTHOGRAPHIC,
+     "Local Orthographic", "ortho", nullptr, paramsNatOriginAzimuthScale},
+
     {PROJ_WKT2_NAME_ORTHOGRAPHIC_SPHERICAL, 0, "Orthographic", "ortho", "f=0",
      paramsNatOrigin},
 

src/proj.h

@@ -2095,6 +2095,12 @@ PJ PROJ_DLL *proj_create_conversion_orthographic(
     double ang_unit_conv_factor, const char *linear_unit_name,
     double linear_unit_conv_factor);
 
+PJ PROJ_DLL *proj_create_conversion_local_orthographic(
+    PJ_CONTEXT *ctx, double center_lat, double center_long, double azimuth,
+    double scale, 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);
+
 PJ PROJ_DLL *proj_create_conversion_american_polyconic(
     PJ_CONTEXT *ctx, double center_lat, double center_long,
     double false_easting, double false_northing, const char *ang_unit_name,

src/proj_constants.h

@@ -214,6 +214,9 @@
 #define EPSG_NAME_METHOD_ORTHOGRAPHIC "Orthographic"
 #define EPSG_CODE_METHOD_ORTHOGRAPHIC 9840
 
+#define EPSG_NAME_METHOD_LOCAL_ORTHOGRAPHIC "Local Orthographic"
+#define EPSG_CODE_METHOD_LOCAL_ORTHOGRAPHIC 1130
+
 #define PROJ_WKT2_NAME_ORTHOGRAPHIC_SPHERICAL "Orthographic (Spherical)"
 
 #define PROJ_WKT2_NAME_METHOD_PATTERSON "Patterson"

src/projections/ortho.cpp

@@ -18,6 +18,8 @@ struct pj_ortho_data {
     double y_shift;
     double y_scale;
     enum pj_ortho_ns::Mode mode;
+    double sinalpha;
+    double cosalpha;
 };
 } // anonymous namespace
 
@@ -72,6 +74,11 @@ static PJ_XY ortho_s_forward(PJ_LP lp, PJ *P) { /* Spheroidal, forward */
         break;
     }
     xy.x = cosphi * sin(lp.lam);
+
+    const double xp = xy.x;
+    const double yp = xy.y;
+    xy.x = (xp * Q->cosalpha - yp * Q->sinalpha) * P->k0;
+    xy.y = (xp * Q->sinalpha + yp * Q->cosalpha) * P->k0;
     return xy;
 }
 
@@ -83,6 +90,11 @@ static PJ_LP ortho_s_inverse(PJ_XY xy, PJ *P) { /* Spheroidal, inverse */
     lp.lam = HUGE_VAL;
     lp.phi = HUGE_VAL;
 
+    const double xf = xy.x;
+    const double yf = xy.y;
+    xy.x = (Q->cosalpha * xf + Q->sinalpha * yf) / P->k0;
+    xy.y = (-Q->sinalpha * xf + Q->cosalpha * yf) / P->k0;
+
     const double rh = hypot(xy.x, xy.y);
     sinc = rh;
     if (sinc > 1.) {
@@ -150,9 +162,11 @@ static PJ_XY ortho_e_forward(PJ_LP lp, PJ *P) { /* Ellipsoidal, forward */
     }
 
     const double nu = 1.0 / sqrt(1.0 - P->es * sinphi * sinphi);
-    xy.x = nu * cosphi * sinlam;
-    xy.y = nu * (sinphi * Q->cosph0 - cosphi * Q->sinph0 * coslam) +
-           P->es * (Q->nu0 * Q->sinph0 - nu * sinphi) * Q->cosph0;
+    const double xp = nu * cosphi * sinlam;
+    const double yp = nu * (sinphi * Q->cosph0 - cosphi * Q->sinph0 * coslam) +
+                      P->es * (Q->nu0 * Q->sinph0 - nu * sinphi) * Q->cosph0;
+    xy.x = (Q->cosalpha * xp - Q->sinalpha * yp) * P->k0;
+    xy.y = (Q->sinalpha * xp + Q->cosalpha * yp) * P->k0;
 
     return xy;
 }
@@ -163,6 +177,11 @@ static PJ_LP ortho_e_inverse(PJ_XY xy, PJ *P) { /* Ellipsoidal, inverse */
 
     const auto SQ = [](double x) { return x * x; };
 
+    const double xf = xy.x;
+    const double yf = xy.y;
+    xy.x = (Q->cosalpha * xf + Q->sinalpha * yf) / P->k0;
+    xy.y = (-Q->sinalpha * xf + Q->cosalpha * yf) / P->k0;
+
     if (Q->mode == pj_ortho_ns::N_POLE || Q->mode == pj_ortho_ns::S_POLE) {
         // Polar case. Forward case equations can be simplified as:
         // xy.x = nu * cosphi * sinlam
@@ -309,6 +328,10 @@ PJ *PJ_PROJECTION(ortho) {
         P->fwd = ortho_e_forward;
     }
 
+    const double alpha = pj_param(P->ctx, P->params, "ralpha").f;
+    Q->sinalpha = sin(alpha);
+    Q->cosalpha = cos(alpha);
+
     return P;
 }
 

test/gie/builtins.gie

@@ -5415,6 +5415,20 @@ direction   inverse
 accept      0               6378137.1
 expect      failure         errno coord_transfm_outside_projection_domain
 
+-------------------------------------------------------------------------------
+# Test the Local Orthographic formulation
+
+# Test data from Guidance Note 7 part 2, August 2024,  p. 101
+# Note: The parameters of the CRS in the example are not
+#       exactly the parameters of the CRS in the database for EPSG:10622
+-------------------------------------------------------------------------------
+operation   +proj=ortho +lat_0=37.628969166666664 +lon_0=-122.39394166666668 +k_0=0.9999968 +alpha=27.7927777777777 +x_0=0 +y_0=0 +ellps=GRS80
+-------------------------------------------------------------------------------
+tolerance   0.1 mm
+accept      -122.3846388888889      37.62607694444444
+expect      876.13676       98.97406
+roundtrip   100
+
 
 ===============================================================================
 # Perspective Conic