Let us have the Spilhaus projection accessible from pscoast.c (#8745)

* Let us have the Spilhaus projection accessible from pscoast.c

For this, it was necessary to add the -g option to pscoast and make a gap detection in gmt_geo_to_xy_line().

We can not yet do the interrupted graticule nor plot a frame.

* Use a default -g in pscoast when proj is Spilhaus

* Fix issues with scale and frame drawing. Now we can use -B

* Document -g in pscoast. Make Spilhaus parameters as #defines.

Shut up (avoid them) PROJ warnings.

Author: joa-quim

doc/rst/source/coast.rst

@@ -36,6 +36,7 @@ Synopsis
 [ |SYN_OPT-Y| ]
 [ |SYN_OPT-bo| ]
 [ |SYN_OPT-d| ]
+[ |SYN_OPT-g| ]
 [ |SYN_OPT-p| ]
 [ |SYN_OPT-t| ]
 [ |SYN_OPT--| ]
@@ -290,6 +291,15 @@ Optional Arguments
 .. |Add_-d| unicode:: 0x20 .. just an invisible code
 .. include:: explain_-d.rst_
 
+.. _-g:
+
+**-gD**\ *dist*
+    Short version of the global **-g** option. Here it is used to set the _dist_ distance, in map units,
+		higher then which we consider to have a gap. Useful for the Spilhaus projection (though we automatically set it)
+		and when line wrapping on dateline was not correctly detected. As an example of this, see for instance
+    the old issue `#3667 <https://github.com/GenericMappingTools/gmt/issues/3667>`_ that can be fixed
+    by adding **-gD**\ 1 to the command line.
+
 .. |Add_perspective| unicode:: 0x20 .. just an invisible code
 .. include:: explain_perspective.rst_
 

doc/rst/source/pscoast.rst

@@ -37,6 +37,7 @@ Synopsis
 [ |SYN_OPT-X| ]
 [ |SYN_OPT-Y| ]
 [ |SYN_OPT-bo| ]
+[ |SYN_OPT-g| ]
 [ |SYN_OPT-p| ]
 [ |SYN_OPT-t| ]
 [ |SYN_OPT--| ]

src/gmt_init.c

@@ -18546,7 +18546,7 @@ 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      = GMT_PROJ4_PROJS;		/* This now make it use the proj4 lib */
+				GMT->current.proj.projection = strstr(dest, "spilhaus") ? GMT_PROJ4_SPILHAUS : GMT_PROJ4_PROJS;	/* Special case for spilhaus */
 				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

@@ -3029,7 +3029,8 @@ GMT_LOCAL void gmtplot_map_boundary (struct GMT_CTRL *GMT) {
 	if ((way = gmtlib_adjust_we_if_central_lon_set (GMT, &w, &e)))
 		GMT_Report (GMT->parent, GMT_MSG_DEBUG, "W/E boundaries shifted by %d\n", way * 360);
 
-	PSL_comment (PSL, "Start of map frame\n");
+	if (!PSL->internal.comments) PSL_command(PSL, "\n%%Start of map frame\n");
+	PSL_comment(PSL, "Start of map frame\n");
 
 	gmt_setpen (GMT, &GMT->current.setting.map_frame_pen);
 	PSL_setcolor (PSL, GMT->current.setting.map_frame_pen.rgb, PSL_IS_STROKE);
@@ -3088,8 +3089,13 @@ GMT_LOCAL void gmtplot_map_boundary (struct GMT_CTRL *GMT) {
 		case GMT_VANGRINTEN:
 			gmtplot_basic_map_boundary (GMT, PSL, w, e, s, n);
 			break;
+		case GMT_PROJ4_SPILHAUS:
+			gmtplot_linear_map_boundary(GMT, PSL, w, e, s, n);
+			break;
 	}
-	PSL_comment (PSL, "End of map frame\n");
+
+	PSL_comment(PSL, "End of map frame\n");
+	if (!PSL->internal.comments) PSL_command(PSL, "\n%%End of map fram\n");
 }
 
 /* gmt_map_basemap will create a basemap for the given area.

src/gmt_plot.c

@@ -90,6 +90,20 @@
 #define GMT_PROJ_CONV_LIMIT	1e-9
 #define gmt_M_proj_is_zero(x) (fabs (x) < GMT_PROJ_CONV_LIMIT)
 
+/* The Spilhaus proj is f. The points bellow were obtained by projecting a global grid at 1
+   deg resolution with gdalwarp and inverting the UL and LR coordinates. But even those
+   had to be moved a bit to inside (by trial en error). Couldn't do the same to the other
+   two corners because permanently got Inf's. This could probably still be refined a tinny bit.
+   They are good for global maps, but will be in excess for some local -R. However, given
+   the complexity of the projection, it is not worth the effort to determine them for all
+   cases. Anyway, the beauty of the Spilhaus projection is when applyied to the whole world,
+   so it is not a big deal.
+*/
+#define GMT_PROJ_SPILH_LON_UL	-113.0447807067042
+#define GMT_PROJ_SPILH_LAT_UL	49.56674656682158
+#define GMT_PROJ_SPILH_LON_LR	-113.06804976730443
+#define GMT_PROJ_SPILH_LAT_LR	49.553963054590234
+
 GMT_LOCAL double gmtproj_robinson_spline (struct GMT_CTRL *GMT, double xp, double *x, double *y, double *c) {
 	/* Returns the interpolated value y(xp) from the Robinson coefficients */
 

src/gmt_proj.c

@@ -120,7 +120,10 @@ enum gmt_enum_misc {GMT_MOLLWEIDE = 400,
 	GMT_WINKEL};
 
 /* All projections from proj.4 lib */
-enum gmt_enum_allprojs {GMT_PROJ4_PROJS = 500};
+#define gmt_M_is_proj4(C) (C->current.proj.projection / 100 == 5)
+enum gmt_enum_allprojs {GMT_PROJ4_PROJS = 500,
+	GMT_PROJ4_SPILHAUS,
+};
 
 /*! The various GMT measurement units */
 enum gmt_enum_units {GMT_IS_METER = 0,

src/gmt_project.h

@@ -1518,7 +1518,7 @@ EXTERN_MSC int GMT_grdimage(void *V_API, int mode, void *args) {
 
 	/* Determine if grid/image is to be projected */
 	need_to_project = (gmt_M_is_nonlinear_graticule (GMT) || Ctrl->E.dpi > 0);
-	if (need_to_project && GMT->current.proj.projection == GMT_PROJ4_PROJS) {
+	if (need_to_project && gmt_M_is_proj4(GMT)) {
 		if (GMT->current.proj.is_proj4)
 			if (strstr(GMT->common.J.proj4string, "longlat") || strstr(GMT->common.J.proj4string, "lonlat") || strstr(GMT->common.J.proj4string, "latlon"))
 				need_to_project = false;

src/grdimage.c

@@ -59,7 +59,7 @@
 #define THIS_MODULE_PURPOSE	"Plot continents, countries, shorelines, rivers, and borders"
 #define THIS_MODULE_KEYS	">?},>DE-lL"
 #define THIS_MODULE_NEEDS	"JR"
-#define THIS_MODULE_OPTIONS "->BJKOPRUVXYbdptxy" GMT_OPT("Zc")
+#define THIS_MODULE_OPTIONS "->BJKOPRUVXYbdgptxy" GMT_OPT("Zc")
 
 #define LAKE	0
 #define RIVER	1
@@ -199,10 +199,10 @@ static int usage (struct GMTAPI_CTRL *API, int level) {
 	if (level == GMT_MODULE_PURPOSE) return (GMT_NOERROR);
 	GMT_Usage (API, 0, "usage: %s %s %s [%s] [%s] [-C<fill>[+l|r]] [-D<resolution>[+f]] [-E%s] "
 		"[-G[<fill>]] [-F%s] [-I<feature>[/<pen>]] %s [-L%s] [-M] [-N<feature>[/<pen>]] %s%s[-Q] [-S[<fill>]] "
-		"[-Td%s] [-Tm%s] [%s] [%s] [-W[<feature>/][<pen>]] [%s] [%s] [%s] %s[%s] [%s] [%s] [%s]%s [%s]\n",
+		"[-Td%s] [-Tm%s] [%s] [%s] [-W[<feature>/][<pen>]] [%s] [%s] [%s] [%s] %s[%s] [%s] [%s] [%s]%s [%s]\n",
 		name, GMT_J_OPT, GMT_Rgeoz_OPT, GMT_A_OPT, GMT_B_OPT, DCW_OPT, GMT_PANEL, API->K_OPT, GMT_SCALE, API->O_OPT,
 		API->P_OPT, GMT_TROSE_DIR, GMT_TROSE_MAG, GMT_U_OPT, GMT_V_OPT, GMT_X_OPT, GMT_Y_OPT, GMT_bo_OPT, API->c_OPT,
-		GMT_do_OPT, GMT_p_OPT, GMT_t_OPT, GMT_colon_OPT, dbg, GMT_PAR_OPT);
+		GMT_do_OPT, GMT_g_OPT, GMT_p_OPT, GMT_t_OPT, GMT_colon_OPT, dbg, GMT_PAR_OPT);
 
 	if (level == GMT_SYNOPSIS) return (GMT_MODULE_SYNOPSIS);
 
@@ -268,6 +268,9 @@ static int usage (struct GMTAPI_CTRL *API, int level) {
 	GMT_Usage (API, 3, "4: Lake in island in lake shores.");
 	GMT_Usage (API, -2, "Note: When feature-specific pens are used, those not set are deactivated.");
 	GMT_Option (API, "X,bo,c,do,p,t");
+	GMT_Usage (API, 1, "\n-gD<dist> Short version of the global -g option. Here it is used to set the distance, in map "
+		"units, above which we consider to have a gap. Useful for the Spilhaus projection (though we automatically set it) "
+		"and when line wrapping on dateline was not correctly detected.");
 #ifdef DEBUG
 	GMT_Usage (API, 1, "\n-+<bin> (repeatable up to 16 times)");
 	GMT_Usage (API, -2, "Plot only the specified bins (debug option).");
@@ -533,6 +536,14 @@ static int parse (struct GMT_CTRL *GMT, struct PSCOAST_CTRL *Ctrl, struct GMT_OP
 		}
 	}
 
+	/* If Spilhaus projection and no -g used, set a default value. */
+	if (!GMT->common.g.active && GMT->current.proj.projection == GMT_PROJ4_SPILHAUS) {
+		gmt_parse_g_option(GMT, "D1");
+		GMT->common.g.active = true;
+		if (!((GMT->common.R.wesn[XHI] - GMT->common.R.wesn[XLO]) == 360 && GMT->common.R.wesn[YLO] == -90 && GMT->common.R.wesn[YHI] == 90))
+			GMT_Report(API, GMT_MSG_WARNING, "Using a non-global region with Spilhaus projection has unknown effects.\n");
+	}
+
 	if ((error = gmt_DCW_list (GMT, &(Ctrl->E.info)))) {	/* This is either success or failure... */
 		if (error != GMT_DCW_LIST)
 			return (1);	/* Not good */