Skip to content

Fix issues in Mercator projection #4120

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
merged 6 commits into from
Aug 30, 2020
Merged

Fix issues in Mercator projection #4120

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
c06f0b8
Requires special -C +m modifier to report Mercator units relative to …
PaulWessel Aug 29, 2020
6a42fe8
Merge branch 'master' into merc-fix
PaulWessel Aug 29, 2020
45af89c
Document -C+m
PaulWessel Aug 29, 2020
9b73638
Update gmt_proj.c
PaulWessel Aug 30, 2020
f9ac5e7
Update gmt_proj.c
PaulWessel Aug 30, 2020
36f044d
Update merc_origin.ps
PaulWessel Aug 30, 2020
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
8 changes: 5 additions & 3 deletions doc/rst/source/mapproject.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,7 @@ Synopsis
**gmt mapproject** [ *tables* ] |-J|\ *parameters*
|SYN_OPT-R|
[ |-A|\ **b**\|\ **B**\|\ **f**\|\ **F**\|\ **o**\|\ **O**\ [*lon0*/*lat0*][**+v**] ]
[ |-C|\ [*dx*/*dy*] ]
[ |-C|\ [*dx*/*dy*][**+m**] ]
[ |-D|\ **c**\|\ **i**\|\ **p** ]
[ |-E|\ [*datum*] ] [ |-F|\ [*unit*] ]
[ |-G|\ [*lon0*/*lat0*][**+a**][**+i**][**+u**\ *unit*][**+v**] ]
Expand Down Expand Up @@ -105,15 +105,17 @@ Optional Arguments

.. _-C:

**-C**\ [*dx*/*dy*]
**-C**\ [*dx*/*dy*][**+m**]
Set center of projected coordinates to be at map projection center
[Default is lower left corner]. Optionally, add offsets in the
projected units to be added (or subtracted when **-I** is set) to
(from) the projected coordinates, such as false eastings and
northings for particular projection zones [0/0]. The unit used for
the offsets is the plot distance unit in effect (see
:term:`PROJ_LENGTH_UNIT`) unless **-F** is used, in which case the
offsets are in meters.
offsets are in meters. Alternatively, for the Mercator projection
only, append **+m** to set the origin of the projected *y* coordinates
to coincide with the standard parallel [Equator].

.. _-D:

