diff --git a/doc/rst/source/supplements/geodesy/psvelo.rst b/doc/rst/source/supplements/geodesy/psvelo.rst
index fa1eede144b..443e45ff094 100644
--- a/doc/rst/source/supplements/geodesy/psvelo.rst
+++ b/doc/rst/source/supplements/geodesy/psvelo.rst
@@ -16,17 +16,19 @@ Synopsis
 |SYN_OPT-R|
 [ |-A|\ *parameters* ]
 [ |SYN_OPT-B| ]
+[ |-C|\ *cpt*]
 [ |-E|\ *fill* ]
 [ |-G|\ *fill* ]
 [ |-K| ]
-[ |-L| ]
+[ |-L|\ [*pen*\ [**+c**\ [**f**\|\ **l**]]] ]
 [ |-N| ] [ |-O| ] [ |-P| ]
 [ |-S|\ *<format><args>*\ [**+f**\ *font*] ]
 [ |SYN_OPT-U| ]
 [ |SYN_OPT-V| ]
-[ |-W|\ *pen* ]
+[ |-W|\ [*pen*][**+c**\ [**f**\|\ **l**]] ]
 [ |SYN_OPT-X| ]
 [ |SYN_OPT-Y| ]
+[ |-Z|\ [**m**\|\ **e**\|\ **n**\|\ **u**\ ][**+e**] ]
 [ |SYN_OPT-di| ]
 [ |SYN_OPT-e| ]
 [ |SYN_OPT-h| ]
diff --git a/doc/rst/source/supplements/geodesy/velo.rst b/doc/rst/source/supplements/geodesy/velo.rst
index bc4f5640608..f311bae58ca 100644
--- a/doc/rst/source/supplements/geodesy/velo.rst
+++ b/doc/rst/source/supplements/geodesy/velo.rst
@@ -16,15 +16,17 @@ Synopsis
 |SYN_OPT-R|
 [ |-A|\ *parameters* ]
 [ |SYN_OPT-B| ]
+[ |-C|\ *cpt*]
 [ |-E|\ *fill* ]
 [ |-G|\ *fill* ]
-[ |-L| ]
+[ |-L|\ [*pen*\ [**+c**\ [**f**\|\ **l**]]] ]
 [ |-N| ]
 [ |-S|\ *<format><args>*\ [**+f**\ *font*] ]
 [ |SYN_OPT-U| ]
 [ |SYN_OPT-V| ]
-[ |-W|\ *pen* ]
-[ |SYN_OPT-X| ]
+[ |-W|\ [*pen*][**+c**\ [**f**\|\ **l**]] ]
+[ |-Z|\ [**m**\|\ **e**\|\ **n**\|\ **u**\ ][**+e**] ]
+|SYN_OPT-X| ]
 [ |SYN_OPT-Y| ]
 [ |SYN_OPT-di| ]
 [ |SYN_OPT-e| ]
diff --git a/doc/rst/source/supplements/geodesy/velo_common.rst_ b/doc/rst/source/supplements/geodesy/velo_common.rst_
index ba1ccdf378e..70031bd1da9 100644
--- a/doc/rst/source/supplements/geodesy/velo_common.rst_
+++ b/doc/rst/source/supplements/geodesy/velo_common.rst_
@@ -4,8 +4,11 @@ Description
 -----------
 
 Reads data values from *files* [or standard input] and
-will plot velocity arrows on a map.
-Most options are the same as for :doc:`plot </plot>`, except **-S**.
+will plot the selected geodesy symbol on a map.
+You may choose from velocity vectors and their uncertainties,
+rotational wedges and their uncertainties, anisotropy bars,
+or strain crosses.  Symbol fills or their outlines may be colored
+based on constant parameters or via color lookup tables.
 
 Required Arguments
 ------------------
@@ -36,9 +39,9 @@ Required Arguments
         velocity arrows. The *confidence* sets the 2-dimensional confidence
         limit for the ellipse, e.g., 0.95 for 95% confidence ellipse. *font*
         sets the font and size of the text [9p,Helvetica,black]. The ellipse will be filled with the
-        color or shade specified by the |-G| option [default transparent]. The
+        color or shade specified by the |-E| option [default is transparent]. The
         arrow and the circumference of the ellipse will be drawn with the pen
-        attributes specified by the |-W| option. Parameters are expected to be
+        attributes specified by the |-W| option and arrow-head can be colored via |-G|. Parameters are expected to be
         in the following columns:
 
             **1**,\ **2**:
@@ -49,7 +52,7 @@ Required Arguments
             uncertainty of eastward, northward velocities (1-sigma) (**-:** option interchanges order)
             **7**:
             correlation between eastward and northward components
-            **8**:
+            **Trailing text**:
             name of station (optional).
 
     **-Sn**\ *barscale*
@@ -68,9 +71,9 @@ Required Arguments
         the velocity arrows. The *confidence* sets the 2-dimensional
         confidence limit for the ellipse, e.g., 0.95 for 95% confidence ellipse.
         *font* sets the font and size of the text [9p,Helvetica,black]. The ellipse will be
