Implementation of Equal Earth projection (#1090)

Implement the Equal Earth projection (closes #1085) 

5 files changed, 99 insertions, 1 deletions

diff --git a/src/Makefile.am b/src/Makefile.am
index 4e1db2e9..3040a403 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -78,6 +78,7 @@ libproj_la_SOURCES = \
 	pj_tsfn.c pj_units.c pj_ctx.c pj_log.c pj_zpoly1.c rtodms.c \
 	vector1.c pj_release.c pj_gauss.c \
 	PJ_healpix.c PJ_natearth.c PJ_natearth2.c PJ_calcofi.c pj_fileapi.c \
+	PJ_eqearth.c \
 	\
 	pj_gc_reader.c pj_gridcatalog.c \
 	nad_cvt.c nad_init.c nad_intr.c emess.c emess.h \
diff --git a/src/PJ_eqearth.c b/src/PJ_eqearth.c
new file mode 100644
index 00000000..2ca5ca20
--- /dev/null
+++ b/src/PJ_eqearth.c
@@ -0,0 +1,94 @@
+/*
+Equal Earth is a projection inspired by the Robinson projection, but unlike
+the Robinson projection retains the relative size of areas. The projection
+was designed in 2018 by Bojan Savric, Tom Patterson and Bernhard Jenny.
+
+Publication:
+Bojan Savric, Tom Patterson & Bernhard Jenny (2018). The Equal Earth map
+projection, International Journal of Geographical Information Science,
+DOI: 10.1080/13658816.2018.1504949
+
+Port to PROJ by Juernjakob Dugge, 16 August 2018
+*/
+#define PJ_LIB__
+
+#include <math.h>
+
+#include "projects.h"
+
+PROJ_HEAD(eqearth, "Equal Earth") "\n\tPCyl., Sph.";
+
+#define A1 1.340264
+#define A2 -0.081106
+#define A3 0.000893
+#define A4 0.003796
+#define M (sqrt(3) / 2.0)
+
+#define MAX_Y 1.3173627591574  /* 90° latitude on a sphere with radius 1 */
+#define EPS 1e-11
+#define MAX_ITER 12
+
+static XY s_forward (LP lp, PJ *P) {           /* Spheroidal, forward */
+    XY xy = {0.0,0.0};
+    double phi, phi2, phi6;
+    (void) P;
+
+    phi = asin(M * sin(lp.phi));
+    phi2 = phi * phi;
+    phi6 = phi2 * phi2 * phi2;
+
+    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 */
+    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;
+        }
+    }
+
+    if( i == 0 ) {
+        pj_ctx_set_errno( P->ctx, PJD_ERR_NON_CONVERGENT );
+        return lp;
+    }
+
+    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;
+}
+
+
+PJ *PROJECTION(eqearth) {
+    P->es = 0;
+    P->inv = s_inverse;
+    P->fwd = s_forward;
+
+    return P;
+}
diff --git a/src/lib_proj.cmake b/src/lib_proj.cmake
index 249e8173..33f37bb0 100644
--- a/src/lib_proj.cmake
+++ b/src/lib_proj.cmake
@@ -67,6 +67,7 @@ SET(SRC_LIBPROJ_PJ
         PJ_eck5.c
         PJ_eqc.c
         PJ_eqdc.c
+        PJ_eqearth.c
         PJ_fahey.c
         PJ_fouc_s.c
         PJ_gall.c
diff --git a/src/makefile.vc b/src/makefile.vc
index ef39b084..d86649dd 100644
--- a/src/makefile.vc
+++ b/src/makefile.vc
@@ -38,7 +38,8 @@ pseudo = \
 	PJ_putp4p.obj PJ_putp5.obj PJ_putp6.obj PJ_robin.obj \
 	PJ_sts.obj PJ_urm5.obj PJ_urmfps.obj PJ_wag2.obj \
 	PJ_wag3.obj PJ_wink1.obj PJ_wink2.obj PJ_isea.obj \
-	PJ_calcofi.obj PJ_natearth.obj PJ_natearth2.obj PJ_times.obj
+	PJ_calcofi.obj PJ_natearth.obj PJ_natearth2.obj PJ_times.obj \
+	PJ_eqearth.obj
 
 support = \
 	aasincos.obj adjlon.obj bch2bps.obj bchgen.obj pj_gauss.obj \
diff --git a/src/pj_list.h b/src/pj_list.h
index 440f7972..0a8d9d1b 100644
--- a/src/pj_list.h
+++ b/src/pj_list.h
@@ -35,6 +35,7 @@ PROJ_HEAD(eck3, "Eckert III")
 PROJ_HEAD(eck4, "Eckert IV")
 PROJ_HEAD(eck5, "Eckert V")
 PROJ_HEAD(eck6, "Eckert VI")
+PROJ_HEAD(eqearth, "Equal Earth")
 PROJ_HEAD(eqc, "Equidistant Cylindrical (Plate Caree)")
 PROJ_HEAD(eqdc, "Equidistant Conic")
 PROJ_HEAD(euler, "Euler")