[Backport 6.1] Add special check for non-rotated global grids

src/gmt_map.c

@@ -7956,13 +7956,24 @@ int gmt_grd_project (struct GMT_CTRL *GMT, struct GMT_GRID *I, struct GMT_GRID *
 
 	O->header->z_min = FLT_MAX; O->header->z_max = -FLT_MAX;	/* Min/max for out */
 	if (GMT->common.n.antialias) {	/* Blockaverage repeat pixels, at least the first ~32767 of them... */
-		int n_columns = O->header->n_columns, n_rows = O->header->n_rows;
+		bool skip_repeat = false, duplicate_east = false;
+		int n_columns = O->header->n_columns, n_rows = O->header->n_rows, kase, duplicate_col;
 		nz = gmt_M_memory (GMT, NULL, O->header->size, short int);
 		/* Cannot do OPENMP yet here since it would require a reduction into an output array (nz) */
+		kase = gmt_whole_earth (GMT, I->header->wesn, GMT->common.R.wesn);
+		if (kase == 1 && I->header->registration == GMT_GRID_NODE_REG) {
+			/* Need to avoid giving the repeated east and west meridian values double weight when they plot inside the image.
+			 * This is only likely to happen when external global grids are passed in via GMT_IS_REFERENCE. */
+			skip_repeat = true;	/* Since both -180 and +180 fall inside the image, we only want to use one of then (-180) */
+			duplicate_east = gmt_M_is_periodic (GMT);	/* Since the meridian corresponding to map west only appears once we may need to duplicate on east */
+			if (duplicate_east)	/* Find the column in I->data that corresponds to the longitude of the left boundary of the map */
+				duplicate_col = gmt_M_grd_x_to_col (GMT, GMT->common.R.wesn[XLO], I->header);
+		}
 
 		gmt_M_row_loop (GMT, I, row_in) {	/* Loop over the input grid row coordinates */
 			if (gmt_M_is_rect_graticule (GMT)) y_proj = y_in_proj[row_in];
 			gmt_M_col_loop (GMT, I, row_in, col_in, ij_in) {	/* Loop over the input grid col coordinates */
+				if (skip_repeat && col_in == 0) continue;
 				if (gmt_M_is_rect_graticule (GMT))
 					x_proj = x_in_proj[col_in];
 				else if (inverse)
@@ -7991,6 +8002,38 @@ int gmt_grd_project (struct GMT_CTRL *GMT, struct GMT_GRID *I, struct GMT_GRID *
 				if (O->data[ij_out] < O->header->z_min) O->header->z_min = O->data[ij_out];
 				else if (O->data[ij_out] > O->header->z_max) O->header->z_max = O->data[ij_out];
 			}
+			if (duplicate_east) {
+				ij_in = ij_in - I->header->n_columns + duplicate_col;	/* Rewind to the col to be duplicated */
+				col_in = duplicate_col;
+				if (gmt_M_is_rect_graticule (GMT))
+					x_proj = GMT->current.proj.rect[XHI];
+				else if (inverse)
+					gmt_xy_to_geo (GMT, &x_proj, &y_proj, x_in[col_in], y_in[row_in]);
+				else {
+					if (GMT->current.map.outside (GMT, x_in[col_in], y_in[row_in])) continue;	/* Quite possible we are beyond the horizon */
+					gmt_geo_to_xy (GMT, GMT->common.R.wesn[XHI], y_in[row_in], &x_proj, &y_proj);
+				}
+
+				/* Here, (x_proj, y_proj) is the projected grid point.  Now find nearest node on the output grid */
+
+				row_out = gmt_M_grd_y_to_row (GMT, y_proj, O->header);
+				if (row_out < 0 || row_out >= n_rows) continue;	/* Outside our grid region */
+				col_out = gmt_M_grd_x_to_col (GMT, x_proj, O->header);
+				if (col_out < 0 || col_out >= n_columns) continue;	/* Outside our grid region */
+
+				/* OK, this projected point falls inside the projected grid's rectangular domain */
+
+				ij_out = gmt_M_ijp (O->header, row_out, col_out);	/* The output node */
+				if (nz[ij_out] == 0) O->data[ij_out] = 0.0f;	/* First time, override the initial value */
+				if (nz[ij_out] < SHRT_MAX) {			/* Avoid overflow */
+					O->data[ij_out] += I->data[ij_in];	/* Add up the z-sum inside this rect... */
+					nz[ij_out]++;				/* ..and how many points there were */
+				}
+				if (gmt_M_is_fnan (O->data[ij_out])) continue;
+				if (O->data[ij_out] < O->header->z_min) O->header->z_min = O->data[ij_out];
+				else if (O->data[ij_out] > O->header->z_max) O->header->z_max = O->data[ij_out];
+
+			}
 		}
 	}
 

src/gmt_project.h

@@ -140,6 +140,9 @@ enum GMT_enum_zdown {GMT_ZDOWN_R = 0,	/* Default: Annotating radius */
 	GMT_ZDOWN_ZP	= 2,	/* Annotating planetary radius - r */
 	GMT_ZDOWN_ZR	= 3};	/* Annotating given radius - r */
 
+/* gmt_M_is_periodic means the east and west meridians of a global map are separated */
+#define gmt_M_is_periodic(C) (gmt_M_is_cylindrical (C) || gmt_M_is_misc (C))
+
 /* gmt_M_is_rect_graticule means parallels and meridians are orthogonal, but does not imply linear spacing */
 #define gmt_M_is_rect_graticule(C) (C->current.proj.projection <= GMT_MILLER)
 

src/psxy.c

@@ -1215,7 +1215,7 @@ EXTERN_MSC int GMT_psxy (void *V_API, int mode, void *args) {
 		/* Determine if we need to worry about repeating periodic symbols */
 		if ((Ctrl->N.mode == PSXY_CLIP_REPEAT || Ctrl->N.mode == PSXY_NO_CLIP_REPEAT) && gmt_M_360_range (GMT->common.R.wesn[XLO], GMT->common.R.wesn[XHI]) && gmt_M_is_geographic (GMT, GMT_IN)) {
 			/* Only do this for projection where west and east are split into two separate repeating boundaries */
-			periodic = (gmt_M_is_cylindrical (GMT) || gmt_M_is_misc (GMT));
+			periodic = gmt_M_is_periodic (GMT);
 		}
 		n_times = (periodic) ? 2 : 1;	/* For periodic boundaries we plot each symbol twice to allow for periodic clipping */
 

src/psxyz.c

@@ -858,7 +858,7 @@ EXTERN_MSC int GMT_psxyz (void *V_API, int mode, void *args) {
 		/* Determine if we need to worry about repeating periodic symbols */
 		if (clip_set && (Ctrl->N.mode == PSXYZ_CLIP_REPEAT || Ctrl->N.mode == PSXYZ_NO_CLIP_REPEAT) && gmt_M_360_range (GMT->common.R.wesn[XLO], GMT->common.R.wesn[XHI]) && gmt_M_is_geographic (GMT, GMT_IN)) {
 			/* Only do this for projection where west and east are split into two separate repeating boundaries */
-			periodic = (gmt_M_is_cylindrical (GMT) || gmt_M_is_misc (GMT));
+			periodic = gmt_M_is_periodic (GMT);
 			if (S.symbol == GMT_SYMBOL_GEOVECTOR) periodic = false;
 		}
 		n_times = (periodic) ? 2 : 1;	/* For periodic boundaries we plot each symbol twice to allow for periodic clipping */