Skip to content

A better fix for: grdimage -Jq -R0/360/-90/90 fails to plot map #9056

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
joa-quim merged 1 commit into from
Jun 2, 2026
+22 −20
Merged

A better fix for: grdimage -Jq -R0/360/-90/90 fails to plot map #9056

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
32 changes: 22 additions & 10 deletions src/gmt_init.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 */
Expand All @@ -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) {
Expand Down
10 changes: 0 additions & 10 deletions src/grdimage.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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) {
Expand Down
Loading