diff --git a/src/gmt_init.c b/src/gmt_init.c index bc0b869a064..f78022483e5 100644 --- a/src/gmt_init.c +++ b/src/gmt_init.c @@ -16570,7 +16570,8 @@ void gmt_detect_oblique_region (struct GMT_CTRL *GMT, char *file) { * We wish to "do no harm" as well, so only some situations will get through this function. */ int k_data; - double d_inc, wesn[4]; + double d_inc, wesn[4], pars[16]; /* pars[] backs up the raw projection parameters since gmt_map_setup below mutates them in place */ + bool units_pr_degree; struct GMTAPI_CTRL *API = GMT->parent; /* Shorthand */ if (gmt_M_is_cartesian (GMT, GMT_IN)) return; /* This check only applies to geographic data */ @@ -16580,19 +16581,30 @@ void gmt_detect_oblique_region (struct GMT_CTRL *GMT, char *file) { if (!(gmt_M_360_range (GMT->common.R.wesn[XLO], GMT->common.R.wesn[XHI]) && gmt_M_180_range (GMT->common.R.wesn[YLO], GMT->common.R.wesn[YHI]))) return; /* Gave -Rd or -Rg so need to probe more*/ if (gmt_M_is_azimuthal (GMT) && doubleAlmostEqual (fabs (GMT->current.proj.lat0), 90.0) && !GMT->common.R.oblique) return; /* Nothing to do */ if (GMT->current.proj.projection == GMT_PROJ4_SPILHAUS) return; /* This is one is square */ - gmt_M_memcpy (wesn, GMT->common.R.wesn, 4, double); /* Save the region we were given */ - - if (gmt_map_setup (GMT, GMT->common.R.wesn)) /* Set up projection */ + gmt_M_memcpy(wesn, GMT->common.R.wesn, 4, double); /* Save the region we were given */ + /* This is a throwaway probe: gmt_map_setup scales the projection parameters (e.g. pars[2] /= M_PR_DEG for -Jq scale), + * so we must restore the raw pars and units_pr_degree flag afterwards or the caller's real gmt_map_setup will scale them + * a second time and collapse the map (https://github.com/GenericMappingTools/gmt issue #8963). */ + gmt_M_memcpy(pars, GMT->current.proj.pars, 16, double); + units_pr_degree = GMT->current.proj.units_pr_degree; + + if (gmt_map_setup (GMT, GMT->common.R.wesn)) { /* Set up projection */ + gmt_M_memcpy(GMT->current.proj.pars, pars, 16, double); GMT->current.proj.units_pr_degree = units_pr_degree; return; /* Something went wrong */ - if (gmt_map_perimeter_search (GMT, GMT->common.R.wesn, false)) /* Refine without 0.1 degree padding */ + } + if (gmt_map_perimeter_search(GMT, GMT->common.R.wesn, false)) { /* Refine without 0.1 degree padding */ + gmt_M_memcpy(GMT->current.proj.pars, pars, 16, double); GMT->current.proj.units_pr_degree = units_pr_degree; return; /* Something went wrong */ + } + gmt_M_memcpy(GMT->current.proj.pars, pars, 16, double); /* Restore raw projection parameters clobbered by the probe map_setup */ + GMT->current.proj.units_pr_degree = units_pr_degree; if (GMT->common.R.wesn[XHI] < GMT->common.R.wesn[XLO] || GMT->common.R.wesn[YHI] < GMT->common.R.wesn[YLO]) - gmt_M_memcpy (GMT->common.R.wesn, wesn, 4, double); /* Reset to give region if junk resulted */ - else if (gmt_M_360_range (wesn[XLO], wesn[XHI]) && gmt_M_360_range (GMT->common.R.wesn[XLO], GMT->common.R.wesn[XHI])) - gmt_M_memcpy (GMT->common.R.wesn, wesn, 2, double); /* Reset to given global w/e in the same format */ - gmt_M_memcpy (API->tile_wesn, GMT->common.R.wesn, 4, double); /* Save the region we found */ + gmt_M_memcpy(GMT->common.R.wesn, wesn, 4, double); /* Reset to give region if junk resulted */ + else if (gmt_M_360_range(wesn[XLO], wesn[XHI]) && gmt_M_360_range (GMT->common.R.wesn[XLO], GMT->common.R.wesn[XHI])) + gmt_M_memcpy(GMT->common.R.wesn, wesn, 2, double); /* Reset to given global w/e in the same format */ + gmt_M_memcpy(API->tile_wesn, GMT->common.R.wesn, 4, double); /* Save the region we found */ d_inc = API->tile_inc; /* Increment in degrees, if set */ - if (d_inc == 0.0 && file && file[0] == '@' && (k_data = gmt_remote_dataset_id (API, file)) != GMT_NOTSET) { /* Got a remote file to work on */ + if (d_inc == 0.0 && file && file[0] == '@' && (k_data = gmt_remote_dataset_id(API, file)) != GMT_NOTSET) { /* Got a remote file to work on */ d_inc = API->remote_info[k_data].d_inc; /* Increment in degrees */ } if (d_inc > 0.0) { diff --git a/src/grdimage.c b/src/grdimage.c index ff0783d5f67..eb39334a9c8 100644 --- a/src/grdimage.c +++ b/src/grdimage.c @@ -1495,16 +1495,6 @@ EXTERN_MSC int GMT_grdimage(void *V_API, int mode, void *args) { Ctrl = New_Ctrl (GMT); /* Allocate and initialize a new control structure */ if ((error = parse (GMT, Ctrl, options)) != 0) Return (error); - if (GMT->common.R.active[RSET] && gmt_M_is_geographic(GMT, GMT_IN) && !gmt_M_is_nonlinear_graticule(GMT) && GMT->common.R.wesn[YHI] == 90.0) { - /* Plotting a global pixel-registered grid with a linear/cylindrical (rectangular graticule) projection and a - * north boundary of exactly +90 makes a later gmt_map_setup collapse the y-scale, producing a degenerate plot - * (https://github.com/GenericMappingTools/gmt issue). Nudge the north boundary just below the pole, which is - * what passing -R.../89.999999999999 does and avoids the problem with no visible difference. - */ - GMT->common.R.wesn[YHI] = 90.0 - GMT_CONV12_LIMIT; - GMT_Report(API, GMT_MSG_DEBUG, "Nudged north boundary from 90 to %.12g to avoid degenerate projection scaling\n", GMT->common.R.wesn[YHI]); - } - /*---------------------------- This is the grdimage main code ----------------------------*/ if ((Conf = gmt_M_memory (GMT, NULL, 1, struct GRDIMAGE_CONF)) == NULL) {