Try to get matrix as grid to work with changing central meridian (#3813)

* Try to get matrix as grid to work with changing central meridian

This effort runs into trouble with the function gmt_grd_setregion that determins what region of the grid we should use.

* Adjust projected limits

* Update globalgrid.ps

* Fix columns

* Update gmt_api.c

Author: PaulWessel

src/gmt_api.c

@@ -5204,7 +5204,7 @@ GMT_LOCAL struct GMT_GRID * gmtapi_import_grid (struct GMTAPI_CTRL *API, int obj
 
 	int item, new_item, new_ID;
 	bool done = true, new = false, row_by_row;
- 	uint64_t row, col, i0, i1, j0, j1, ij, ij_orig;
+ 	uint64_t row, col, kol, i0, i1, j0, j1, ij, ij_orig;
 	size_t size;
 	unsigned int both_set = (GMT_CONTAINER_ONLY | GMT_DATA_ONLY);
 	unsigned int method;
@@ -5231,7 +5231,7 @@ GMT_LOCAL struct GMT_GRID * gmtapi_import_grid (struct GMTAPI_CTRL *API, int obj
 		GMT_Report (API, GMT_MSG_ERROR, "Can only use method GMT_IS_FILE when row-by-row reading of grid is selected\n");
 		return_null (API, GMT_NOT_A_VALID_METHOD);
 	}
-	if ((mode & GMT_CONTAINER_ONLY) && S_obj->region) {
+	if ((mode & GMT_CONTAINER_ONLY) && S_obj->region && S_obj->method == GMT_IS_FILE) {
 		GMT_Report (API, GMT_MSG_ERROR, "Cannot request a subset when just inquiring about the grid header\n");
 		return_null (API, GMT_SUBSET_NOT_ALLOWED);		
 	}
@@ -5466,8 +5466,9 @@ GMT_LOCAL struct GMT_GRID * gmtapi_import_grid (struct GMTAPI_CTRL *API, int obj
 
 			for (row = j0; row <= j1; row++) {
 				for (col = i0; col <= i1; col++, ij++) {
-					ij_orig = GMT_2D_to_index (row, col, M_obj->dim);	/* Position of this (row,col) in input matrix organization */
-					ij = gmt_M_ijp (G_obj->header, row, col);	/* Position of this (row,col) in output grid organization */
+					kol = col % M_obj->n_columns;
+					ij_orig = GMT_2D_to_index (row, kol, M_obj->dim);	/* Position of this (row,col) in input matrix organization */
+					ij = gmt_M_ijp (G_obj->header, row, kol);	/* Position of this (row,col) in output grid organization */
 					api_get_val (&(M_obj->data), ij_orig, &d);	/* Get the next item from the matrix */
 					G_obj->data[ij] = (gmt_grdfloat)d;
 					if (gmt_M_is_dnan (d))
@@ -5478,6 +5479,10 @@ GMT_LOCAL struct GMT_GRID * gmtapi_import_grid (struct GMTAPI_CTRL *API, int obj
 					}
 				}
 			}
+			if (gmt_M_is_geographic (GMT, GMT_IN) && gmt_M_360_range (M_obj->range[XLO], M_obj->range[XHI]) && gmt_M_360_range (G_obj->header->wesn[XLO], G_obj->header->wesn[XHI])) {
+				/* Global grids passed via matrix are not rotated to fit the desired global region, so we need to correct the wesn for this grid to match the matrix */
+				gmt_M_memcpy (G_obj->header->wesn, M_obj->range, 4U, double);
+			}
 			gmt_BC_init (GMT, G_obj->header);	/* Initialize grid interpolation and boundary condition parameters */
 			if (gmt_M_err_pass (GMT, gmt_grd_BC_set (GMT, G_obj, GMT_IN), "Grid memory"))
 				return_null (API, GMT_GRID_BC_ERROR);	/* Set boundary conditions */
@@ -5492,7 +5497,7 @@ GMT_LOCAL struct GMT_GRID * gmtapi_import_grid (struct GMTAPI_CTRL *API, int obj
 			//if (S_obj->region) return_null (API, GMT_SUBSET_NOT_ALLOWED);
 			/* This method requires the input data to be a GMT_GRD_FORMAT matrix - otherwise we should be DUPLICATING */
 			MH = gmt_get_M_hidden (M_obj);
-			if (!(M_obj->shape == GMT_IS_ROW_FORMAT && M_obj->type == GMT_GRDFLOAT && MH->alloc_mode == GMT_ALLOC_EXTERNALLY && (mode & GMT_GRID_IS_COMPLEX_MASK) == 0))
+			if (!(M_obj->shape == GMT_IS_ROW_FORMAT && M_obj->type == GMT_GRDFLOAT && (mode & GMT_GRID_IS_COMPLEX_MASK) == 0))
 				 return_null (API, GMT_NOT_A_VALID_IO_ACCESS);
 			if (grid == NULL) {	/* Only allocate when not already allocated.  Note cannot have pad since input matrix wont have one */
 				uint64_t dim[3] = {M_obj->n_rows, M_obj->n_columns, 1};

src/gmt_grdio.c

@@ -2170,13 +2170,16 @@ int gmt_grd_setregion (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *h, double *
 	/* Periodic grid with 360 degree range is easy */
 
 	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]));
 		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) {
+		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];
 		}
+		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;
 		return (1);
 	}
 	if (GMT->current.map.is_world) {

src/grdimage.c

@@ -567,6 +567,13 @@ GMT_LOCAL void grdimage_set_proj_limits (struct GMT_CTRL *GMT, struct GMT_GRID_H
 			r->wesn[YLO] = GMT->current.proj.rect[YLO];	r->wesn[YHI] = GMT->current.proj.rect[YHI];
 		}
 	}
+	else if (gmt_M_x_is_lon (GMT, GMT_IN)) {	/* Extra check for non-projected longitudes that wrap */
+		double x1;
+		gmt_geo_to_xy (GMT, g->wesn[XHI]-g->inc[GMT_X], g->wesn[YHI], &x1, &y);
+		gmt_geo_to_xy (GMT, g->wesn[XHI], g->wesn[YHI], &x, &y);
+		if (x < x1) /* Wrapped around because end of pixel is outside east; use plot width instead */
+			r->wesn[XHI] = r->wesn[XLO] + GMT->current.proj.rect[XHI];
+	}
 }
 
 GMT_LOCAL int grdimage_set_rgb_three_grids (struct GMT_GRID *Grid_proj[], uint64_t node, double *NaN_rgb, double *rgb) {

src/testapi_matrix_360.c

@@ -13,7 +13,7 @@ int main () {
 
 	API = GMT_Create_Session ("test", 2U, mode, NULL);
 
-	/* Read in earth_relief_01d as a matrix from text file and set correct region, inc, registartion */
+	/* Read in earth_relief_01d as a matrix from text file and set correct region, inc, registration */
 	M = GMT_Read_Data (API, GMT_IS_MATRIX, GMT_IS_FILE, GMT_IS_SURFACE, GMT_READ_NORMAL, NULL, "@earth_relief_01d.txt", NULL);
 	M->range[0] = -180;	M->range[1] = 180; M->range[2] = -90;	M->range[3] = 90.0;
 	M->inc[0] = M->inc[1] = 1.0;
@@ -24,8 +24,12 @@ int main () {
 	sprintf (args, "%s -Rg -JH0/6i -Bg30 -K -Cgeo -P", input);
 	GMT_Call_Module (API, "grdimage", GMT_MODULE_CMD, args);
 	GMT_Init_VirtualFile (API, 0, input);
+	/* Call grdimage with central longitude 98.75W, which is arbitrary */
+	sprintf (args, "%s -R -JH98.7/6i -Bg30 -K -Cgeo -O -Y3.25i", input);
+	GMT_Call_Module (API, "grdimage", GMT_MODULE_CMD, args);
+	GMT_Init_VirtualFile (API, 0, input);
 	/* Call grdimage with central longitude 180 which means grid needs to be rotated 180 */
-	sprintf (args, "%s -R -JH180/6i -Bg30 -O -Cgeo -Y3.5i", input);
+	sprintf (args, "%s -R -JH180/6i -Bg30 -O -Cgeo -Y3.25i", input);
 	GMT_Call_Module (API, "grdimage", GMT_MODULE_CMD, args);
 	GMT_Close_VirtualFile (API, input);