Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Adding ellipsoidal equations for the Equal Earth #1101

Merged
merged 5 commits into from
Aug 22, 2018
Merged

Adding ellipsoidal equations for the Equal Earth #1101

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
73ca7e6
Adding Equal Earth ellipsoidal math
beuan Aug 22, 2018
67b9ea5
Adding tolerance to the tests...
beuan Aug 22, 2018
8197803
Joining ellipsoidal and sphere functions
beuan Aug 22, 2018
670e874
Adding doc changes
beuan Aug 22, 2018
5e5787f
Adding +ellps parm
beuan Aug 22, 2018
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
32 changes: 17 additions & 15 deletions docs/source/operations/projections/eqearth.rst
View file
Open in desktop
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
Open in desktop
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) {
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

A little bit of nitpicking here: For clarity, I would prefer to have block-scoped comments precede the block they refer to:

/* Ellipsoidal case, converting to authalic latitude */
if (P->es != 0.0) {
    sbeta = pj_qsfn(sin(lp.phi), P->e, 1.0 - P->es) / Q->qp;

    /* Rounding error. */
    if (fabs(sbeta) > 1) 
        sbeta = sbeta > 0 ? 1 : -1;
}

/* Spheroidal case, using latitude */
else 
    sbeta = sin(lp.phi);

Or with even less branching:

sbeta = sin(lp.phi);

/* In the ellipsoidal case, we convert sbeta to authalic latitude */
if (P->es != 0.0) {
    sbeta = pj_qsfn(sbeta, P->e, 1.0 - P->es) / Q->qp;

    /* Rounding error. */
    if (fabs(sbeta) > 1) 
        sbeta = sbeta > 0 ? 1 : -1;
}

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sorry Thomas, I merged this while you were writing...

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Don't worry - although I think my suggestion has some merit in terms of clarity, it also is mostly cosmetical and may be simply a matter of personal preference

Copy link
Contributor Author

@beuan beuan Aug 22, 2018 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I happy to do the change and another PR if you would like to...

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@beuan: Well - I wrote the comment because I would prefer to minimize branching and have block-scope comments precede the block.

But it's your code, and it's already merged, so please only make the change if you personally find my (cosmetical) arguments highly convincing, and would like to see the change - in which case I will obviously feel honoured, that you found the comment convincing and/or useful :-)

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@busstoptaktik No problem. Any additional inputs are welcome.
I have changes on the way. I just need to figure out how to make a pull request without any "merge conflicts"

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

#1102

/* 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;
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It would be nice to have a test case that exercises this condition.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I am not sure I can hit this test case. The best testing case is probably phi = +/-90. In this case pj_qsfn() function returns +/-Q->qp and for some rounding error it might hit this case. Those tests are already in and looked they are passing.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Using phi > 90 would not help, because sin function would already "fix" this.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

All right, we'll live without a test for that

}
}
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
Open in desktop
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