Skip to content

Commit

Permalink
Do not apply the gmtgrdio_doctor_geo_increments() small modification …
Browse files Browse the repository at this point in the history
…if it breaks. (#8660)

See issue #8655. The issue is that gmtgrdio_doctor_geo_increments() may change the increments in such a way that a grid with a previous good range as a multiple of inc becomes bad (because only the incs were changed). This PR checks that the _doctor_ works does not break the condition. If it does, do not apply it.
  • Loading branch information
joa-quim authored Jan 4, 2025
1 parent ae8103f commit 5811308
Showing 1 changed file with 10 additions and 4 deletions.
14 changes: 10 additions & 4 deletions src/gmt_grdio.c
Original file line number Diff line number Diff line change
Expand Up @@ -1396,25 +1396,31 @@ int gmtlib_get_grdtype (struct GMT_CTRL *GMT, unsigned int direction, struct GMT
return (GMT_GRID_CARTESIAN);
}

GMT_LOCAL void gmtgrdio_doctor_geo_increments (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header) {
GMT_LOCAL void gmtgrdio_doctor_geo_increments(struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header) {
/* Check for sloppy arc min/sec increments due to divisions by 60 or 3600 */
double round_inc[2], scale[2], inc, slop;
unsigned int side, n_fix = 0;
static char *type[2] = {"longitude", "latitude"};

GMT_Report (GMT->parent, GMT_MSG_DEBUG, "Call gmtgrdio_doctor_geo_increments on a geographic grid\n");
GMT_Report(GMT->parent, GMT_MSG_DEBUG, "Call gmtgrdio_doctor_geo_increments on a geographic grid\n");
for (side = GMT_X; side <= GMT_Y; side++) { /* Check both increments */
scale[side] = (header->inc[side] < GMT_MIN2DEG) ? 3600.0 : 60.0; /* Check for clean multiples of minutes or seconds */
inc = header->inc[side] * scale[side];
round_inc[side] = rint (inc);
slop = fabs (inc - round_inc[side]);
slop = fabs(inc - round_inc[side]);
if (slop > 0 && slop < GMT_CONV4_LIMIT) n_fix++;
}
if (n_fix == 2) {
unsigned int good = gmt_minmaxinc_verify(GMT, header->wesn[XLO], header->wesn[XHI], header->inc[GMT_X], GMT_CONV4_LIMIT);
for (side = GMT_X; side <= GMT_Y; side++) { /* Check both increments */
inc = header->inc[side];
header->inc[side] = round_inc[side] / scale[side];
GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Round-off patrol changed geographic grid increment for %s from %.18g to %.18g\n",
if (good == 0 && gmt_minmaxinc_verify(GMT, header->wesn[XLO], header->wesn[XHI], header->inc[GMT_X], GMT_CONV4_LIMIT) != 0) {
/* We are in a situation where -R -I were right but the slightly modified increments make it wrong. So do nothing & return */
header->inc[side] = inc;
break;
}
GMT_Report(GMT->parent, GMT_MSG_INFORMATION, "Round-off patrol changed geographic grid increment for %s from %.18g to %.18g\n",
type[side], inc, header->inc[side]);
}
}
Expand Down

0 comments on commit 5811308

Please sign in to comment.