Skip to content

Commit

Permalink
Make stricter checks on central longitude (#7295)
Browse files Browse the repository at this point in the history 
There are a few cases where GMT was inconsistent and wrong - we will label these as bugs:

1. When the user specified a specific central longitude, it is to be respected, and for a global map the w/e borders are set at 180 distance from that longitude.
2. If not a global map then if the distance from the specified central meridian to either border exceeds 180 then we yield a map region error.
3. If no central meridian is set then it continues to default to the mid longitude.
  • Loading branch information
@PaulWessel @anbj
PaulWessel authored and anbj committed Mar 7, 2023
1 parent b50845a commit 12a6faa
Show file tree
Hide file tree
Showing 2 changed files with 58 additions and 28 deletions.
2 changes: 1 addition & 1 deletion src/gmt_init.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5249,7 +5249,7 @@ GMT_LOCAL bool gmtinit_parse_J_option (struct GMT_CTRL *GMT, char *args_in) {
gmt_M_memset (l_pos, 3, int); gmt_M_memset (p_pos, 3, int);
gmt_M_memset (t_pos, 3, int); gmt_M_memset (d_pos, 3, int);
if (!GMT->common.J.active) /* Down want to clobber this during -Jz/Z after the horizontal part has been set */
GMT->current.proj.lon0 = GMT->current.proj.lat0 = GMT->session.d_NaN; /* Projection center, to be set via -J */
GMT->current.proj.lon0 = GMT->current.proj.central_meridian = GMT->current.proj.lat0 = GMT->session.d_NaN; /* Projection center, to be set via -J */

project = gmtinit_project_type (args, &i, &width_given);
if (project == GMT_NO_PROJ) return (true); /* No valid projection specified */
Expand Down
84 changes: 57 additions & 27 deletions src/gmt_map.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -321,38 +321,47 @@ GMT_LOCAL void gmtmap_set_default_central_meridian (struct GMT_CTRL *GMT) {
else
GMT->current.proj.pars[0] = 0.5 * (GMT->common.R.wesn[XLO] + GMT->common.R.wesn[XHI]); /* Set to middle lon */
GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Central meridian not given, default to %g\n", GMT->current.proj.pars[0]);
GMT->current.proj.central_meridian = GMT->current.proj.pars[0];
}

/*! . */
GMT_LOCAL void gmtmap_cyl_validate_clon (struct GMT_CTRL *GMT, unsigned int mode) {
GMT_LOCAL int gmtmap_cyl_validate_clon (struct GMT_CTRL *GMT, unsigned int mode) {
/* Make sure that for global (360-range) cylindrical projections, the central meridian is neither west nor east.
* If so then we reset it to the middle value or we change -R:
* mode == 0: <clon> should be reset based on w/e mid-point
* mode == 1: -J<clon> is firm so w/e is centered on <c.lon>
* mode == 2: -J<clon> is firm and e-w < 360 so must give an error instead.
*/
if (gmtmap_central_meridian_not_set (GMT))
int error = GMT_NOERROR;

if (gmtmap_central_meridian_not_set (GMT)) /* If not set then we pick halfway between w and e */
gmtmap_set_default_central_meridian (GMT);
else if (GMT->current.map.is_world && (GMT->current.proj.pars[0] == GMT->common.R.wesn[XLO] || GMT->current.proj.pars[0] == GMT->common.R.wesn[XHI])) {
/* Reset central meridian since cannot be 360 away from one of the boundaries since that gives xmin == xmax below */
if (mode == 1) { /* Change -R to fit central meridian */
double w = GMT->current.proj.pars[0] - 180.0, e = GMT->current.proj.pars[0] + 180.0;
if (GMT->current.map.is_world) { /* For full 360 range the central meridian must be in the middle */
double w = GMT->current.proj.pars[0] - 180.0, e = GMT->current.proj.pars[0] + 180.0;
if (!doubleAlmostEqualZero (GMT->common.R.wesn[XLO], w)) { /* Not yet aligned */
GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Region for global cylindrical projection had to be reset from %g/%g to %g/%g\n",
GMT->common.R.wesn[XLO], GMT->common.R.wesn[XHI], w, e);
GMT->common.R.wesn[XLO] = w; GMT->common.R.wesn[XHI] = e;
}
else { /* Change central meridian to fit -R */
double new_lon = 0.5 * (GMT->common.R.wesn[XLO] + GMT->common.R.wesn[XHI]);
GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Central meridian for global cylindrical projection had to be reset from %g to %g\n", GMT->current.proj.pars[0], new_lon);
GMT->current.proj.pars[0] = new_lon;
}
}
else if (!GMT->current.map.is_world) { /* For reginal areas we cannot have clon > 180 away from either boundary */
if (fabs (GMT->current.proj.pars[0] - GMT->common.R.wesn[XLO]) > 180.0 || fabs (GMT->current.proj.pars[0] - GMT->common.R.wesn[XHI]) > 180.0) {
double new_lon = 0.5 * (GMT->common.R.wesn[XLO] + GMT->common.R.wesn[XHI]);
GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Central meridian for cylindrical projection had to be reset from %g to %g\n", GMT->current.proj.pars[0], new_lon);
GMT->current.proj.pars[0] = new_lon;
else if (!GMT->common.R.oblique) { /* For regional (<360) areas we cannot have clon > 180 away from either boundary */
double dw = fabs (GMT->current.proj.pars[0] - GMT->common.R.wesn[XLO]);
double de = fabs (GMT->current.proj.pars[0] - GMT->common.R.wesn[XHI]);
if (dw > 180.0 || de > 180.0) {
if (mode == 2) { /* Yield an error if fixed central longitude, range < 360, and exceed 180 to the border from central longitude */
static char *border[2] = {"Western", "Eastern"};
unsigned int kase = (dw > 180.0) ? 0 : 1;
GMT_Report (GMT->parent, GMT_MSG_ERROR, "%s boundary is > 180 degrees from specified central meridian and thus your region is invalid\n", border[kase]);
error = GMT_MAP_EXCEEDS_360;
}
else { /* Else we just adapt to the situation */
double new_lon = 0.5 * (GMT->common.R.wesn[XLO] + GMT->common.R.wesn[XHI]);
GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Central meridian for cylindrical projection had to be reset from %g to %g\n", GMT->current.proj.pars[0], new_lon);
GMT->current.proj.pars[0] = new_lon;
}
}
}
return (error);
}

/*! . */
Expand Down Expand Up @@ -2888,7 +2897,8 @@ GMT_LOCAL int gmtmap_init_merc (struct GMT_CTRL *GMT, bool *search) {
return GMT_PROJECTION_ERROR;
}
GMT->current.map.is_world = gmt_M_360_range (GMT->common.R.wesn[XLO], GMT->common.R.wesn[XHI]);
gmtmap_cyl_validate_clon (GMT, 0); /* Make sure the central longitude is valid */
if (gmtmap_cyl_validate_clon (GMT, 0)) /* Make sure the central longitude is valid */
return GMT_PROJECTION_ERROR;
gmtproj_vmerc (GMT, GMT->current.proj.pars[0], GMT->current.proj.pars[1]);
GMT->current.proj.j_x *= D;
GMT->current.proj.j_yc *= D;
Expand Down Expand Up @@ -2940,7 +2950,8 @@ GMT_LOCAL int gmtmap_init_cyleq (struct GMT_CTRL *GMT, bool *search) {
GMT->current.proj.iDx = 1.0 / GMT->current.proj.Dx;
GMT->current.proj.iDy = 1.0 / GMT->current.proj.Dy;
GMT->current.map.is_world = gmt_M_360_range (GMT->common.R.wesn[XLO], GMT->common.R.wesn[XHI]);
gmtmap_cyl_validate_clon (GMT, 1); /* Make sure the central longitude is valid */
if (gmtmap_cyl_validate_clon (GMT, 1)) /* Make sure the central longitude is valid */
return GMT_PROJECTION_ERROR;
gmtproj_vcyleq (GMT, GMT->current.proj.pars[0], GMT->current.proj.pars[1]);
gmtproj_cyleq (GMT, GMT->common.R.wesn[XLO], GMT->common.R.wesn[YLO], &xmin, &ymin);
gmtproj_cyleq (GMT, GMT->common.R.wesn[XHI], GMT->common.R.wesn[YHI], &xmax, &ymax);
Expand Down Expand Up @@ -2976,7 +2987,8 @@ GMT_LOCAL bool gmtmap_init_cyleqdist (struct GMT_CTRL *GMT, bool *search) {
gmtmap_set_spherical (GMT, true); /* Force spherical for now */

GMT->current.map.is_world = gmt_M_360_range (GMT->common.R.wesn[XLO], GMT->common.R.wesn[XHI]);
gmtmap_cyl_validate_clon (GMT, 1); /* Make sure the central longitude is valid */
if (gmtmap_cyl_validate_clon (GMT, 1)) /* Make sure the central longitude is valid */
return GMT_PROJECTION_ERROR;
gmtproj_vcyleqdist (GMT, GMT->current.proj.pars[0], GMT->current.proj.pars[1]);
gmtproj_cyleqdist (GMT, GMT->common.R.wesn[XLO], GMT->common.R.wesn[YLO], &xmin, &ymin);
gmtproj_cyleqdist (GMT, GMT->common.R.wesn[XHI], GMT->common.R.wesn[YHI], &xmax, &ymax);
Expand Down Expand Up @@ -3012,7 +3024,8 @@ GMT_LOCAL int gmtmap_init_miller (struct GMT_CTRL *GMT, bool *search) {
gmtmap_set_spherical (GMT, true); /* Force spherical for now */

GMT->current.map.is_world = gmt_M_360_range (GMT->common.R.wesn[XLO], GMT->common.R.wesn[XHI]);
gmtmap_cyl_validate_clon (GMT, 1); /* Make sure the central longitude is valid */
if (gmtmap_cyl_validate_clon (GMT, 1)) /* Make sure the central longitude is valid */
return GMT_PROJECTION_ERROR;
gmtproj_vmiller (GMT, GMT->current.proj.pars[0]);
gmtproj_miller (GMT, GMT->common.R.wesn[XLO], GMT->common.R.wesn[YLO], &xmin, &ymin);
#ifdef CHRISTMAS
Expand Down Expand Up @@ -3058,7 +3071,8 @@ GMT_LOCAL int gmtmap_init_cylstereo (struct GMT_CTRL *GMT, bool *search) {
gmtmap_set_spherical (GMT, true); /* Force spherical for now */

GMT->current.map.is_world = gmt_M_360_range (GMT->common.R.wesn[XLO], GMT->common.R.wesn[XHI]);
gmtmap_cyl_validate_clon (GMT, 1); /* Make sure the central longitude is valid */
if (gmtmap_cyl_validate_clon (GMT, 1)) /* Make sure the central longitude is valid */
return GMT_PROJECTION_ERROR;
gmtproj_vcylstereo (GMT, GMT->current.proj.pars[0], GMT->current.proj.pars[1]);
gmtproj_cylstereo (GMT, GMT->common.R.wesn[XLO], GMT->common.R.wesn[YLO], &xmin, &ymin);
gmtproj_cylstereo (GMT, GMT->common.R.wesn[XHI], GMT->common.R.wesn[YHI], &xmax, &ymax);
Expand Down Expand Up @@ -4416,6 +4430,8 @@ GMT_LOCAL int gmtmap_init_mollweide (struct GMT_CTRL *GMT, bool *search) {
gmtmap_set_default_central_meridian (GMT);
if (GMT->current.proj.pars[0] < 0.0) GMT->current.proj.pars[0] += 360.0;
GMT->current.map.is_world = gmt_M_360_range (GMT->common.R.wesn[XLO], GMT->common.R.wesn[XHI]);
if (gmtmap_cyl_validate_clon (GMT, 2)) /* Make sure the central longitude is valid */
return GMT_PROJECTION_ERROR;
if (GMT->current.proj.units_pr_degree) GMT->current.proj.pars[1] /= GMT->current.proj.M_PR_DEG;
GMT->current.proj.scale[GMT_X] = GMT->current.proj.scale[GMT_Y] = M_PI * GMT->current.proj.pars[1] / sqrt (8.0);
gmtproj_vmollweide (GMT, GMT->current.proj.pars[0], GMT->current.proj.pars[1]);
Expand Down Expand Up @@ -4475,6 +4491,8 @@ GMT_LOCAL int gmtmap_init_hammer (struct GMT_CTRL *GMT, bool *search) {
gmtmap_set_default_central_meridian (GMT);
if (GMT->current.proj.pars[0] < 0.0) GMT->current.proj.pars[0] += 360.0;
GMT->current.map.is_world = gmt_M_360_range (GMT->common.R.wesn[XLO], GMT->common.R.wesn[XHI]);
if (gmtmap_cyl_validate_clon (GMT, 2)) /* Make sure the central longitude is valid */
return GMT_PROJECTION_ERROR;
if (GMT->current.proj.units_pr_degree) GMT->current.proj.pars[1] /= GMT->current.proj.M_PR_DEG;
GMT->current.proj.scale[GMT_X] = GMT->current.proj.scale[GMT_Y] = 0.5 * M_PI * GMT->current.proj.pars[1] / M_SQRT2;
gmtproj_vhammer (GMT, GMT->current.proj.pars[0], GMT->current.proj.pars[1]);
Expand Down Expand Up @@ -4530,8 +4548,10 @@ GMT_LOCAL int gmtmap_init_grinten (struct GMT_CTRL *GMT, bool *search) {

if (gmtmap_central_meridian_not_set (GMT))
gmtmap_set_default_central_meridian (GMT);
if (GMT->current.proj.pars[0] < 0.0) GMT->current.proj.pars[0] += 360.0;
GMT->current.map.is_world = gmt_M_360_range (GMT->common.R.wesn[XLO], GMT->common.R.wesn[XHI]);
if (gmtmap_cyl_validate_clon (GMT, 2)) /* Make sure the central longitude is valid */
return GMT_PROJECTION_ERROR;
if (GMT->current.proj.pars[0] < 0.0) GMT->current.proj.pars[0] += 360.0;
if (GMT->current.proj.units_pr_degree) GMT->current.proj.pars[1] /= GMT->current.proj.M_PR_DEG;
GMT->current.proj.scale[GMT_X] = GMT->current.proj.scale[GMT_Y] = GMT->current.proj.pars[1];
gmtproj_vgrinten (GMT, GMT->current.proj.pars[0], GMT->current.proj.pars[1]);
Expand Down Expand Up @@ -4594,8 +4614,10 @@ GMT_LOCAL int gmtmap_init_winkel (struct GMT_CTRL *GMT, bool *search) {

if (gmtmap_central_meridian_not_set (GMT))
gmtmap_set_default_central_meridian (GMT);
if (GMT->current.proj.pars[0] < 0.0) GMT->current.proj.pars[0] += 360.0;
GMT->current.map.is_world = gmt_M_360_range (GMT->common.R.wesn[XLO], GMT->common.R.wesn[XHI]);
if (gmtmap_cyl_validate_clon (GMT, 2)) /* Make sure the central longitude is valid */
return GMT_PROJECTION_ERROR;
if (GMT->current.proj.pars[0] < 0.0) GMT->current.proj.pars[0] += 360.0;
if (GMT->current.proj.units_pr_degree) GMT->current.proj.pars[1] /= GMT->current.proj.M_PR_DEG;
gmtproj_vwinkel (GMT, GMT->current.proj.pars[0], GMT->current.proj.pars[1]);
GMT->current.proj.scale[GMT_X] = GMT->current.proj.scale[GMT_Y] = 2.0 * GMT->current.proj.pars[1] / (1.0 + GMT->current.proj.r_cosphi1);
Expand Down Expand Up @@ -4651,8 +4673,10 @@ GMT_LOCAL int gmtmap_init_eckert4 (struct GMT_CTRL *GMT, bool *search) {

if (gmtmap_central_meridian_not_set (GMT))
gmtmap_set_default_central_meridian (GMT);
if (GMT->current.proj.pars[0] < 0.0) GMT->current.proj.pars[0] += 360.0;
GMT->current.map.is_world = gmt_M_360_range (GMT->common.R.wesn[XLO], GMT->common.R.wesn[XHI]);
if (gmtmap_cyl_validate_clon (GMT, 2)) /* Make sure the central longitude is valid */
return GMT_PROJECTION_ERROR;
if (GMT->current.proj.pars[0] < 0.0) GMT->current.proj.pars[0] += 360.0;
if (GMT->current.proj.units_pr_degree) GMT->current.proj.pars[1] /= GMT->current.proj.M_PR_DEG;
gmtproj_veckert4 (GMT, GMT->current.proj.pars[0]);
GMT->current.proj.scale[GMT_X] = GMT->current.proj.scale[GMT_Y] = GMT->current.proj.pars[1];
Expand Down Expand Up @@ -4709,8 +4733,10 @@ GMT_LOCAL int gmtmap_init_eckert6 (struct GMT_CTRL *GMT, bool *search) {

if (gmtmap_central_meridian_not_set (GMT))
gmtmap_set_default_central_meridian (GMT);
if (GMT->current.proj.pars[0] < 0.0) GMT->current.proj.pars[0] += 360.0;
GMT->current.map.is_world = gmt_M_360_range (GMT->common.R.wesn[XLO], GMT->common.R.wesn[XHI]);
if (gmtmap_cyl_validate_clon (GMT, 2)) /* Make sure the central longitude is valid */
return GMT_PROJECTION_ERROR;
if (GMT->current.proj.pars[0] < 0.0) GMT->current.proj.pars[0] += 360.0;
if (GMT->current.proj.units_pr_degree) GMT->current.proj.pars[1] /= GMT->current.proj.M_PR_DEG;
gmtproj_veckert6 (GMT, GMT->current.proj.pars[0]);
GMT->current.proj.scale[GMT_X] = GMT->current.proj.scale[GMT_Y] = 0.5 * GMT->current.proj.pars[1] * sqrt (2.0 + M_PI);
Expand Down Expand Up @@ -4766,8 +4792,10 @@ GMT_LOCAL int gmtmap_init_robinson (struct GMT_CTRL *GMT, bool *search) {

if (gmtmap_central_meridian_not_set (GMT))
gmtmap_set_default_central_meridian (GMT);
if (GMT->current.proj.pars[0] < 0.0) GMT->current.proj.pars[0] += 360.0;
GMT->current.map.is_world = gmt_M_360_range (GMT->common.R.wesn[XLO], GMT->common.R.wesn[XHI]);
if (gmtmap_cyl_validate_clon (GMT, 2)) /* Make sure the central longitude is valid */
return GMT_PROJECTION_ERROR;
if (GMT->current.proj.pars[0] < 0.0) GMT->current.proj.pars[0] += 360.0;
if (GMT->current.proj.units_pr_degree) GMT->current.proj.pars[1] /= GMT->current.proj.M_PR_DEG;
gmtproj_vrobinson (GMT, GMT->current.proj.pars[0]);
GMT->current.proj.scale[GMT_X] = GMT->current.proj.scale[GMT_Y] = GMT->current.proj.pars[1] / 0.8487;
Expand Down Expand Up @@ -4824,8 +4852,10 @@ GMT_LOCAL int gmtmap_init_sinusoidal (struct GMT_CTRL *GMT, bool *search) {

if (gmtmap_central_meridian_not_set (GMT))
gmtmap_set_default_central_meridian (GMT);
if (GMT->current.proj.pars[0] < 0.0) GMT->current.proj.pars[0] += 360.0;
GMT->current.map.is_world = gmt_M_360_range (GMT->common.R.wesn[XLO], GMT->common.R.wesn[XHI]);
if (gmtmap_cyl_validate_clon (GMT, 2)) /* Make sure the central longitude is valid */
return GMT_PROJECTION_ERROR;
if (GMT->current.proj.pars[0] < 0.0) GMT->current.proj.pars[0] += 360.0;
if (GMT->common.R.wesn[YLO] <= -90.0) GMT->current.proj.edge[0] = false;
if (GMT->common.R.wesn[YHI] >= 90.0) GMT->current.proj.edge[2] = false;
gmtproj_vsinusoidal (GMT, GMT->current.proj.pars[0]);
Expand Down

0 comments on commit 12a6faa

Please sign in to comment.