Skip to content

Commit 08a50c8

Browse files
PaulWesselseisman
authored andcommitted
[backport-6.1] Fix issues in Mercator projection (#4120)
* Requires special -C +m modifier to report Mercator units relative to standard parallel Normallly, there is no such thing as Mercator should report relative to Equator. However, when you wish to make a Mercator map centered on a point and give spacial distances then we temporarily need to ability to measure Mercator distances from that point. * Document -C+m * Update gmt_proj.c * Update gmt_proj.c * Update merc_origin.ps
1 parent b82cfb9 commit 08a50c8

File tree

5 files changed

+40
-16
lines changed

5 files changed

+40
-16
lines changed

doc/rst/source/mapproject.rst

Lines changed: 5 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -15,7 +15,7 @@ Synopsis
1515
**gmt mapproject** [ *tables* ] |-J|\ *parameters*
1616
|SYN_OPT-R|
1717
[ |-A|\ **b**\|\ **B**\|\ **f**\|\ **F**\|\ **o**\|\ **O**\ [*lon0*/*lat0*][**+v**] ]
18-
[ |-C|\ [*dx*/*dy*] ]
18+
[ |-C|\ [*dx*/*dy*][**+m**] ]
1919
[ |-D|\ **c**\|\ **i**\|\ **p** ]
2020
[ |-E|\ [*datum*] ] [ |-F|\ [*unit*] ]
2121
[ |-G|\ [*lon0*/*lat0*][**+a**][**+i**][**+u**\ *unit*][**+v**] ]
@@ -105,15 +105,17 @@ Optional Arguments
105105

106106
.. _-C:
107107

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

118120
.. _-D:
119121

src/gmt_init.c

Lines changed: 5 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -853,7 +853,7 @@ GMT_LOCAL int gmtinit_rectR_to_geoR (struct GMT_CTRL *GMT, char unit, double rec
853853
uint64_t dim[GMT_DIM_SIZE] = {1, 1, 2, 2}; /* Just a single data table with one segment with two 2-column records */
854854
bool was_R, was_J;
855855
double wesn[4];
856-
char buffer[GMT_LEN256] = {""}, Jstring[GMT_LEN128] = {""}, in_string[GMT_VF_LEN] = {""}, out_string[GMT_VF_LEN] = {""}, *v = NULL;
856+
char buffer[GMT_LEN256] = {""}, Jstring[GMT_LEN128] = {""}, in_string[GMT_VF_LEN] = {""}, out_string[GMT_VF_LEN] = {""}, origin_flag[4] = {""}, *v = NULL;
857857
struct GMT_DATASET *In = NULL, *Out = NULL;
858858

859859
GMT_Report (GMT->parent, GMT_MSG_DEBUG, "Call gmtinit_rectR_to_geoR to convert projected -R to geo -R\n");
@@ -899,6 +899,8 @@ GMT_LOCAL int gmtinit_rectR_to_geoR (struct GMT_CTRL *GMT, char unit, double rec
899899
GMT_Report (GMT->parent, GMT_MSG_ERROR, "UTM projection insufficiently specified to auto-determine geographic region\n");
900900
return (GMT_MAP_NO_PROJECTION);
901901
}
902+
if (GMT->current.proj.projection_GMT == GMT_MERCATOR) /* Special use of Mercator units relative to stated origin */
903+
strcpy (origin_flag, "+m");
902904
break;
903905
case 2: /* Conical: Use default patch */
904906
break;
@@ -936,8 +938,8 @@ GMT_LOCAL int gmtinit_rectR_to_geoR (struct GMT_CTRL *GMT, char unit, double rec
936938
v[0] = '\0';
937939
}
938940
}
939-
snprintf (buffer, GMT_LEN256, "-R%g/%g/%g/%g -J%s -I -F%c -C -bi2d -bo2d -<%s ->%s --GMT_HISTORY=false",
940-
wesn[XLO], wesn[XHI], wesn[YLO], wesn[YHI], Jstring, unit, in_string, out_string);
941+
snprintf (buffer, GMT_LEN256, "-R%g/%g/%g/%g -J%s -I -F%c -C%s -bi2d -bo2d -<%s ->%s --GMT_HISTORY=false",
942+
wesn[XLO], wesn[XHI], wesn[YLO], wesn[YHI], Jstring, unit, origin_flag, in_string, out_string);
941943
if (get_R) GMT_Report (GMT->parent, GMT_MSG_DEBUG, "Obtaining geographic corner coordinates via mapproject %s\n", buffer);
942944
if (GMT_Call_Module (GMT->parent, "mapproject", GMT_MODULE_CMD, buffer) != GMT_OK) /* Get the corners in degrees via mapproject */
943945
return (GMT->parent->error);

src/gmt_proj.c

Lines changed: 8 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -652,15 +652,17 @@ GMT_LOCAL void gmtproj_ipolar (struct GMT_CTRL *GMT, double *x, double *y, doubl
652652

653653
/* -JM MERCATOR PROJECTION */
654654

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

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

660660
GMT->current.proj.central_meridian = lon0;
661-
GMT->current.proj.j_x = cosd (slat) / d_sqrt (1.0 - GMT->current.proj.ECC2 * sind (slat) * sind (slat)) * GMT->current.proj.EQ_RAD;
661+
/* Need geodetic latitude in this expression: */
662+
GMT->current.proj.j_x = cosd (lat0) / d_sqrt (1.0 - GMT->current.proj.ECC2 * sind (lat0) * sind (lat0)) * GMT->current.proj.EQ_RAD;
662663
GMT->current.proj.j_ix = 1.0 / GMT->current.proj.j_x;
663-
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;
664+
/* Need conformal latitude in this expression (same as in gmtproj_merc_sph) */
665+
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;
664666
}
665667

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

674676
*x = GMT->current.proj.j_x * D2R * lon;
675-
*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);
677+
*y = (fabs (lat) < 90.0) ? GMT->current.proj.j_x * d_log (GMT, tand (45.0 + 0.5 * lat)) : copysign (DBL_MAX, lat);
676678
}
677679

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

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

