diff --git a/src/gmt_api.c b/src/gmt_api.c index 29c85c588c4..d1dce11242a 100644 --- a/src/gmt_api.c +++ b/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}; diff --git a/src/gmt_grdio.c b/src/gmt_grdio.c index d9f8822a194..8de8b7ae32d 100644 --- a/src/gmt_grdio.c +++ b/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) { diff --git a/src/grdimage.c b/src/grdimage.c index 1bbbadff636..eeb04be91b2 100644 --- a/src/grdimage.c +++ b/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) { diff --git a/src/testapi_matrix_360.c b/src/testapi_matrix_360.c index 7419b61c072..130cdb56a18 100644 --- a/src/testapi_matrix_360.c +++ b/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); diff --git a/test/api/apimat_360.ps b/test/api/apimat_360.ps new file mode 100644 index 00000000000..a490aef33fc Binary files /dev/null and b/test/api/apimat_360.ps differ diff --git a/test/grdimage/globalgrid.ps b/test/grdimage/globalgrid.ps index 7698103bc59..18c06651fcf 100644 Binary files a/test/grdimage/globalgrid.ps and b/test/grdimage/globalgrid.ps differ