From 58113087ddc917d92710bf7d6af3a04a0a32d637 Mon Sep 17 00:00:00 2001 From: Joaquim Date: Sat, 4 Jan 2025 16:58:46 +0000 Subject: [PATCH] Do not apply the gmtgrdio_doctor_geo_increments() small modification 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. --- src/gmt_grdio.c | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) diff --git a/src/gmt_grdio.c b/src/gmt_grdio.c index 3fc055da148..129db7cff7a 100644 --- a/src/gmt_grdio.c +++ b/src/gmt_grdio.c @@ -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]); } }