Add polar resampling via -Ar|t in plot, plot3d, and clip

doc/rst/source/clip.rst

@@ -14,7 +14,7 @@ Synopsis
 
 **gmt clip** [ *table* ] |-J|\ *parameters* |-C|\ [*n*]
 |SYN_OPT-Rz|
-[ |-A|\ [**m**\|\ **p**\|\ **x**\|\ **y**] ]
+[ |-A|\ [**m**\|\ **p**\|\ **x**\|\ **y**\|\ **r**\|\ **t**] ]
 [ |SYN_OPT-B| ]
 |-J|\ **z**\|\ **Z**\ *parameters* ]
 [ |-N| ]

doc/rst/source/clip_common.rst_

@@ -52,15 +52,16 @@ Optional Arguments
 
 .. _-A:
 
-**-A**\ [**m**\|\ **p**\|\ **x**\|\ **y**]
+**-A**\ [**m**\|\ **p**\|\ **x**\|\ **y**\|\ **r**\|\ **t**]
     By default, geographic line segments are connected as great circle arcs. To connect them as
     straight lines, use the **-A** flag. Alternatively, add **m** to connect
     the line by first following a meridian, then a parallel. Or append **p**
     to start following a parallel, then a meridian. (This can be practical
     to connect lines along parallels, for example).
     For Cartesian data, points are simply connected, unless you append
     **x** or **y** to construct stair-case paths whose first move is along
-    *x* or *y*, respectively.
+    *x* or *y*, respectively. For polar projection, append **r** or **t** to
+    draw stair-case curves that whose first move is along *r* or *theta*, respectively.
 
 .. |Add_-B| replace:: |Add_-B_links|
 .. include:: explain_-B.rst_

doc/rst/source/plot.rst

@@ -14,7 +14,7 @@ Synopsis
 
 **gmt plot** [ *table* ] |-J|\ *parameters*
 |SYN_OPT-Rz|
