Special handling of global pixel grids plotted with grdimage -JQ (#4451)

* Addres limitation with global pixel grids and lack of projection

See #4440 for background.  This implements the solution of adjusting -R in the rare case we cannot plot a rotated grid.

* Update grdimage.c

* Update PS files

src/gmt_grdio.c

@@ -2182,12 +2182,13 @@ int gmt_grd_setregion (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *h, double *
 
 	if (grid_global) {
 		bool true_global_region = (gmt_M_360_range (h->wesn[XLO], h->wesn[XHI]) && gmt_M_360_range (GMT->common.R.wesn[XLO], GMT->common.R.wesn[XHI]));
+		double x_noise = GMT_CONV4_LIMIT * h->inc[GMT_X];
 		wesn[XLO] = h->wesn[XLO] + floor ((wesn[XLO] - h->wesn[XLO]) * HH->r_inc[GMT_X] + GMT_CONV4_LIMIT) * h->inc[GMT_X];
 		wesn[XHI] = h->wesn[XLO] + ceil  ((wesn[XHI] - h->wesn[XLO]) * HH->r_inc[GMT_X] - GMT_CONV4_LIMIT) * h->inc[GMT_X];
 		/* For the odd chance that xmin or xmax are outside the region: bring them in */
 		if ((wesn[XHI] - wesn[XLO]) >= 360.0) {
-			while ((wesn[XLO]) < GMT->common.R.wesn[XLO]) wesn[XLO] += h->inc[GMT_X];
-			while ((wesn[XHI]) > GMT->common.R.wesn[XHI]) wesn[XHI] -= h->inc[GMT_X];
+			while ((wesn[XLO] + x_noise) < GMT->common.R.wesn[XLO]) wesn[XLO] += h->inc[GMT_X];
+			while ((wesn[XHI]- x_noise) > GMT->common.R.wesn[XHI]) wesn[XHI] -= h->inc[GMT_X];
 		}
 		if (true_global_region && (wesn[XHI] - wesn[XLO]) < 360.0)	/* Need to enforce 360 range since that is what it should be */
 			wesn[XHI] = wesn[XLO] + 360.0;

src/grdimage.c

@@ -1075,6 +1075,34 @@ GMT_LOCAL void grdimage_img_color_with_intensity (struct GMT_CTRL *GMT, struct G
 	}
 }
 
+GMT_LOCAL bool grdimage_adjust_R_consideration (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *h) {
+	/* As per https://github.com/GenericMappingTools/gmt/issues/4440, when the user wants
+	 * to plot a pixel-registered global grid using a lon-lat scaling with periodic boundaries
+	 * and a central meridian that is not a multiple of the grid increment, we must actually
+	 * adjust the plot domain -R to be such a multiple.  That, or require projection. */
+
+	double delta;
+
+	if (!gmt_M_is_geographic (GMT, GMT_IN)) return false;	/* No geographic */
+	if (h->registration == GMT_GRID_NODE_REG) return false;	/* gridline-registration has the repeated column needed */
+	if (gmt_M_is_nonlinear_graticule (GMT)) return true;	/* Always have to project when given most projections except -JQ */
+	if (!gmt_M_360_range (GMT->common.R.wesn[XLO], GMT->common.R.wesn[XHI])) return false;	/* No repeating columns would be visible */
+	delta = remainder (h->wesn[XLO] - GMT->common.R.wesn[XLO], h->inc[GMT_X]);
+	if (gmt_M_is_zero (delta)) return false;	/* No need to project if it is lining up */
+	/* Here we need to adjust plot region */
+	GMT_Report (GMT->parent, GMT_MSG_WARNING, "Your grid is pixel-registered, the projection is simply longlat, and plot region is a full 360 degrees in longitude.\n");
+	GMT_Report (GMT->parent, GMT_MSG_WARNING, "Due to lack of repeated boundary nodes, -Rw/e must be given in multiples of the grid increment (%g)\n", h->inc[GMT_X]);
+	GMT_Report (GMT->parent, GMT_MSG_WARNING, "Current region (-R) setting in longitude is %g to %g\n", GMT->common.R.wesn[XLO], GMT->common.R.wesn[XHI]);
+	GMT->common.R.wesn[XLO] += delta;
+	GMT->common.R.wesn[XHI] += delta;
+	if (GMT->common.R.wesn[XHI] <= 0.0) {
+		GMT->common.R.wesn[XLO] += 360.0;
+		GMT->common.R.wesn[XHI] += 360.0;
+	}
+	GMT_Report (GMT->parent, GMT_MSG_WARNING, "Adjusted region (-R) setting in longitude is %g to %g\n", GMT->common.R.wesn[XLO], GMT->common.R.wesn[XHI]);
+	return false;
+}
+
 #define bailout(code) {gmt_M_free_options (mode); return (code);}
 #define Return(code) {Free_Ctrl (GMT, Ctrl); gmt_M_free (GMT, Conf); gmt_end_module (GMT, GMT_cpy); bailout (code);}
 
@@ -1212,7 +1240,7 @@ EXTERN_MSC int GMT_grdimage (void *V_API, int mode, void *args) {
 
 		if (use_intensity_grid && GMT->common.R.active[RSET]) {
 			/* Make sure the region of the intensity grid and -R are in agreement within a noise threshold */
-            double xnoise = Intens_orig->header->inc[GMT_X]*GMT_CONV4_LIMIT, ynoise = Intens_orig->header->inc[GMT_Y]*GMT_CONV4_LIMIT;
+			double xnoise = Intens_orig->header->inc[GMT_X]*GMT_CONV4_LIMIT, ynoise = Intens_orig->header->inc[GMT_Y]*GMT_CONV4_LIMIT;
 			if (GMT->common.R.wesn[XLO] < (Intens_orig->header->wesn[XLO]-xnoise) || GMT->common.R.wesn[XHI] > (Intens_orig->header->wesn[XHI]+xnoise) ||
 			    GMT->common.R.wesn[YLO] < (Intens_orig->header->wesn[YLO]-ynoise) || GMT->common.R.wesn[YHI] > (Intens_orig->header->wesn[YHI]+ynoise)) {
 				GMT_Report (API, GMT_MSG_ERROR, "Requested region exceeds illumination extent\n");
@@ -1302,6 +1330,8 @@ EXTERN_MSC int GMT_grdimage (void *V_API, int mode, void *args) {
 	if (!GMT->common.R.active[RSET] && got_z_grid)	/* -R was not set so we use the grid domain */
 		gmt_set_R_from_grd (GMT, Grid_orig->header);
 
+	if (Ctrl->E.dpi == 0) grdimage_adjust_R_consideration (GMT, header_work);	/* Special check for global pixel-registered plots */
+
 	/* Initialize the projection for the selected -R -J */
 	if (gmt_map_setup (GMT, GMT->common.R.wesn)) Return (GMT_PROJECTION_ERROR);
 

test/grdimage/globalgrid.sh

@@ -1,6 +1,8 @@
 #!/usr/bin/env bash
 #
 # Check if we can wrap global grids over longitude
+# Note: This recently changed due to https://github.com/GenericMappingTools/gmt/pull/4451
+# so the PS has been updated.
 #
 
 ps=globalgrid.ps
@@ -29,4 +31,3 @@ $plot $R1 -Y6i -X5i -O -K >> $ps
 $plot $R2 -Y-2i -O -K >> $ps
 $plot $R3 -Y-2i -O -K >> $ps
 $plot $R4 -Y-2i -O >> $ps
-