Skip to content

Commit 0fac95e

Browse files
PaulWesselgithub-actions[bot]
authored andcommitted
Make sure intensity subset matches tiled region if given (#4130)
* Make sure intensity subset matches tiled region if given If an intensity grid is given, and a remote earth_relief grid is given, then when the intensity grid is read in and we specify the subregion wesn, make sure this region matches the possibly adjusted tile_wesn due to rounding in gmt_remote.c See #4124 for background. * Update grdimage.c * Update grdimage.c Actually, similar problem cased the issue to begin with.
1 parent c49a896 commit 0fac95e

File tree

1 file changed

+5
-4
lines changed

1 file changed

+5
-4
lines changed

src/grdimage.c

Lines changed: 5 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -901,8 +901,9 @@ EXTERN_MSC int GMT_grdimage (void *V_API, int mode, void *args) {
901901
}
902902

903903
if (n_grids) { /* Get grid dimensions */
904-
n_columns = gmt_M_get_n (GMT, wesn[XLO], wesn[XHI], Grid_orig[0]->header->inc[GMT_X], Grid_orig[0]->header->registration);
905-
n_rows = gmt_M_get_n (GMT, wesn[YLO], wesn[YHI], Grid_orig[0]->header->inc[GMT_Y], Grid_orig[0]->header->registration);
904+
double *region = (gmt_file_is_tiled_list (API, Ctrl->In.file[0], NULL, NULL, NULL)) ? API->tile_wesn : wesn; /* Make sure we get correct dimensions if tiled grids are used */
905+
n_columns = gmt_M_get_n (GMT, region[XLO], region[XHI], Grid_orig[0]->header->inc[GMT_X], Grid_orig[0]->header->registration);
906+
n_rows = gmt_M_get_n (GMT, region[YLO], region[YHI], Grid_orig[0]->header->inc[GMT_Y], Grid_orig[0]->header->registration);
906907
}
907908

908909
if (Ctrl->D.active) { /* Trust the info from gdal to make it more stable against pixel vs grid registration troubles */
@@ -932,11 +933,11 @@ EXTERN_MSC int GMT_grdimage (void *V_API, int mode, void *args) {
932933
/* If given, get intensity grid or compute intensities (for a constant intensity) */
933934

934935
if (use_intensity_grid) { /* Illumination wanted */
935-
936+
double *region = (gmt_file_is_tiled_list (API, Ctrl->In.file[0], NULL, NULL, NULL)) ? API->tile_wesn : wesn; /* Subset to pass to GMT_Read_Data if data set is tiled */
936937
GMT_Report (API, GMT_MSG_INFORMATION, "Allocates memory and read intensity file\n");
937938

938939
/* Remember, the illumination header was already read at the start of grdimage */
939-
if (!Ctrl->I.derive && GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_DATA_ONLY, wesn, Ctrl->I.file, Intens_orig) == NULL) {
940+
if (!Ctrl->I.derive && GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_DATA_ONLY, region, Ctrl->I.file, Intens_orig) == NULL) {
940941
Return (API->error); /* Failed to read the intensity grid data */
941942
}
942943
mixed = grdimage_clean_global_headers (GMT, Intens_orig->header);

0 commit comments

Comments
 (0)