Add -l to grdvector for legend entry (#6199)

* Add -l to grdvector for legend entry

The magnitude of the legend vector (in data units) can be set via -S modifier +sref.

* Add test script

* Truncate vector size if head width exceeds H

* Update docs and how we set default

* Add test PS

* Update grdvector.rst

* Update grdvector.c

* Minor fixes and PS update

* Update grdvector.rst

* Add support for geovector legend

* Oops

* More oops

doc/rst/source/grdvector.rst

@@ -20,7 +20,7 @@ Synopsis
 [ |-I|\ [**x**]\ *dx*\ [/*dy*] ]
 [ |-N| ] [ |-Q|\ *parameters* ]
 [ |SYN_OPT-R| ]
-[ |-S|\ [**i**\|\ **l**]\ *scale* ]
+[ |-S|\ [**i**\|\ **l**]\ *scale*\ [**+c**\ [[*slon*/]\ *slat*]][**+s**\ *refsize*] ]
 [ |-T| ]
 [ |SYN_OPT-U| ]
 [ |SYN_OPT-V| ]
@@ -29,6 +29,7 @@ Synopsis
 [ |SYN_OPT-Y| ]
 [ |-Z| ]
 [ |SYN_OPT-f| ]
+[ |SYN_OPT-l| ]
 [ |SYN_OPT-p| ]
 [ |SYN_OPT-t| ]
 [ |SYN_OPT--| ]
@@ -111,22 +112,33 @@ Optional Arguments
 
 .. _-S:
 
-**-S**\ [**i**\|\ **l**]\ *scale*
+**-S**\ [**i**\|\ **l**]\ *scale*\ [**+c**\ [[*slon*/]\ *slat*]][**+s**\ *refsize*]
     Sets scale for vector plot lengths in data units per plot distance measurement unit.
     Append **c**, **i**, or **p** to indicate the desired plot distance measurement
     unit (cm, inch, or point); if no unit is given we use the default value that
     is controlled by :term:`PROJ_LENGTH_UNIT`.  Vector lengths converted via plot unit
     scaling will plot as straight Cartesian vectors and their lengths are not
-    affected by map projection and coordinate locations.
+    affected by map projections and coordinate locations.
     For geographic data you may alternatively give *scale* in data units per map distance
     unit (see `Units`_). Then, your vector magnitudes (in data units) are scaled to map
     *distances* in the given distance unit, and finally projected onto the Earth to give
     *plot* dimensions.  These are geo-vectors that follow great circle paths and their
     lengths may be affected by the map projection and their coordinates.  Finally, use
     **-Si** if it is simpler to give the reciprocal scale in plot length or distance units
-    per data unit.  To report the minimum, maximum, and mean data and plot vector lengths
-    of all vectors plotted, use **-V**.  Alternatively, use **-Sl**\ *length* to set a fixed
-    plot length for all vectors.
+    per data unit.  Alternatively, use **-Sl**\ *length* to set a fixed plot length for all
+    vectors. To report the minimum, maximum, and mean data and plot vector lengths
+    of all vectors plotted, use **-V**. If an automatic legend entry is desired via **-l**,
+    or or two modifiers will be required:
+
+    - **+c**\ [[*slon*/]\ *slat*] controls where on a geographic map a geovector's *refsize*
+      length applies. The modifier is neither needed nor available when plotting Cartesian vectors.
+      The length is calculated for latitude *slat* (optionally supply longitude *slon* for
+      oblique projections [default is central meridian]). If **+c** is given with no arguments
+      then we select the reference length origin to be the middle of the map.
+    - **+s**\ *refsize* sets the desired reference vector magnitude in data units. E.g., for a
+      reference length of 25 mm/yr for plate motions, use modifier **+s**\ 25 with a corresponding
+      option **-l**\ "Velocity (25 mm/yr)".  If *refsize* is not specified we default to the *scale*
+      given above.
 
 .. _-T:
 
@@ -165,6 +177,9 @@ Optional Arguments
 .. |Add_-f| unicode:: 0x20 .. just an invisible code
 .. include:: explain_-f.rst_
 
+.. |Add_-l| unicode:: 0x20 .. just an invisible code
+.. include:: explain_-l.rst_
+
 .. |Add_perspective| unicode:: 0x20 .. just an invisible code
 .. include:: explain_perspective.rst_
 

src/gmt_init.c

@@ -15936,56 +15936,64 @@ int gmt_init_vector_param (struct GMT_CTRL *GMT, struct GMT_SYMBOL *S, bool set,
 	return 0;
 }
 
