Skip to content

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

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 6 commits into from
Aug 3, 2020
Merged

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

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
80974d0
Try to get matrix as grid to work with changing central meridian
PaulWessel Aug 2, 2020
1294a40
Adjust projected limits
PaulWessel Aug 2, 2020
509eeac
Update globalgrid.ps
PaulWessel Aug 2, 2020
65fea0a
Fix columns
PaulWessel Aug 2, 2020
858102e
Merge branch 'master' into matrix-grid-fix
PaulWessel Aug 2, 2020
b9e5ebd
Update gmt_api.c
PaulWessel Aug 2, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
15 changes: 10 additions & 5 deletions src/gmt_api.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand All @@ -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);
}
Expand Down Expand Up @@ -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))
Expand All @@ -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 */
Expand All @@ -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};
Expand Down
5 changes: 4 additions & 1 deletion src/gmt_grdio.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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) {
Expand Down
7 changes: 7 additions & 0 deletions src/grdimage.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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) {
Expand Down
8 changes: 6 additions & 2 deletions src/testapi_matrix_360.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand All @@ -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);

Expand Down
Binary file added test/api/apimat_360.ps
View file
Open in desktop
Binary file not shown.
Binary file modified test/grdimage/globalgrid.ps
View file
Open in desktop
Binary file not shown.