Implementation of Equal Earth projection (#1090)

Implement the Equal Earth projection (closes #1085) 

11 files changed, 247 insertions, 1 deletions

diff --git a/docs/plot/plotdefs.json b/docs/plot/plotdefs.json
index b61466f0..c494d6a7 100644
--- a/docs/plot/plotdefs.json
+++ b/docs/plot/plotdefs.json
@@ -308,6 +308,17 @@
         "type": "poly"
     },
     {
+        "filename": "eqearth.png",
+        "latmax": 90,
+        "latmin": -90,
+        "lonmax": 180,
+        "lonmin": -180,
+        "name": "eqearth",
+        "projstring": "+proj=eqearth",
+        "res": "low",
+        "type": "poly"
+    },
+    {
         "filename": "euler.png",
         "latmax": 90,
         "latmin": -90,
diff --git a/docs/source/operations/projections/eqearth.rst b/docs/source/operations/projections/eqearth.rst
new file mode 100644
index 00000000..e8b9417d
--- /dev/null
+++ b/docs/source/operations/projections/eqearth.rst
@@ -0,0 +1,59 @@
+.. _eqearth:
+
+********************************************************************************
+Equal Earth
+********************************************************************************
++---------------------+--------------------------------------------------------+
+| **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                                  |
++---------------------+--------------------------------------------------------+
+
+
+
+.. image:: ./images/eqearth.png
+   :alt:   Equal Earth
+
+
+The Equal Earth projection is intended for making world maps. 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 :cite:`Savric2018`.
+
+
+Usage
+###############################################################################
+
+The Equal Earth projection has no special options. Here is
+an example of an forward projection on a sphere with a radius of 1 m::
+
+    $ echo 122 47 | src/proj +proj=eqearth +R=1
+    1.55	0.89
+
+Parameters
+################################################################################
+
+.. note:: All parameters for the projection are optional.
+
+.. include:: ../options/lon_0.rst
+
+.. include:: ../options/R.rst
+
+.. include:: ../options/x_0.rst
+
+.. include:: ../options/y_0.rst
+
+Further reading
+################################################################################
+
+#. Bojan Savric, Tom Patterson & Bernhard Jenny (2018). `The Equal Earth map projection <https://www.researchgate.net/profile/Bojan_Savric2/publication/326879978_The_Equal_Earth_map_projection/links/5b69d0ae299bf14c6d951b77/The-Equal-Earth-map-projection.pdf>`_, International Journal of Geographical Information Science
diff --git a/docs/source/operations/projections/images/eqearth.png b/docs/source/operations/projections/images/eqearth.png
new file mode 100644
index 00000000..04d23c2d
--- /dev/null
+++ b/docs/source/operations/projections/images/eqearth.png
diff --git a/docs/source/operations/projections/index.rst b/docs/source/operations/projections/index.rst
index 77decd6c..a48eabe6 100644
--- a/docs/source/operations/projections/index.rst
+++ b/docs/source/operations/projections/index.rst
@@ -40,6 +40,7 @@ Projections map the spherical 3D space to a flat 2D space.
    eck6
    eqc
    eqdc
+   eqearth
    euler
    etmerc
    fahey
diff --git a/docs/source/references.bib b/docs/source/references.bib
index aceeb6cc..bd3089ee 100644
--- a/docs/source/references.bib
+++ b/docs/source/references.bib
@@ -196,6 +196,15 @@
   Year                     = {1993}
 }
 
+@Article{Savric2018,
+  Author                   = {B. Šavrič and T. Patterson and B. Jenny},
+  Title                    = {The Equal Earth map projection},
+  Journal                  = {International Journal of Geographical Information Science},
+  Year                     = {2018},
+  Publisher                = {Taylor & Francis},
+  Doi                      = {10.1080/13658816.2018.1504949},
+  Url                      = {https://www.researchgate.net/profile/Bojan_Savric2/publication/326879978_The_Equal_Earth_map_projection/links/5b69d0ae299bf14c6d951b77/The-Equal-Earth-map-projection.pdf}
+}
 
 @Article{Snyder1988,
   Author                   = {J. P. Snyder },
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")
diff --git a/test/gie/more_builtins.gie b/test/gie/more_builtins.gie
index 3a6c1650..2908fd63 100644
--- a/test/gie/more_builtins.gie
+++ b/test/gie/more_builtins.gie
@@ -408,6 +408,74 @@ accept       -81d00'00.000 34d32'30.000 0.0
 expect       2_000_000.000 561_019.077 0.0
 -------------------------------------------------------------------------------
 
+-------------------------------------------------------------------------------
+Some tests from PJ_eqearth.c
+-------------------------------------------------------------------------------
+operation +proj=eqearth +ellps=WGS84
+-------------------------------------------------------------------------------
+direction forward
+tolerance 1cm
+
+accept    0 0
+expect    0 0
+
+accept    -180 90
+expect    -10227908.09 8402320.16
+
+accept    0 90
+expect    0.00 8402320.16
+
+accept    180 90
+expect    10227908.09 8402320.16
+
+accept    180 45
+expect    14795421.79 5486671.72
+
+accept    180 0
+expect    17263256.84 0.00
+
+accept    -70 -31.2
+expect    -6244707.88 -3924893.29
+
+direction inverse
+
+accept    -6244707.88 -3924893.29
+expect    -70 -31.2
+
+accept    17263256.84 0.00
+expect    180 0
+
+accept    14795421.79 5486671.72
+expect    180 45
+
+accept    0 0
+expect    0 0
+
+accept    -10227908.09 8402320.16
+expect    -180 90
+
+accept    0.00 8402320.16
+expect    0 90
+
+accept    10227908.09 8402320.16
+expect    180 90
+
+operation +proj=eqearth +R=1
+direction inverse
+
+# coordinate in valid region
+accept    0 -1.3
+expect    0 -82.318
+
+# coordinate on edge
+accept    0 -1.3173627591574
+expect    0 -90
+
+# coordinate outside valid region, should be clamped
+accept    0 -1.4
+expect    0 -90
+
+-------------------------------------------------------------------------------
 
 -------------------------------------------------------------------------------
 run the few gie-builtin tests, which are currently either awkward or impossible