Let grdtrack -C also accept a fixed azimuth for all profiles

doc/rst/source/grdtrack.rst

@@ -14,7 +14,7 @@ Synopsis
 
 **gmt grdtrack** [ *table* ] |-G|\ *grd1* |-G|\ *grd2* ...
 [ |-A|\ [**f**\|\ **p**\|\ **m**\|\ **r**\|\ **R**][**+l**] ]
-[ |-C|\ *length*/\ *ds*\ [*/spacing*][**+a**\|\ **v**][**d**\ *deviation*][**l**\|\ **r**] ]
+[ |-C|\ *length*/\ *ds*\ [*/spacing*][**+a**\|\ **v**][**d**\|\**f**\ *value*][**l**\|\ **r**] ]
 [ |-D|\ *dfile* ]
 [ |-E|\ *line* ]
 [ |-F|\ [**+b**][**+n**][**+r**][**+z**\ *z0*] ]
@@ -105,7 +105,7 @@ Optional Arguments
 
 .. _-C:
 
-**-C**\ *length*/\ *ds*\ [*/spacing*][**+a**\|\ **v**][**d**\ *deviation*][**l**\|\ **r**]
+**-C**\ *length*/\ *ds*\ [*/spacing*][**+a**\|\ **v**][**d**\|\**f**\ *value*][**l**\|\ **r**]
     Use input line segments to create an equidistant and (optionally)
     equally-spaced set of crossing profiles along which we sample the
     grid(s) [Default simply samples the grid(s) at the input locations].
@@ -125,11 +125,12 @@ Optional Arguments
     `Units`_ below). The default unit for geographic grids is meter while
     Cartesian grids implies the user unit.  The output columns will be
     *lon*, *lat*, *dist*, *azimuth*, *z1*, *z2*, ..., *zn* (The *zi* are
-    the sampled values for each of the *n* grids). Finally, use **+d** to
+    the sampled values for each of the *n* grids). Use **+d** to
     change the profiles from being orthogonal to the line by the given
     *deviation* [0]. Looking in the direction of the line, a positive *deviation*
     will rotate the crosslines clockwise and a negative one will rotate them
-    counter-clockwise.
+    counter-clockwise.  Finally, you can use **+f** to set a fixed azimuth
+    for all profiles.
 
 .. _-D:
 

src/gmt_constants.h

