diff --git a/src/gmt_init.c b/src/gmt_init.c index 7c4db821161..d44517e8967 100644 --- a/src/gmt_init.c +++ b/src/gmt_init.c @@ -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 */ diff --git a/src/gmt_map.c b/src/gmt_map.c index 5e8e3eb4bf4..b9e2b0d99ec 100644 --- a/src/gmt_map.c +++ b/src/gmt_map.c @@ -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: should be reset based on w/e mid-point * mode == 1: -J is firm so w/e is centered on + * mode == 2: -J 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); } /*! . */ @@ -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; @@ -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); @@ -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); @@ -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 @@ -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); @@ -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]); @@ -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]); @@ -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]); @@ -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); @@ -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]; @@ -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); @@ -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; @@ -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]);