Skip to content

Commit

Permalink
Adding ellipsoidal equations for the Equal Earth (#1101)
Browse files Browse the repository at this point in the history
  • Loading branch information
@beuan @kbevers
beuan authored and kbevers committed Aug 22, 2018
1 parent 3e9ae6f commit b2fe227
Show file tree
Hide file tree
Showing 3 changed files with 155 additions and 32 deletions.
32 changes: 17 additions & 15 deletions docs/source/operations/projections/eqearth.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 |
+---------------------+-----------------------------------------------------------+



Expand Down Expand Up @@ -50,6 +50,8 @@ Parameters

.. include:: ../options/lon_0.rst

.. include:: ../options/ellps.rst

.. include:: ../options/R.rst

.. include:: ../options/x_0.rst
Expand Down
105 changes: 88 additions & 17 deletions src/PJ_eqearth.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 <errno.h>
#include <math.h>

#include "projects.h"

PROJ_HEAD(eqearth, "Equal Earth") "\n\tPCyl., Sph.";
PROJ_HEAD(eqearth, "Equal Earth") "\n\tPCyl., Sph&Ell";

#define A1 1.340264
#define A2 -0.081106
Expand All @@ -28,28 +30,59 @@ PROJ_HEAD(eqearth, "Equal Earth") "\n\tPCyl., Sph.";
#define EPS 1e-11
#define MAX_ITER 12

static XY s_forward (LP lp, PJ *P) { /* Spheroidal, forward */
struct pj_opaque {
double qp;
double rqda;
double *apa;
};

static XY e_forward (LP lp, PJ *P) { /* Ellipsoidal/Spheroidal, forward */
XY xy = {0.0,0.0};
double phi, phi2, phi6;
(void) P;
struct pj_opaque *Q = P->opaque;
double sbeta;
double psi, psi2, psi6;

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);
}

phi = asin(M * sin(lp.phi));
phi2 = phi * phi;
phi6 = phi2 * phi2 * phi2;
/* Equal Earth projection */
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;

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));
return xy;
}


static LP s_inverse (XY xy, PJ *P) { /* Spheroidal, 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;
double beta;
int i;
(void) P;

/* make sure y is inside valid range */
/* 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) {
Expand All @@ -75,20 +108,58 @@ static LP s_inverse (XY xy, PJ *P) { /* Spheroidal, inverse */
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);
lp.phi = asin(sin(yc) / M);


/* Latitude */
beta = asin(sin(yc) / M);
if (P->es != 0.0) {
/* Ellipsoidal case, converting authalic latitude */
lp.phi = pj_authlat(beta, Q->apa);
}
else {
/* Spheroidal case, using latitude */
lp.phi = beta;
}

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 */
}
else {
Q->rqda = 1.0;
}

P->fwd = e_forward;
P->inv = e_inverse;

return P;
}
50 changes: 50 additions & 0 deletions test/gie/more_builtins.gie
View file
Original file line number Diff line number Diff line change
Expand Up @@ -419,6 +419,56 @@ 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
tolerance 1cm

accept 0 0
expect 0 0

accept -180 90
expect -10227908.09 8402320.16

Expand Down

0 comments on commit b2fe227

Please sign in to comment.