Skip to content
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

Special attention when tile grids are passed to grdimage #3799

Merged
merged 1 commit into from
Aug 1, 2020
Merged

Special attention when tile grids are passed to grdimage #3799

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
5 changes: 3 additions & 2 deletions src/gmt_map.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8029,8 +8029,9 @@ int gmt_grd_project (struct GMT_CTRL *GMT, struct GMT_GRID *I, struct GMT_GRID *

if (!GMT->common.n.antialias || nz[ij_out] < 2) /* Just use the interpolated value */
O->data[ij_out] = (gmt_grdfloat)z_int;
else if (gmt_M_is_dnan (z_int)) /* Take the average of what we accumulated */
O->data[ij_out] /= nz[ij_out]; /* Plain average */
else if (gmt_M_is_dnan (z_int)) { /* Take the average of what we accumulated */
if (nz[ij_out]) O->data[ij_out] /= nz[ij_out]; /* Plain average */
}
else { /* Weighted average between blockmean'ed and interpolated values */
inv_nz = 1.0 / nz[ij_out];
O->data[ij_out] = (gmt_grdfloat) ((O->data[ij_out] + z_int * inv_nz) / (nz[ij_out] + inv_nz));
Expand Down
22 changes: 13 additions & 9 deletions src/grdimage.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -598,7 +598,7 @@ EXTERN_MSC int gmtlib_read_grd_info (struct GMT_CTRL *GMT, char *file, struct GM

EXTERN_MSC int GMT_grdimage (void *V_API, int mode, void *args) {
bool done, need_to_project, normal_x, normal_y, resampled = false, gray_only = false;
bool nothing_inside = false, use_intensity_grid = false, got_data_grid = false, set_gray, rgb_from_z, rgb_cube_scan;
bool nothing_inside = false, use_intensity_grid = false, got_data_tiles = false, set_gray, rgb_from_z, rgb_cube_scan;
bool has_content = false, mem_G[3] = {false, false, false}, mem_I = false, mem_D = false;
unsigned int n_columns = 0, n_rows = 0, grid_registration = GMT_GRID_NODE_REG, n_grids, intensity_mode;
unsigned int colormask_offset = 0, try, row, col, mixed = 0, *actual_row = NULL, *actual_col = NULL;
Expand Down Expand Up @@ -681,7 +681,7 @@ EXTERN_MSC int GMT_grdimage (void *V_API, int mode, void *args) {
Return (API->error);
if (GMT_Open_VirtualFile (API, GMT_IS_GRID, GMT_IS_SURFACE, GMT_IN|GMT_IS_REFERENCE, Grid_orig[0], data_grd))
Return (API->error);
got_data_grid = true;
got_data_tiles = true;
}
}
else { /* Illumination grid present and can be read */
Expand Down Expand Up @@ -772,7 +772,7 @@ EXTERN_MSC int GMT_grdimage (void *V_API, int mode, void *args) {
if (!Ctrl->D.active) { /* Read the headers of 1 or 3 grids */
for (k = 0; k < n_grids; k++) {
mem_G[k] = gmt_M_file_is_memory (Ctrl->In.file[k]);
if (got_data_grid && k == 0) continue; /* Only true if we already read a SRTM tile bunch earlier under I.derive = true */
if (got_data_tiles && k == 0) continue; /* Only true if we already read a SRTM tile bunch earlier under I.derive = true */
GMT_Report (API, GMT_MSG_INFORMATION, "Read header from file %s\n", Ctrl->In.file[k]);
if ((Grid_orig[k] = GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_ONLY, NULL, Ctrl->In.file[k], NULL)) == NULL) { /* Get header only */
Return (API->error);
Expand All @@ -795,7 +795,7 @@ EXTERN_MSC int GMT_grdimage (void *V_API, int mode, void *args) {
Ctrl->C.active = true; /* Use default CPT (GMT_DEFAULT_CPT_NAME) and autostretch or under modern reuse current CPT */
}

if (n_grids) header_work = Grid_orig[0]->header; /* OK, we are in GRID mode and this was not set previuosly. Do it now */
if (n_grids) header_work = Grid_orig[0]->header; /* OK, we are in GRID mode and this was not set previously. Do it now */

if (n_grids && Ctrl->In.do_rgb) { /* Must ensure all three grids are coregistered */
if (!gmt_M_grd_same_region (GMT, Grid_orig[0], Grid_orig[1])) error++;
Expand Down Expand Up @@ -852,6 +852,7 @@ EXTERN_MSC int GMT_grdimage (void *V_API, int mode, void *args) {
if (Ctrl->I.derive) { /* Auto-create intensity grid from data grid using the now determined data region */
bool got_int4_grid = false;
char int_grd[GMT_VF_LEN] = {""}, int4_grd[GMT_VF_LEN] = {""};
double *region = (got_data_tiles) ? header_work->wesn : wesn; /* Region to pass to grdgradient */
struct GMT_GRID *I_data = NULL;

GMT_Report (API, GMT_MSG_INFORMATION, "Derive intensity grid from data grid\n");
Expand All @@ -865,10 +866,13 @@ EXTERN_MSC int GMT_grdimage (void *V_API, int mode, void *args) {
Return (API->error);
got_int4_grid = true;
}
/* Prepare the grdgradient arguments using selected -A -N and the data region in effect */
/* Prepare the grdgradient arguments using selected -A -N and the data region in effect. If we read in a tiled grid
* then it was made with 0/360 so we must use that region in grdgradient. For non-tiles, we read the actual grid
* AFTER calling gmt_mapsetup which changes the region. The tiled region remains unchanged because it is read in
* all at once as it is being assembled. This fixes https://github.com/GenericMappingTools/gmt/issues/3694 */
sprintf (cmd, "-G%s -A%s -N%s+a%s -R%.16g/%.16g/%.16g/%.16g --GMT_HISTORY=readonly ",
int_grd, Ctrl->I.azimuth, Ctrl->I.method, Ctrl->I.ambient, wesn[XLO], wesn[XHI], wesn[YLO], wesn[YHI]);
if (got_data_grid) /* Use the virtual file we made earlier */
int_grd, Ctrl->I.azimuth, Ctrl->I.method, Ctrl->I.ambient, region[XLO], region[XHI], region[YLO], region[YHI]);
if (got_data_tiles) /* Use the virtual file we made earlier */
strcat (cmd, data_grd);
else if (got_int4_grid) /* Use the virtual file we made earlier */
strcat (cmd, int4_grd);
Expand All @@ -883,7 +887,7 @@ EXTERN_MSC int GMT_grdimage (void *V_API, int mode, void *args) {
/* Obtain the data from the virtual file */
if ((Intens_orig = GMT_Read_VirtualFile (API, int_grd)) == NULL)
Return (API->error);
if (got_data_grid)
if (got_data_tiles)
GMT_Close_VirtualFile (API, data_grd);
else if (got_int4_grid) /* Use the virtual file we made earlier */
GMT_Close_VirtualFile (API, int4_grd);
Expand All @@ -910,7 +914,7 @@ EXTERN_MSC int GMT_grdimage (void *V_API, int mode, void *args) {

for (k = 0; k < n_grids; k++) {
GMT_Report (API, GMT_MSG_INFORMATION, "Allocate and read data from file %s\n", Ctrl->In.file[k]);
if (got_data_grid && k == 0) continue; /* Only true if we already read a SRTM tile bunch earlier under I.derive = true */
if (got_data_tiles && k == 0) continue; /* Only true if we already read a SRTM tile bunch earlier under I.derive = true */
if (GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_DATA_ONLY, wesn, Ctrl->In.file[k], Grid_orig[k]) == NULL) { /* Get grid data */
Return (API->error);
}
Expand Down