From 73ca7e68c6d0acd1e00597557f7f83b8246d49f6 Mon Sep 17 00:00:00 2001 From: boja8320 Date: Tue, 21 Aug 2018 20:31:57 -0700 Subject: [PATCH 1/5] Adding Equal Earth ellipsoidal math --- src/PJ_eqearth.c | 135 ++++++++++++++++++++++++++++++++++--- test/gie/more_builtins.gie | 49 ++++++++++++++ 2 files changed, 173 insertions(+), 11 deletions(-) diff --git a/src/PJ_eqearth.c b/src/PJ_eqearth.c index 2ca5ca208a..df7380eec7 100644 --- a/src/PJ_eqearth.c +++ b/src/PJ_eqearth.c @@ -9,14 +9,16 @@ projection, International Journal of Geographical Information Science, DOI: 10.1080/13658816.2018.1504949 Port to PROJ by Juernjakob Dugge, 16 August 2018 +Added ellipsoidal equations by Bojan Savric, 22 August 2018 */ #define PJ_LIB__ +#include #include #include "projects.h" -PROJ_HEAD(eqearth, "Equal Earth") "\n\tPCyl., Sph."; +PROJ_HEAD(eqearth, "Equal Earth") "\n\tPCyl., Sph&Ell\n\th="; #define A1 1.340264 #define A2 -0.081106 @@ -28,28 +30,112 @@ PROJ_HEAD(eqearth, "Equal Earth") "\n\tPCyl., Sph."; #define EPS 1e-11 #define MAX_ITER 12 +struct pj_opaque { + double qp; + double rqda; + double *apa; +}; + +static XY e_forward (LP lp, PJ *P) { /* Ellipsoidal, forward */ + XY xy = {0.0,0.0}; + struct pj_opaque *Q = P->opaque; + double sbeta; + double psi, psi2, psi6; + + /* Authalic latitude */ + sbeta = pj_qsfn(sin(lp.phi), P->e, 1.0 - P->es) / Q->qp; + if (fabs(sbeta) > 1) { + /* Rounding error. */ + sbeta = sbeta > 0 ? 1 : (sbeta < 0 ? -1 : 0); + } + + /* Equal Earth on an authalic sphere */ + psi = asin(M * sbeta); + psi2 = psi * psi; + psi6 = psi2 * psi2 * psi2; + + xy.x = lp.lam * cos(psi) / (M * (A1 + 3 * A2 * psi2 + psi6 * (7 * A3 + 9 * A4 * psi2))); + xy.y = psi * (A1 + A2 * psi2 + psi6 * (A3 + A4 * psi2)); + + /* Adjusting x and y for authalic radius */ + xy.x *= Q->rqda; + xy.y *= Q->rqda; + + return xy; +} + static XY s_forward (LP lp, PJ *P) { /* Spheroidal, forward */ XY xy = {0.0,0.0}; - double phi, phi2, phi6; + double psi, psi2, psi6; (void) P; - phi = asin(M * sin(lp.phi)); - phi2 = phi * phi; - phi6 = phi2 * phi2 * phi2; + psi = asin(M * sin(lp.phi)); + psi2 = psi * psi; + psi6 = psi2 * psi2 * psi2; - xy.x = lp.lam * cos(phi) / (M * (A1 + 3 * A2 * phi2 + phi6 * (7 * A3 + 9 * A4 * phi2))); - xy.y = phi * (A1 + A2 * phi2 + phi6 * (A3 + A4 * phi2)); + xy.x = lp.lam * cos(psi) / (M * (A1 + 3 * A2 * psi2 + psi6 * (7 * A3 + 9 * A4 * psi2))); + xy.y = psi * (A1 + A2 * psi2 + psi6 * (A3 + A4 * psi2)); return xy; } +static LP e_inverse (XY xy, PJ *P) { /* Ellipsoidal, inverse */ + LP lp = {0.0,0.0}; + struct pj_opaque *Q = P->opaque; + double yc, tol, y2, y6, f, fder; + double beta; + int i; + + /* Adjusting x and y for authalic radius */ + xy.x /= Q->rqda; + xy.y /= Q->rqda; + + /* Make sure y is inside valid range */ + if (xy.y > MAX_Y) { + xy.y = MAX_Y; + } else if (xy.y < -MAX_Y) { + xy.y = -MAX_Y; + } + + yc = xy.y; + + for (i = MAX_ITER; i ; --i) { /* Newton-Raphson */ + y2 = yc * yc; + y6 = y2 * y2 * y2; + f = yc * (A1 + A2 * y2 + y6 * (A3 + A4 * y2)) - xy.y; + fder = A1 + 3 * A2 * y2 + y6 * (7 * A3 + 9 * A4 * y2); + tol = f / fder; + yc -= tol; + if (fabs(tol) < EPS) { + break; + } + } + + if( i == 0 ) { + pj_ctx_set_errno( P->ctx, PJD_ERR_NON_CONVERGENT ); + return lp; + } + + /* Longitude */ + y2 = yc * yc; + y6 = y2 * y2 * y2; + + lp.lam = M * xy.x * (A1 + 3 * A2 * y2 + y6 * (7 * A3 + 9 * A4 * y2)) / cos(yc); + + /* Latitude */ + beta = asin(sin(yc) / M); + lp.phi = pj_authlat(beta, Q->apa); + + return lp; +} + static LP s_inverse (XY xy, PJ *P) { /* Spheroidal, inverse */ LP lp = {0.0,0.0}; double yc, tol, y2, y6, f, fder; int i; (void) P; - /* make sure y is inside valid range */ + /* Make sure y is inside valid range */ if (xy.y > MAX_Y) { xy.y = MAX_Y; } else if (xy.y < -MAX_Y) { @@ -84,11 +170,38 @@ static LP s_inverse (XY xy, PJ *P) { /* Spheroidal, inverse */ return lp; } +static void *destructor (PJ *P, int errlev) { /* Destructor */ + if (0==P) + return 0; + + if (0==P->opaque) + return pj_default_destructor (P, errlev); + + pj_dealloc (P->opaque->apa); + return pj_default_destructor (P, errlev); +} + PJ *PROJECTION(eqearth) { - P->es = 0; - P->inv = s_inverse; - P->fwd = s_forward; + struct pj_opaque *Q = pj_calloc (1, sizeof (struct pj_opaque)); + if (0==Q) + return pj_default_destructor (P, ENOMEM); + P->opaque = Q; + P->destructor = destructor; + + if (P->es != 0.0) { + Q->apa = pj_authset(P->es); /* For auth_lat(). */ + if (0 == Q->apa) + return destructor(P, ENOMEM); + Q->qp = pj_qsfn(1.0, P->e, P->one_es); /* For auth_lat(). */ + Q->rqda = sqrt(0.5*Q->qp); /* Authalic radius devided by major axis */ + P->fwd = e_forward; + P->inv = e_inverse; + } + else { + P->fwd = s_forward; + P->inv = s_inverse; + } return P; } diff --git a/test/gie/more_builtins.gie b/test/gie/more_builtins.gie index 2908fd633e..48139b1bd0 100644 --- a/test/gie/more_builtins.gie +++ b/test/gie/more_builtins.gie @@ -419,6 +419,55 @@ tolerance 1cm accept 0 0 expect 0 0 +accept -180 90 +expect -10216474.79 8392927.6 + +accept 0 90 +expect 0 8392927.6 + +accept 180 90 +expect 10216474.79 8392927.6 + +accept 180 45 +expect 14792474.75 5466867.76 + +accept 180 0 +expect 17243959.06 0 + +accept -70 -31.2 +expect -6241081.64 -3907019.16 + + +direction inverse + +accept -6241081.64 -3907019.16 +expect -70 -31.2 + +accept 17243959.06 0 +expect 180 0 + +accept 14792474.75 5466867.76 +expect 180 45 + +accept 0 0 +expect 0 0 + +accept -10216474.79 8392927.6 +expect -180 90 + +accept 0 8392927.6 +expect 0 90 + +accept 10216474.79 8392927.6 +expect 180 90 + + +operation +proj=eqearth +R=6378137 +direction forward + +accept 0 0 +expect 0 0 + accept -180 90 expect -10227908.09 8402320.16 From 67b9ea5d99472238657857aa1a93f50ba585b98d Mon Sep 17 00:00:00 2001 From: boja8320 Date: Wed, 22 Aug 2018 07:11:24 -0700 Subject: [PATCH 2/5] Adding tolerance to the tests... --- test/gie/more_builtins.gie | 1 + 1 file changed, 1 insertion(+) diff --git a/test/gie/more_builtins.gie b/test/gie/more_builtins.gie index 48139b1bd0..6fa081eb93 100644 --- a/test/gie/more_builtins.gie +++ b/test/gie/more_builtins.gie @@ -464,6 +464,7 @@ expect 180 90 operation +proj=eqearth +R=6378137 direction forward +tolerance 1cm accept 0 0 expect 0 0 From 8197803bd40d9861de642445baa2687d5f4c825a Mon Sep 17 00:00:00 2001 From: boja8320 Date: Wed, 22 Aug 2018 09:21:16 -0700 Subject: [PATCH 3/5] Joining ellipsoidal and sphere functions --- src/PJ_eqearth.c | 92 +++++++++++++----------------------------------- 1 file changed, 25 insertions(+), 67 deletions(-) diff --git a/src/PJ_eqearth.c b/src/PJ_eqearth.c index df7380eec7..35795bbc5c 100644 --- a/src/PJ_eqearth.c +++ b/src/PJ_eqearth.c @@ -36,20 +36,26 @@ struct pj_opaque { double *apa; }; -static XY e_forward (LP lp, PJ *P) { /* Ellipsoidal, forward */ +static XY e_forward (LP lp, PJ *P) { /* Ellipsoidal/Spheroidal, forward */ XY xy = {0.0,0.0}; struct pj_opaque *Q = P->opaque; double sbeta; double psi, psi2, psi6; - /* Authalic latitude */ - sbeta = pj_qsfn(sin(lp.phi), P->e, 1.0 - P->es) / Q->qp; - if (fabs(sbeta) > 1) { - /* Rounding error. */ - sbeta = sbeta > 0 ? 1 : (sbeta < 0 ? -1 : 0); + if (P->es != 0.0) { + /* Ellipsoidal case, converting to authalic latitude */ + sbeta = pj_qsfn(sin(lp.phi), P->e, 1.0 - P->es) / Q->qp; + if (fabs(sbeta) > 1) { + /* Rounding error. */ + sbeta = sbeta > 0 ? 1 : -1; + } + } + else { + /* Spheroidal case, using latitude */ + sbeta = sin(lp.phi); } - /* Equal Earth on an authalic sphere */ + /* Equal Earth projection */ psi = asin(M * sbeta); psi2 = psi * psi; psi6 = psi2 * psi2 * psi2; @@ -64,22 +70,8 @@ static XY e_forward (LP lp, PJ *P) { /* Ellipsoidal, forward */ return xy; } -static XY s_forward (LP lp, PJ *P) { /* Spheroidal, forward */ - XY xy = {0.0,0.0}; - double psi, psi2, psi6; - (void) P; - - psi = asin(M * sin(lp.phi)); - psi2 = psi * psi; - psi6 = psi2 * psi2 * psi2; - - xy.x = lp.lam * cos(psi) / (M * (A1 + 3 * A2 * psi2 + psi6 * (7 * A3 + 9 * A4 * psi2))); - xy.y = psi * (A1 + A2 * psi2 + psi6 * (A3 + A4 * psi2)); - return xy; -} - -static LP e_inverse (XY xy, PJ *P) { /* Ellipsoidal, inverse */ +static LP e_inverse (XY xy, PJ *P) { /* Ellipsoidal/Spheroidal, inverse */ LP lp = {0.0,0.0}; struct pj_opaque *Q = P->opaque; double yc, tol, y2, y6, f, fder; @@ -123,50 +115,16 @@ static LP e_inverse (XY xy, PJ *P) { /* Ellipsoidal, inverse */ lp.lam = M * xy.x * (A1 + 3 * A2 * y2 + y6 * (7 * A3 + 9 * A4 * y2)) / cos(yc); /* Latitude */ - beta = asin(sin(yc) / M); - lp.phi = pj_authlat(beta, Q->apa); - - return lp; -} - -static LP s_inverse (XY xy, PJ *P) { /* Spheroidal, inverse */ - LP lp = {0.0,0.0}; - double yc, tol, y2, y6, f, fder; - int i; - (void) P; - - /* Make sure y is inside valid range */ - if (xy.y > MAX_Y) { - xy.y = MAX_Y; - } else if (xy.y < -MAX_Y) { - xy.y = -MAX_Y; - } - - yc = xy.y; - - for (i = MAX_ITER; i ; --i) { /* Newton-Raphson */ - y2 = yc * yc; - y6 = y2 * y2 * y2; - f = yc * (A1 + A2 * y2 + y6 * (A3 + A4 * y2)) - xy.y; - fder = A1 + 3 * A2 * y2 + y6 * (7 * A3 + 9 * A4 * y2); - tol = f / fder; - yc -= tol; - if (fabs(tol) < EPS) { - break; - } + beta = asin(sin(yc) / M); + if (P->es != 0.0) { + /* Ellipsoidal case, converting authalic latitude */ + lp.phi = pj_authlat(beta, Q->apa); } - - if( i == 0 ) { - pj_ctx_set_errno( P->ctx, PJD_ERR_NON_CONVERGENT ); - return lp; + else { + /* Spheroidal case, using latitude */ + lp.phi = beta; } - y2 = yc * yc; - y6 = y2 * y2 * y2; - - lp.lam = M * xy.x * (A1 + 3 * A2 * y2 + y6 * (7 * A3 + 9 * A4 * y2)) / cos(yc); - lp.phi = asin(sin(yc) / M); - return lp; } @@ -195,13 +153,13 @@ PJ *PROJECTION(eqearth) { return destructor(P, ENOMEM); Q->qp = pj_qsfn(1.0, P->e, P->one_es); /* For auth_lat(). */ Q->rqda = sqrt(0.5*Q->qp); /* Authalic radius devided by major axis */ - P->fwd = e_forward; - P->inv = e_inverse; } else { - P->fwd = s_forward; - P->inv = s_inverse; + Q->rqda = 1.0; } + P->fwd = e_forward; + P->inv = e_inverse; + return P; } From 670e874e6adb420e6e8bf6f94b31fbcf12b59994 Mon Sep 17 00:00:00 2001 From: boja8320 Date: Wed, 22 Aug 2018 09:59:18 -0700 Subject: [PATCH 4/5] Adding doc changes --- .../source/operations/projections/eqearth.rst | 30 +++++++++---------- src/PJ_eqearth.c | 2 +- 2 files changed, 16 insertions(+), 16 deletions(-) diff --git a/docs/source/operations/projections/eqearth.rst b/docs/source/operations/projections/eqearth.rst index dc5dd8f9e6..cb0d607d25 100644 --- a/docs/source/operations/projections/eqearth.rst +++ b/docs/source/operations/projections/eqearth.rst @@ -6,21 +6,21 @@ Equal Earth .. versionadded:: 5.2.0 -+---------------------+--------------------------------------------------------+ -| **Classification** | Pseudo cylindrical | -+---------------------+--------------------------------------------------------+ -| **Available forms** | Forward and inverse, spherical projection | -+---------------------+--------------------------------------------------------+ -| **Defined area** | Global | -+---------------------+--------------------------------------------------------+ -| **Alias** | eqearth | -+---------------------+--------------------------------------------------------+ -| **Domain** | 2D | -+---------------------+--------------------------------------------------------+ -| **Input type** | Geodetic coordinates | -+---------------------+--------------------------------------------------------+ -| **Output type** | Projected coordinates | -+---------------------+--------------------------------------------------------+ ++---------------------+-----------------------------------------------------------+ +| **Classification** | Pseudo cylindrical | ++---------------------+-----------------------------------------------------------+ +| **Available forms** | Forward and inverse, spherical and ellipsoidal projection | ++---------------------+-----------------------------------------------------------+ +| **Defined area** | Global | ++---------------------+-----------------------------------------------------------+ +| **Alias** | eqearth | ++---------------------+-----------------------------------------------------------+ +| **Domain** | 2D | ++---------------------+-----------------------------------------------------------+ +| **Input type** | Geodetic coordinates | ++---------------------+-----------------------------------------------------------+ +| **Output type** | Projected coordinates | ++---------------------+-----------------------------------------------------------+ diff --git a/src/PJ_eqearth.c b/src/PJ_eqearth.c index 35795bbc5c..2633cbc870 100644 --- a/src/PJ_eqearth.c +++ b/src/PJ_eqearth.c @@ -18,7 +18,7 @@ Added ellipsoidal equations by Bojan Savric, 22 August 2018 #include "projects.h" -PROJ_HEAD(eqearth, "Equal Earth") "\n\tPCyl., Sph&Ell\n\th="; +PROJ_HEAD(eqearth, "Equal Earth") "\n\tPCyl., Sph&Ell"; #define A1 1.340264 #define A2 -0.081106 From 5e5787fe36430bbf24f3cc7a0bed48d3fe79f8d5 Mon Sep 17 00:00:00 2001 From: boja8320 Date: Wed, 22 Aug 2018 10:20:53 -0700 Subject: [PATCH 5/5] Adding +ellps parm --- docs/source/operations/projections/eqearth.rst | 2 ++ 1 file changed, 2 insertions(+) diff --git a/docs/source/operations/projections/eqearth.rst b/docs/source/operations/projections/eqearth.rst index cb0d607d25..cbc2ab154f 100644 --- a/docs/source/operations/projections/eqearth.rst +++ b/docs/source/operations/projections/eqearth.rst @@ -50,6 +50,8 @@ Parameters .. include:: ../options/lon_0.rst +.. include:: ../options/ellps.rst + .. include:: ../options/R.rst .. include:: ../options/x_0.rst