diff --git a/src/gmt_init.c b/src/gmt_init.c index 7d7fb9fb4bf..50368f0265b 100644 --- a/src/gmt_init.c +++ b/src/gmt_init.c @@ -14830,28 +14830,49 @@ GMT_LOCAL bool gmtinit_is_PS_module (struct GMTAPI_CTRL *API, const char *name, return true; /* Remaining PostScript producing modules always write PostScript */ } -void gmt_round_wesn (double wesn[], bool geo) { /* Use data range to round to nearest reasonable multiples */ +GMT_LOCAL void gmtinit_expand_to_decades (double *lo, double *hi) { + /* Expand a positive range outwards to the enclosing powers of ten */ + double p_lo = floor (log10 (*lo)), p_hi = ceil (log10 (*hi)); + if (pow (10.0, p_lo + 1.0) <= *lo) p_lo += 1.0; /* Avoid rounding error on an exact decade */ + if (pow (10.0, p_hi - 1.0) >= *hi) p_hi -= 1.0; + *lo = pow (10.0, p_lo); + *hi = pow (10.0, p_hi); +} + +GMT_LOCAL void gmtinit_round_wesn_log (double wesn[], bool geo, bool log_axis[2]) { + /* Use data range to round to nearest reasonable multiples. + * log_axis[GMT_X]/[GMT_Y]: if true, round that side to enclosing powers of ten + * instead, so the region never has a bound <= 0 on a log10 axis (issue #6614). */ bool set[2] = {false, false}; unsigned int side, item; double mag, inc, range[2] = {0.0, 0.0}; + if (log_axis[GMT_X] && wesn[XLO] > 0.0 && wesn[XHI] > 0.0) { + gmtinit_expand_to_decades (&wesn[XLO], &wesn[XHI]); + set[GMT_X] = true; + } + if (log_axis[GMT_Y] && wesn[YLO] > 0.0 && wesn[YHI] > 0.0) { + gmtinit_expand_to_decades (&wesn[YLO], &wesn[YHI]); + set[GMT_Y] = true; + } range[GMT_X] = wesn[XHI] - wesn[XLO]; range[GMT_Y] = wesn[YHI] - wesn[YLO]; if (geo) { /* Special checks due to periodicity */ - if (range[GMT_X] > 306.0) { /* If within 15% of a full 360 we promote to 360 */ + if (!set[GMT_X] && range[GMT_X] > 306.0) { /* If within 15% of a full 360 we promote to 360 */ wesn[XLO] = 0.0; wesn[XHI] = 360.0; set[GMT_X] = true; } - if (range[GMT_Y] > 153.0) { /* If within 15% of a full 180 we promote to 180 */ + if (!set[GMT_Y] && range[GMT_Y] > 153.0) { /* If within 15% of a full 180 we promote to 180 */ wesn[YLO] = -90.0; wesn[YHI] = 90.0; set[GMT_Y] = true; } } else { /* Add a tinny pad so that the rounding algorithm always round to next fifth of decade */ double dx, dy; - dx = range[GMT_X] * 0.001; dy = range[GMT_Y] * 0.001; + dx = set[GMT_X] ? 0.0 : range[GMT_X] * 0.001; dy = set[GMT_Y] ? 0.0 : range[GMT_Y] * 0.001; wesn[0] -= dx; wesn[1] += dx; wesn[2] -= dy; wesn[3] += dy; } - for (side = GMT_X, item = XLO; side <= GMT_Y; side++) { + for (side = GMT_X; side <= GMT_Y; side++) { + item = 2 * side; /* XLO|YLO; a skipped side below must not throw off the other side's index */ if (set[side]) continue; /* Done above */ mag = rint (log10 (range[side])) - 1.0; inc = pow (10.0, mag); @@ -14883,6 +14904,37 @@ void gmt_round_wesn (double wesn[], bool geo) { /* Use data range to round to ne } } +/*! . */ +void gmt_round_wesn (double wesn[], bool geo) { /* Use data range to round to nearest reasonable multiples */ + bool no_log[2] = {false, false}; + gmtinit_round_wesn_log (wesn, geo, no_log); +} + +GMT_LOCAL void gmtinit_get_log_axes (struct GMT_OPTION *options, bool log_axis[2]) { + /* Check the raw -J string for a log10 x- and/or y-axis (-Jx|X only). + * -J is not parsed for real until later, so we must not call gmtinit_parse_J_option + * here: that would reset column types and complain about unresolved ?-widths. */ + struct GMT_OPTION *opt = NULL; + char arg[GMT_LEN128] = {""}, *slash = NULL, *plus = NULL, *ypart = NULL; + + log_axis[GMT_X] = log_axis[GMT_Y] = false; + for (opt = options; opt; opt = opt->next) /* Find the horizontal -J option, if any */ + if (opt->option == 'J' && !(opt->arg[0] == 'z' || opt->arg[0] == 'Z')) break; + if (opt == NULL || opt->arg[0] == '\0') return; /* No -J argument to examine */ + if (!(opt->arg[0] == 'x' || opt->arg[0] == 'X')) return; /* Only -Jx|X scale the axes individually */ + strncpy (arg, &opt->arg[1], GMT_LEN128-1); /* Skip the projection code */ + if ((plus = strchr (arg, '+'))) plus[0] = '\0'; /* Chop off any +d modifier */ + if ((slash = strchr (arg, '/'))) { /* Separate y-axis specification given */ + slash[0] = '\0'; + ypart = &slash[1]; + } + else /* Just one specification, which then applies to both axes */ + ypart = arg; + /* Only l|L can appear here for log10; the length units are c|i|p and d means degrees */ + log_axis[GMT_X] = (strpbrk (arg, "lL") != NULL); + log_axis[GMT_Y] = (strpbrk (ypart, "lL") != NULL); +} + GMT_LOCAL int gmtinit_get_region_from_data(struct GMTAPI_CTRL *API, int family, bool exact, struct GMT_OPTION **options, double wesn[], int *aspect) { /* Determines the data region by examining the input data. This could be a grid or datasets. The * latter may be one or many files, or none, meaning we must capture stdin. If we do, then we must @@ -15114,7 +15166,11 @@ GMT_LOCAL int gmtinit_get_region_from_data(struct GMTAPI_CTRL *API, int family, else wesn[YLO] *= 0.9, wesn[YHI] *= 1.1; /* +/- 10% of values */ } - if (!exact) gmt_round_wesn (wesn, geo); /* Use data range to round to nearest reasonable multiples */ + if (!exact) { /* Use data range to round to nearest reasonable multiples */ + bool log_axis[2]; + gmtinit_get_log_axes (*options, log_axis); /* Check -J for log10 axes so we do not round a bound to <= 0 there */ + gmtinit_round_wesn_log (wesn, geo, log_axis); + } break; default: GMT_Report (API, GMT_MSG_DEBUG, "gmtinit_get_region_from_data: Family %d not supported", family); diff --git a/test/common/log_axis_auto_region.sh b/test/common/log_axis_auto_region.sh new file mode 100755 index 00000000000..c791290f8e9 --- /dev/null +++ b/test/common/log_axis_auto_region.sh @@ -0,0 +1,17 @@ +#!/usr/bin/env bash +# Test that -Ra (auto-determine region from data) never rounds a bound on a +# log10 axis to a non-positive value, which is invalid on that axis (#6614). + +gmt begin +gmt math -T0/10240/1 T 10240 DIV 360 MUL 400 MUL COSD = t.txt +gmt spectrum1d t.txt -S256 -W --GMT_FFT=brenner -N -i1 > pow5.txt +gmt plot pow5.txt -Ra -JX-15cl/4c -Bxa2f3 -Bya -Vd > plot.log 2>&1 +gmt end + +region=$(sed -n 's/.*Modern: Adding -R\([^ ]*\) to options\..*/\1/p' plot.log) +xlo=$(cut -d/ -f1 <<< "$region") +if [ -z "$region" ]; then + echo "No -R was determined for -Ra" > fail +elif ! awk -v xlo="$xlo" 'BEGIN {exit !(xlo+0 > 0)}'; then + echo "-Ra gave a non-positive x-min ($xlo) for a log10 x-axis: -R$region" > fail +fi