-GMT_LOCAL unsigned int gmtinit_get_length (struct GMT_CTRL *GMT, char symbol, char *string, float *value, bool *user_unit) {
-	/* Used by gmt_parse_vector to set length scales or limits related to vectors */
+GMT_LOCAL unsigned int gmtinit_get_length (struct GMT_CTRL *GMT, char symbol, char *string, bool norm, bool is_grdvector, float *value, bool *user_unit) {
+	/* Used by gmt_parse_vector to set length scales or limits related to vectors.
+	 * Note: is_grdvector is true when being called from grdvector because we need to honor normalizations +n<lengtH> without unit to mean unit = q */
 	unsigned int error = GMT_NOERROR;
-	size_t len = strlen (string) - 1;	/* Position of last character in string */
+	size_t len = (string) ? strlen (string) : 0;	/* Length of string */
 
+	if (len == 0) {
+		GMT_Report (GMT->parent, GMT_MSG_ERROR, "No vector attributes given for symbol: %c\n", symbol);
+		return GMT_NOERROR;
+	}
+	len -= 1;	/* Position of last character in string */
 	if (string[0] == 'i')	/* Get reciprocal value first, then deal with unit */
 		*value = 1.0f / (float)atof (&string[1]);
 	else
 		*value = (float)atof (string);
 	*user_unit = false;
 	/* value is now the normalizing length in given units, not (yet) converted to inches or degrees (but see next lines) */
 	if (symbol == '=') {	/* Since we have map distance units for geovectors we convert to km */
-		if (len) {	/* Examine if a unit was given */
-			if (string[len] == 'q') {	/* Got data units so we assume scale is in km/unit */
-				*user_unit = true;	/* This turns off unit conversion during reading */
-				*value = (float)gmtlib_conv_distance (GMT, *value, 'k', 'k');	/* Convert to km/unit */
-			}
-			else if (strchr (GMT_DIM_UNITS, string[len])) {	/* Confused user gave geovector length in c|i|p, this is an error */
-				GMT_Report (GMT->parent, GMT_MSG_ERROR, "Vector scale or length limit for geovectors must be given in units of %s or q (data quantity unit) [k]!\n", GMT_LEN_UNITS);
-				error++;
-			}
-			else if (strchr (GMT_LEN_UNITS, string[len]))	/* Got length with valid unit, otherwise we assume it was given in km */
-				*value = (float)gmtlib_conv_distance (GMT, *value, string[len], 'k');	/* Convert to km first */
-			else if (strchr (".0123456789", string[len]) == NULL) {	/* Unless this triggers we got km */
-				GMT_Report (GMT->parent, GMT_MSG_ERROR, "Unrecognized map distance unit: %c\n", string[len]);
-				error++;
-			}
+		if (string[len] == 'q') {	/* Got data units so we assume scale is in km/unit */
+			*user_unit = true;	/* This turns off unit conversion during reading */
+			if (!norm) *value = (float)gmtlib_conv_distance (GMT, *value, 'k', 'k');	/* Convert to km/unit */
+		}
+		else if (strchr (GMT_DIM_UNITS, string[len])) {	/* Confused user gave geovector length in c|i|p, this is an error */
+			GMT_Report (GMT->parent, GMT_MSG_ERROR, "Vector scale or length limit for geovectors must be given in units of %s or q (data quantity unit) [k]!\n", GMT_LEN_UNITS_DISPLAY);
+			error++;
+		}
+		else if (strchr (GMT_LEN_UNITS, string[len]))	/* Got length with valid unit, otherwise we assume it was given in km */
+			*value = (float)gmtlib_conv_distance (GMT, *value, string[len], 'k');	/* Convert to km first */
+		else if (strchr (".0123456789", string[len]) == NULL) {	/* Unless this triggers we got km */
+			GMT_Report (GMT->parent, GMT_MSG_ERROR, "Unrecognized map distance unit: %c\n", string[len]);
+			error++;
 		}
 	}
 	else {	/* Cartesian vector */
-		if (len) {	/* Examine if a unit was given */
-			int j = GMT->current.setting.proj_length_unit;	/* Default unit index */
-			if (string[len] == 'q') {	/* Got data units so we assume scale is in c|i|p per user quantity unit */
-				*user_unit = true;	/* This turns off unit conversion during reading */
-				*value *= (float)GMT->session.u2u[j][GMT_INCH];	/* This places unit conversion into the scale instead */
-			}
-			else if (strchr (GMT_LEN_UNITS, string[len])) {	/* Confused user gave Cartesian vector length in d|m|s|e|f|k|M|n|u, this is an error */
-				GMT_Report (GMT->parent, GMT_MSG_ERROR, "Vector scale or length limit for Caretesian vectors must be given in units of %s or q (data quantity unit) [k]!\n",
-					GMT_DIM_UNITS, GMT->session.unit_name[GMT->current.setting.proj_length_unit][0]);
-				error++;
-			}
-			else if (strchr (GMT_DIM_UNITS, string[len])) {	/* Got length with valid unit c|i|p, convert to inch */
-				j = gmt_get_dim_unit (GMT, string[len]);	/* Override j */
-				*value *= (float)GMT->session.u2u[j][GMT_INCH];
-			}
-			else if (strchr (".0123456789", string[len]) == NULL) {	/* Unless this triggers we got default plot units */
-				GMT_Report (GMT->parent, GMT_MSG_ERROR, "Unrecognized plot length unit: %c\n", string[len]);
-				error++;
+		int j = GMT->current.setting.proj_length_unit;	/* Default unit index */
+		if (string[len] == 'q') {	/* Got data units so we assume scale is in c|i|p per user quantity unit */
+			*user_unit = true;	/* This turns off unit conversion during reading */
+			if (!norm) *value *= (float)GMT->session.u2u[j][GMT_INCH];	/* This places unit conversion into the scale instead */
+		}
+		else if (strchr (GMT_LEN_UNITS, string[len])) {	/* Confused user gave Cartesian vector length in d|m|s|e|f|k|M|n|u, this is an error */
+			GMT_Report (GMT->parent, GMT_MSG_ERROR, "Vector scale or length limit for Cartesian vectors must be given in units of %s or q (data quantity unit) [k]!\n",
+				GMT_DIM_UNITS_DISPLAY, GMT->session.unit_name[GMT->current.setting.proj_length_unit][0]);
+			error++;
+		}
+		else if (strchr (GMT_DIM_UNITS, string[len])) {	/* Got length with valid unit c|i|p, convert to inch */
+			j = gmt_get_dim_unit (GMT, string[len]);	/* Override j */
+			*value *= (float)GMT->session.u2u[j][GMT_INCH];
+		}
+		else if (strchr (".0123456789", string[len]) == NULL) {	/* Unless this triggers we got default plot units */
+			GMT_Report (GMT->parent, GMT_MSG_ERROR, "Unrecognized plot length unit: %c\n", string[len]);
+			error++;
+		}
+		else {	/* Here we got a plot distance with or without unit; if the latter we use current default unit to convert to inch */
+			if (is_grdvector) {
+				if (!norm) *value *= (float)GMT->session.u2u[j][GMT_INCH];
+				*user_unit = norm;
 			}
-			else	/* Here we got a plot distance with or without unit; if the latter we use current default unit to convert to inch */
+			else
 				*value *= (float)GMT->session.u2u[j][GMT_INCH];
 		}
 	}
@@ -15997,7 +16005,7 @@ GMT_LOCAL unsigned int gmtinit_get_length (struct GMT_CTRL *GMT, char symbol, ch
 int gmt_parse_vector (struct GMT_CTRL *GMT, char symbol, char *text, struct GMT_SYMBOL *S) {
 
 	unsigned int pos = 0, k, f = 0, end, error = 0;
-	bool p_opt = false, g_opt = false;
+	bool p_opt = false, g_opt = false, is_grdvector = !strncmp (GMT->init.module_name, "grdvector", 9U);
 	int j;
 	char p[GMT_BUFSIZ], *c = NULL;
 	double value[2];
@@ -16144,7 +16152,7 @@ int gmt_parse_vector (struct GMT_CTRL *GMT, char symbol, char *text, struct GMT_
 				if (p[1] == '\0')	/* No shrink, and no skipping heads regardless of vector length */
 					S->v.status |= PSL_VEC_LINE;
 				else {	/* Parse the cutoff size */
-					error += gmtinit_get_length (GMT, symbol, &p[1], &(S->v.v_norm), &(S->v.v_norm_d));
+					error += gmtinit_get_length (GMT, symbol, &p[1], true, is_grdvector, &(S->v.v_norm), &(S->v.v_norm_d));
 					if (symbol == '=' && !S->v.v_norm_d)	/* Since norm distance is now in km we convert to spherical degrees */
 						S->v.v_norm /= (float)GMT->current.proj.DIST_KM_PR_DEG;	/* Finally, convert km to degrees */
 					/* Here, v_norm is either in inches (if Cartesian vector), spherical degrees (if geovector), or unitless (if q was used) */
@@ -16217,7 +16225,7 @@ int gmt_parse_vector (struct GMT_CTRL *GMT, char symbol, char *text, struct GMT_
 			case 'z':	/* Input (angle,length) are vector components (dx,dy) instead */
 				S->v.status |= PSL_VEC_COMPONENTS;
 				S->v.status |= PSL_VEC_MAGNIFY;
-				if (p[1]) error += gmtinit_get_length (GMT, symbol, &p[1], &(S->v.comp_scale), &(S->v.v_unit_d));
+				if (p[1]) error += gmtinit_get_length (GMT, symbol, &p[1], false, is_grdvector, &(S->v.comp_scale), &(S->v.v_unit_d));
 				break;
 			case 'v':	/* Scale vector polar length component, or get inverse scale, or a fixed magnitude */
 				S->v.status |= PSL_VEC_MAGNIFY;
@@ -16227,7 +16235,7 @@ int gmt_parse_vector (struct GMT_CTRL *GMT, char symbol, char *text, struct GMT_
 				}
 				else
 					j = 1;
-				error += gmtinit_get_length (GMT, symbol, &p[j], &(S->v.comp_scale), &(S->v.v_unit_d));
+				error += gmtinit_get_length (GMT, symbol, &p[j], false, is_grdvector, &(S->v.comp_scale), &(S->v.v_unit_d));
 				break;
 			default:
 				GMT_Report (GMT->parent, GMT_MSG_ERROR, "Bad modifier +%c\n", p[0]);

src/gmt_plot.c

@@ -4566,24 +4566,23 @@ GMT_LOCAL uint64_t gmtplot_geo_polygon_segment (struct GMT_CTRL *GMT, struct GMT
 	return (k);	/* Number of points plotted */
 }
 
-GMT_LOCAL float gmtplot_inch_to_degree_scale (struct GMT_CTRL *GMT, double lon0, double lat0, double azimuth) {
+double gmt_inch_to_degree_scale (struct GMT_CTRL *GMT, double lon0, double lat0, double azimuth) {
 	/* Used to determine delta radius in degrees for thickness of vector lines to be drawn as small or
 	 * great circles.  We must convert pen thickness to some sense of spherical degrees.
 	 * Determine the map scale at (lon, lat) in a direction normal to the vector and use that
 	 * to scale items in inches to spherical degrees.  We repeat this for all locations except
 	 * if the projection is perspective, when we pick the center as the map distortion will be the least here.
 	 * This scaling is approximate only but needed to convert geovector head lengths to degrees. */
 
-	double tlon, tlat, x0, y0, dx, x1, y1, length;
-	float scale;
+	double tlon, tlat, x0, y0, dx, x1, y1, length, scale;
 
 	length = 0.001 * (GMT->common.R.wesn[YHI] - GMT->common.R.wesn[YLO]);		/* 0.1 percent of latitude extent is fairly small */
 	gmt_geo_to_xy (GMT, lon0, lat0, &x0, &y0);					/* Get map position in inches for close point */
 	gmtlib_get_point_from_r_az (GMT, lon0, lat0, length, azimuth-90.0, &tlon, &tlat);	/* ANearby arbitrary 2nd point in direction normal to vector at (lon0,lat0) */
 	gmt_geo_to_xy (GMT, tlon, tlat, &x1, &y1);					/* Get map position in inches for close point */
 	dx = fabs (x1 - x0);
 	if (dx > (0.25 * GMT->current.map.half_width)) dx = GMT->current.map.width - dx;
-	scale = (float) (length / hypot (dx, y1 - y0));				/* This scales a length in inches to degrees, approximately */
+	scale = length / hypot (dx, y1 - y0);				/* This scales a length in inches to degrees, approximately */
 	return (scale);
 }
 
@@ -4956,21 +4955,20 @@ GMT_LOCAL unsigned int gmtplot_geo_vector_smallcircle (struct GMT_CTRL *GMT, dou
 		max_length = MAX (dr[0], dr[1]);
 		if (max_length > h_length_limit) {
 			s2 = h_length_limit / max_length;
+			if (s2 < S->v.v_norm_limit) s2 = S->v.v_norm_limit;
 			warn = 2;
 		}
 	}
 
 	/* Might have to shrink things */
-	if (C.r0 < S->v.v_norm)
-		s1 = C.r0 / S->v.v_norm;
+	s1 = gmt_get_vector_shrinking (GMT, &(S->v), S->v.value, C.r0); /* Vector attribute shrinking factor or 1 */
+
 	if (s1 < s2) {	/* Pick the smallest scale required for head shrinking */
 		s = s1;
 		warn = 0;	/* When requesting +n there is no warning unless s2 is the smaller scale */
 	}
 	else
 		s = s2;
-	if (s < S->v.v_norm_limit) s = S->v.v_norm_limit;
-	if (s1 < S->v.v_norm_limit) s1 = S->v.v_norm_limit;
 	if (s1 > s) s1 = s;
 	head_length = s * S->size_x;
 	arc_width   = s1 * S->v.v_width;	/* Use scale s1 for pen shrinking */
@@ -5209,21 +5207,21 @@ GMT_LOCAL unsigned int gmtplot_geo_vector_greatcircle (struct GMT_CTRL *GMT, dou
 		max_length = MAX (dr[0], dr[1]);
 		if (max_length > h_length_limit) {
 			s2 = h_length_limit / max_length;
+			if (s2 < S->v.v_norm_limit) s2 = S->v.v_norm_limit;
 			warn = 2;
 		}
 	}
 
 	/* Might have to shrink things */
-	if (C.r0 < S->v.v_norm)
-		s1 = C.r0 / S->v.v_norm;
+
+	s1 = gmt_get_vector_shrinking (GMT, &(S->v), S->v.value, C.r0); /* Vector attribute shrinking factor or 1 */
+
 	if (s1 < s2) {	/* Pick the smallest scale required for head shrinking */
 		s = s1;
 		warn = 0;	/* When requesting +n there is no warning unless s2 is the smaller scale */
 	}
 	else
 		s = s2;
-	if (s < S->v.v_norm_limit) s = S->v.v_norm_limit;
-	if (s1 < S->v.v_norm_limit) s1 = S->v.v_norm_limit;
 	if (s1 > s) s1 = s;
 	head_length = s * S->size_x;
 	arc_width   = s1 * S->v.v_width;	/* Use scale s1 for pen shrinking */
@@ -9797,12 +9795,12 @@ unsigned int gmt_geo_vector (struct GMT_CTRL *GMT, double lon0, double lat0, dou
 		if (gmt_M_is_perspective (GMT)) {
 			double clon, clat;
 			gmt_xy_to_geo (GMT, &clon, &clat, GMT->current.map.half_width, GMT->current.map.half_height);	/* Geographic coordinates of middle map point */
-			S->v.scale = gmtplot_inch_to_degree_scale (GMT, clon, clat, azimuth);
+			S->v.scale = (float)gmt_inch_to_degree_scale (GMT, clon, clat, azimuth);
 			S->v.status |= PSL_VEC_SCALE;
 			GMT_Report (GMT->parent, GMT_MSG_DEBUG, "Vector stem scale is %g degrees/inch at (%g, %g) for az = %g\n", S->v.scale, clon, clat, azimuth);
 		}
 		else {	/* Set scale each time locally */
-			S->v.scale = gmtplot_inch_to_degree_scale (GMT, lon0, lat0, azimuth);
+			S->v.scale = (float)gmt_inch_to_degree_scale (GMT, lon0, lat0, azimuth);
 			GMT_Report (GMT->parent, GMT_MSG_DEBUG, "Vector stem scale is %g degrees/inch at (%g, %g) for az = %g\n", S->v.scale, lon0, lat0, azimuth);
 		}
 	}

src/gmt_plot.h

@@ -123,6 +123,7 @@ struct GMT_VECT_ATTR {
 	float h_width;		/* Width of vector head in inches */
 	float pole[2];		/* Longitude and latitude of geovector pole */
 	float scale;		/* Converts inches to spherical degrees */
+	float value;		/* Original data quantity */
 	float comp_scale;	/* Converts hypot (dx, dy) to inches */
 	float v_trim[2];	/* Offsets from begin/end point in inches */
 	struct GMT_PEN pen;	/* Pen for outline of head */

src/gmt_prototypes.h

@@ -254,6 +254,7 @@ EXTERN_MSC struct GMT_GRID * gmt_duplicate_grid (struct GMT_CTRL *GMT, struct GM
 #ifdef _POSTSCRIPTLIGHT_H
 /* gmt_plot.c prototypes only included if postscriptlight has been included */
 
+EXTERN_MSC double gmt_inch_to_degree_scale (struct GMT_CTRL *GMT, double lon0, double lat0, double azimuth);
 EXTERN_MSC bool gmt_text_is_latex (struct GMT_CTRL *GMT, const char *string);
 EXTERN_MSC void gmt_map_text (struct GMT_CTRL *GMT, double x, double y, struct GMT_FONT *font, char *label, double angle, int just, unsigned int form);
 EXTERN_MSC void gmt_map_title (struct GMT_CTRL *GMT, double x, double y);