Merge pull request #2637 from rouault/hyperbolic_cass

Add support for hyperbolic Cassini-Soldner

7 files changed, 94 insertions, 41 deletions

diff --git a/docs/source/operations/projections/cass.rst b/docs/source/operations/projections/cass.rst
index 95684e51..5428260a 100644
--- a/docs/source/operations/projections/cass.rst
+++ b/docs/source/operations/projections/cass.rst
@@ -65,6 +65,10 @@ Options
 
 .. include:: ../options/R.rst
 
+.. option:: +hyperbolic
+
+    Use modified form of the standard Cassini-Soldner projection known as the Hyperbolic Cassini-Soldner.
+    This is used in particular for the "Vanua Levu Grid" of the the island of Vanua Levu, Fiji (EPSG:3139)
 
 
 Mathematical definition
diff --git a/src/generic_inverse.cpp b/src/generic_inverse.cpp
index 41cf6c50..28d939fb 100644
--- a/src/generic_inverse.cpp
+++ b/src/generic_inverse.cpp
@@ -55,7 +55,7 @@ PJ_LP pj_generic_inverse_2d(PJ_XY xy, PJ *P, PJ_LP lpInitial) {
             return lp;
         }
 
-        if (fabs(deltaX) > 1e-6 || fabs(deltaY) > 1e-6) {
+        if (i == 0 || fabs(deltaX) > 1e-6 || fabs(deltaY) > 1e-6) {
             // Compute Jacobian matrix (only if we aren't close to the final
             // result to speed things a bit)
             PJ_LP lp2;
diff --git a/src/iso19111/operation/parammappings.cpp b/src/iso19111/operation/parammappings.cpp
index 03200468..15ecb24f 100644
--- a/src/iso19111/operation/parammappings.cpp
+++ b/src/iso19111/operation/parammappings.cpp
@@ -642,6 +642,10 @@ static const MethodMapping projectionMethodMappings[] = {
     {EPSG_NAME_METHOD_CASSINI_SOLDNER, EPSG_CODE_METHOD_CASSINI_SOLDNER,
      "Cassini_Soldner", "cass", nullptr, paramsNatOrigin},
 
+    {EPSG_NAME_METHOD_HYPERBOLIC_CASSINI_SOLDNER,
+     EPSG_CODE_METHOD_HYPERBOLIC_CASSINI_SOLDNER, nullptr, "cass", "hyperbolic",
+     paramsNatOrigin},
+
     {PROJ_WKT2_NAME_METHOD_EQUIDISTANT_CONIC, 0, "Equidistant_Conic", "eqdc",
      nullptr, paramsEQDC},
 
diff --git a/src/proj_constants.h b/src/proj_constants.h
index 0976dc2a..baea05da 100644
--- a/src/proj_constants.h
+++ b/src/proj_constants.h
@@ -88,6 +88,9 @@
 #define EPSG_NAME_METHOD_CASSINI_SOLDNER "Cassini-Soldner"
 #define EPSG_CODE_METHOD_CASSINI_SOLDNER 9806
 
+#define EPSG_NAME_METHOD_HYPERBOLIC_CASSINI_SOLDNER "Hyperbolic Cassini-Soldner"
+#define EPSG_CODE_METHOD_HYPERBOLIC_CASSINI_SOLDNER 9833
+
 #define PROJ_WKT2_NAME_METHOD_EQUIDISTANT_CONIC "Equidistant Conic"
 
 #define PROJ_WKT2_NAME_METHOD_ECKERT_I "Eckert I"
diff --git a/src/projections/cass.cpp b/src/projections/cass.cpp
index a7fb48be..d45532ba 100644
--- a/src/projections/cass.cpp
+++ b/src/projections/cass.cpp
@@ -17,34 +17,40 @@ PROJ_HEAD(cass, "Cassini") "\n\tCyl, Sph&Ell";
 
 
 namespace { // anonymous namespace
-struct pj_opaque {
+struct cass_data {
     double *en;
     double m0;
+    bool hyperbolic;
 };
 } // anonymous namespace
 
 
 
 static PJ_XY cass_e_forward (PJ_LP lp, PJ *P) {          /* Ellipsoidal, forward */
-    double n, t, a1, c, a2, tn;
     PJ_XY xy = {0.0, 0.0};
-    struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque);
-
-    n = sin (lp.phi);
-    c = cos (lp.phi);
-    xy.y = pj_mlfn (lp.phi, n, c, Q->en);
-
-    n  = 1./sqrt(1. - P->es * n*n);
-    tn = tan(lp.phi);
-    t = tn * tn;
-    a1 = lp.lam * c;
-    c *= P->es * c / (1 - P->es);
-    a2 = a1 * a1;
-
-    xy.x = n * a1 * (1. - a2 * t *
-        (C1 - (8. - t + 8. * c) * a2 * C2));
-    xy.y -= Q->m0 - n * tn * a2 *
-        (.5 + (5. - t + 6. * c) * a2 * C3);
+    struct cass_data *Q = static_cast<struct cass_data*>(P->opaque);
+
+    const double sinphi = sin (lp.phi);
+    const double cosphi = cos (lp.phi);
+    const double M = pj_mlfn (lp.phi, sinphi, cosphi, Q->en);
+
+    const double nu_square = 1./(1. - P->es * sinphi*sinphi);
+    const double nu  = sqrt(nu_square);
+    const double tanphi = tan(lp.phi);
+    const double T = tanphi * tanphi;
+    const double A = lp.lam * cosphi;
+    const double C = P->es * (cosphi * cosphi) / (1 - P->es);
+    const double A2 = A * A;
+
+    xy.x = nu * A * (1. - A2 * T *
+        (C1 - (8. - T + 8. * C) * A2 * C2));
+    xy.y = M - Q->m0 + nu * tanphi * A2 *
+        (.5 + (5. - T + 6. * C) * A2 * C3);
+    if( Q->hyperbolic )
+    {
+        const double rho = nu_square * (1. - P->es) * nu;
+        xy.y -= xy.y * xy.y * xy.y / (6 * rho * nu);
+    }
 
     return xy;
 }
@@ -59,23 +65,31 @@ static PJ_XY cass_s_forward (PJ_LP lp, PJ *P) {           /* Spheroidal, forward
 
 
 static PJ_LP cass_e_inverse (PJ_XY xy, PJ *P) {          /* Ellipsoidal, inverse */
-    double n, t, r, dd, d2, tn, ph1;
     PJ_LP lp = {0.0, 0.0};
-    struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque);
-
-    ph1 = pj_inv_mlfn (P->ctx, Q->m0 + xy.y, P->es, Q->en);
-    tn  = tan (ph1);
-    t   = tn*tn;
-    n   = sin (ph1);
-    r   = 1. / (1. - P->es * n * n);
-    n   = sqrt (r);
-    r  *= (1. - P->es) * n;
-    dd  = xy.x / n;
-    d2  = dd * dd;
-    lp.phi = ph1 - (n * tn / r) * d2 *
-        (.5 - (1. + 3. * t) * d2 * C3);
-    lp.lam = dd * (1. + t * d2 *
-        (-C4 + (1. + 3. * t) * d2 * C5)) / cos (ph1);
+    struct cass_data *Q = static_cast<struct cass_data*>(P->opaque);
+
+    const double phi1 = pj_inv_mlfn (P->ctx, Q->m0 + xy.y, P->es, Q->en);
+    const double tanphi1 = tan (phi1);
+    const double T1 = tanphi1*tanphi1;
+    const double sinphi1 = sin (phi1);
+    const double nu1_square = 1. / (1. - P->es * sinphi1 * sinphi1);
+    const double nu1 = sqrt (nu1_square);
+    const double rho1 = nu1_square * (1. - P->es) * nu1;
+    const double D  = xy.x / nu1;
+    const double D2 = D * D;
+    lp.phi = phi1 - (nu1 * tanphi1 / rho1) * D2 *
+        (.5 - (1. + 3. * T1) * D2 * C3);
+    lp.lam = D * (1. + T1 * D2 *
+        (-C4 + (1. + 3. * T1) * D2 * C5)) / cos (phi1);
+
+    if( Q->hyperbolic )
+    {
+        // EPSG guidance note 7-2 suggests a custom approximation for the
+        // 'Vanua Levu 1915 / Vanua Levu Grid' case, but better use the
+        // generic inversion method
+        lp = pj_generic_inverse_2d(xy, P, lp);
+    }
+
     return lp;
 }
 
@@ -95,7 +109,7 @@ static PJ *destructor (PJ *P, int errlev) {                        /* Destructor
     if (nullptr==P->opaque)
         return pj_default_destructor (P, errlev);
 
-    free (static_cast<struct pj_opaque*>(P->opaque)->en);
+    free (static_cast<struct cass_data*>(P->opaque)->en);
     return pj_default_destructor (P, errlev);
 }
 
@@ -111,16 +125,19 @@ PJ *PROJECTION(cass) {
     }
 
     /* otherwise it's ellipsoidal */
-    P->opaque = static_cast<struct pj_opaque*>(calloc (1, sizeof (struct pj_opaque)));
+    auto Q = static_cast<struct cass_data*>(calloc (1, sizeof (struct cass_data)));
+    P->opaque = Q;
     if (nullptr==P->opaque)
         return pj_default_destructor (P, PROJ_ERR_OTHER /*ENOMEM*/);
     P->destructor = destructor;
 
-    static_cast<struct pj_opaque*>(P->opaque)->en = pj_enfn (P->es);
-    if (nullptr==static_cast<struct pj_opaque*>(P->opaque)->en)
+    Q->en = pj_enfn (P->es);
+    if (nullptr==Q->en)
         return pj_default_destructor (P, PROJ_ERR_OTHER /*ENOMEM*/);
 
-    static_cast<struct pj_opaque*>(P->opaque)->m0 = pj_mlfn (P->phi0,  sin (P->phi0),  cos (P->phi0),  static_cast<struct pj_opaque*>(P->opaque)->en);
+    Q->m0 = pj_mlfn (P->phi0,  sin (P->phi0),  cos (P->phi0), Q->en);
+    if (pj_param_exists(P->params, "hyperbolic"))
+        Q->hyperbolic = true;
     P->inv = cass_e_inverse;
     P->fwd = cass_e_forward;
 
diff --git a/test/gie/builtins.gie b/test/gie/builtins.gie
index b3c84f74..8c0cf6d2 100644
--- a/test/gie/builtins.gie
+++ b/test/gie/builtins.gie
@@ -809,6 +809,7 @@ operation +proj=cass   +ellps=GRS80
 tolerance 0.1 mm
 accept  2 1
 expect  222605.285776991 110642.229253999
+roundtrip 1
 accept  2 -1
 expect  222605.285776991 -110642.229253999
 accept  -2 1
@@ -849,6 +850,18 @@ expect  -0.001790493 0.000895247
 accept  -200 -100
 expect  -0.001790493 -0.000895247
 
+-------------------------------------------------------------------------------
+# Hyperbolic variant: test point from EPSG Guidance Note 7.2
+operation +proj=cass +hyperbolic +a=6378306.376305601 +rf=293.466307 \
+          +lat_0=-16.25 +lon_0=179.33333333333333 +to_meter=20.1168 \
+          +x_0=251727.9155424 +y_0=334519.953768
+-------------------------------------------------------------------------------
+
+tolerance 0.1 mm
+accept  179.99433652777776 -16.841456527777776
+expect  16015.289017555102 13369.660053668682
+roundtrip 1
+
 ===============================================================================
 # Central Conic
 # 	Sph
diff --git a/test/unit/test_operation.cpp b/test/unit/test_operation.cpp
index 6ddbd9f4..2830ea94 100644
--- a/test/unit/test_operation.cpp
+++ b/test/unit/test_operation.cpp
@@ -4223,6 +4223,18 @@ TEST(operation, adams_ws2_export_failure) {
 
 // ---------------------------------------------------------------------------
 
+TEST(operation, hyperbolic_cassini_soldner) {
+    auto dbContext = DatabaseContext::create();
+    auto crs =
+        AuthorityFactory::create(dbContext, "EPSG")->createProjectedCRS("3139");
+    EXPECT_EQ(crs->exportToPROJString(PROJStringFormatter::create().get()),
+              "+proj=cass +hyperbolic +lat_0=-16.25 +lon_0=179.333333333333 "
+              "+x_0=251727.9155424 +y_0=334519.953768 "
+              "+a=6378306.3696 +b=6356571.996 +units=link +no_defs +type=crs");
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(operation, PROJ_based) {
     auto conv = SingleOperation::createPROJBased(PropertyMap(), "+proj=merc",
                                                  nullptr, nullptr);