Add the Equal Earth projection support following what was done for Spilhaus

doc/rst/source/changes.rst

@@ -6,10 +6,17 @@
 Changelog
 =========
 
-New Features in GMT 6.7
-=======================
+New Features in GMT 6.7.1
+=========================
+
+* New Equal Earth map projection, selected with the **PROJ** syntax **-J+proj=eqearth**. GMT now draws the proper
+  pseudo-cylindrical frame, annotations, ticks and gridlines for it, as it does for the Robinson projection.
+  See :ref:`Equal Earth <-Jeqearth>`
 
-GMT 6.7 brings two new programs (a tsunami modeling supplement and a stand-alone updater), several new
+New Features in GMT 6.7.0
+=========================
+
+GMT 6.7.0 brings two new programs (a tsunami modeling supplement and a stand-alone updater), several new
 plotting capabilities built on smooth color gradients (Gouraud shading), Bezier curves in the custom symbol
 language, and a large number of bug fixes across the whole code base.
 
@@ -32,16 +39,16 @@ language, and a large number of bug fixes across the whole code base.
 * Change the default GMT data server to ``remotedata.generic-mapping-tools.org`` (`#9102 <https://github.com/GenericMappingTools/gmt/pull/9102>`_)
 * New documentation chapter on how to create custom supplement libraries (`#9029 <https://github.com/GenericMappingTools/gmt/pull/9029>`_)
 
-New Modules and Programs in GMT 6.7
+New Modules and Programs in GMT 6.7.0
 -----------------------------------
 
 * New supplement :doc:`/supplements/nswing/nswing`: a finite-differences tsunami (shallow-water) modeling
   module (`#9093 <https://github.com/GenericMappingTools/gmt/pull/9093>`_)
 * New stand-alone program :doc:`gmt_update` that updates a GMT installation between release cycles.  It
   downloads and verifies (SHA-256) the new files.  Currently tested on Windows only (`#9032 <https://github.com/GenericMappingTools/gmt/pull/9032>`_)
 
-Enhancements in GMT 6.7
------------------------
+Enhancements in GMT 6.7.0
+-------------------------
 
 * :doc:`plot`, :doc:`plot3d`: Add **-G+g** for vertex-based color gradients (Gouraud shading) in polygons.
   Vertex colors may be given as *r g b*, as color names/hex/HSV/CMYK strings, or as *z* values mapped through

doc/rst/source/reference/map-projections.rst

@@ -1103,7 +1103,8 @@ Miscellaneous projections
 GMT supports eight common projections for global presentation of data or models. These are the :ref:`Hammer <-Jh>`,
 :ref:`Mollweide <-Jw>`, :ref:`Winkel Tripel <-Jr>`, :ref:`Robinson <-Jn>`, :ref:`Eckert IV and VI <-Jk>`,
 :ref:`Sinusoidal <-Ji>`, and :ref:`Van der Grinten <-Jv>` projections. Due to the small scale used for global maps these
-projections all use the spherical approximation rather than more elaborate elliptical formulae.
+projections all use the spherical approximation rather than more elaborate elliptical formulae. In addition, the
+:ref:`Equal Earth <-Jeqearth>` projection is available through the **PROJ** string syntax of **-J**.
 
 In all cases, the specification of the central meridian can be skipped. The default is the middle of the longitude
 range of the plot, specified by the (**-R**) option.
@@ -1389,6 +1390,42 @@ Centered on the Dateline, the example below was created by this command:
 
    World map using the Van der Grinten projection.
 
