Skip to content
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

Make stricter checks on central longitude #7295

Merged
merged 1 commit into from
Mar 5, 2023
Merged

Make stricter checks on central longitude #7295

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
2 changes: 1 addition & 1 deletion src/gmt_init.c
View file
Open in desktop
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
Open in desktop
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