Expand Down
8 changes: 5 additions & 3 deletions src/gmt_init.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -865,7 +865,7 @@ GMT_LOCAL int gmtinit_rectR_to_geoR (struct GMT_CTRL *GMT, char unit, double rec
uint64_t dim[GMT_DIM_SIZE] = {1, 1, 2, 2}; /* Just a single data table with one segment with two 2-column records */
bool was_R, was_J;
double wesn[4];
char buffer[GMT_LEN256] = {""}, Jstring[GMT_LEN128] = {""}, in_string[GMT_VF_LEN] = {""}, out_string[GMT_VF_LEN] = {""}, *v = NULL;
char buffer[GMT_LEN256] = {""}, Jstring[GMT_LEN128] = {""}, in_string[GMT_VF_LEN] = {""}, out_string[GMT_VF_LEN] = {""}, origin_flag[4] = {""}, *v = NULL;
struct GMT_DATASET *In = NULL, *Out = NULL;

GMT_Report (GMT->parent, GMT_MSG_DEBUG, "Call gmtinit_rectR_to_geoR to convert projected -R to geo -R\n");
Expand Down Expand Up @@ -911,6 +911,8 @@ GMT_LOCAL int gmtinit_rectR_to_geoR (struct GMT_CTRL *GMT, char unit, double rec
GMT_Report (GMT->parent, GMT_MSG_ERROR, "UTM projection insufficiently specified to auto-determine geographic region\n");
return (GMT_MAP_NO_PROJECTION);
}
if (GMT->current.proj.projection_GMT == GMT_MERCATOR) /* Special use of Mercator units relative to stated origin */
strcpy (origin_flag, "+m");
break;
case 2: /* Conical: Use default patch */
break;
Expand Down Expand Up @@ -948,8 +950,8 @@ GMT_LOCAL int gmtinit_rectR_to_geoR (struct GMT_CTRL *GMT, char unit, double rec
v[0] = '\0';
}
}
snprintf (buffer, GMT_LEN256, "-R%g/%g/%g/%g -J%s -I -F%c -C -bi2d -bo2d -<%s ->%s --GMT_HISTORY=readonly",
wesn[XLO], wesn[XHI], wesn[YLO], wesn[YHI], Jstring, unit, in_string, out_string);
snprintf (buffer, GMT_LEN256, "-R%g/%g/%g/%g -J%s -I -F%c -C%s -bi2d -bo2d -<%s ->%s --GMT_HISTORY=readonly",
wesn[XLO], wesn[XHI], wesn[YLO], wesn[YHI], Jstring, unit, origin_flag, in_string, out_string);
if (get_R) GMT_Report (GMT->parent, GMT_MSG_DEBUG, "Obtaining geographic corner coordinates via mapproject %s\n", buffer);
if (GMT_Call_Module (GMT->parent, "mapproject", GMT_MODULE_CMD, buffer) != GMT_OK) /* Get the corners in degrees via mapproject */
return (GMT->parent->error);
Expand Down
14 changes: 8 additions & 6 deletions src/gmt_proj.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -652,15 +652,17 @@ GMT_LOCAL void gmtproj_ipolar (struct GMT_CTRL *GMT, double *x, double *y, doubl

/* -JM MERCATOR PROJECTION */

GMT_LOCAL void gmtproj_vmerc (struct GMT_CTRL *GMT, double lon0, double slat) {
GMT_LOCAL void gmtproj_vmerc (struct GMT_CTRL *GMT, double lon0, double lat0) {
/* Set up a Mercator transformation with origin at (lon0, lat0) */

if (GMT->current.proj.GMT_convert_latitudes) slat = gmt_M_latg_to_latc (GMT, slat);
double aux_lat0 = (GMT->current.proj.GMT_convert_latitudes) ? gmt_M_latg_to_latc (GMT, lat0) : lat0;

GMT->current.proj.central_meridian = lon0;
GMT->current.proj.j_x = cosd (slat) / d_sqrt (1.0 - GMT->current.proj.ECC2 * sind (slat) * sind (slat)) * GMT->current.proj.EQ_RAD;
/* Need geodetic latitude in this expression: */
GMT->current.proj.j_x = cosd (lat0) / d_sqrt (1.0 - GMT->current.proj.ECC2 * sind (lat0) * sind (lat0)) * GMT->current.proj.EQ_RAD;
GMT->current.proj.j_ix = 1.0 / GMT->current.proj.j_x;
GMT->current.proj.j_yc = (fabs (slat) > 0.0) ? GMT->current.proj.j_x * d_log (GMT, tand (45.0 + 0.5 * slat)) : 0.0;
/* Need conformal latitude in this expression (same as in gmtproj_merc_sph) */
GMT->current.proj.j_yc = (fabs (lat0) > 0.0) ? GMT->current.proj.j_x * d_log (GMT, tand (45.0 + 0.5 * aux_lat0)) : 0.0;
}

/* Mercator projection for the sphere */
Expand All @@ -672,14 +674,14 @@ GMT_LOCAL void gmtproj_merc_sph (struct GMT_CTRL *GMT, double lon, double lat, d
if (GMT->current.proj.GMT_convert_latitudes) lat = gmt_M_latg_to_latc (GMT, lat);

*x = GMT->current.proj.j_x * D2R * lon;
*y = (fabs (lat) < 90.0) ? GMT->current.proj.j_x * d_log (GMT, tand (45.0 + 0.5 * lat)) - GMT->current.proj.j_yc : copysign (DBL_MAX, lat);
*y = (fabs (lat) < 90.0) ? GMT->current.proj.j_x * d_log (GMT, tand (45.0 + 0.5 * lat)) : copysign (DBL_MAX, lat);
}

GMT_LOCAL void gmtproj_imerc_sph (struct GMT_CTRL *GMT, double *lon, double *lat, double x, double y) {
/* Convert Mercator x/y to lon/lat (GMT->current.proj.EQ_RAD in GMT->current.proj.j_ix) */

*lon = x * GMT->current.proj.j_ix * R2D + GMT->current.proj.central_meridian;
*lat = atand (sinh ((y + GMT->current.proj.j_yc) * GMT->current.proj.j_ix));
*lat = atand (sinh (y * GMT->current.proj.j_ix));
if (GMT->current.proj.GMT_convert_latitudes) *lat = gmt_M_latc_to_latg (GMT, *lat);
}

Expand Down
26 changes: 22 additions & 4 deletions src/mapproject.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -97,9 +97,10 @@ struct MAPPROJECT_CTRL { /* All control options for this program (except common
unsigned int mode;
double lon, lat; /* Fixed point of reference */
} A;
struct MAPPROJECT_C { /* -C[<false_easting>/<false_northing>] */
struct MAPPROJECT_C { /* -C[<false_easting>/<false_northing>][+m] */
bool active;
bool shift;
bool m_origin; /* True if we want projected Mercator y-values relative to standard latitude [Equator] */
double easting, northing; /* Shifts */
} C;
struct MAPPROJECT_D { /* -D<c|i|p> */
Expand Down Expand Up @@ -189,7 +190,7 @@ static int usage (struct GMTAPI_CTRL *API, int level) {
const char *name = gmt_show_name_and_purpose (API, THIS_MODULE_LIB, THIS_MODULE_CLASSIC_NAME, THIS_MODULE_PURPOSE);
if (level == GMT_MODULE_PURPOSE) return (GMT_NOERROR);
GMT_Message (API, GMT_TIME_NONE, "usage: %s <table> %s %s\n", name, GMT_J_OPT, GMT_Rgeo_OPT);
GMT_Message (API, GMT_TIME_NONE, "\t[-Ab|B|f|F|o|O[<lon0>/<lat0>][+v]] [-C[<dx></dy>]] [-D%s] [-E[<datum>]] [-F[<unit>]]\n\t[-G[<lon0>/<lat0>][+a][+i][+u<unit>][+v]]", GMT_DIM_UNITS_DISPLAY);
GMT_Message (API, GMT_TIME_NONE, "\t[-Ab|B|f|F|o|O[<lon0>/<lat0>][+v]] [-C[<dx></dy>][+m]] [-D%s] [-E[<datum>]] [-F[<unit>]]\n\t[-G[<lon0>/<lat0>][+a][+i][+u<unit>][+v]]", GMT_DIM_UNITS_DISPLAY);
GMT_Message (API, GMT_TIME_NONE, " [-I] [-L<table>[+u<unit>][+p]] [-N[a|c|g|m]]\n\t[-Q[e|d]] [-S] [-T[h]<from>[/<to>] [%s] [-W[g|h|j|n|w|x]] [-Z[<speed>][+a][+i][+f][+t<epoch>]]\n", GMT_V_OPT);
GMT_Message (API, GMT_TIME_NONE, "\t[%s] [%s] [%s] [%s]\n\t[%s] [%s]\n\t[%s] [%s] [%s]\n\t[%s]\n\t[%s] [%s] [%s] [%s]\n\n",
GMT_b_OPT, GMT_d_OPT, GMT_e_OPT, GMT_f_OPT, GMT_g_OPT, GMT_h_OPT, GMT_i_OPT, GMT_j_OPT, GMT_o_OPT, GMT_p_OPT, GMT_q_OPT, GMT_s_OPT, GMT_colon_OPT, GMT_PAR_OPT);
Expand All @@ -209,6 +210,7 @@ static int usage (struct GMTAPI_CTRL *API, int level) {
GMT_Message (API, GMT_TIME_NONE, "\t-C Return x/y relative to projection center [Default is relative to lower left corner].\n");
GMT_Message (API, GMT_TIME_NONE, "\t Optionally append <dx></dy> to add (or subtract if -I) (i.e., false easting & northing) [0/0].\n");
GMT_Message (API, GMT_TIME_NONE, "\t Units are plot units unless -F is set in which case the unit is meters.\n");
GMT_Message (API, GMT_TIME_NONE, "\t Mercator only: Append +m to set origin for projected y-values to the standard parallel [Equator].\n");
GMT_Message (API, GMT_TIME_NONE, "\t-D Temporarily reset PROJ_LENGTH_UNIT to be c (cm), i (inch), or p (point).\n");
GMT_Message (API, GMT_TIME_NONE, "\t Cannot be used if -F is set.\n");
GMT_Message (API, GMT_TIME_NONE, "\t-E Convert (lon, lat, h) to Earth Centered Earth Fixed (ECEF) coordinates [-I for inverse].\n");
Expand Down Expand Up @@ -482,12 +484,17 @@ static int parse (struct GMT_CTRL *GMT, struct MAPPROJECT_CTRL *Ctrl, struct GMT
break;
case 'C':
Ctrl->C.active = true;
if (opt->arg[0] && (p = strstr (opt->arg, "+m"))) { /* Gave +m for reset relative offsets to Mercator origin */
Ctrl->C.m_origin = true;
p[0] = '\0'; /* Temporarily chop off modifier */
}
if (opt->arg[0]) { /* Also gave shifts */
n_errors += gmt_M_check_condition (GMT, sscanf (opt->arg, "%lf/%lf", &Ctrl->C.easting, &Ctrl->C.northing) != 2,
"Option -C: Expected -C[<false_easting>/<false_northing>]\n");
Ctrl->C.shift = true;
}
will_need_RJ = true; /* Since -C is used with projections only */
if (p) p[0] = '+'; /* Restore modifier */
break;
case 'D':
Ctrl->D.active = true;
Expand Down Expand Up @@ -752,6 +759,8 @@ static int parse (struct GMT_CTRL *GMT, struct MAPPROJECT_CTRL *Ctrl, struct GMT
n_errors += gmt_M_check_condition (GMT, ((Ctrl->T.active && GMT->current.proj.datum.h_given) || Ctrl->E.active) &&
GMT->common.b.active[GMT_IN] && gmt_get_cols (GMT, GMT_IN) < 3,
"For -E or -T, binary input data (-bi) must have at least 3 columns\n");
n_errors += gmt_M_check_condition (GMT, Ctrl->C.m_origin && GMT->current.proj.projection_GMT != GMT_MERCATOR, "Option -C: Can only give +m for Mercator projection\n");
n_errors += gmt_M_check_condition (GMT, Ctrl->C.m_origin && Ctrl->C.northing != 0.0, "Option -C: Cannot mix +m and the \"northing\" setting\n");

if (!(n_errors || GMT->common.R.active[RSET])) {
GMT->common.R.wesn[XLO] = 0.0; GMT->common.R.wesn[XHI] = 360.0;
Expand All @@ -777,7 +786,7 @@ EXTERN_MSC int GMT_mapproject (void *V_API, int mode, void *args) {
uint64_t row, n_read_in_seg, seg, n_read = 0, n = 0, k, n_output = 0;

double x_in = 0.0, y_in = 0.0, d = 0.0, fwd_scale, inv_scale, xtmp, ytmp, *out = NULL;
double xmin, xmax, ymin, ymax, inch_to_unit, unit_to_inch, u_scale, y_out_min;
double xmin, xmax, ymin, ymax, inch_to_unit, unit_to_inch, u_scale, y_out_min, m_standard_y_value = 0.0;
double x_in_min, x_in_max, y_in_min, y_in_max, x_out_min, x_out_max, y_out_max;
double xnear = 0.0, ynear = 0.0, lon_prev = 0, lat_prev = 0, **data = NULL, *in = NULL;
double speed = 0, last_speed = -1.0, extra[MP_COL_N]; /* Max possible extra output columns from -A -G -L -Z */
Expand Down Expand Up @@ -1163,9 +1172,14 @@ EXTERN_MSC int GMT_mapproject (void *V_API, int mode, void *args) {
two = (Ctrl->E.active || (Ctrl->T.active && GMT->current.proj.datum.h_given)) ? 3 : 2; /* # of output points from conversion */

if (Ctrl->C.shift && Ctrl->F.active) { /* Use same units in -C and -F */
Ctrl->C.easting *= u_scale;
Ctrl->C.easting *= u_scale;
Ctrl->C.northing *= u_scale;
}
if (Ctrl->C.m_origin) { /* Use same units in -C and -F */
m_standard_y_value = GMT->current.proj.j_yc; /* Special Mercator adjustment for measuring distance relative to standard latitude */
if (Ctrl->F.active) /* Use same units in -C and -F */
m_standard_y_value *= fwd_scale;
}

if (Ctrl->L.active) { /* Possibly adjust output types */
if (Ctrl->L.mode == GMT_MP_GIVE_FRAC) /* Want fractional point locations */
Expand Down Expand Up @@ -1258,6 +1272,8 @@ EXTERN_MSC int GMT_mapproject (void *V_API, int mode, void *args) {
in[GMT_X] -= Ctrl->C.easting;
in[GMT_Y] -= Ctrl->C.northing;
}
if (Ctrl->C.m_origin) /* Special Mercator adjustment for measuring distance relative to standard latitude */
in[GMT_Y] += m_standard_y_value;
if (Ctrl->N.active) {
out[GMT_X] = in[GMT_X];
out[GMT_Y] = gmt_lat_swap (GMT, in[GMT_Y], lat_mode);
Expand Down Expand Up @@ -1366,6 +1382,8 @@ EXTERN_MSC int GMT_mapproject (void *V_API, int mode, void *args) {
out[GMT_X] += Ctrl->C.easting;
out[GMT_Y] += Ctrl->C.northing;
}
else if (Ctrl->C.m_origin) /* Special Mercator adjustment for measuring distance relative to standard latitude */
out[GMT_Y] -= m_standard_y_value;
if (GMT->current.proj.three_D) {
double xx = out[GMT_X], yy = out[GMT_Y];
gmt_xyz_to_xy (GMT, xx, yy, gmt_z_to_zz (GMT, in[GMT_Z]), &out[GMT_X], &out[GMT_Y]);
Expand Down
Binary file modified test/psbasemap/merc_origin.ps
View file
Open in desktop
Binary file not shown.