+.. _-Jeqearth:
+
+Equal Earth projection (**-J+proj=eqearth**)
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+**Syntax**
+
+    **-J+proj=eqearth**\ [**+lon_0=**\ *lon0*]\ **+width=**\ *width*\|\ **+scale=**\ *scale*
+
+**Parameters**
+
+- The central meridian (*lon0*) [0].
+- The map *width* or *scale* in :ref:`plot-units <plt-units>`.
+
+**Description**
+
+The Equal Earth projection, presented by Bojan Šavrič, Tom Patterson and Bernhard Jenny in 2018, is an equal-area
+pseudo-cylindrical projection whose overall look was designed to resemble the (not equal-area) Robinson projection.
+Parallels are straight lines, the central meridian is a straight line half as long as the Equator, and all other
+meridians are curved. It has no GMT 1-letter code: it is done by `PROJ <https://proj.org/en/stable/operations/projections/eqearth.html>`_
+and is selected with the **PROJ** string syntax of the **-J** option. Because **PROJ** does the projection, the
+parameters are the **PROJ** ones (e.g., **+lon_0**, **+ellps**), to which GMT adds **+width** or **+scale** to set
+the map dimension. Note that the whole string must be quoted if you separate the **+**\ *keys* with spaces.
+
+**Example**
+
+Centered on Greenwich, the example below was created by this command:
+
+.. literalinclude:: /_verbatim/GMT_eqearth.txt
+
+.. figure:: /_images/GMT_eqearth.*
+   :width: 500 px
+   :align: center
+
+   World map using the Equal Earth projection.
+
 Footnotes
 ---------
 

doc/scripts/GMT_eqearth.sh

@@ -0,0 +1,2 @@
+#!/usr/bin/env bash
+gmt coast -Rd -J+proj=eqearth+width=12c -Bg -Dc -A10000 -Ggoldenrod -Ssnow2 -view GMT_eqearth

src/gmt_init.c