-[ |-A|\ [**m**\|\ **p**\|\ **x**\|\ **y**] ]
+[ |-A|\ [**m**\|\ **p**\|\ **x**\|\ **y**\|\ **r**\|\ **t**] ]
 [ |SYN_OPT-B| ]
 [ |-C|\ *cpt* ]
 [ |-D|\ *dx*/*dy* ]

doc/rst/source/plot3d.rst

@@ -15,7 +15,7 @@ Synopsis
 **gmt plot3d** [ *table* ] |-J|\ *parameters*
 |-Jz|\ \|\ **Z**\ *parameters*
 |SYN_OPT-Rz|
-[ |-A|\ [**m**\|\ **p**\|\ **x**\|\ **y**] ]
+[ |-A|\ [**m**\|\ **p**\|\ **x**\|\ **y**\|\ **r**\|\ **t**] ]
 [ |SYN_OPT-B| ]
 [ |-C|\ *cpt* ]
 [ |-D|\ *dx*/*dy*\ [/*dz*] ]

doc/rst/source/plot3d_common.rst_

@@ -44,15 +44,17 @@ Optional Arguments
 
 .. _-A:
 
-**-A**\ [**m**\|\ **p**\|\ **x**\|\ **y**]
+**-A**\ [**m**\|\ **p**\|\ **x**\|\ **y**\|\ **r**\|\ **t**]
     By default, geographic line segments are drawn as great circle arcs by resampling
     coarse input data along such arcs. To disable this sampling and draw them as
     straight lines, use the **-A** flag.  Alternatively, add **m** to draw
     the line by first following a meridian, then a parallel. Or append **p**
     to start following a parallel, then a meridian. (This can be practical
     to draw a line along parallels, for example).  For Cartesian data, points
     are simply connected, unless you append **x** or **y** to draw stair-case
-    curves that whose first move is along *x* or *y*, respectively. **Note**:
+    curves that whose first move is along *x* or *y*, respectively. For polar
+    projection, append **r** or **t** to draw stair-case curves that whose first
+    move is along *r* or *theta*, respectively. **Note**:
     The **-A** option requires constant *z*-coordinates.
 
 .. |Add_-B| replace:: |Add_-B_links|

doc/rst/source/plot_common.rst_

@@ -34,15 +34,17 @@ Optional Arguments
 
 .. _-A:
 
-**-A**\ [**m**\|\ **p**\|\ **x**\|\ **y**]
+**-A**\ [**m**\|\ **p**\|\ **x**\|\ **y**\|\ **r**\|\ **t**]
     By default, geographic line segments are drawn as great circle arcs by resampling
     coarse input data along such arcs. To disable this sampling and draw them as
     straight lines, use the **-A** flag.  Alternatively, add **m** to draw
     the line by first following a meridian, then a parallel. Or append **p**
     to start following a parallel, then a meridian. (This can be practical
     to draw a line along parallels, for example).  For Cartesian data, points
     are simply connected, unless you append **x** or **y** to draw stair-case
-    curves that whose first move is along *x* or *y*, respectively.
+    curves that whose first move is along *x* or *y*, respectively. For polar
+    projection, append **r** or **t** to draw stair-case curves that whose first
+    move is along *r* or *theta*, respectively.
 
 .. |Add_-B| replace:: |Add_-B_links|
 .. include:: explain_-B.rst_

doc/rst/source/psclip.rst

@@ -14,7 +14,7 @@ Synopsis
 
 **gmt psclip** [ *table* ] |-J|\ *parameters* |-C|\ [*n*]
 |SYN_OPT-Rz|
-[ |-A|\ [**m**\|\ **p**\|\ **x**\|\ **y**] ]
+[ |-A|\ [**m**\|\ **p**\|\ **x**\|\ **y**\|\ **r**\|\ **t**] ]
 [ |SYN_OPT-B| ]
 [ |-K| ] [ |-N| ] [ |-O| ]
 [ |-P| ] [ |-T| ]

doc/rst/source/psxy.rst

@@ -14,7 +14,7 @@ Synopsis
 
 **gmt psxy** [ *table* ] |-J|\ *parameters*
 |SYN_OPT-Rz|
-[ |-A|\ [**m**\|\ **p**\|\ **x**\|\ **y**] ]
+[ |-A|\ [**m**\|\ **p**\|\ **x**\|\ **y**\|\ **r**\|\ **t**] ]
 [ |SYN_OPT-B| ]
 [ |-C|\ *cpt* ]
 [ |-D|\ *dx*/*dy* ]

doc/rst/source/psxyz.rst

@@ -15,7 +15,7 @@ Synopsis
 **gmt psxyz** [ *table* ] |-J|\ *parameters*
 |-J|\ **z**\|\ **Z**\ *parameters*
 |SYN_OPT-Rz|
-[ |-A|\ [**m**\|\ **p**\|\ **x**\|\ **y**] ]
+[ |-A|\ [**m**\|\ **p**\|\ **x**\|\ **y**\|\ **r**\|\ **t**] ]
 [ |SYN_OPT-B| ]
 [ |-C|\ *cpt* ]
 [ |-D|\ *dx*/*dy*\ [/*dz*] ]

src/gmt_vector.c

@@ -340,6 +340,61 @@ GMT_LOCAL uint64_t gmtvector_fix_up_path_cartonly (struct GMT_CTRL *GMT, double
 	return (n_new);
 }
 
+GMT_LOCAL uint64_t gmtvector_fix_up_path_polar (struct GMT_CTRL *GMT, double **a_x, double **a_y, uint64_t n, double step, unsigned int mode) {
+	/* Takes pointers to a list of <n> theta/r pairs (in user units) and adds
+	 * auxiliary points if the distance between two given points exceeds
+	 * <step> units.
+	 * If mode=0: returns points along a straight line
+	 * If mode=1: staircase; first follows r, then theta
+	 * If mode=2: staircase; first follows theta, then r
+	 * Returns the new number of points (original plus auxiliary).
+	 */
+
+	uint64_t i, j, n_new, n_step = 0;
+	double *x = NULL, *y = NULL, c;
+
+	if (mode == GMT_STAIRS_OFF) return n;	/* Nothing to do */
+
+	x = *a_x;	y = *a_y;
+
+	gmt_prep_tmp_arrays (GMT, GMT_NOTSET, 1, 2);	/* Init or reallocate tmp vectors */
+	GMT->hidden.mem_coord[GMT_X][0] = x[0];	GMT->hidden.mem_coord[GMT_Y][0] = y[0];	n_new = 1;
+	if (step <= 0.0) step = 1.0;	/* Sanity valve; if step not given we set it to 1 */
+
+	for (i = 1; i < n; i++) {
+		if (mode == GMT_STAIRS_X) {	/* First follow theta, then r */
+			n_step = lrint (fabs (x[i] - x[i-1]) / step);
+			for (j = 1; j <= n_step; j++) {
+				c = j / (double)n_step;
+				gmt_prep_tmp_arrays (GMT, GMT_NOTSET, n_new, 2);	/* Init or reallocate tmp read vectors */
+				GMT->hidden.mem_coord[GMT_X][n_new] = x[i-1] * (1 - c) + x[i] * c;
+				GMT->hidden.mem_coord[GMT_Y][n_new] = y[i-1];
+				n_new++;
+			}
+			GMT->hidden.mem_coord[GMT_X][n_new] = x[i];	GMT->hidden.mem_coord[GMT_Y][n_new] = y[i];	n_new++;
+		}
+		else if (mode == GMT_STAIRS_Y) {	/* First follow r, then theta */
+			GMT->hidden.mem_coord[GMT_X][n_new] = x[i-1];	GMT->hidden.mem_coord[GMT_Y][n_new] = y[i];	n_new++;
+			n_step = lrint (fabs (x[i] - x[i-1]) / step);
+			for (j = 1; j <= n_step; j++) {
+				c = j / (double)n_step;
+				gmt_prep_tmp_arrays (GMT, GMT_NOTSET, n_new, 2);	/* Init or reallocate tmp read vectors */
+				GMT->hidden.mem_coord[GMT_X][n_new] = x[i-1] * (1 - c) + x[i] * c;
+				GMT->hidden.mem_coord[GMT_Y][n_new] = y[i];
+				n_new++;
+			}
+		}
+		gmt_prep_tmp_arrays (GMT, GMT_NOTSET, n_new, 2);	/* Init or reallocate tmp read vectors */
+	}
+
+	/* Destroy old allocated memory and put the new one in place */
+	gmt_M_free (GMT, x);	gmt_M_free (GMT, y);
+	*a_x = gmtlib_assign_vector (GMT, n_new, GMT_X);
+	*a_y = gmtlib_assign_vector (GMT, n_new, GMT_Y);
+
+	return (n_new);
+}
+
 GMT_LOCAL uint64_t gmtvector_fix_up_path_cartesian (struct GMT_CTRL *GMT, double **a_x, double **a_y, uint64_t n, double step, unsigned int mode) {
 	/* Takes pointers to a list of <n> x/y pairs (in user units) and adds
 	 * auxiliary points if the distance between two given points exceeds
@@ -1429,6 +1484,8 @@ uint64_t gmt_fix_up_path (struct GMT_CTRL *GMT, double **a_lon, double **a_lat,
 	double c, d, fraction, theta, minlon, maxlon;
 	double dlon, lon_i, boost, f_lat_a, f_lat_b;
 
+	if (GMT->current.proj.projection_GMT == GMT_POLAR) return (gmtvector_fix_up_path_polar (GMT, a_lon, a_lat, n, 1.0, mode));	/* r-theta stepping */
+
 	if (gmt_M_is_cartesian (GMT, GMT_IN)) return (gmtvector_fix_up_path_cartonly (GMT, a_lon, a_lat, n, mode));	/* Stair case only */
 
 	lon = *a_lon;	lat = *a_lat;	/* Input arrays */

src/psclip.c

@@ -36,7 +36,7 @@
 #define THIS_MODULE_OPTIONS "-:>BJKOPRUVXYbdefghipqstxy" GMT_OPT("EZMmc")
 
 struct PSCLIP_CTRL {
-	struct PSCLIP_A {	/* -A[m|p|step] */
+	struct PSCLIP_A {	/* -A[m|y|p|x|r|t<step>] */
 		bool active;
 		unsigned int mode;
 		double step;
@@ -76,7 +76,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 -C[a|<n>] [-A[m|p|x|y]] [-K] [-O]  OR\n", name);
+	GMT_Message (API, GMT_TIME_NONE, "usage: %s -C[a|<n>] [-A[m|p|x|y|r|t]] [-K] [-O]  OR\n", name);
 	GMT_Message (API, GMT_TIME_NONE, "\t%s <table> %s %s [%s]\n", name, GMT_J_OPT, GMT_Rgeoz_OPT, GMT_B_OPT);
 	GMT_Message (API, GMT_TIME_NONE, "\t%s[-N] %s%s[-T] [%s] [%s] [-W[<pen>]\n", API->K_OPT, API->O_OPT, API->P_OPT, GMT_U_OPT, GMT_V_OPT);
 	GMT_Message (API, GMT_TIME_NONE, "\t[%s] [%s] [%s] %s[%s]\n", GMT_X_OPT, GMT_Y_OPT, GMT_bi_OPT, API->c_OPT, GMT_di_OPT);
@@ -92,6 +92,7 @@ static int usage (struct GMTAPI_CTRL *API, int level) {
 	GMT_Message (API, GMT_TIME_NONE, "\t-A Suppress connecting geographic points using great circle arcs, i.e., connect by straight lines,\n");
 	GMT_Message (API, GMT_TIME_NONE, "\t   unless m or p is appended to first follow meridian then parallel, or vice versa.\n");
 	GMT_Message (API, GMT_TIME_NONE, "\t   For Cartesian data, use -Ax or -Ay to connect first in x, then y, or vice versa.\n");
+	GMT_Message (API, GMT_TIME_NONE, "\t   For Polar projections, use -At or -Ar connect first in theta, then r, or vice versa.\n");
 	GMT_Option (API, "<,B-,K");
 	GMT_Message (API, GMT_TIME_NONE, "\t-N Use the outside of the polygons and the map boundary as clip paths.\n");
 	GMT_Option (API, "O,P");
@@ -129,8 +130,8 @@ static int parse (struct GMT_CTRL *GMT, struct PSCLIP_CTRL *Ctrl, struct GMT_OPT
 			case 'A':	/* Turn off draw_arc mode */
 				Ctrl->A.active = true;
 				switch (opt->arg[0]) {
-					case 'm': case 'y': Ctrl->A.mode = GMT_STAIRS_Y; break;
-					case 'p': case 'x': Ctrl->A.mode = GMT_STAIRS_X; break;
+					case 'm': case 'y': case 'r': Ctrl->A.mode = GMT_STAIRS_Y; break;
+					case 'p': case 'x': case 't': Ctrl->A.mode = GMT_STAIRS_X; break;
 #ifdef DEBUG
 					default: Ctrl->A.step = atof (opt->arg); break; /* Undocumented test feature */
 #endif

src/psxy.c

@@ -37,7 +37,7 @@
 
 struct PSXY_CTRL {
 	bool no_RJ_needed;	/* Special case of -T and no -B when -R -J is not required */
-	struct PSXY_A {	/* -A[m|y|p|x|step] */
+	struct PSXY_A {	/* -A[m|y|p|x|r|t<step>] */
 		bool active;
 		unsigned int mode;
 		double step;
@@ -473,7 +473,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);
 	const char *mod_name = &name[4];	/* To skip the leading gmt for usage messages */
 	if (level == GMT_MODULE_PURPOSE) return (GMT_NOERROR);
-	GMT_Message (API, GMT_TIME_NONE, "usage: %s [<table>] %s %s [-A[m|p|x|y]]\n", name, GMT_J_OPT, GMT_Rgeoz_OPT);
+	GMT_Message (API, GMT_TIME_NONE, "usage: %s [<table>] %s %s [-A[m|p|x|y|r|t]]\n", name, GMT_J_OPT, GMT_Rgeoz_OPT);
 	GMT_Message (API, GMT_TIME_NONE, "\t[%s] [-C<cpt>] [-D<dx>/<dy>] [-E[x|y|X|Y][+a][+c[l|f]][+n][+p<pen>][+w<width>]] [-F<arg>] [-G<fill>|+z]\n", GMT_B_OPT);
 	GMT_Message (API, GMT_TIME_NONE, "\t[-H[<scale>]] [-I[<intens>]] %s[-L[+b|d|D][+xl|r|x0][+yb|t|y0][+p<pen>]] [-N[c|r]] %s%s\n", API->K_OPT, API->O_OPT, API->P_OPT);
 	if (API->GMT->current.setting.run_mode == GMT_CLASSIC)	/* -T has no purpose in modern mode */
@@ -492,6 +492,7 @@ static int usage (struct GMTAPI_CTRL *API, int level) {
 	GMT_Message (API, GMT_TIME_NONE, "\t   straight lines unless m or p is appended to first follow meridian\n");
 	GMT_Message (API, GMT_TIME_NONE, "\t   then parallel, or vice versa.\n");
 	GMT_Message (API, GMT_TIME_NONE, "\t   For Cartesian data, use -Ax or -Ay to draw x- or y-staircase curves.\n");
+	GMT_Message (API, GMT_TIME_NONE, "\t   For Polar projections, use -At or -Ar to draw theta- or r-staircase curves.\n");
 	GMT_Option (API, "B-");
 	GMT_Message (API, GMT_TIME_NONE, "\t-C Use CPT (or specify -Ccolor1,color2[,color3,...]) to assign symbol\n");
 	GMT_Message (API, GMT_TIME_NONE, "\t   colors based on z-value in 3rd column.\n");
@@ -725,8 +726,8 @@ static int parse (struct GMT_CTRL *GMT, struct PSXY_CTRL *Ctrl, struct GMT_OPTIO
 			case 'A':	/* Turn off draw_arc mode */
 				Ctrl->A.active = true;
 				switch (opt->arg[0]) {
-					case 'm': case 'y': Ctrl->A.mode = GMT_STAIRS_Y; break;
-					case 'p': case 'x': Ctrl->A.mode = GMT_STAIRS_X; break;
+					case 'm': case 'y': case 'r': Ctrl->A.mode = GMT_STAIRS_Y; break;
+					case 'p': case 'x': case 't': Ctrl->A.mode = GMT_STAIRS_X; break;
 
 #ifdef DEBUG
 					default: Ctrl->A.step = atof (opt->arg); break; /* Undocumented test feature */

src/psxyz.c

@@ -36,7 +36,7 @@
 /* Control structure for psxyz */
 
 struct PSXYZ_CTRL {
-	struct PSXYZ_A {	/* -A[m|y|p|x|step] */
+	struct PSXYZ_A {	/* -A[m|y|p|x|r|t|step] */
 		bool active;
 		unsigned int mode;
 		double step;
@@ -159,7 +159,7 @@ static int usage (struct GMTAPI_CTRL *API, int level) {
 	const char *mod_name = &name[4];	/* To skip the leading gmt for usage messages */
 	if (level == GMT_MODULE_PURPOSE) return (GMT_NOERROR);
 	GMT_Message (API, GMT_TIME_NONE, "usage: %s [<table>] %s %s [%s]\n", name, GMT_J_OPT, GMT_Rgeoz_OPT, GMT_B_OPT);
-	GMT_Message (API, GMT_TIME_NONE, "\t[%s] [-A[m|p|x|y]] [-C<cpt>] [-D<dx>/<dy>[/<dz>]] [-G<fill>] [-H[<scale>]] [-I[<intens>]] %s\n\t[-L[+b|d|D][+xl|r|x0][+yb|t|y0][+p<pen>]] [-N[c|r]] %s\n", GMT_Jz_OPT, API->K_OPT, API->O_OPT);
+	GMT_Message (API, GMT_TIME_NONE, "\t[%s] [-A[m|p|x|y|r|t]] [-C<cpt>] [-D<dx>/<dy>[/<dz>]] [-G<fill>] [-H[<scale>]] [-I[<intens>]] %s\n\t[-L[+b|d|D][+xl|r|x0][+yb|t|y0][+p<pen>]] [-N[c|r]] %s\n", GMT_Jz_OPT, API->K_OPT, API->O_OPT);
 	if (API->GMT->current.setting.run_mode == GMT_CLASSIC)	/* -T has no purpose in modern mode */
 		GMT_Message (API, GMT_TIME_NONE, "\t%s[-Q] [-S[<symbol>][<size>][/size_y]] [-T]\n\t[%s] [%s] [-W[<pen>][<attr>]]\n", API->P_OPT, GMT_U_OPT, GMT_V_OPT);
 	else
@@ -176,6 +176,7 @@ static int usage (struct GMTAPI_CTRL *API, int level) {
 	GMT_Message (API, GMT_TIME_NONE, "\t   straight lines unless m or p is appended to first follow meridian\n");
 	GMT_Message (API, GMT_TIME_NONE, "\t   then parallel, or vice versa. Note: -A requires constant z-coordinates.\n");
 	GMT_Message (API, GMT_TIME_NONE, "\t   For Cartesian data, use -Ax or -Ay to draw x- or y-staircase curves.\n");
+	GMT_Message (API, GMT_TIME_NONE, "\t   For Polar projections, use -At or -Ar to draw theta- or r-staircase curves.\n");
 	GMT_Message (API, GMT_TIME_NONE, "\t-C Use CPT (or specify -Ccolor1,color2[,color3,...]) to assign symbol\n");
 	GMT_Message (API, GMT_TIME_NONE, "\t   colors based on t-value in 4rd column.\n");
 	GMT_Message (API, GMT_TIME_NONE, "\t   Note: requires -S.  Without -S, %s excepts lines/polygons\n", mod_name);
@@ -355,8 +356,8 @@ static int parse (struct GMT_CTRL *GMT, struct PSXYZ_CTRL *Ctrl, struct GMT_OPTI
 			case 'A':	/* Turn off draw_arc mode */
 				Ctrl->A.active = true;
 				switch (opt->arg[0]) {
-					case 'm': case 'y': Ctrl->A.mode = GMT_STAIRS_Y; break;
-					case 'p': case 'x': Ctrl->A.mode = GMT_STAIRS_X; break;
+					case 'm': case 'y': case 'r': Ctrl->A.mode = GMT_STAIRS_Y; break;
+					case 'p': case 'x': case 't': Ctrl->A.mode = GMT_STAIRS_X; break;
 
 #ifdef DEBUG
 					default: Ctrl->A.step = atof (opt->arg); break; /* Undocumented test feature */

test/psxy/polarcurves.sh

@@ -0,0 +1,20 @@
+#!/usr/bin/env bash
+# test the -Ar|t line-resampling for polar projections
+cat << EOF > t.txt
+10  5
+30 4
+50 7
+EOF
+gmt begin polarcurves ps
+	gmt subplot begin 3x1 -R0/60/0/10 -Fs8c -JP8c -A
+		gmt subplot set 0 -A"No -A"
+		gmt plot -W2p,red t.txt
+		gmt plot t.txt -Sc0.2c -Wthin
+		gmt subplot set 1 -A"-At: @~q@~, then r"
+		gmt plot -W2p,red -At t.txt
+		gmt plot t.txt -Sc0.2c -Wthin
+		gmt subplot set 2 -A"-Ar: r, then @~q@~"
+		gmt plot -W2p,red -Ar t.txt
+		gmt plot t.txt -Sc0.2c -Wthin
+	gmt subplot end
+gmt end show