Add support for Equidistant Cylindrical ellipsoidal method (EPSG:1028)

docs/source/operations/projections/eqc.rst

@@ -6,10 +6,13 @@ Equidistant Cylindrical (Plate Carrée)
 
 The simplest of all projections. Standard parallels (0° when omitted) may be specified that determine latitude of true scale (k=h=1).
 
+This projection is available in both spherical (EPSG:1029) and ellipsoidal (EPSG:1028) forms.
+When using an ellipsoid, PROJ uses the ellipsoidal formulas from IOGP Guidance Note 7-2.
+
 +---------------------+----------------------------------------------------------+
 | **Classification**  | Conformal cylindrical                                    |
 +---------------------+----------------------------------------------------------+
-| **Available forms** | Forward and inverse                                      |
+| **Available forms** | Forward and inverse, spherical and ellipsoidal           |
 +---------------------+----------------------------------------------------------+
 | **Defined area**    | Global, but best used near the equator                   |
 +---------------------+----------------------------------------------------------+
@@ -92,16 +95,19 @@ Parameters
 Mathematical definition
 #######################
 
-The formulas describing the Equidistant Cylindrical projection are all taken from :cite:`Snyder1987`.
-
 :math:`\phi_{ts}` is the latitude of true scale, i.e., the standard parallel where the scale of the projection is true. It can be set with ``+lat_ts``.
 
 :math:`\phi_0` is the latitude of origin that match the center of the map. It can be set with ``+lat_0``.
 
 
-Forward projection
+Spherical formulas
 ==================
 
+The spherical formulas are taken from :cite:`Snyder1987` and correspond to EPSG method 1029.
+
+Forward projection
+------------------
+
 .. math::
 
    x = \lambda \cos \phi_{ts}
@@ -111,21 +117,58 @@ Forward projection
    y = \phi - \phi_0
 
 Inverse projection
-==================
+------------------
 
 .. math::
 
-   \lambda = x / cos \phi_{ts}
+   \lambda = x / \cos \phi_{ts}
 
 .. math::
 
    \phi = y + \phi_0
 
 
+Ellipsoidal formulas
+====================
+
+The ellipsoidal formulas are from IOGP Guidance Note 7-2, Section 3.2.5, and correspond to EPSG method 1028.
+
+Forward projection
+------------------
+
+.. math::
+
+   x = \nu_1 \cos \phi_{ts} \cdot \lambda
+
+.. math::
+
+   y = M(\phi) - M(\phi_0)
+
+where :math:`\nu_1` is the radius of curvature in the prime vertical at the standard parallel:
+
+.. math::
+
+   \nu_1 = \frac{a}{\sqrt{1 - e^2 \sin^2 \phi_{ts}}}
+
+and :math:`M(\phi)` is the meridional arc distance from the equator to latitude :math:`\phi`.
+
+Inverse projection
+------------------
+
+.. math::
+
+   \lambda = x / (\nu_1 \cos \phi_{ts})
+
+.. math::
+
+   \phi = M^{-1}(y + M(\phi_0))
+
+
 Further reading
 ###############
 
 #. `Wikipedia <https://en.wikipedia.org/wiki/Equirectangular_projection>`_
 #. `Wolfram Mathworld <http://mathworld.wolfram.com/CylindricalEquidistantProjection.html>`_
+#. `IOGP Guidance Note 7-2 <https://www.iogp.org/wp-content/uploads/2019/09/373-07-02.pdf>`_ - Section 3.2.5
 
 

src/projections/eqc.cpp

@@ -1,19 +1,42 @@
+/******************************************************************************
+ * Project:  PROJ
+ * Purpose:  Implementation of the eqc (Equidistant Cylindrical) projection.
+ *           Also known as Plate Carree when lat_ts=0.
+ *
+ * This file implements both the spherical (EPSG:1029) and ellipsoidal
+ * (EPSG:1028) forms of the Equidistant Cylindrical projection.
+ *
+ * Spherical formulas (EPSG:1029):
+ *   E = FE + R cos(lat_ts) (lon - lon_0)
+ *   N = FN + R lat
+ *
+ * Ellipsoidal formulas (EPSG:1028) from IOGP Guidance Note 7-2:
+ *   E = FE + nu1 cos(lat_ts) (lon - lon_0)
+ *   N = FN + M
+ *
+ *
+ * Reference: IOGP Publication 373-7-2, Section 3.2.5
+ *
 
+ *****************************************************************************/
 
 #include <errno.h>
 #include <math.h>
+#include <stdlib.h>
 
 #include "proj.h"
 #include "proj_internal.h"
 
 namespace { // anonymous namespace
 struct pj_eqc_data {
-    double rc;
+    double rc;  // Spherical: cos(lat_ts); Ellipsoidal: nu1 * cos(lat_ts)
+    double M0;  // Meridional arc at latitude of origin (lat_0)
+    double *en; // Coefficients for meridional arc computation
 };
 } // anonymous namespace
 
 PROJ_HEAD(eqc, "Equidistant Cylindrical (Plate Carree)")
