From d436d4fbbfa030b12d6d051a809a6be6c220d747 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Sat, 3 Apr 2021 15:34:46 +0200 Subject: [PATCH 1/4] cass: rewrite ellipsoidal formulas in a clearer way using EPSG guidance note names --- src/projections/cass.cpp | 73 ++++++++++++++++++++-------------------- test/gie/builtins.gie | 1 + 2 files changed, 37 insertions(+), 37 deletions(-) diff --git a/src/projections/cass.cpp b/src/projections/cass.cpp index a7fb48be72..aeddb27cb8 100644 --- a/src/projections/cass.cpp +++ b/src/projections/cass.cpp @@ -17,7 +17,7 @@ PROJ_HEAD(cass, "Cassini") "\n\tCyl, Sph&Ell"; namespace { // anonymous namespace -struct pj_opaque { +struct cass_data { double *en; double m0; }; @@ -26,25 +26,24 @@ struct pj_opaque { 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(P->opaque); + struct cass_data *Q = static_cast(P->opaque); - n = sin (lp.phi); - c = cos (lp.phi); - xy.y = pj_mlfn (lp.phi, n, c, Q->en); + const double sinphi = sin (lp.phi); + const double cosphi = cos (lp.phi); + const double M = pj_mlfn (lp.phi, sinphi, cosphi, 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; + const double nu = 1./sqrt(1. - P->es * sinphi*sinphi); + 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 = 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); + 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); return xy; } @@ -59,23 +58,22 @@ 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(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(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); return lp; } @@ -95,7 +93,7 @@ static PJ *destructor (PJ *P, int errlev) { /* Destructor if (nullptr==P->opaque) return pj_default_destructor (P, errlev); - free (static_cast(P->opaque)->en); + free (static_cast(P->opaque)->en); return pj_default_destructor (P, errlev); } @@ -111,16 +109,17 @@ PJ *PROJECTION(cass) { } /* otherwise it's ellipsoidal */ - P->opaque = static_cast(calloc (1, sizeof (struct pj_opaque))); + auto Q = static_cast(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(P->opaque)->en = pj_enfn (P->es); - if (nullptr==static_cast(P->opaque)->en) + Q->en = pj_enfn (P->es); + if (nullptr==Q->en) return pj_default_destructor (P, PROJ_ERR_OTHER /*ENOMEM*/); - static_cast(P->opaque)->m0 = pj_mlfn (P->phi0, sin (P->phi0), cos (P->phi0), static_cast(P->opaque)->en); + Q->m0 = pj_mlfn (P->phi0, sin (P->phi0), cos (P->phi0), Q->en); P->inv = cass_e_inverse; P->fwd = cass_e_forward; diff --git a/test/gie/builtins.gie b/test/gie/builtins.gie index b3c84f7475..3b4b098be1 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 From bf79113d1de05d139e5a517c705e61eed20a12d8 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Sat, 3 Apr 2021 20:17:15 +0200 Subject: [PATCH 2/4] Generic inverse: fix when initial guess if super close to result --- src/generic_inverse.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/generic_inverse.cpp b/src/generic_inverse.cpp index 41cf6c506c..28d939fb52 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; From 970c849cc70ccb8a016c9834daef70d28e224ae9 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Sat, 3 Apr 2021 20:18:08 +0200 Subject: [PATCH 3/4] cass: add +hyperbolic switch for variant used by EPSG:3139 'Vanua Levu 1915 / Vanua Levu Grid' --- docs/source/operations/projections/cass.rst | 4 ++++ src/projections/cass.cpp | 20 +++++++++++++++++++- test/gie/builtins.gie | 12 ++++++++++++ 3 files changed, 35 insertions(+), 1 deletion(-) diff --git a/docs/source/operations/projections/cass.rst b/docs/source/operations/projections/cass.rst index 95684e51dd..5428260a6f 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/projections/cass.cpp b/src/projections/cass.cpp index aeddb27cb8..d45532baa7 100644 --- a/src/projections/cass.cpp +++ b/src/projections/cass.cpp @@ -20,6 +20,7 @@ namespace { // anonymous namespace struct cass_data { double *en; double m0; + bool hyperbolic; }; } // anonymous namespace @@ -33,7 +34,8 @@ static PJ_XY cass_e_forward (PJ_LP lp, PJ *P) { /* Ellipsoidal, forward const double cosphi = cos (lp.phi); const double M = pj_mlfn (lp.phi, sinphi, cosphi, Q->en); - const double nu = 1./sqrt(1. - P->es * sinphi*sinphi); + 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; @@ -44,6 +46,11 @@ static PJ_XY cass_e_forward (PJ_LP lp, PJ *P) { /* Ellipsoidal, forward (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; } @@ -74,6 +81,15 @@ static PJ_LP cass_e_inverse (PJ_XY xy, PJ *P) { /* Ellipsoidal, inverse (.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; } @@ -120,6 +136,8 @@ PJ *PROJECTION(cass) { return pj_default_destructor (P, PROJ_ERR_OTHER /*ENOMEM*/); 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 3b4b098be1..8c0cf6d28d 100644 --- a/test/gie/builtins.gie +++ b/test/gie/builtins.gie @@ -850,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 From 2594b3a2240db7d7007e1c98f08762613b7f7e73 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Sat, 3 Apr 2021 16:16:22 +0200 Subject: [PATCH 4/4] Add mapping between EPSG method 'Hyperbolic Cassini-Soldner' and +proj=cass +hyperbolic --- src/iso19111/operation/parammappings.cpp | 4 ++++ src/proj_constants.h | 3 +++ test/unit/test_operation.cpp | 12 ++++++++++++ 3 files changed, 19 insertions(+) diff --git a/src/iso19111/operation/parammappings.cpp b/src/iso19111/operation/parammappings.cpp index 03200468f8..15ecb24f3f 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 0976dc2a32..baea05dad8 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/test/unit/test_operation.cpp b/test/unit/test_operation.cpp index 6ddbd9f448..2830ea94a1 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);