-        filled with the color or shade specified by the |-G| option [default
+        filled with the color or shade specified by the |-E| option [default
         transparent]. The arrow and the circumference of the ellipse will be
-        drawn with the pen attributes specified by the |-W| option. Parameters
+        drawn with the pen attributes specified by the |-W| option and arrow-head can be colored via |-G|. Parameters
         are expected to be in the following columns:
 
             **1**,\ **2**:
@@ -81,16 +84,16 @@ Required Arguments
             semi-major, semi-minor axes
             **7**:
             counter-clockwise angle, in degrees, from horizontal axis to major axis of ellipse.
-            **8**:
+            **Trailing text**:
             name of station (optional)
 
     **-Sw**\ *wedgescale/wedgemag*
 
-        Rotational wedges. *wedgescale* sets the size of the wedges. Values are
+        Rotational wedges. *wedgescale* sets the size of the wedges. Rotation values are
         multiplied by *wedgemag* before plotting. For example, setting
         *wedgemag* to 1.e7 works well for rotations of the order of 100
-        nanoradians/yr. Use **-G** to set the fill color or shade for the wedge,
-        and **-E** to set the color or shade for the uncertainty. Parameters are
+        nanoradians/yr. Use |-G| to set the fill color or shade for the wedge,
+        and |-E| to set the color or shade for the uncertainty. Parameters are
         expected to be in the following columns:
 
             **1**,\ **2**:
@@ -128,6 +131,12 @@ Optional Arguments
 
 .. include:: ../../explain_-B.rst_
 
+.. _-C:
+
+**-C**\ *cpt*
+    Give a CPT and let symbol color normally set by |-G| be
+    determined from the magnitude.  See **-Z** for other selections.
+
 .. _-D:
 
 **-D**\ *Sigma\_scale*
@@ -138,19 +147,28 @@ Optional Arguments
 
 **-E**\ *fill* :ref:`(more ...) <-Gfill_attrib>`
     Sets the color or shade used for filling uncertainty wedges
-    (**-Sw**) or velocity error ellipses (**-Se** or **-Sr**). [If
-    **-E** is not specified, the uncertainty regions will be transparent.]
+    (**-Sw**) or velocity error ellipses (**-Se** or **-Sr**). If
+    **-E** is not specified, the uncertainty regions will be transparent.
+    **Note**: Using **-C** and **-Z+e** will update the uncertainty fill
+    color based on the selected measure in **-Z** [magnitude error].
 
 .. _-G:
 
 **-G**\ *fill* :ref:`(more ...) <-Gfill_attrib>`
-    Select color or pattern for filling of symbols or polygons [Default is no fill].
+    Select color or pattern for filling of symbols [Default is no fill].
+    **Note**: Using **-C** (and optionally **-Z**) will update the symbol fill
+    color based on the selected measure in **-Z** [magnitude].
 
 .. _-L:
 
-**-L**
-    Draw lines. Ellipses and fault planes will have their outlines drawn
-    using current pen (see **-W**).
+**-L**\ [*pen*\ [**+c**\ [**f**\|\ **l**]]]
+    Draw lines. Ellipses and rotational wedges will have their outlines drawn
+    using current pen (see **-W**).  Alternatively, append a separate pen to
+    use for the error outlines.
+    If the modifier **+cl** is appended then the color of the pen are updated from the CPT (see
+    **-C**). If instead modifier **+cf** is appended then the color from the cpt
+    file is applied to error fill only [Default].  Use just **+c** to set both
+    pen and fill color.
 
 .. _-N:
 
@@ -169,14 +187,29 @@ Optional Arguments
 
 .. _-W:
 
-**-W**
+**-W**\ [*pen*][**+c**\ [**f**\|\ **l**]]
     Set pen attributes for velocity arrows, ellipse circumference and
     fault plane edges. [Defaults: width = default, color = black, style = solid].
+    If the modifier **+cl** is appended then the color of the pen are updated from the CPT (see
+    **-C**). If instead modifier **+cf** is appended then the color from the cpt
+    file is applied to symbol fill only [Default].  Use just **+c** to set both
+    pen and fill color.
 
 .. _-X:
 
 .. include:: ../../explain_-XY.rst_
 
+.. _-Z:
+
+**Z**\ [**m**\|\ **e**\|\ **n**\|\ **u**\ ][**+e**]
+    Select the quantity that will be used with the CPT given via **-C** to
+    set the fill color.  Choose from **m**\ agnitude (vector magnitude or
+    rotation magnitude), **e**\ ast-west velocity, **n**\ orth-south velocity,
+    or **u**\ ser-supplied data column (supplied after the required columns).
+    To instead use the corresponding error estimates (i.e., vector or rotation
+    uncertainty) to lookup the color and paint the error ellipse or wedge
+    instead, append **+e**.
+
 .. |Add_-di| unicode:: 0x20 .. just an invisible code
 .. include:: ../../explain_-di.rst_
 
diff --git a/src/geodesy/psvelo.c b/src/geodesy/psvelo.c
index b4d901e1a47..e3398562a3d 100644
--- a/src/geodesy/psvelo.c
+++ b/src/geodesy/psvelo.c
@@ -1,35 +1,41 @@
 /*--------------------------------------------------------------------
  *
- *    Copyright (c) 1996-2012 by G. Patau
- *    Copyright (c) 2013-2021 by the GMT Team (https://www.generic-mapping-tools.org/team.html)
- *    Donated to the GMT project by G. Patau upon her retirement from IGPG
- *    Distributed under the Lesser GNU Public Licence
- *    See README file for copying and redistribution conditions.
+ *	Copyright (c) 1991-2021 by the GMT Team (https://www.generic-mapping-tools.org/team.html)
+ *	See LICENSE.TXT file for copying and redistribution conditions.
+ *
+ *	This program is free software; you can redistribute it and/or modify
+ *	it under the terms of the GNU Lesser General Public License as published by
+ *	the Free Software Foundation; version 3 or any later version.
+ *
+ *	This program is distributed in the hope that it will be useful,
+ *	but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *	GNU Lesser General Public License for more details.
+ *
+ *	Contact info: www.generic-mapping-tools.org
  *--------------------------------------------------------------------*/
 /*
-
-psvelo will read <x,y> pairs (or <lon,lat>) from inputfile and
-plot symbols on a map. Velocity ellipses, strain
-crosses, or strain wedges, may be specified, some of which require
-additional columns of data.  Only one symbol may be plotted at a time.
-PostScript code is written to stdout.
-
-
- Author:	Kurt Feigl
- Date:		7 July 1998
- Version:	5
- Roots:		based on psxy.c
- Adapted to version 3.3 by Genevieve Patau (25 June 1999)
- Last modified : 18 February 2000.  Ported to GMT 5 by P. Wessel
-
-*/
+ * Author:	G. Patau, IGPG, w/ Kurt Feigl
+ * Date:	1-JAN-2010
+ * Version:	6 API
+ *    Copyright (c) 1996-2012 by G. Patau, then donated to the GMT project
+ *    by G. Patau upon her retirement from IGPG
+ *
+ * Roots:		based on psxy.c
+ * Adapted to version 3.3 by Genevieve Patau (25 June 1999)
+ * Last modified : 18 February 2000.  Ported to GMT 5 by P. Wessel in 2013.
+ * Updated 3 March 2021 to add color enhancements by P. Wessel
+ *
+ * Brief synopsis: psvelo will read coordinates and plot geodetic symbols
+ * such as velocity ellipses, strain crosses, or strain wedges on maps.
+ */
 
 #include "gmt_dev.h"
 
 #define THIS_MODULE_CLASSIC_NAME	"psvelo"
 #define THIS_MODULE_MODERN_NAME	"velo"
 #define THIS_MODULE_LIB		"geodesy"
-#define THIS_MODULE_PURPOSE	"Plot velocity vectors, crosses, and wedges"
+#define THIS_MODULE_PURPOSE	"Plot velocity vectors, crosses, anisotropy bars and wedges"
 #define THIS_MODULE_KEYS	"<D{,>X}"
 #define THIS_MODULE_NEEDS	"Jd"
 #define THIS_MODULE_OPTIONS "-:>BHJKOPRUVXYdehipqt" GMT_OPT("c")
@@ -58,6 +64,10 @@ struct PSVELO_CTRL {
 		bool active;
 		struct GMT_SYMBOL S;
 	} A;
+	struct PSVELO_C {	/* -C<cpt> */
+		bool active;
+		char *file;
+	} C;
 	struct PSVELO_D {	/* -D */
 		bool active;
 		double scale;
@@ -70,8 +80,10 @@ struct PSVELO_CTRL {
 		bool active;
 		struct GMT_FILL fill;
 	} G;
-	struct PSVELO_L {	/* -L */
+	struct PSVELO_L {	/* -L[<pen>] */
 		bool active;
+		bool error_pen;
+		struct GMT_PEN pen;
 	} L;
 	struct PSVELO_N {	/* -N */
 		bool active;
@@ -90,9 +102,24 @@ struct PSVELO_CTRL {
 		bool active;
 		struct GMT_PEN pen;
 	} W;
+	struct PSVELO_Z {	/* -Z */
+		bool active;
+		unsigned int mode;
+		unsigned int item;
+	} Z;
+};
+
+enum psvelo_types {
+	PSVELO_G_FILL = 0,
+	PSVELO_E_FILL = 1,
+	PSVELO_V_MAG	= 0,
+	PSVELO_V_EAST,
+	PSVELO_V_NORTH,
+	PSVELO_R_MAG,
+	PSVELO_V_USER
 };
 
-/* COntent of old utilvelo.c is here */
+/* Content of old utilvelo.c is here */
 
 #define squared(x) ((x) * (x))
 #define EPSIL 0.0001
@@ -622,10 +649,10 @@ 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 [-A<vecpar>] [%s] [-D<sigscale>]\n", name, GMT_J_OPT, GMT_Rgeo_OPT, GMT_B_OPT);
-	GMT_Message (API, GMT_TIME_NONE, "\t[-G<fill>] %s[-L] [-N] %s%s[-S<symbol><args>[+f<font>]]\n", API->K_OPT, API->O_OPT, API->P_OPT);
-	GMT_Message (API, GMT_TIME_NONE, "\t[%s] [-V] [-W<pen>] [%s]\n", GMT_U_OPT, GMT_X_OPT);
-	GMT_Message (API, GMT_TIME_NONE, "\t[%s] %s[%s] [%s] [%s]\n\t[%s] [%s]\n\t[%s] [%s] [%s] [%s]\n\n", GMT_Y_OPT, API->c_OPT, GMT_di_OPT, GMT_e_OPT, GMT_h_OPT, GMT_i_OPT, GMT_p_OPT, GMT_qi_OPT, GMT_t_OPT, GMT_colon_OPT, GMT_PAR_OPT);
+	GMT_Message (API, GMT_TIME_NONE, "usage: %s [<table>] %s %s [-A<vecpar>] [%s]\n", name, GMT_J_OPT, GMT_Rgeo_OPT, GMT_B_OPT);
+	GMT_Message (API, GMT_TIME_NONE, "\t[-C<cpt>] [-D<sigscale>] [-G<fill>] %s[-L[<pen>][+c[f|l]]] [-N] %s%s[-S<symbol><args>[+f<font>]]\n", API->K_OPT, API->O_OPT, API->P_OPT);
+	GMT_Message (API, GMT_TIME_NONE, "\t[%s] [-V] [-W[<pen>][+c[f|l]]] [%s] [%s]\n", GMT_U_OPT, GMT_X_OPT, GMT_Y_OPT);
+	GMT_Message (API, GMT_TIME_NONE, "\t[-Z[m|e|n|u][+e] %s[%s] [%s] [%s]\n\t[%s] [%s]\n\t[%s] [%s] [%s] [%s]\n\n", API->c_OPT, GMT_di_OPT, GMT_e_OPT, GMT_h_OPT, GMT_i_OPT, GMT_p_OPT, GMT_qi_OPT, GMT_t_OPT, GMT_colon_OPT, GMT_PAR_OPT);
 
 	if (level == GMT_SYNOPSIS) return (GMT_MODULE_SYNOPSIS);
 
@@ -635,13 +662,15 @@ static int usage (struct GMTAPI_CTRL *API, int level) {
 	GMT_Message (API, GMT_TIME_NONE, "\t-A Specify arrow head attributes:\n");
 	gmt_vector_syntax (API->GMT, 15);
 	GMT_Message (API, GMT_TIME_NONE, "\t   Default is %gp+gblack+p1p\n", VECTOR_HEAD_LENGTH);
-	GMT_Message (API, GMT_TIME_NONE, "\t-D Multiply uncertainties by <sigscale>. (Se and Sw only)i\n");
-	GMT_Message (API, GMT_TIME_NONE, "\t-E Set color used for uncertainty wedges in -Sw option.\n");
-	GMT_Message (API, GMT_TIME_NONE, "\t-G Specify color (for symbols/polygons) or pattern (for polygons). fill can be either\n");
-	GMT_Message (API, GMT_TIME_NONE, "\t   1) <r/g/b> (each 0-255) for color or <gray> (0-255) for gray-shade [0].\n");
-	GMT_Message (API, GMT_TIME_NONE, "\t   2) p[or P]<iconsize>/<pattern> for predefined patterns (0-90).\n");
+	GMT_Message (API, GMT_TIME_NONE, "\t-C Use CPT to assign colors based on vector magnitude.\n");
+	GMT_Message (API, GMT_TIME_NONE, "\t   For other coloring options, see -W and -Z.\n");
+	GMT_Message (API, GMT_TIME_NONE, "\t-D Multiply uncertainties by <sigscale> (Se and Sw only).\n");
+	GMT_Message (API, GMT_TIME_NONE, "\t-E Set color used for uncertainty ellipses and wedges [no fill].\n");
+	GMT_Message (API, GMT_TIME_NONE, "\t   For other coloring options, see -L and -Z.\n");
+	gmt_fill_syntax (API->GMT, 'G', NULL, "Specify color or pattern for symbol fill [no fill].");
 	GMT_Option (API, "K");
 	GMT_Message (API, GMT_TIME_NONE, "\t-L Draw line or symbol outline using the current pen (see -W).\n");
+	GMT_Message (API, GMT_TIME_NONE, "\t   Optionally, append separate pen for error outlines [Same as -W].\n");
 	GMT_Message (API, GMT_TIME_NONE, "\t-N Do Not skip/clip symbols that fall outside map border [Default will ignore those outside].\n");
 	GMT_Option (API, "O,P");
 	GMT_Message (API, GMT_TIME_NONE, "\t-S Select symbol type and scale (plus optional font; see documentation). Choose between:\n");
@@ -652,7 +681,14 @@ static int usage (struct GMTAPI_CTRL *API, int level) {
 	GMT_Message (API, GMT_TIME_NONE, "\t     x  Strain crosses : in X,Y,Eps1,Eps2,Theta.\n");
 	GMT_Option (API, "U,V");
 	GMT_Message (API, GMT_TIME_NONE, "\t-W Set pen attributes [%s].\n", gmt_putpen (API->GMT, &API->GMT->current.setting.map_default_pen));
-	GMT_Option (API, "X,c,di,e,h,i,p,qi,t,:,.");
+	GMT_Option (API, "X");
+	GMT_Message (API, GMT_TIME_NONE, "\t-Z Select quantity to use with -C to look-up colors.  Choose among:\n");
+	GMT_Message (API, GMT_TIME_NONE, "\t   m: Magnitude of vector or rotation [Default].\n");
+	GMT_Message (API, GMT_TIME_NONE, "\t   e: East velocity component.\n");
+	GMT_Message (API, GMT_TIME_NONE, "\t   n: North velocity component.\n");
+	GMT_Message (API, GMT_TIME_NONE, "\t   u: User column (given after required columns).\n");
+	GMT_Message (API, GMT_TIME_NONE, "\t   Color selected will replace -G<fill>.  Append +e to instead act as -E<fill>.\n");
+	GMT_Option (API, "c,di,e,h,i,p,qi,t,:,.");
 
 	return (GMT_MODULE_USAGE);
 }
@@ -707,6 +743,10 @@ static int parse (struct GMT_CTRL *GMT, struct PSVELO_CTRL *Ctrl, struct GMT_OPT
 					Ctrl->A.S.symbol = PSL_VECTOR;
 				}
 				break;
+			case 'C':	/* Select CPT for coloring */
+				Ctrl->C.active = true;
+				if (opt->arg[0]) Ctrl->C.file = strdup (opt->arg);
+				break;
 			case 'D':	/* Rescale Sigmas */
 				Ctrl->D.active = true;
 				sscanf (opt->arg, "%lf",&Ctrl->D.scale);
@@ -727,6 +767,13 @@ static int parse (struct GMT_CTRL *GMT, struct PSVELO_CTRL *Ctrl, struct GMT_OPT
 				break;
 			case 'L':	/* Draw the outline */
 				Ctrl->L.active = true;
+				if (opt->arg[0]) {
+					Ctrl->L.error_pen = true;
+					if (gmt_getpen (GMT, opt->arg, &Ctrl->L.pen)) {
+						gmt_pen_syntax (GMT, 'L', NULL, " ", 0);
+						n_errors++;
+					}
+				}
 				break;
 			case 'N':	/* Do not skip points outside border */
 				Ctrl->N.active = true;
@@ -739,21 +786,23 @@ static int parse (struct GMT_CTRL *GMT, struct PSVELO_CTRL *Ctrl, struct GMT_OPT
 				}
  				if (opt->arg[0] == 'e' || opt->arg[0] == 'r') {
 					strncpy (txt, &opt->arg[1], GMT_LEN256);
-					n = 0; while (txt[n] && txt[n] != '/') n++; txt[n] = 0;
-					Ctrl->S.scale = gmt_M_to_inch (GMT, txt);
-					sscanf (strchr(&opt->arg[1],'/')+1, "%lf/%s", &Ctrl->S.confidence, txt_b);
+					n = 0; while (txt[n] && txt[n] != '/') n++; txt[n] = 0;	/* Hide the /confidence part */
+					Ctrl->S.scale = gmt_M_to_inch (GMT, txt);	/* Get symbol size */
+					sscanf (strchr (&opt->arg[1],'/')+1, "%lf/%s", &Ctrl->S.confidence, txt_b);
 					/* confidence scaling */
 					Ctrl->S.conrad = sqrt (-2.0 * log (1.0 - Ctrl->S.confidence));
+					/* Check for deprecated font syntax */
 					if (txt_b[0]) Ctrl->S.font.size = gmt_convert_units (GMT, txt_b, GMT_PT, GMT_PT);
 				}
-				if (opt->arg[0] == 'n' || opt->arg[0] == 'x' ) Ctrl->S.scale = gmt_M_to_inch (GMT, &opt->arg[1]);
+				if (opt->arg[0] == 'n' || opt->arg[0] == 'x')	/* Simple one-parameter argument */
+					Ctrl->S.scale = gmt_M_to_inch (GMT, &opt->arg[1]);
 				if (opt->arg[0] == 'w' && strlen(opt->arg) > 3) {
-					strncpy(txt, &opt->arg[1], GMT_LEN256);
-					n=0; while (txt[n] && txt[n] != '/') n++; txt[n]=0;
-					Ctrl->S.scale = gmt_M_to_inch (GMT, txt);
-					sscanf(strchr(&opt->arg[1],'/')+1, "%lf", &Ctrl->S.wedge_amp);
+					strncpy (txt, &opt->arg[1], GMT_LEN256);
+					n = 0; while (txt[n] && txt[n] != '/') n++; txt[n] = 0;	/* Hide the /wedgemag part */
+					Ctrl->S.scale = gmt_M_to_inch (GMT, txt);	/* Get symbol size */
+					sscanf (strchr (&opt->arg[1],'/')+1, "%lf", &Ctrl->S.wedge_amp);
 				}
-				switch (opt->arg[0]) {
+				switch (opt->arg[0]) {	/* Set modes and expected input columns */
 					case 'e':
 						Ctrl->S.symbol = CINE;	Ctrl->S.n_cols = 7;
 						Ctrl->S.readmode = READ_ELLIPSE;
@@ -782,17 +831,38 @@ static int parse (struct GMT_CTRL *GMT, struct PSVELO_CTRL *Ctrl, struct GMT_OPT
 				break;
 			case 'W':	/* Set line attributes */
 				Ctrl->W.active = true;
-				if (opt->arg && gmt_getpen (GMT, opt->arg, &Ctrl->W.pen)) {
+				if (opt->arg[0] && gmt_getpen (GMT, opt->arg, &Ctrl->W.pen)) {
 					gmt_pen_syntax (GMT, 'W', NULL, " ", 0);
 					n_errors++;
 				}
 				break;
+			case 'Z':	/* Set items to control CPT coloring */
+				Ctrl->Z.active = true;
+				if (opt->arg[0] && (c = strstr (opt->arg, "+e"))) {	/* Paint error part of symbol instead (-E) */
+					Ctrl->Z.item = PSVELO_E_FILL;
+					c[0] = '\0';	/* Temporarily chop off the modifier */
+				}
+				switch (opt->arg[0]) {
+					case 'm':	case '\0': Ctrl->Z.mode = PSVELO_V_MAG;	break;
+					case 'e':	Ctrl->Z.mode = PSVELO_V_EAST;	break;
+					case 'n':	Ctrl->Z.mode = PSVELO_V_NORTH;	break;
+					case 'r':	Ctrl->Z.mode = PSVELO_R_MAG;	break;
+					case 'u':	Ctrl->Z.mode = PSVELO_V_USER;	break;
+					default:
+						GMT_Report (GMT->parent, GMT_MSG_ERROR, "Option -Z: Unrecognized mode %s\n", opt->arg[0]);
+						n_errors++;
+						break;
+				}
+				if (c) c[0] = '+';	/* Restore modifier */
+				break;
 
 			/* Illegal options */
 
 		}
 	}
 
+	gmt_consider_current_cpt (GMT->parent, &Ctrl->C.active, &(Ctrl->C.file));
+
 	no_size_needed = (Ctrl->S.readmode == READ_ELLIPSE || Ctrl->S.readmode == READ_ROTELLIPSE || Ctrl->S.readmode == READ_ANISOTROPY || Ctrl->S.readmode == READ_CROSS || Ctrl->S.readmode == READ_WEDGE );
         /* Only one allowed */
 	n_set = (Ctrl->S.readmode == READ_ELLIPSE) + (Ctrl->S.readmode == READ_ROTELLIPSE) + (Ctrl->S.readmode == READ_ANISOTROPY) + (Ctrl->S.readmode == READ_CROSS) + (Ctrl->S.readmode == READ_WEDGE);
@@ -800,6 +870,9 @@ static int parse (struct GMT_CTRL *GMT, struct PSVELO_CTRL *Ctrl, struct GMT_OPT
 	n_errors += gmt_M_check_condition (GMT, n_set > 1, "Only one -S setting is allowed.\n");
 	n_errors += gmt_M_check_condition (GMT, !no_size_needed && (Ctrl->S.symbol > 1 && Ctrl->S.scale <= 0.0), "Option -S: Must specify symbol size.\n");
 	n_errors += gmt_M_check_condition (GMT, Ctrl->D.active && ! (Ctrl->S.readmode == READ_ELLIPSE || Ctrl->S.readmode == READ_WEDGE), "Option -D requires -Se|w.\n");
+	n_errors += gmt_M_check_condition (GMT, Ctrl->Z.active && !Ctrl->C.active, "Option -Z requires -C.\n");
+	n_errors += gmt_M_check_condition (GMT, Ctrl->C.active && Ctrl->Z.item == PSVELO_E_FILL && Ctrl->E.active, "Options -C -Z+e cannot be combined with -E.\n");
+	n_errors += gmt_M_check_condition (GMT, Ctrl->C.active && Ctrl->Z.item == PSVELO_G_FILL && Ctrl->G.active, "Options -C -Z cannot be combined with -G.\n");
 
 	if (!got_A && Ctrl->W.active) Ctrl->A.S.v.pen = Ctrl->W.pen;	/* Set vector pen to that given by -W  */
 	if (Ctrl->A.S.v.status & PSL_VEC_OUTLINE2 && Ctrl->W.active) gmt_M_rgb_copy (Ctrl->A.S.v.pen.rgb, Ctrl->W.pen.rgb);	/* Set vector pen color from -W but not thickness */
@@ -809,18 +882,37 @@ static int parse (struct GMT_CTRL *GMT, struct PSVELO_CTRL *Ctrl, struct GMT_OPT
 #define bailout(code) {gmt_M_free_options (mode); return (code);}
 #define Return(code) {Free_Ctrl (GMT, Ctrl); gmt_end_module (GMT, GMT_cpy); bailout (code);}
 
+GMT_LOCAL void psvelo_set_colorfill (struct GMT_CTRL *GMT, struct PSVELO_CTRL *Ctrl, struct GMT_PALETTE *P, double value) {
+	/* Called if -C was given.  Selects and updates color fills and possibly pen colors */
+	struct GMT_FILL *F = (Ctrl->Z.item == PSVELO_G_FILL) ? &Ctrl->G.fill : &Ctrl->E.fill;
+	gmt_get_fill_from_z (GMT, P, value, F);
+	if (Ctrl->L.pen.cptmode & 1) {	/* Also change error pen color via CPT */
+		gmt_M_rgb_copy (Ctrl->L.pen.rgb, F->rgb);
+		gmt_setpen (GMT, &Ctrl->L.pen);
+	}
+	if (Ctrl->W.pen.cptmode & 1) {	/* Also change pen color via CPT */
+		gmt_M_rgb_copy (Ctrl->W.pen.rgb, F->rgb);
+		gmt_setpen (GMT, &Ctrl->W.pen);
+		if (!Ctrl->L.error_pen)	/* No -L pen so duplicate the change */
+			gmt_M_rgb_copy (Ctrl->L.pen.rgb, Ctrl->W.pen.rgb);
+	}
+}
+
 EXTERN_MSC int GMT_psvelo (void *V_API, int mode, void *args) {
 	int ix = 0, iy = 1, n_rec = 0, justify;
-	int des_ellipse = true, des_arrow = true, error = false;
+	int plot_ellipse = true, plot_vector = true, error = false;
+	bool set_g_fill, set_e_fill;
 
 	double plot_x, plot_y, vxy[2], plot_vx, plot_vy, length, s, dim[PSL_MAX_DIMS];
 	double eps1 = 0.0, eps2 = 0.0, spin = 0.0, spinsig = 0.0, theta = 0.0, *in = NULL;
 	double direction = 0, small_axis = 0, great_axis = 0, sigma_x, sigma_y, corr_xy;
 	double t11 = 1.0, t12 = 0.0, t21 = 0.0, t22 = 1.0, hl, hw, vw, ssize, headpen_width = 0.0;
+	double z_val, e_val, value;
 
 	char *station_name = NULL;
 
 	struct GMT_RECORD *In = NULL;
+	struct GMT_PALETTE *CPT = NULL;
 	struct PSVELO_CTRL *Ctrl = NULL;
 	struct GMT_CTRL *GMT = NULL, *GMT_cpy = NULL;		/* General GMT internal parameters */
 	struct GMT_OPTION *options = NULL;
@@ -844,6 +936,17 @@ EXTERN_MSC int GMT_psvelo (void *V_API, int mode, void *args) {
 
 	/*---------------------------- This is the psvelo main code ----------------------------*/
 
+	set_e_fill = Ctrl->E.active;	set_g_fill = Ctrl->G.active;
+	if (Ctrl->C.active) {
+		if ((CPT = GMT_Read_Data (API, GMT_IS_PALETTE, GMT_IS_FILE, GMT_IS_NONE, GMT_READ_NORMAL, NULL, Ctrl->C.file, NULL)) == NULL) {
+			Return (API->error);
+		}
+		if (Ctrl->Z.item == PSVELO_G_FILL) set_g_fill = true;	/* Since we will set it via CPT lookup */
+		if (Ctrl->Z.item == PSVELO_E_FILL) set_e_fill = true;	/* Since we will set it via CPT lookup */
+	}
+	if (!Ctrl->L.error_pen)	/* Duplicate -W to -L */
+		gmt_M_memcpy (&Ctrl->L.pen, &Ctrl->W.pen, 1, struct GMT_PEN);
+
 	if (gmt_map_setup (GMT, GMT->common.R.wesn)) Return (GMT_PROJECTION_ERROR);
 
 	if ((PSL = gmt_plotinit (GMT, options)) == NULL) Return (GMT_RUNTIME_ERROR);
@@ -863,6 +966,8 @@ EXTERN_MSC int GMT_psvelo (void *V_API, int mode, void *args) {
 
 	ix = (GMT->current.setting.io_lonlat_toggle[0]);	iy = 1 - ix;
 
+	if (Ctrl->Z.mode == PSVELO_V_USER) Ctrl->S.n_cols++;	/* Need to read one extra column */
+
 	GMT_Set_Columns (API, GMT_IN, Ctrl->S.n_cols, GMT_COL_FIX);
 
 	if (GMT_Init_IO (API, GMT_IS_DATASET, GMT_IS_POINT, GMT_IN, GMT_ADD_DEFAULT, 0, options) != GMT_NOERROR) {	/* Register data input */
@@ -874,6 +979,8 @@ EXTERN_MSC int GMT_psvelo (void *V_API, int mode, void *args) {
 
 	if (Ctrl->S.readmode == READ_ELLIPSE || Ctrl->S.readmode == READ_ROTELLIPSE) GMT_Report (API, GMT_MSG_INFORMATION, "psvelo: 2-D confidence interval and scaling factor %f %f\n", Ctrl->S.confidence, Ctrl->S.conrad);
 
+	if (Ctrl->D.active)  GMT_Report (API, GMT_MSG_INFORMATION, "Rescaling uncertainties by a factor of %f\n", Ctrl->D.scale);
+
 	if (Ctrl->S.symbol == CINE || Ctrl->S.symbol == CROSS) {
 		if (Ctrl->A.S.v.status & PSL_VEC_OUTLINE2) {	/* Vector head outline pen specified separately */
 			PSL_defpen (PSL, "PSL_vecheadpen", Ctrl->A.S.v.pen.width, Ctrl->A.S.v.pen.style, Ctrl->A.S.v.pen.offset, Ctrl->A.S.v.pen.rgb);
@@ -914,15 +1021,23 @@ EXTERN_MSC int GMT_psvelo (void *V_API, int mode, void *args) {
 			sigma_x = in[4];
 			sigma_y = in[5];
 			corr_xy = in[6];
+			if (Ctrl->C.active) {
+				switch (Ctrl->Z.mode) {
+					case PSVELO_V_MAG:   z_val = hypot (vxy[0], vxy[1]);	e_val = hypot (sigma_x, sigma_y); break;
+					case PSVELO_V_EAST:  z_val = vxy[0]; e_val = sigma_x;	break;
+					case PSVELO_V_NORTH: z_val = vxy[1]; e_val = sigma_y;	break;
+					case PSVELO_V_USER:  z_val = e_val = in[7];	break;
+				}
+			}
 			/* rescale uncertainties if necessary */
 			if (Ctrl->D.active) {
 				sigma_x = Ctrl->D.scale * sigma_x;
 				sigma_y = Ctrl->D.scale * sigma_y;
 			}
 			if (fabs (sigma_x) < EPSIL && fabs (sigma_y) < EPSIL)
-				des_ellipse = false;
+				plot_ellipse = false;
 			else {
-				des_ellipse = true;
+				plot_ellipse = true;
 				psvelo_ellipse_convert (sigma_x, sigma_y, corr_xy, Ctrl->S.conrad, &small_axis, &great_axis, &direction);
 
 				/* convert to degrees */
@@ -935,10 +1050,18 @@ EXTERN_MSC int GMT_psvelo (void *V_API, int mode, void *args) {
 			great_axis = Ctrl->S.conrad*in[4];
 			small_axis = Ctrl->S.conrad*in[5];
 			direction = in[6];
+			if (Ctrl->C.active) {
+				switch (Ctrl->Z.mode) {
+					case PSVELO_V_MAG:   z_val = hypot (vxy[0], vxy[1]);	e_val = hypot (great_axis, small_axis); break;
+					case PSVELO_V_EAST:  z_val = vxy[0]; e_val = great_axis;	break;
+					case PSVELO_V_NORTH: z_val = vxy[1]; e_val = small_axis;	break;
+					case PSVELO_V_USER:  z_val = e_val = in[7];	break;
+				}
+			}
 			if (fabs (great_axis) < EPSIL && fabs (small_axis) < EPSIL)
-				des_ellipse = false;
+				plot_ellipse = false;
 			else
-				des_ellipse = true;
+				plot_ellipse = true;
 		}
 		else if (Ctrl->S.readmode == READ_ANISOTROPY) {
 			vxy[ix] = in[2];
@@ -952,6 +1075,12 @@ EXTERN_MSC int GMT_psvelo (void *V_API, int mode, void *args) {
 		else if (Ctrl->S.readmode == READ_WEDGE) {
 			spin    = in[2];
 			spinsig = in[3];
+			if (Ctrl->C.active) {
+				switch (Ctrl->Z.mode) {
+					case PSVELO_V_MAG:  z_val = spin;	e_val = spinsig; break;
+					case PSVELO_V_USER: z_val = e_val = in[4];	break;
+				}
+			}
 			if (Ctrl->D.active) spinsig = spinsig * Ctrl->D.scale;
 		}
 
@@ -962,12 +1091,19 @@ EXTERN_MSC int GMT_psvelo (void *V_API, int mode, void *args) {
 
 		gmt_geo_to_xy (GMT, in[GMT_X], in[GMT_Y], &plot_x, &plot_y);
 
+		value = (Ctrl->Z.item == PSVELO_E_FILL) ? e_val : z_val;	/* Select which value for color lookup - if active */
+		if (Ctrl->C.active) {	/* Possibly update E or G fills based on value, then set in PS */
+			psvelo_set_colorfill (GMT, Ctrl, CPT, value);
+			gmt_init_vector_param (GMT, &Ctrl->A.S, true, Ctrl->W.active, &Ctrl->W.pen, Ctrl->G.active, &Ctrl->G.fill);
+		}
+
 		switch (Ctrl->S.symbol) {
 			case CINE:
-				des_arrow = hypot (vxy[0], vxy[1]) < 1.e-8 ? false : true;
+				plot_vector = (hypot (vxy[0], vxy[1]) < 1.e-8) ? false : true;
 				psvelo_trace_arrow (GMT, in[GMT_X], in[GMT_Y], vxy[0], vxy[1], Ctrl->S.scale, &plot_x, &plot_y, &plot_vx, &plot_vy);
 				psvelo_get_trans (GMT, in[GMT_X], in[GMT_Y], &t11, &t12, &t21, &t22);
-				if (des_ellipse) {
+				if (plot_ellipse) {
+					if (Ctrl->L.active) gmt_setpen (GMT, &Ctrl->L.pen);
 					if (Ctrl->E.active)
 						psvelo_paint_ellipse (GMT, plot_vx, plot_vy, direction, great_axis, small_axis, Ctrl->S.scale,
 							t11,t12,t21,t22, Ctrl->E.active, &Ctrl->E.fill, Ctrl->L.active);
@@ -975,7 +1111,7 @@ EXTERN_MSC int GMT_psvelo (void *V_API, int mode, void *args) {
 						psvelo_paint_ellipse (GMT, plot_vx, plot_vy, direction, great_axis, small_axis, Ctrl->S.scale,
 							t11,t12,t21,t22, Ctrl->E.active, &Ctrl->G.fill, Ctrl->L.active);
 				}
-				if (des_arrow) {	/* verify that arrow is not ridiculously small */
+				if (plot_vector) {	/* verify that vector length is not ridiculously small */
 					length = hypot (plot_x-plot_vx, plot_y-plot_vy);	/* Length of arrow */
 					if (length < Ctrl->A.S.v.h_length && Ctrl->A.S.v.v_norm < 0.0)	/* No shrink requested yet head length exceeds total vector length */
 						GMT_Report (API, GMT_MSG_INFORMATION, "Vector head length exceeds overall vector length near line %d. Consider adding +n<norm> to -A\n", n_rec);
@@ -988,13 +1124,9 @@ EXTERN_MSC int GMT_psvelo (void *V_API, int mode, void *args) {
 					dim[0] = plot_vx, dim[1] = plot_vy;
 					dim[2] = vw, dim[3] = hl, dim[4] = hw;
 					dim[5] = Ctrl->A.S.v.v_shape;
-					if (Ctrl->A.S.symbol == GMT_SYMBOL_VECTOR_V4) {
-						double *this_rgb = NULL;
-						if (Ctrl->G.active)
-							this_rgb = Ctrl->G.fill.rgb;
-						else
-							this_rgb = GMT->session.no_rgb;
-						if (Ctrl->L.active) gmt_setpen (GMT, &Ctrl->W.pen);
+					if (Ctrl->L.active) gmt_setpen (GMT, &Ctrl->W.pen);
+					if (Ctrl->A.S.symbol == GMT_SYMBOL_VECTOR_V4) {	/* Old GMT4 vector selected */
+						double *this_rgb = (set_g_fill) ? Ctrl->G.fill.rgb : GMT->session.no_rgb;
 						psl_vector_v4 (PSL, plot_x, plot_y, dim, this_rgb, Ctrl->L.active);
 					}
 					else {
@@ -1003,29 +1135,29 @@ EXTERN_MSC int GMT_psvelo (void *V_API, int mode, void *args) {
 						dim[11] = (headpen_width > 0.0) ? headpen_width : 0.5 * Ctrl->W.pen.width;
 						if (Ctrl->A.S.v.status & PSL_VEC_FILL2)
 							gmt_setfill (GMT, &Ctrl->A.S.v.fill, Ctrl->L.active);
-						else if (Ctrl->G.active)
+						else if (set_g_fill)
 							gmt_setfill (GMT, &Ctrl->G.fill, Ctrl->L.active);
 						PSL_plotsymbol (PSL, plot_x, plot_y, dim, PSL_VECTOR);
 					}
 					if (Ctrl->A.S.v.status & PSL_VEC_OUTLINE2) gmt_setpen (GMT, &Ctrl->W.pen);
 
-					justify = plot_vx - plot_x > 0. ? PSL_MR : PSL_ML;
-					if (Ctrl->S.font.size > 0.0 && station_name)	/* 1 inch = 2.54 cm */
-						PSL_plottext (PSL, plot_x + (6 - justify) / 25.4 , plot_y, Ctrl->S.font.size, station_name, ANGLE, justify, FORM);
+					justify = ((plot_vx - plot_x) > 0.0) ? PSL_MR : PSL_ML;
+					if (Ctrl->S.font.size > 0.0 && station_name && station_name[0])	/* 1 inch = 2.54 cm */
+						PSL_plottext (PSL, plot_x + (PSL_MC - justify) / 25.4 , plot_y, Ctrl->S.font.size, station_name, ANGLE, justify, FORM);
 				}
-				else {
-					gmt_setfill (GMT, &Ctrl->G.fill, 1);
+				else {	/* vector too small, just place an circle there instead */
+					if (set_g_fill)
+						gmt_setfill (GMT, &Ctrl->G.fill, 1);
 					ssize = GMT_DOT_SIZE;
 					PSL_plotsymbol (PSL, plot_x, plot_y, &ssize, PSL_CIRCLE);
 					justify = PSL_TC;
-					if (Ctrl->S.font.size > 0.0 && station_name) {
+					if (Ctrl->S.font.size > 0.0 && station_name && station_name[0])	/* Place station name */
 						PSL_plottext (PSL, plot_x, plot_y - 1. / 25.4, Ctrl->S.font.size, station_name, ANGLE, justify, FORM);
-					}
-					/*  1 inch = 2.54 cm */
 				}
 				break;
 			case ANISO:
 				psvelo_trace_arrow (GMT, in[GMT_X], in[GMT_Y], vxy[0], vxy[1], Ctrl->S.scale, &plot_x, &plot_y, &plot_vx, &plot_vy);
+				gmt_setpen (GMT, &Ctrl->W.pen);
 				PSL_plotsegment (PSL, plot_x, plot_y, plot_vx, plot_vy);
 				break;
 			case CROSS:
@@ -1037,8 +1169,8 @@ EXTERN_MSC int GMT_psvelo (void *V_API, int mode, void *args) {
 				PSL_comment (PSL, "begin wedge number %li", n_rec);
 				gmt_geo_to_xy (GMT, in[GMT_X], in[GMT_Y], &plot_x, &plot_y);
 				psvelo_get_trans (GMT, in[GMT_X], in[GMT_Y], &t11, &t12, &t21, &t22);
-				psvelo_paint_wedge (PSL, plot_x, plot_y, spin, spinsig, Ctrl->S.scale, Ctrl->S.wedge_amp, t11,t12,t21,t22,
-					Ctrl->G.active, Ctrl->G.fill.rgb, Ctrl->E.active, Ctrl->E.fill.rgb, Ctrl->L.active);
+				psvelo_paint_wedge (PSL, plot_x, plot_y, spin, spinsig, Ctrl->S.scale, Ctrl->S.wedge_amp, t11, t12, t21, t22,
+					set_g_fill, Ctrl->G.fill.rgb, set_e_fill, Ctrl->E.fill.rgb, Ctrl->L.active);
 				break;
 		}
 	} while (true);
@@ -1049,8 +1181,6 @@ EXTERN_MSC int GMT_psvelo (void *V_API, int mode, void *args) {
 
 	GMT_Report (API, GMT_MSG_INFORMATION, "Number of records read: %li\n", n_rec);
 
-	if (Ctrl->D.active)  GMT_Report (API, GMT_MSG_INFORMATION, "Rescaling uncertainties by a factor of %f\n", Ctrl->D.scale);
-
 	if (!Ctrl->N.active) gmt_map_clip_off (GMT);
 
 	PSL_setdash (PSL, NULL, 0);
diff --git a/test/geodesy/geodesy_01.sh b/test/geodesy/geodesy_01.sh
index bd686ec5a1c..df3d068e326 100755
--- a/test/geodesy/geodesy_01.sh
+++ b/test/geodesy/geodesy_01.sh
@@ -52,7 +52,7 @@ EOF
 
 # simpler colors, labeled with following font
 gmt set FONT_ANNOT_PRIMARY Helvetica
-gmt psvelo -Se0.2/0.39+f18p -R -J -A0.25c+p0.25p+e -O -Umeca_4 << EOF >> $ps
+gmt psvelo -Se0.2/0.39+f18p -R -J -A0.25c+p0.25p+e+gblack -O -Umeca_4 << EOF >> $ps
 # Long.   Lat.   Evel   Nvel   Esig   Nsig  CorEN SITE
 # (deg)  (deg)    (mm/yr)        (mm/yr)
    0.    -8.     0.0    0.0     4.0    6.0  0.100  4x6
diff --git a/test/geodesy/geodesy_07.ps b/test/geodesy/geodesy_07.ps
new file mode 100644
index 00000000000..8f4bcaadf7c
Binary files /dev/null and b/test/geodesy/geodesy_07.ps differ
diff --git a/test/geodesy/geodesy_07.sh b/test/geodesy/geodesy_07.sh
new file mode 100755
index 00000000000..ec0d1e0446a
--- /dev/null
+++ b/test/geodesy/geodesy_07.sh
@@ -0,0 +1,44 @@
+#!/usr/bin/env bash
+# Test various combinations of pen/fill/cpt color in velo, using a small dataset
+# provided by Katherine Guns
+cat << EOF > GPS_test.txt
+235.6036	40.4414	-0.598	3.184	0.026	0.012	1e-05	CME6
+235.6036	40.4417	-0.537	3.142	0.024	0.013	1e-05	CME5
+235.6036	40.4417	-0.798	3.109	0.026	0.007	1e-05	CME1
+235.6919	40.2475	-1.89	3.287	0.014	0.005	1e-05	P157
+235.7172	40.5047	-0.231	2.678	0.023	0.008	1e-05	P159
+235.7629	40.691	0.267	2.143	0.026	0.008	1e-05	P162
+235.7869	40.6373	0.105	2.147	0.02	0.018	1e-05	P161
+235.8667	40.5512	-0.229	2.137	0.025	0.006	1e-05	P160
+235.8927	40.4224	-0.527	2.307	0.013	0.034	1e-05	P158
+235.9246	40.8763	0.18	1.724	0.017	0.008	1e-05	P058
+235.9427	40.2195	-1.437	2.486	0.01	0.008	1e-05	P163
+236.0323	40.7911	0.071	1.605	0.018	0.005	1e-05	P169
+236.0938	40.0244	-1.974	2.438	0.009	0.006	1e-05	P156
+236.1185	40.6686	-0.168	1.493	0.01	0.018	1e-05	P168
+236.1198	40.5437	-0.362	1.616	0.011	0.006	1e-05	P167
+236.1367	40.8802	-0.051	1.184	0.052	0.045	1e-05	P170
+236.1371	40.4351	-0.565	1.685	0.011	0.006	1e-05	P166
+236.1467	40.2455	-1.147	1.827	0.02	0.008	1e-05	P165
+236.3010	40.5753	-0.592	0.994	0.02	0.005	1e-05	P326
+236.3066	40.1192	-1.466	1.589	0.009	0.004	1e-05	P164
+236.3442	40.2568	-1.093	1.454	0.021	0.007	1e-05	P324
+236.4268	40.4787	-0.801	0.974	0.018	0.011	1e-05	P793
+236.4269	40.4788	-0.769	0.95	0.031	0.015	1e-05	P327
+236.5479	40.0822	-1.425	1.283	0.013	0.01	1e-05	P329
+EOF
+gmt begin -C geodesy_07 ps
+	gmt makecpt -Cturbo -T0/5
+	gmt set MAP_FRAME_TYPE plain
+	gmt subplot begin 4x2 -Fs5c -JM5c -R-124.6/-123.3/40/41 -A+jTR+gwhite -SCb -SRl+p -BWSrt -M6p -Ba1
+		gmt velo GPS_test.txt -A10p+e+a40+n6k+h0.5   -Se0.50/0/0 -Gorange -W1.5p -c
+		gmt velo GPS_test.txt -A10p+e+a40+n6k+h0.5   -Se0.50/0/0 -Gorange -W1.5p,orange -c
+		gmt velo GPS_test.txt -A10p+e+a40+n6k+h0.5   -Se0.50/0/0 -C -W1.5p -c
+		gmt velo GPS_test.txt -A10p+e+a40+n6k+h0.5   -Se0.50/0/0 -C -W1.5p+c -c
+		gmt velo GPS_test.txt -A10p+e+a40+n6k+h0.5+p -Se0.50/0/0 -C -W1.5p+c -c
+		gmt velo GPS_test.txt -A10p+e+a40+n6k+h0.5+p -Se0.50/0.95/0 -D10 -C -W1.5p+c -c
+		gmt velo GPS_test.txt -A10p+e+a40+n6k+h0.5+p -Se0.50/0.95/0 -D10 -C -W1.5p+c -L0.25p -Elightgray -c
+		gmt velo GPS_test.txt -A10p+e+a40+n6k+h0.5+p -Se0.50/0.0/0 -D10 -C -W1.5p+c
+		gmt velo GPS_test.txt -A10p+e+a40+n6k+h0.5+p -Se0.50/0.95/0 -D10 -C -W1.5p+c -L0.25p+c -c -Zn
+	gmt subplot end
+gmt end show