-"\n\tCyl, Sph\n\tlat_ts=[, lat_0=0]";
+"\n\tCyl, Sph&Ell\n\tlat_ts=[, lat_0=0]";
 
 static PJ_XY eqc_s_forward(PJ_LP lp, PJ *P) { /* Spheroidal, forward */
     PJ_XY xy = {0.0, 0.0};
@@ -35,21 +58,88 @@ static PJ_LP eqc_s_inverse(PJ_XY xy, PJ *P) { /* Spheroidal, inverse */
     return lp;
 }
 
+// Ellipsoidal forward (EPSG method 1028)
+// E = FE + nu1 cos(lat_ts) (lon - lon_0)
+// N = FN + M - M0
+// where M is the meridional arc from equator to latitude
+static PJ_XY eqc_e_forward(PJ_LP lp, PJ *P) { /* Ellipsoidal, forward */
+    PJ_XY xy = {0.0, 0.0};
+    struct pj_eqc_data *Q = static_cast<struct pj_eqc_data *>(P->opaque);
+
+    const double sinphi = sin(lp.phi);
+    const double cosphi = cos(lp.phi);
+
+    xy.x = Q->rc * lp.lam;
+    xy.y = pj_mlfn(lp.phi, sinphi, cosphi, Q->en) - Q->M0;
+
+    return xy;
+}
+
+// Ellipsoidal inverse (EPSG method 1028)
+// lon = lon_0 + (E - FE) / (nu1 cos lat_ts)
+// lat is found by inverting M = N - FN + M0
+static PJ_LP eqc_e_inverse(PJ_XY xy, PJ *P) { /* Ellipsoidal, inverse */
+    PJ_LP lp = {0.0, 0.0};
+    struct pj_eqc_data *Q = static_cast<struct pj_eqc_data *>(P->opaque);
+
+    lp.lam = xy.x / Q->rc;
+    lp.phi = pj_inv_mlfn(xy.y + Q->M0, Q->en);
+
+    return lp;
+}
+
+static PJ *pj_eqc_destructor(PJ *P, int errlev) { /* Destructor */
+    if (nullptr == P)
+        return nullptr;
+
+    if (nullptr == P->opaque)
+        return pj_default_destructor(P, errlev);
+
+    free(static_cast<struct pj_eqc_data *>(P->opaque)->en);
+    return pj_default_destructor(P, errlev);
+}
+
 PJ *PJ_PROJECTION(eqc) {
     struct pj_eqc_data *Q = static_cast<struct pj_eqc_data *>(
         calloc(1, sizeof(struct pj_eqc_data)));
     if (nullptr == Q)
         return pj_default_destructor(P, PROJ_ERR_OTHER /*ENOMEM*/);
     P->opaque = Q;
+    P->destructor = pj_eqc_destructor;
 
-    if ((Q->rc = cos(pj_param(P->ctx, P->params, "rlat_ts").f)) <= 0.) {
+    const double phi1 = pj_param(P->ctx, P->params, "rlat_ts").f;
+    const double cos_phi1 = cos(phi1);
+
+    if (cos_phi1 <= 0.) {
         proj_log_error(
             P, _("Invalid value for lat_ts: |lat_ts| should be <= 90°"));
         return pj_default_destructor(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
     }
-    P->inv = eqc_s_inverse;
-    P->fwd = eqc_s_forward;
-    P->es = 0.;
+
+    if (P->es != 0.0) {
+        // Ellipsoidal case (EPSG:1028)
+        const double sin_phi1 = sin(phi1);
+        // nu1 = a / sqrt(1 - e^2 sin^2(lat_ts)), normalized by a
+        const double nu1 = 1.0 / sqrt(1.0 - P->es * sin_phi1 * sin_phi1);
+        Q->rc = nu1 * cos_phi1;
+
+        Q->en = pj_enfn(P->n);
+        if (nullptr == Q->en)
+            return pj_default_destructor(P, PROJ_ERR_OTHER /*ENOMEM*/);
+
+        Q->M0 = pj_mlfn(P->phi0, sin(P->phi0), cos(P->phi0), Q->en);
+
+        P->inv = eqc_e_inverse;
+        P->fwd = eqc_e_forward;
+    } else {
+        // Spheroidal case (EPSG:1029)
+        Q->rc = cos_phi1;
+        Q->en = nullptr;
+        Q->M0 = 0.0;
+
+        P->inv = eqc_s_inverse;
+        P->fwd = eqc_s_forward;
+    }
 
     return P;
 }

test/cli/test_cs2cs_ignf.yaml

@@ -35,52 +35,52 @@ tests:
   out: "1238837.253\t5057451.037 0.000"
 - args: +init=IGNF:LAMBE +to +init=IGNF:GEOPORTALFXX -f %.3f
   in: 600000.0000 2600545.4523  0.0000
-  out: "179040.148\t5610495.275 0.000"
+  out: "179356.306\t5585333.153 0.000"
 - args: +init=IGNF:LAMBE +to +init=IGNF:GEOPORTALFXX -f %.3f
   in: 135638.3592 2418760.4094  0.0000
-  out: "-303729.363\t5410118.356 0.000"
+  out: "-304265.704\t5385136.401 0.000"
 - args: +init=IGNF:LAMBE +to +init=IGNF:GEOPORTALFXX -f %.3f
   in: 998137.3947 2413822.2844  0.0000
-  out: "592842.792\t5410120.554 0.000"
+  out: "593889.664\t5385138.596 0.000"
 - args: +init=IGNF:LAMBE +to +init=IGNF:GEOPORTALFXX -f %.3f
   in: 600000.0000 2200000.0000  0.0000
-  out: "179041.670\t5209746.080 0.000"
+  out: "179357.831\t5185007.149 0.000"
 - args: +init=IGNF:LAMBE +to +init=IGNF:GEOPORTALFXX -f %.3f
   in: 311552.5340 1906457.4840  0.0000
-  out: "-96825.465\t4909184.136 0.000"
+  out: "-96996.444\t4884928.293 0.000"
 - args: +init=IGNF:LAMBE +to +init=IGNF:GEOPORTALFXX -f %.3f
   in: 960488.4138 1910172.8812  0.0000
-  out: "523880.019\t4909191.141 0.000"
+  out: "524805.113\t4884935.286 0.000"
 - args: +init=IGNF:LAMBE +to +init=IGNF:GEOPORTALFXX -f %.3f
   in: 600000.0000 1699510.8340  0.0000
-  out: "179047.633\t4708817.007 0.000"
+  out: "179363.804\t4684962.279 0.000"
 - args: +init=IGNF:LAMBE +to +init=IGNF:GEOPORTALFXX -f %.3f
   in: 1203792.5981 626873.17210  0.0000
-  out: "658259.467\t3623786.764 0.000"
+  out: "659421.855\t3603178.891 0.000"
 - args: +init=IGNF:RGF93G +to +init=IGNF:GEOPORTALFXX -f %.3f
   in: 2d20'11.4239243" 50d23'59.7718445" 0.0
-  out: "179040.151\t5610495.281 0.000"
+  out: "179356.309\t5585333.158 0.000"
 - args: +init=IGNF:RGF93G +to +init=IGNF:GEOPORTALFXX -f %.3f
   in: -3d57'49.4051448" 48d35'59.7121716" 0.0
-  out: "-303729.365\t5410118.352 0.000"
+  out: "-304265.706\t5385136.397 0.000"
 - args: +init=IGNF:RGF93G +to +init=IGNF:GEOPORTALFXX -f %.3f
   in: 7d44'12.1439796" 48d35'59.7832558" 0.0
-  out: "592842.794\t5410120.550 0.000"
+  out: "593889.666\t5385138.593 0.000"
 - args: +init=IGNF:RGF93G +to +init=IGNF:GEOPORTALFXX -f %.3f
   in: 2d20'11.4951975" 46d47'59.8029841" 0.0
-  out: "179041.668\t5209746.077 0.000"
+  out: "179357.829\t5185007.147 0.000"
 - args: +init=IGNF:RGF93G +to +init=IGNF:GEOPORTALFXX -f %.3f
   in: -1d15'48.9240599" 44d05'59.8251878" 0.0
-  out: "-96825.467\t4909184.138 0.000"
+  out: "-96996.446\t4884928.296 0.000"
 - args: +init=IGNF:RGF93G +to +init=IGNF:GEOPORTALFXX -f %.3f
   in: 6d50'12.2276489" 44d06'00.0517019" 0.0
-  out: "523880.022\t4909191.143 0.000"
+  out: "524805.116\t4884935.287 0.000"
 - args: +init=IGNF:RGF93G +to +init=IGNF:GEOPORTALFXX -f %.3f
   in: 2d20'11.7754730" 42d18'00.0824436" 0.0
-  out: "179047.634\t4708817.010 0.000"
+  out: "179363.805\t4684962.283 0.000"
 - args: +init=IGNF:RGF93G +to +init=IGNF:GEOPORTALFXX -f %.3f
   in: 9d32'12.6680218" 41d24'00.3542556" 0.0
-  out: "730783.054\t4608637.873 0.000"
+  out: "732073.508\t4585007.331 0.000"
 - args: +init=IGNF:RGF93G +to +init=IGNF:MILLER -f %.3f
   in: 2d20'11.4239243" 50d23'59.7718445" 0.0
   out: "260098.730\t6140682.441 0.000"