@@ -18784,7 +18784,16 @@ int gmt_parse_common_options (struct GMT_CTRL *GMT, char *list, char option, cha
 				if (GMT->current.gdal_read_in.hCT_inv) OCTDestroyCoordinateTransformation(GMT->current.gdal_read_in.hCT_inv);
 				GMT->current.gdal_read_in.hCT_fwd = gmt_OGRCoordinateTransformation(GMT, source, dest);
 				GMT->current.gdal_read_in.hCT_inv = gmt_OGRCoordinateTransformation(GMT, dest, source);
-				GMT->current.proj.projection = strstr(dest, "spilhaus") ? GMT_PROJ4_SPILHAUS : GMT_PROJ4_PROJS;	/* Special case for spilhaus */
+				if (strstr(dest, "spilhaus"))		/* Special case for spilhaus */
+					GMT->current.proj.projection = GMT_PROJ4_SPILHAUS;
+				else if (strstr(dest, "eqearth")) {	/* Special case for the Equal Earth projection */
+					GMT->current.proj.projection = GMT_PROJ4_EQEARTH;
+					/* gmtinit_parse_J_option does this for the projections GMT knows, but it is not called for
+					 * the PROJ-only ones, so a non-global -R would not be recognized as geographic */
+					gmt_set_geographic(GMT, GMT_IN);	/* This may be overridden by mapproject -I */
+				}
+				else
+					GMT->current.proj.projection = GMT_PROJ4_PROJS;
 				GMT->common.J.active = true;
 				if (GMT->current.gdal_read_in.hCT_fwd == NULL || GMT->current.gdal_read_in.hCT_inv == NULL)
 					error = 1;

src/gmt_map.c

@@ -5330,6 +5330,14 @@ GMT_LOCAL int gmtmap_init_polyconic (struct GMT_CTRL *GMT, bool *search) {
 			GMT->current.proj.inv = &gmt_proj4_inv;
 			GMT->current.proj.scale[GMT_X] = GMT->current.proj.scale[GMT_Y] = GMT->current.proj.proj4_scl;
 			GMT->current.map.n_lon_nodes = 360;	GMT->current.map.n_lat_nodes = 180;
+			if (GMT->current.proj.projection == GMT_PROJ4_EQEARTH) {	/* Need the central meridian and world status as for the GMT pseudo-cylindricals */
+				char *pch = strstr(GMT->common.J.proj4string, "+lon_0=");
+				GMT->current.proj.central_meridian = (pch) ? atof(&pch[7]) : 0.0;
+				GMT->current.map.is_world = gmt_M_360_range(GMT->common.R.wesn[XLO], GMT->common.R.wesn[XHI]);
+				/* Must be in place before the xy search below since those wrap the GDAL transforms to honor the central meridian */
+				GMT->current.proj.fwd = &gmtproj_proj4_pcyl_fwd;
+				GMT->current.proj.inv = &gmtproj_proj4_pcyl_inv;
+			}
 			if (GMT->common.R.oblique) {
 				gmt_proj4_fwd(GMT, GMT->common.R.wesn[XLO], GMT->common.R.wesn[YLO], &xmin, &ymin);
 				gmt_proj4_fwd(GMT, GMT->common.R.wesn[XHI], GMT->common.R.wesn[YHI], &xmax, &ymax);
@@ -5347,8 +5355,17 @@ GMT_LOCAL int gmtmap_init_polyconic (struct GMT_CTRL *GMT, bool *search) {
 				GMT->current.map.crossing = &gmtmap_wesn_crossing;
 				GMT->current.map.overlap = &gmtmap_wesn_overlap;
 				GMT->current.map.clip = &gmt_map_wesn_clip;
-				GMT->current.map.left_edge = &gmtmap_left_rect;
-				GMT->current.map.right_edge = &gmtmap_right_rect;
+				if (GMT->current.proj.projection == GMT_PROJ4_EQEARTH) {	/* Pseudo-cylindrical, so curved W/E boundaries as in Robinson */
+					GMT->current.map.left_edge = &gmtproj_left_proj4_pcyl;
+					GMT->current.map.right_edge = &gmtproj_right_proj4_pcyl;
+					GMT->current.map.parallel_straight = 1;	/* Parallels are straight horizontal lines */
+					/* Unless MAP_ANNOT_OBLIQUE is set manually, we must add normal ticks for this pole-is-line projection */
+					if (!GMT->current.setting.map_annot_oblique_set) GMT->current.setting.map_annot_oblique |= GMT_OBL_ANNOT_NORMAL_TICKS;
+				}
+				else {
+					GMT->current.map.left_edge = &gmtmap_left_rect;
+					GMT->current.map.right_edge = &gmtmap_right_rect;
+				}
 				GMT->current.map.frame.horizontal = 2;
 			}
 			gmtmap_setinfo(GMT, xmin, xmax, ymin, ymax, GMT->current.proj.proj4_scl);
@@ -5367,6 +5384,11 @@ GMT_LOCAL int gmtmap_init_polyconic (struct GMT_CTRL *GMT, bool *search) {
 	/* Now we only have to replace the pointers to the FWD and INV transform functions */
 	GMT->current.proj.fwd = &gmt_proj4_fwd;
 	GMT->current.proj.inv = &gmt_proj4_inv;
+	if (GMT->current.proj.projection == GMT_PROJ4_EQEARTH) {	/* So that the frame/annotation/gridline machinery recognizes this projection */
+		GMT->current.proj.projection_GMT = GMT->current.proj.projection;
+		GMT->current.proj.fwd = &gmtproj_proj4_pcyl_fwd;	/* These wrap the GDAL transforms to honor the central meridian */
+		GMT->current.proj.inv = &gmtproj_proj4_pcyl_inv;
+	}
 	return error;
 }
 
@@ -7657,6 +7679,7 @@ double gmt_half_map_width (struct GMT_CTRL *GMT, double y) {
 		case GMT_ROBINSON:
 		case GMT_ECKERT4:
 		case GMT_ECKERT6:
+		case GMT_PROJ4_EQEARTH:
 			if (!GMT->common.R.oblique && GMT->current.map.is_world)
 				half_width = gmtlib_right_boundary (GMT, y) - GMT->current.map.half_width;
 			else
@@ -9102,6 +9125,9 @@ uint64_t gmt_map_clip_path(struct GMT_CTRL *GMT, double **x, double **y, bool *d
 			case GMT_PROJ4_SPILHAUS:
 				np = 4;
 				break;
+			case GMT_PROJ4_EQEARTH:	/* Pseudo-cylindrical: curved W/E boundaries, straight N/S, as for Robinson */
+				np = 2 * GMT->current.map.n_lat_nodes + 2;
+				break;
 			default:
 				GMT_Report (GMT->parent, GMT_MSG_ERROR, "Bad case in gmt_map_clip_path (%d)\n", GMT->current.proj.projection_GMT);
 				np = 0;
@@ -9235,6 +9261,7 @@ uint64_t gmt_map_clip_path(struct GMT_CTRL *GMT, double **x, double **y, bool *d
 			case GMT_MOLLWEIDE:
 			case GMT_SINUSOIDAL:
 			case GMT_ROBINSON:
+			case GMT_PROJ4_EQEARTH:
 				for (i = j = 0; i <= GMT->current.map.n_lat_nodes; i++, j++) {	/* Right */
 					lat = (i == GMT->current.map.n_lat_nodes) ? GMT->common.R.wesn[YHI] : GMT->common.R.wesn[YLO] + i * GMT->current.map.dlat;
 					gmt_geo_to_xy (GMT, GMT->common.R.wesn[XHI], lat, &work_x[j], &work_y[j]);

src/gmt_plot.c

@@ -2800,7 +2800,8 @@ GMT_LOCAL void gmtplot_map_annotate (struct GMT_CTRL *GMT, struct PSL_CTRL *PSL,
 		proj_A = (GMT->current.proj.projection_GMT == GMT_MERCATOR || GMT->current.proj.projection_GMT == GMT_OBLIQUE_MERC ||
 			GMT->current.proj.projection_GMT == GMT_WINKEL || GMT->current.proj.projection_GMT == GMT_ECKERT4 || GMT->current.proj.projection_GMT == GMT_ECKERT6 ||
 			GMT->current.proj.projection_GMT == GMT_ROBINSON || GMT->current.proj.projection_GMT == GMT_CYL_EQ || GMT->current.proj.projection_GMT == GMT_CYL_STEREO ||
-			GMT->current.proj.projection_GMT == GMT_CYL_EQDIST || GMT->current.proj.projection_GMT == GMT_MILLER || GMT->current.proj.projection_GMT == GMT_LINEAR);
+			GMT->current.proj.projection_GMT == GMT_CYL_EQDIST || GMT->current.proj.projection_GMT == GMT_MILLER || GMT->current.proj.projection_GMT == GMT_LINEAR ||
+			GMT->current.proj.projection_GMT == GMT_PROJ4_EQEARTH);
 		proj_B = (GMT->current.proj.projection_GMT == GMT_HAMMER || GMT->current.proj.projection_GMT == GMT_MOLLWEIDE ||
 			GMT->current.proj.projection_GMT == GMT_SINUSOIDAL);
 		if (gmt_M_is_conical (GMT) && gmt_M_360_range (w, e))	/* Special case since 360 longitudes do not form a circle but a pacman shape */
@@ -3091,6 +3092,7 @@ GMT_LOCAL void gmtplot_map_boundary (struct GMT_CTRL *GMT) {
 		case GMT_ECKERT6:
 		case GMT_ROBINSON:
 		case GMT_VANGRINTEN:
+		case GMT_PROJ4_EQEARTH:
 			gmtplot_basic_map_boundary (GMT, PSL, w, e, s, n);
 			break;
 		case GMT_PROJ4_SPILHAUS:

src/gmt_proj.c

@@ -2646,6 +2646,58 @@ double gmtproj_right_robinson (struct GMT_CTRL *GMT, double y) {
 }
 #endif
 
+/* -J+proj=eqearth EQUAL EARTH PROJECTION (and, in principle, any other PROJ pseudo-cylindrical projection)
+ *
+ * The projection itself is done by PROJ (via GDAL), so we have no analytical expressions here.  What GMT needs
+ * on top of the transforms are the left and right map boundaries as functions of the plot y-coordinate so that
+ * the frame, annotations, ticks and gridlines come out curved as they should (like the -Jn Robinson case) and
+ * not as a plain rectangle.  We get them by inverting (x = central meridian, y) to a latitude and then project
+ * the west/east boundary at that latitude.  Both steps go through the same PROJ transforms used for the data,
+ * so the border is by construction consistent with what is plotted.
+ */
+
+GMT_LOCAL void gmtproj_proj4_pcyl_fwd(struct GMT_CTRL *GMT, double lon, double lat, double *x, double *y) {
+	/* PROJ first normalizes the geographic longitude to [-180,180] and only then subtracts lon_0.  For a global
+	   map whose west and east boundaries are the same meridian (e.g., -Rg with +lon_0=180) both boundaries then
+	   collapse onto the western edge of the map.  We therefore wind the longitude relative to the central
+	   meridian ourselves (GMT keeps -180 and +180 apart) and pull the antimeridian a hair inside the domain. */
+	double dlon = lon;
+
+	gmt_M_wind_lon(GMT, dlon);	/* Now dlon is in [-180,180] relative to the central meridian */
+	if (dlon >= 180.0) dlon = 180.0 - GMT_PROJ_CONV_LIMIT;
+	else if (dlon <= -180.0) dlon = -180.0 + GMT_PROJ_CONV_LIMIT;
+	gmt_proj4_fwd(GMT, GMT->current.proj.central_meridian + dlon, lat, x, y);
+}
+
+GMT_LOCAL void gmtproj_proj4_pcyl_inv(struct GMT_CTRL *GMT, double *lon, double *lat, double x, double y) {
+	/* The inverse: PROJ hands us a longitude in [-180,180] so put it back in the map's own longitude range */
+	gmt_proj4_inv(GMT, lon, lat, x, y);
+	if (gmt_M_is_dnan(*lon) || gmt_M_is_dinf(*lon)) return;
+	while (*lon < GMT->current.proj.central_meridian - 180.0) *lon += 360.0;
+	while (*lon > GMT->current.proj.central_meridian + 180.0) *lon -= 360.0;
+}
+
+GMT_LOCAL double gmtproj_proj4_pcyl_edge(struct GMT_CTRL *GMT, double y, double lon) {
+	/* Return the plot x-coordinate of meridian lon at plot y-coordinate y */
+	double lon_c, lat, x, y_dummy;
+
+	/* proj.origin[GMT_X] is the plot x where the projected x is zero, i.e., the central meridian */
+	gmt_xy_to_geo (GMT, &lon_c, &lat, GMT->current.proj.origin[GMT_X], y);
+	if (gmt_M_is_dnan(lat) || gmt_M_is_dinf(lat)) return GMT->session.d_NaN;
+	if (lat >  90.0) lat =  90.0;	/* Round-off may place us just outside the domain */
+	if (lat < -90.0) lat = -90.0;
+	gmt_geo_to_xy(GMT, lon, lat, &x, &y_dummy);
+	return (x);
+}
+
+GMT_LOCAL double gmtproj_left_proj4_pcyl(struct GMT_CTRL *GMT, double y) {
+	return (gmtproj_proj4_pcyl_edge(GMT, y, GMT->common.R.wesn[XLO]));
+}
+
+GMT_LOCAL double gmtproj_right_proj4_pcyl(struct GMT_CTRL *GMT, double y) {
+	return (gmtproj_proj4_pcyl_edge(GMT, y, GMT->common.R.wesn[XHI]));
+}
+
 /* -JI SINUSOIDAL EQUAL AREA PROJECTION */
 
 GMT_LOCAL void gmtproj_vsinusoidal (struct GMT_CTRL *GMT, double lon0) {

src/gmt_project.h

@@ -123,6 +123,7 @@ enum gmt_enum_misc {GMT_MOLLWEIDE = 400,
 #define gmt_M_is_proj4(C) (C->current.proj.projection / 100 == 5)
 enum gmt_enum_allprojs {GMT_PROJ4_PROJS = 500,
 	GMT_PROJ4_SPILHAUS,
+	GMT_PROJ4_EQEARTH,
 };
 
 /*! The various GMT measurement units */

test/baseline/psbasemap.dvc

@@ -1,6 +1,6 @@
 outs:
-- md5: 499cd79959131cb42466afb2ba0bda90.dir
-  nfiles: 58
+- md5: 8272438084e659d6dee76d6112f6842c.dir
+  nfiles: 59
   path: psbasemap
   hash: md5
-  size: 2981097
+  size: 3023834

test/psbasemap/eqearth.sh

@@ -0,0 +1,8 @@
+#!/usr/bin/env bash
+# Test frame, annotations, ticks and gridlines for the PROJ Equal Earth projection.
+# Global with the default central meridian, global centered on the Dateline, and a subregion.
+gmt begin eqearth ps
+	gmt basemap -Rd -J+proj=eqearth+width=11c -Bxa60g30 -Bya30g30 -BWSen
+	gmt basemap -Rg -J"+proj=eqearth +lon_0=180 +width=11c" -Bxa60g30 -Bya30g30 -BWSen -Y7c
+	gmt basemap -R-100/40/-20/70 -J+proj=eqearth+width=11c -Bafg -BWSen -Y7c
+gmt end show