@@ -352,7 +352,8 @@ enum GMT_enum_tracklayout {
 	GMT_EW_SN = 2,
 	GMT_LEFT_ONLY = 4,
 	GMT_RIGHT_ONLY = 8,
-	GMT_ALTERNATE = 16};
+	GMT_ALTERNATE = 16,
+	GMT_FIXED_AZIM = 32};
 
 /*! Codes for first segment header */
 enum GMT_enum_firstseg {

src/gmt_support.c

@@ -5621,11 +5621,11 @@ GMT_LOCAL struct GMT_DATASET * gmtsupport_crosstracks_spherical (struct GMT_CTRL
 	uint64_t tbl, row, left, right, seg, seg_no, n_tot_cols;
 	size_t n_x_seg = 0, n_x_seg_alloc = 0;
 
-	bool skip_seg_no;
+	bool skip_seg_no, set_fixed_azim = false;
 
 	char buffer[GMT_BUFSIZ] = {""}, seg_name[GMT_BUFSIZ] = {""}, ID[GMT_BUFSIZ] = {""};
 
-	double dist_inc, d_shift, orientation, sign = 1.0, az_cross, x, y;
+	double dist_inc, d_shift, orientation, sign = 1.0, az_cross, x, y, fixed_azim, dist_to_end;
 	double dist_across_seg, angle_radians, across_ds_radians, ds_phase = 0.0, n_cross_float;
 	double Rot[3][3], Rot0[3][3], E[3], P[3], L[3], R[3], T[3], X[3];
 
@@ -5639,6 +5639,11 @@ GMT_LOCAL struct GMT_DATASET * gmtsupport_crosstracks_spherical (struct GMT_CTRL
 	}
 
 	if (mode & GMT_ALTERNATE) sign = -1.0;	/* Survey layout */
+	if (mode & GMT_FIXED_AZIM) {	/* Got fixed azimuth */
+		fixed_azim = deviation;
+		deviation = 0.0;	/* To prevent any deviation below */
+		set_fixed_azim = true;
+	}
 
 	/* Get resampling step size and zone width in degrees */
 
@@ -5649,10 +5654,14 @@ GMT_LOCAL struct GMT_DATASET * gmtsupport_crosstracks_spherical (struct GMT_CTRL
 	cross_length = cross_length / GMT->current.map.dist[GMT_MAP_DIST].scale;	/* Now in meters [or degrees] */
 	n_cross++;	/* Since one more node than increments */
 	n_half_cross = (n_cross % 2) ? (n_cross - 1) / 2 : n_cross / 2;	/* Half-width of points in a cross profile depending on odd/even */
-	if (unit && strchr (GMT_ARC_UNITS, unit))	/* Gave increments in arc units (already in degrees at this point) */
+	if (unit && strchr (GMT_ARC_UNITS, unit)) {	/* Gave increments in arc units (already in degrees at this point) */
 		across_ds_radians = D2R * cross_length / (n_cross - 1);	/* Angular change from point to point */
-	else	/* Must convert distances to degrees */
+		dist_to_end = 0.5 * cross_length;	/* Distance from center to end of profile in degrees */
+	}
+	else {	/* Must convert distances to degrees */
+		dist_to_end = 0.5 * cross_length / GMT->current.proj.DIST_M_PR_DEG;		/* Distance from center to end of profile in degrees */
 		across_ds_radians = D2R * (cross_length / GMT->current.proj.DIST_M_PR_DEG) / (n_cross - 1);	/* Angular change from point to point */
+	}
 	if ((n_cross % 2) == 0) ds_phase = 0.5;
 	k_start = -n_half_cross;
 	k_stop = k_start + n_cross - 1;
@@ -5686,15 +5695,23 @@ GMT_LOCAL struct GMT_DATASET * gmtsupport_crosstracks_spherical (struct GMT_CTRL
 
 				x = Tin->segment[seg]->data[GMT_X][row];	y = Tin->segment[seg]->data[GMT_Y][row];	/* Reset since now we want lon/lat regardless of grid format */
 				gmt_geo_to_cart (GMT, y, x, P, true);		/* 3-D vector of current point P */
-				left = (row) ? row - 1 : row;			/* Left point (if there is one) */
-				x = Tin->segment[seg]->data[GMT_X][left];	y = Tin->segment[seg]->data[GMT_Y][left];
-				gmt_geo_to_cart (GMT, y, x, L, true);		/* 3-D vector of left point L */
-				right = (row < (Tin->segment[seg]->n_rows-1)) ? row + 1 : row;	/* Right point (if there is one) */
-				x = Tin->segment[seg]->data[GMT_X][right];	y = Tin->segment[seg]->data[GMT_Y][right];
-				gmt_geo_to_cart (GMT, y, x, R, true);		/* 3-D vector of right point R */
-				gmt_cross3v (GMT, L, R, T);			/* Get pole T of plane trough L and R (and center of Earth) */
-				gmt_normalize3v (GMT, T);			/* Make sure T has unit length */
-				gmt_cross3v (GMT, T, P, E);			/* Get pole E to plane trough P normal to L,R (hence going trough P) */
+				if (set_fixed_azim) {	/* Need 2nd point in azim direction to cross with P */
+					double b_x, b_y;
+					gmt_translate_point (GMT, x, y, fixed_azim, dist_to_end, &b_x, &b_y, NULL);
+					gmt_geo_to_cart (GMT, b_y, b_x, R, true);		/* 3-D vector of end point R */
+					gmt_cross3v (GMT, P, R, E);			/* Get pole E to plane trough P and R */
+				}
+				else {	/* Erect an orthogonal (- deviation) profile */
+					left = (row) ? row - 1 : row;			/* Left point (if there is one) */
+					x = Tin->segment[seg]->data[GMT_X][left];	y = Tin->segment[seg]->data[GMT_Y][left];
+					gmt_geo_to_cart (GMT, y, x, L, true);		/* 3-D vector of left point L */
+					right = (row < (Tin->segment[seg]->n_rows-1)) ? row + 1 : row;	/* Right point (if there is one) */
+					x = Tin->segment[seg]->data[GMT_X][right];	y = Tin->segment[seg]->data[GMT_Y][right];
+					gmt_geo_to_cart (GMT, y, x, R, true);		/* 3-D vector of right point R */
+					gmt_cross3v (GMT, L, R, T);			/* Get pole T of plane trough L and R (and center of Earth) */
+					gmt_normalize3v (GMT, T);			/* Make sure T has unit length */
+					gmt_cross3v (GMT, T, P, E);			/* Get pole E to plane trough P normal to L,R (hence going trough P) */
+				}
 				gmt_normalize3v (GMT, E);			/* Make sure E has unit length */
 				if (!gmt_M_is_zero (deviation)) {	/* Must now rotate E about P by the deviation first */
 					gmtlib_load_rot_matrix (-D2R * deviation, Rot0, P);	/* Negative so that rotation goes clockwise */
@@ -5779,6 +5796,8 @@ GMT_LOCAL struct GMT_DATASET * gmtsupport_crosstracks_cartesian (struct GMT_CTRL
 	 * Dout is the new data set with all the crossing profiles;
 	*/
 
+	bool set_fixed_azim = false;
+
 	int k, ndig, sdig, n_half_cross, k_start, k_stop;
 
 	unsigned int ii, np_cross;
@@ -5803,6 +5822,11 @@ GMT_LOCAL struct GMT_DATASET * gmtsupport_crosstracks_cartesian (struct GMT_CTRL
 	}
 
 	if (mode & GMT_ALTERNATE) sign = -1.0;	/* Survey mode */
+	if (mode & GMT_FIXED_AZIM) {	/* Got fixed azimuth */
+		az_cross = deviation;
+		deviation = 0.0;	/* To prevent any deviation below */
+		set_fixed_azim = true;
+	}
 
 	/* Get resampling step size and zone width in degrees */
 
@@ -5837,7 +5861,7 @@ GMT_LOCAL struct GMT_DATASET * gmtsupport_crosstracks_cartesian (struct GMT_CTRL
 				GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Working on cross profile %" PRIu64 " [local line orientation = %06.1f]\n", row, orientation);
 
 				x = Tin->segment[seg]->data[GMT_X][row];	y = Tin->segment[seg]->data[GMT_Y][row];	/* Reset since now we want lon/lat regardless of grid format */
-				az_cross = fmod (Tin->segment[seg]->data[SEG_AZIM][row] + 270.0, 360.0);	/* Azimuth of cross-profile in 0-360 range */
+				if (!set_fixed_azim) az_cross = fmod (Tin->segment[seg]->data[SEG_AZIM][row] + 270.0, 360.0);	/* Azimuth of cross-profile in 0-360 range */
 				if (mode & GMT_ALTERNATE)
 					sign = -sign;
 				else if (mode & GMT_EW_SN)

src/grdtrack.c

@@ -164,7 +164,7 @@ static void Free_Ctrl (struct GMT_CTRL *GMT, struct GRDTRACK_CTRL *C) {	/* Deall
 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_Usage (API, 0, "usage: %s -G%s -G<grid2> ... [<table>] [-A[f|m|p|r|R][+l]] [-C<length>/<ds>[/<spacing>][+a|v][+d<deviation>][+l|r]] "
+	GMT_Usage (API, 0, "usage: %s -G%s -G<grid2> ... [<table>] [-A[f|m|p|r|R][+l]] [-C<length>/<ds>[/<spacing>][+a|v][+d|f<value>][+l|r]] "
 		"[-D<dfile>] [-E<line1>[,<line2>,...][+a<az>][+c][+d][+g][+i<step>][+l<length>][+n<np][+o<az>][+r<radius>]] "
 		"[-F[+b][+n][+r][+z<z0>]] [-N] [%s] [-S[a|l|L|m|p|u|U][+a][+c][+d][+r][+s[<file>]]] [-T<radius>>[+e|p]] [%s] [-Z] [%s] [%s] [%s] "
 		"[%s] [%s] [%s] [%s] [%s] [%s] [%s] %s] [%s] [%s] [%s] [%s]\n",
@@ -190,7 +190,7 @@ static int usage (struct GMTAPI_CTRL *API, int level) {
 	GMT_Usage (API, 3, "r: Resample at equidistant locations; input points not necessarily included.");
 	GMT_Usage (API, 3, "R: Same, but adjust given spacing to fit the track length exactly.");
 	GMT_Usage (API, -2, "Append +l to compute distances along rhumblines (loxodromes) [no].");
-	GMT_Usage (API, 1, "\n-C<length>/<ds>[/<spacing>][+a|v][+d<deviation>][+l|r]");
+	GMT_Usage (API, 1, "\n-C<length>/<ds>[/<spacing>][+a|v][+d|f<value>][+l|r]");
 	GMT_Usage (API, -2, "Create equidistant cross-profiles from input line segments. Append 2-3 parameters:");
 	GMT_Usage (API, 3, "1. <length>: The full-length of each cross profile. "
 		"Append distance unit u (%s); this unit also applies to <ds> and <spacing>. "
@@ -200,6 +200,7 @@ static int usage (struct GMTAPI_CTRL *API, int level) {
 	GMT_Usage (API, -2, "Several modifiers controls the creation of the profiles:");
 	GMT_Usage (API, 3, "+a Alternate the direction of cross-profiles [Default orients all the same way].");
 	GMT_Usage (API, 3, "+d Set deviation (-90/+90 range) from orthogonal cross-profiles [0 (no deviation)].");
+	GMT_Usage (API, 3, "+f Set fixed azimuth for all cross-profiles.");
 	GMT_Usage (API, 3, "+l Only use the left half of the profiles [entire profile].");
 	GMT_Usage (API, 3, "+r Only use the right half of the profiles [entire profile].");
 	GMT_Usage (API, 3, "+v Adjust direction of cross-profiles for E-W or S-N viewing [Default orients all the same way].");
@@ -338,11 +339,12 @@ static int parse (struct GMT_CTRL *GMT, struct GRDTRACK_CTRL *Ctrl, struct GMT_O
 			case 'C':	/* Create cross profiles */
 				n_errors += gmt_M_repeated_module_option (API, Ctrl->C.active);
 				Ctrl->C.active = true;
-				if ((c = gmt_first_modifier (GMT, opt->arg, "adlrv"))) {	/* Process any modifiers */
+				if ((c = gmt_first_modifier (GMT, opt->arg, "adflrv"))) {	/* Process any modifiers */
 					pos = 0;	/* Reset to start of new word */
-					while (gmt_getmodopt (GMT, 'C', c, "adlrv", &pos, p, &n_errors) && n_errors == 0) {
+					while (gmt_getmodopt (GMT, 'C', c, "adflrv", &pos, p, &n_errors) && n_errors == 0) {
 						switch (p[0]) {
 							case 'a': Ctrl->C.mode |= GMT_ALTERNATE; break;		/* Select alternating direction of cross-profiles */
+							case 'f': Ctrl->C.mode |= GMT_FIXED_AZIM;	/* Deliberatily fall through to get fixed azimuth */
 							case 'd': Ctrl->C.deviation = atof (&p[1]); break;		/* Less than orthogonal vy given deviation */
 							case 'l': Ctrl->C.mode |= GMT_LEFT_ONLY; break;		/* cross-profile starts at line and go to left side only */
 							case 'r': Ctrl->C.mode |= GMT_RIGHT_ONLY; break;	/* cross-profile starts at line and go to right side only */

test/grdtrack/layout_fixed.sh

@@ -0,0 +1,24 @@
+#!/usr/bin/env bash
+ps=layout_fixed.ps
+function plot_it {
+	gmt grdtrack @Matthews_2016_subduction_subset.txt -Gdummy.nc -C500k/100/100${1} > sz_pol_left.gmt
+	awk '{ if ( $1 == ">") print $0 ; else if ($3 <= 0) print $0 }' sz_pol_left.gmt > Lhalfxprofile.gmt
+	awk '{ if ( $1 == ">") print $0 ; else if ($3 >= 0) print $0 }' sz_pol_left.gmt > Rhalfxprofile.gmt
+	gmt psxy -R -J -O -K -W1p @Matthews_2016_subduction_subset.txt
+	gmt psxy -R -J -O Lhalfxprofile.gmt -W1p,red -K
+	gmt psxy -R -J Rhalfxprofile.gmt -W1p,blue -O -K
+	gmt psxy -R -J Rhalfxprofile.gmt -Sc0.4c -Gblue -O -K
+	gmt psxy -R -J Lhalfxprofile.gmt -Sc0.4c -Gred -O -K
+	gmt pstext -R -J -O -K sz_pol_left.gmt -F+f9p,Helvetica,white+r0
+	gmt psxy -R -J -O -K -Ss0.25i -Gblack @Matthews_2016_subduction_subset.txt
+	gmt pstext -R -J -O -K @Matthews_2016_subduction_subset.txt -F+f12p,Helvetica,white+r0
+	echo $1 | gmt pstext -R -J -O -K -F+cTL+jTL+f18p -Dj0.2i -Gwhite
+}
+gmt grdmath -R135/162/42/51 -I1 X = dummy.nc
+gmt psbasemap -R135/162/42/51 -JM6i -Bafg -BWSnE -P -K -Xc > $ps
+plot_it +f130 >> $ps
+gmt psbasemap -R -J -O -Bafg -BWsnE -K -Y3.1i >> $ps
+plot_it +f90 >> $ps
+gmt psbasemap -R -J -O -Bafg -BWsnE -K -Y3.1i >> $ps
+plot_it >> $ps
+gmt psxy -R -J -O -T >> $ps