src/mapproject.c

Lines changed: 22 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -97,9 +97,10 @@ struct MAPPROJECT_CTRL { /* All control options for this program (except common
9797
unsigned int mode;
9898
double lon, lat; /* Fixed point of reference */
9999
} A;
100-
struct MAPPROJECT_C { /* -C[<false_easting>/<false_northing>] */
100+
struct MAPPROJECT_C { /* -C[<false_easting>/<false_northing>][+m] */
101101
bool active;
102102
bool shift;
103+
bool m_origin; /* True if we want projected Mercator y-values relative to standard latitude [Equator] */
103104
double easting, northing; /* Shifts */
104105
} C;
105106
struct MAPPROJECT_D { /* -D<c|i|p> */
@@ -189,7 +190,7 @@ static int usage (struct GMTAPI_CTRL *API, int level) {
189190
const char *name = gmt_show_name_and_purpose (API, THIS_MODULE_LIB, THIS_MODULE_CLASSIC_NAME, THIS_MODULE_PURPOSE);
190191
if (level == GMT_MODULE_PURPOSE) return (GMT_NOERROR);
191192
GMT_Message (API, GMT_TIME_NONE, "usage: %s <table> %s %s\n", name, GMT_J_OPT, GMT_Rgeo_OPT);
192-
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);
193+
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);
193194
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);
194195
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",
195196
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);
@@ -209,6 +210,7 @@ static int usage (struct GMTAPI_CTRL *API, int level) {
209210
GMT_Message (API, GMT_TIME_NONE, "\t-C Return x/y relative to projection center [Default is relative to lower left corner].\n");
210211
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");
211212
GMT_Message (API, GMT_TIME_NONE, "\t Units are plot units unless -F is set in which case the unit is meters.\n");
213+
GMT_Message (API, GMT_TIME_NONE, "\t Mercator only: Append +m to set origin for projected y-values to the standard parallel [Equator].\n");
212214
GMT_Message (API, GMT_TIME_NONE, "\t-D Temporarily reset PROJ_LENGTH_UNIT to be c (cm), i (inch), or p (point).\n");
213215
GMT_Message (API, GMT_TIME_NONE, "\t Cannot be used if -F is set.\n");
214216
GMT_Message (API, GMT_TIME_NONE, "\t-E Convert (lon, lat, h) to Earth Centered Earth Fixed (ECEF) coordinates [-I for inverse].\n");
@@ -482,12 +484,17 @@ static int parse (struct GMT_CTRL *GMT, struct MAPPROJECT_CTRL *Ctrl, struct GMT
482484
break;
483485
case 'C':
484486
Ctrl->C.active = true;
487+
if (opt->arg[0] && (p = strstr (opt->arg, "+m"))) { /* Gave +m for reset relative offsets to Mercator origin */
488+
Ctrl->C.m_origin = true;
489+
p[0] = '\0'; /* Temporarily chop off modifier */
490+
}
485491
if (opt->arg[0]) { /* Also gave shifts */
486492
n_errors += gmt_M_check_condition (GMT, sscanf (opt->arg, "%lf/%lf", &Ctrl->C.easting, &Ctrl->C.northing) != 2,
487493
"Option -C: Expected -C[<false_easting>/<false_northing>]\n");
488494
Ctrl->C.shift = true;
489495
}
490496
will_need_RJ = true; /* Since -C is used with projections only */
497+
if (p) p[0] = '+'; /* Restore modifier */
491498
break;
492499
case 'D':
493500
Ctrl->D.active = true;
@@ -752,6 +759,8 @@ static int parse (struct GMT_CTRL *GMT, struct MAPPROJECT_CTRL *Ctrl, struct GMT
752759
n_errors += gmt_M_check_condition (GMT, ((Ctrl->T.active && GMT->current.proj.datum.h_given) || Ctrl->E.active) &&
753760
GMT->common.b.active[GMT_IN] && gmt_get_cols (GMT, GMT_IN) < 3,
754761
"For -E or -T, binary input data (-bi) must have at least 3 columns\n");
762+
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");
763+
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");
755764

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

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

11651174
if (Ctrl->C.shift && Ctrl->F.active) { /* Use same units in -C and -F */
1166-
Ctrl->C.easting *= u_scale;
1175+
Ctrl->C.easting *= u_scale;
11671176
Ctrl->C.northing *= u_scale;
11681177
}
1178+
if (Ctrl->C.m_origin) { /* Use same units in -C and -F */
1179+
m_standard_y_value = GMT->current.proj.j_yc; /* Special Mercator adjustment for measuring distance relative to standard latitude */
1180+
if (Ctrl->F.active) /* Use same units in -C and -F */
1181+
m_standard_y_value *= fwd_scale;
1182+
}
11691183

11701184
if (Ctrl->L.active) { /* Possibly adjust output types */
11711185
if (Ctrl->L.mode == GMT_MP_GIVE_FRAC) /* Want fractional point locations */
@@ -1258,6 +1272,8 @@ EXTERN_MSC int GMT_mapproject (void *V_API, int mode, void *args) {
12581272
in[GMT_X] -= Ctrl->C.easting;
12591273
in[GMT_Y] -= Ctrl->C.northing;
12601274
}
1275+
if (Ctrl->C.m_origin) /* Special Mercator adjustment for measuring distance relative to standard latitude */
1276+
in[GMT_Y] += m_standard_y_value;
12611277
if (Ctrl->N.active) {
12621278
out[GMT_X] = in[GMT_X];
12631279
out[GMT_Y] = gmt_lat_swap (GMT, in[GMT_Y], lat_mode);
@@ -1366,6 +1382,8 @@ EXTERN_MSC int GMT_mapproject (void *V_API, int mode, void *args) {
13661382
out[GMT_X] += Ctrl->C.easting;
13671383
out[GMT_Y] += Ctrl->C.northing;
13681384
}
1385+
else if (Ctrl->C.m_origin) /* Special Mercator adjustment for measuring distance relative to standard latitude */
1386+
out[GMT_Y] -= m_standard_y_value;
13691387
if (GMT->current.proj.three_D) {
13701388
double xx = out[GMT_X], yy = out[GMT_Y];
13711389
gmt_xyz_to_xy (GMT, xx, yy, gmt_z_to_zz (GMT, in[GMT_Z]), &out[GMT_X], &out[GMT_Y]);

test/psbasemap/merc_origin.ps

-5 Bytes
Binary file not shown.

0 commit comments

Comments
 (0)