1 /*--------------------------------------------------------------------
2 *
3 * Copyright (c) 1991-2021 by the GMT Team (https://www.generic-mapping-tools.org/team.html)
4 * See LICENSE.TXT file for copying and redistribution conditions.
5 *
6 * This program is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU Lesser General Public License as published by
8 * the Free Software Foundation; version 3 or any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU Lesser General Public License for more details.
14 *
15 * Contact info: www.generic-mapping-tools.org
16 *--------------------------------------------------------------------*/
17 /*
18 * Brief synopsis: grdvector reads 2 grid files that contains the 2 components of a vector
19 * field (cartesian or polar) and plots vectors at the grid positions.
20 *
21 * Author: Paul Wessel
22 * Date: 1-JAN-2010
23 * Version: 6 API
24 */
25
26 #define THIS_MODULE_CLASSIC_NAME "grdvector"
27 #define THIS_MODULE_MODERN_NAME "grdvector"
28 #define THIS_MODULE_LIB "core"
29 #define THIS_MODULE_PURPOSE "Plot vector field from two component grids"
30 #define THIS_MODULE_KEYS "<G{2,CC(,>X}"
31 #define THIS_MODULE_NEEDS "Jg"
32
33 #include "gmt_dev.h"
34
35 #define THIS_MODULE_OPTIONS "->BJKOPRUVXYfptxy" GMT_OPT("c")
36
37 struct GRDVECTOR_CTRL {
38 struct GRDVECTOR_In {
39 bool active;
40 char *file[2];
41 } In;
42 struct GRDVECTOR_A { /* -A */
43 bool active;
44 } A;
45 struct GRDVECTOR_C { /* -C<cpt>[+i<dz>] */
46 bool active;
47 double dz;
48 char *file;
49 } C;
50 struct GRDVECTOR_G { /* -G<fill> */
51 bool active;
52 struct GMT_FILL fill;
53 } G;
54 struct GRDVECTOR_I { /* -I[x]<dx>[/<dy>] */
55 bool active;
56 unsigned int mode;
57 } I;
58 struct GRDVECTOR_N { /* -N */
59 bool active;
60 } N;
61 struct GRDVECTOR_Q { /* -Q<size>[+<mods>] */
62 bool active;
63 struct GMT_SYMBOL S;
64 } Q;
65 struct GRDVECTOR_S { /* -S[l|i]<length|scale>[<unit>] */
66 bool active;
67 bool constant;
68 bool invert;
69 char unit;
70 double factor;
71 } S;
72 struct GRDVECTOR_T { /* -T */
73 bool active;
74 } T;
75 struct GRDVECTOR_W { /* -W<pen> */
76 bool active;
77 bool cpt_effect;
78 struct GMT_PEN pen;
79 } W;
80 struct GRDVECTOR_Z { /* -Z */
81 bool active;
82 } Z;
83 };
84
New_Ctrl(struct GMT_CTRL * GMT)85 static void *New_Ctrl (struct GMT_CTRL *GMT) { /* Allocate and initialize a new control structure */
86 struct GRDVECTOR_CTRL *C = NULL;
87 char unit[5] = "cimp";
88
89 C = gmt_M_memory (GMT, NULL, 1, struct GRDVECTOR_CTRL);
90
91 /* Initialize values whose defaults are not 0/false/NULL */
92 gmt_init_fill (GMT, &C->G.fill, -1.0, -1.0, -1.0);
93 C->Q.S.symbol = PSL_VECTOR;
94 C->W.pen = GMT->current.setting.map_default_pen;
95 C->S.factor = 1.0;
96 C->S.unit = unit[GMT->current.setting.proj_length_unit];
97 return (C);
98 }
99
Free_Ctrl(struct GMT_CTRL * GMT,struct GRDVECTOR_CTRL * C)100 static void Free_Ctrl (struct GMT_CTRL *GMT, struct GRDVECTOR_CTRL *C) { /* Deallocate control structure */
101 if (!C) return;
102 gmt_M_str_free (C->In.file[GMT_IN]);
103 gmt_M_str_free (C->In.file[GMT_OUT]);
104 gmt_M_str_free (C->C.file);
105 gmt_M_free (GMT, C);
106 }
107
usage(struct GMTAPI_CTRL * API,int level)108 static int usage (struct GMTAPI_CTRL *API, int level) {
109 const char *name = gmt_show_name_and_purpose (API, THIS_MODULE_LIB, THIS_MODULE_CLASSIC_NAME, THIS_MODULE_PURPOSE);
110 if (level == GMT_MODULE_PURPOSE) return (GMT_NOERROR);
111 GMT_Usage (API, 0, "usage: %s <gridx> <gridy> %s [-A] [%s] [-C[<cpt>]] [-G<fill>] [-I[x]<dx>/<dy>] "
112 "%s[-N] %s%s[-Q<params>] [%s] [-S[i|l]<length>|<scale>] [-T] [%s] [%s] [-W<pen>] [%s] [%s] [-Z] "
113 "%s [%s] [%s] [%s] [%s]\n", name, GMT_J_OPT, GMT_B_OPT, API->K_OPT, API->O_OPT, API->P_OPT,
114 GMT_Rgeo_OPT, GMT_U_OPT, GMT_V_OPT, GMT_X_OPT, GMT_Y_OPT, API->c_OPT, GMT_f_OPT, GMT_p_OPT, GMT_t_OPT, GMT_PAR_OPT);
115
116 if (level == GMT_SYNOPSIS) return (GMT_MODULE_SYNOPSIS);
117
118 GMT_Message (API, GMT_TIME_NONE, " REQUIRED ARGUMENTS:\n");
119 GMT_Usage (API, 1, "\n<gridx> <gridy> are grid files with the two vector components.");
120 GMT_Option (API, "J-");
121 GMT_Message (API, GMT_TIME_NONE, "\n OPTIONAL ARGUMENTS:\n");
122 GMT_Usage (API, 1, "\n-A Grids have polar (r, theta) components [Default is Cartesian (x, y) components].");
123 GMT_Option (API, "B-");
124 GMT_Usage (API, 1, "\n-C[<cpt>]");
125 GMT_Usage (API, -2, "Color palette file to convert vector length to colors. Optionally, name a master cpt "
126 "to automatically assign continuous colors over the data range [%s]; if so, "
127 "optionally append +i<dz> to quantize the range [the exact grid range]. "
128 "Another option is to specify -Ccolor1,color2[,color3,...] to build a linear "
129 "continuous cpt from those colors automatically.", API->GMT->current.setting.cpt);
130 gmt_fill_syntax (API->GMT, 'G', NULL, "Select vector fill [Default is outlines only].");
131 GMT_Usage (API, 1, "\n-I[x]<dx>/<dy>");
132 GMT_Usage (API, -2, "Plot only those nodes that are <dx>/<dy> apart [Default is all nodes]. "
133 "Optionally, use -Ix<fact>[/<yfact>] to give multiples of grid spacing.");
134 GMT_Option (API, "K");
135 GMT_Usage (API, 1, "\n-N Do Not clip vectors that exceed the map boundaries [Default will clip].");
136 GMT_Option (API, "O,P");
137 GMT_Usage (API, 1, "\n-Q<params>");
138 GMT_Usage (API, -2, "Modify vector attributes [Default gives stick-plot].");
139 gmt_vector_syntax (API->GMT, 15, 3);
140 GMT_Option (API, "R");
141 GMT_Usage (API, 1, "\n-S[i|l]<length>|<scale>");
142 GMT_Usage (API, -2, "Set lengths for vectors in <data-units> per length unit (e.g., 10 nTesla/yr per cm).");
143 GMT_Usage (API, 2, "%s Cartesian vectors: Append %s to indicate cm, inch, or point as the desired plot length unit [%s]. "
144 "These vectors are straight and plot lengths are independent of projection.",
145 GMT_LINE_BULLET, GMT_DIM_UNITS_DISPLAY, API->GMT->session.unit_name[API->GMT->current.setting.proj_length_unit]);
146 GMT_Usage (API, 2, "%s Geographic vectors: Alternatively give <data-units> per map distance unit "
147 "by appending any of the distance units in %s to the length. "
148 "These vectors may curve and plot lengths may depend on the projection.",
149 GMT_LINE_BULLET, GMT_LEN_UNITS_DISPLAY);
150 GMT_Usage (API, -2, "Optional directives:");
151 GMT_Usage (API, 3, "i: The given <scale> is the reciprocal scale, e.g., in %s or km per <data-unit>.",
152 API->GMT->session.unit_name[API->GMT->current.setting.proj_length_unit]);
153 GMT_Usage (API, 3, "l: Fixed length (in given unit) for all vectors.");
154 GMT_Usage (API, -2, "Note: Use -V to see the min, max, and mean vector length of plotted vectors.");
155 GMT_Usage (API, 1, "\n-T Transform angles for Cartesian grids when x- and y-scales differ [Leave alone].");
156 GMT_Option (API, "U,V");
157 gmt_pen_syntax (API->GMT, 'W', NULL, "Set pen attributes.", NULL, 0);
158 GMT_Usage (API, -2, "Default pen attributes [%s].", gmt_putpen(API->GMT, &API->GMT->current.setting.map_default_pen));
159 GMT_Option (API, "X");
160 GMT_Usage (API, 1, "\n-Z The theta grid provided has azimuths rather than directions (implies -A).");
161 GMT_Option (API, "c,f,p,t,.");
162
163 return (GMT_MODULE_USAGE);
164 }
165
parse(struct GMT_CTRL * GMT,struct GRDVECTOR_CTRL * Ctrl,struct GMT_OPTION * options)166 static int parse (struct GMT_CTRL *GMT, struct GRDVECTOR_CTRL *Ctrl, struct GMT_OPTION *options) {
167 /* This parses the options provided to grdvector and sets parameters in Ctrl.
168 * Note Ctrl has already been initialized and non-zero default values set.
169 * Any GMT common options will override values set previously by other commands.
170 * It also replaces any file names specified as input or output with the data ID
171 * returned when registering these sources/destinations with the API.
172 */
173
174 unsigned int n_errors = 0, n_files = 0;
175 int j;
176 size_t len;
177 char txt_a[GMT_LEN256] = {""}, txt_b[GMT_LEN256] = {""}, txt_c[GMT_LEN256] = {""}, symbol;
178 struct GMT_OPTION *opt = NULL;
179 struct GMTAPI_CTRL *API = GMT->parent;
180
181 symbol = (gmt_M_is_geographic (GMT, GMT_IN)) ? '=' : 'v'; /* Type of vector */
182
183 for (opt = options; opt; opt = opt->next) { /* Process all the options given */
184
185 switch (opt->option) {
186 case '<': /* Input file (only one is accepted) */
187 Ctrl->In.active = true;
188 if (n_files >= 2) {n_errors++; continue; }
189 if (opt->arg[0]) Ctrl->In.file[n_files] = strdup (opt->arg);
190 if (GMT_Get_FilePath (API, GMT_IS_DATASET, GMT_IN, GMT_FILE_REMOTE, &(Ctrl->In.file[n_files]))) n_errors++;
191 n_files++;
192 break;
193
194 /* Processes program-specific parameters */
195
196 case 'A': /* Polar data grids */
197 n_errors += gmt_M_repeated_module_option (API, Ctrl->A.active);
198 Ctrl->A.active = true;
199 break;
200 case 'C': /* Vary symbol color with z */
201 n_errors += gmt_M_repeated_module_option (API, Ctrl->C.active);
202 Ctrl->C.active = true;
203 gmt_M_str_free (Ctrl->C.file);
204 if (opt->arg[0]) Ctrl->C.file = strdup (opt->arg);
205 gmt_cpt_interval_modifier (GMT, &(Ctrl->C.file), &(Ctrl->C.dz));
206 break;
207 case 'E': /* Center vectors [OBSOLETE; use modifier +jc in -Q ] */
208 if (gmt_M_compat_check (GMT, 4)) {
209 GMT_Report (API, GMT_MSG_COMPAT, "Option -E is deprecated; use modifier +jc in -Q instead.\n");
210 Ctrl->Q.S.v.status |= PSL_VEC_JUST_C;
211 }
212 else
213 n_errors += gmt_default_error (GMT, opt->option);
214 break;
215 case 'G': /* Set fill for vectors */
216 n_errors += gmt_M_repeated_module_option (API, Ctrl->G.active);
217 Ctrl->G.active = true;
218 if (gmt_getfill (GMT, opt->arg, &Ctrl->G.fill)) {
219 gmt_fill_syntax (GMT, 'G', NULL, " ");
220 n_errors++;
221 }
222 break;
223 case 'I': /* Only use gridnodes GMT->common.R.inc[GMT_X],GMT->common.R.inc[GMT_Y] apart */
224 n_errors += gmt_M_repeated_module_option (API, Ctrl->I.active);
225 Ctrl->I.active = true;
226 if (opt->arg[0] == 'x') Ctrl->I.mode = 1;
227 n_errors += gmt_parse_inc_option (GMT, 'I', &opt->arg[Ctrl->I.mode]);
228 break;
229 case 'N': /* Do not clip at border */
230 n_errors += gmt_M_repeated_module_option (API, Ctrl->N.active);
231 Ctrl->N.active = true;
232 break;
233 case 'Q': /* Vector plots, with parameters */
234 n_errors += gmt_M_repeated_module_option (API, Ctrl->Q.active);
235 Ctrl->Q.active = true;
236 if (gmt_M_compat_check (GMT, 4) && (strchr (opt->arg, '/') && !strchr (opt->arg, '+'))) { /* Old-style args */
237 if (gmt_M_is_geographic (GMT, GMT_IN))
238 GMT_Report (API, GMT_MSG_COMPAT, "Vector arrowwidth/headlength/headwidth is deprecated for geo-vectors; see -Q documentation.\n");
239 Ctrl->Q.S.v.status = PSL_VEC_END + PSL_VEC_FILL + PSL_VEC_OUTLINE;
240 Ctrl->Q.S.size_x = VECTOR_HEAD_LENGTH * GMT->session.u2u[GMT_PT][GMT_INCH]; /* 9p */
241 Ctrl->Q.S.v.h_length = (float)Ctrl->Q.S.size_x; /* 9p */
242 Ctrl->Q.S.v.v_angle = 60.0f;
243 Ctrl->Q.S.v.pen = GMT->current.setting.map_default_pen;
244 for (j = 0; opt->arg[j] && opt->arg[j] != 'n'; j++);
245 if (opt->arg[j]) { /* Normalize option used */
246 Ctrl->Q.S.v.v_norm = (float)gmt_M_to_inch (GMT, &opt->arg[j+1]); /* Getting inches directly here */
247 n_errors += gmt_M_check_condition (GMT, Ctrl->Q.S.v.v_norm <= 0.0, "Option -Qn: No reference length given\n");
248 opt->arg[j] = '\0'; /* Temporarily chop of the n<norm> string */
249 }
250 if (opt->arg[0] && opt->arg[1] != 'n') { /* We specified the three parameters */
251 if (sscanf (opt->arg, "%[^/]/%[^/]/%s", txt_a, txt_b, txt_c) != 3) {
252 GMT_Report (API, GMT_MSG_ERROR, "Option -Q: Could not decode arrowwidth/headlength/headwidth\n");
253 n_errors++;
254 }
255 else { /* Turn the old args into new +a<angle> and pen width */
256 Ctrl->Q.S.v.v_width = (float)gmt_M_to_inch (GMT, txt_a);
257 Ctrl->Q.S.v.pen.width = gmt_M_to_points (GMT, txt_a);
258 Ctrl->Q.S.v.h_length = (float)gmt_M_to_inch (GMT, txt_b);
259 Ctrl->Q.S.v.h_width = (float)gmt_M_to_inch (GMT, txt_c);
260 }
261 }
262 if (Ctrl->Q.S.v.v_norm > 0.0) opt->arg[j] = 'n'; /* Restore the n<norm> string */
263 Ctrl->Q.S.v.status |= (PSL_VEC_JUST_B + PSL_VEC_FILL); /* Start filled vector at node location */
264 Ctrl->Q.S.symbol = GMT_SYMBOL_VECTOR_V4;
265 }
266 else {
267 if (opt->arg[0] == '+') { /* No size (use default), just attributes */
268 Ctrl->Q.S.size_x = VECTOR_HEAD_LENGTH * GMT->session.u2u[GMT_PT][GMT_INCH]; /* 9p */
269 n_errors += gmt_parse_vector (GMT, symbol, opt->arg, &Ctrl->Q.S);
270 }
271 else { /* Size, plus possible attributes */
272 j = sscanf (opt->arg, "%[^+]%s", txt_a, txt_b); /* txt_a should be symbols size with any +<modifiers> in txt_b */
273 if (j == 1) txt_b[0] = 0; /* No modifiers present, set txt_b to empty */
274 Ctrl->Q.S.size_x = gmt_M_to_inch (GMT, txt_a); /* Length of vector */
275 n_errors += gmt_parse_vector (GMT, symbol, txt_b, &Ctrl->Q.S);
276 }
277 /* Must possibly change v_norm to inches if given in another Cartesian unit */
278 if (Ctrl->Q.S.u_set && Ctrl->Q.S.u != GMT_INCH) {
279 Ctrl->Q.S.v.v_norm *= GMT->session.u2u[Ctrl->Q.S.u][GMT_INCH]; /* Since we are not reading this again we change to inches */
280 Ctrl->Q.S.u = GMT_INCH;
281 }
282 }
283 break;
284 case 'S': /* Scale -S[l|i]<length|scale>[<unit>] */
285 n_errors += gmt_M_repeated_module_option (API, Ctrl->S.active);
286 Ctrl->S.active = true;
287 len = strlen (opt->arg) - 1; /* Location of expected unit */
288 j = (opt->arg[0] == 'i') ? 1 : 0; /* j = 1 if -Si was selected */
289 if (strchr (GMT_DIM_UNITS GMT_LEN_UNITS, (int)opt->arg[len])) /* Recognized plot length or map distance unit character */
290 Ctrl->S.unit = opt->arg[len];
291 else if (! (opt->arg[len] == '.' || isdigit ((int)opt->arg[len]))) { /* Not decimal point or digit at the end means trouble */
292 GMT_Report (API, GMT_MSG_ERROR, "Option -S: Unrecognized length unit %c\n", opt->arg[len]);
293 n_errors++;
294 }
295 if (opt->arg[0] == 'l') { /* Want a fixed length for all vectors (ignore magnitudes) */
296 Ctrl->S.constant = true;
297 Ctrl->S.factor = atof (&opt->arg[1]);
298 }
299 else /* Get the length or scale */
300 Ctrl->S.factor = atof (&opt->arg[j]);
301 if (j == 1) Ctrl->S.invert = true;
302 break;
303 case 'T': /* Rescale Cartesian angles */
304 n_errors += gmt_M_repeated_module_option (API, Ctrl->T.active);
305 Ctrl->T.active = true;
306 break;
307 case 'W': /* Set line attributes */
308 n_errors += gmt_M_repeated_module_option (API, Ctrl->W.active);
309 Ctrl->W.active = true;
310 if (gmt_getpen (GMT, opt->arg, &Ctrl->W.pen)) {
311 gmt_pen_syntax (GMT, 'W', NULL, " ", NULL, 0);
312 n_errors++;
313 }
314 if (Ctrl->W.pen.cptmode) Ctrl->W.cpt_effect = true;
315 break;
316 case 'Z': /* Azimuths given */
317 n_errors += gmt_M_repeated_module_option (API, Ctrl->Z.active);
318 Ctrl->A.active = true;
319 Ctrl->Z.active = true;
320 break;
321
322 default: /* Report bad options */
323 n_errors += gmt_default_error (GMT, opt->option);
324 break;
325 }
326 }
327
328 if (!Ctrl->W.active) { /* Accept -W default pen for stem */
329 GMT_Report (API, GMT_MSG_DEBUG, "Option -W: Not given so we accept default pen\n");
330 Ctrl->W.active = true;
331 }
332 if (!Ctrl->G.active && (Ctrl->Q.S.v.status & PSL_VEC_FILL2)) { /* Gave fill via +g instead */
333 GMT_Report (API, GMT_MSG_DEBUG, "Option -G: Not given but -Q+g was set so we use it to fill head\n");
334 gmt_M_rgb_copy (Ctrl->G.fill.rgb, Ctrl->Q.S.v.fill.rgb);
335 Ctrl->G.active = true;
336 }
337 gmt_consider_current_cpt (API, &Ctrl->C.active, &(Ctrl->C.file));
338
339 n_errors += gmt_M_check_condition (GMT, !GMT->common.J.active, "Must specify a map projection with the -J option\n");
340 n_errors += gmt_M_check_condition (GMT, GMT->common.R.active[ISET] && (GMT->common.R.inc[GMT_X] <= 0.0 || GMT->common.R.inc[GMT_Y] <= 0.0),
341 "Option -I: Must specify positive increments\n");
342 n_errors += gmt_M_check_condition (GMT, Ctrl->S.factor == 0.0 && !Ctrl->S.constant, "Option -S: Scale must be nonzero\n");
343 n_errors += gmt_M_check_condition (GMT, Ctrl->S.factor <= 0.0 && Ctrl->S.constant, "Option -Sl: Length must be positive\n");
344 n_errors += gmt_M_check_condition (GMT, Ctrl->S.constant && Ctrl->Q.S.v.v_norm > 0.0,
345 "Option -Sl, -Q options: Cannot use -Q..n<size> with -Sl\n");
346 n_errors += gmt_M_check_condition (GMT, !(Ctrl->G.active || Ctrl->W.active || Ctrl->C.active),
347 "Must specify at least one of -G, -W, -C\n");
348 n_errors += gmt_M_check_condition (GMT, n_files != 2, "Must specify two input grid files\n");
349 n_errors += gmt_M_check_condition (GMT, Ctrl->W.cpt_effect && !Ctrl->C.active, "Option -W: modifier +c only makes sense if -C is given\n");
350
351 return (n_errors ? GMT_PARSE_ERROR : GMT_NOERROR);
352 }
353
354 #define bailout(code) {gmt_M_free_options (mode); return (code);}
355 #define Return(code) {Free_Ctrl (GMT, Ctrl); gmt_end_module (GMT, GMT_cpy); bailout (code);}
356
GMT_grdvector(void * V_API,int mode,void * args)357 EXTERN_MSC int GMT_grdvector (void *V_API, int mode, void *args) {
358 unsigned int justify, row, col, col_0, row_0, d_col, d_row, k, n_warn[3] = {0, 0, 0}, warn;
359 int error = 0;
360 bool Geographic;
361
362 uint64_t ij;
363
364 double tmp, x, y, plot_x, plot_y, x_off, y_off, f, headpen_width = 0.0;
365 double x2, y2, wesn[4], value, vec_data_length, vec_azim, scaled_vec_length, c, s, dim[PSL_MAX_DIMS];
366
367 struct GMT_GRID *Grid[2] = {NULL, NULL};
368 struct GMT_PALETTE *P = NULL;
369 struct GMT_PEN last_headpen;
370 struct GRDVECTOR_CTRL *Ctrl = NULL;
371 struct GMT_CTRL *GMT = NULL, *GMT_cpy = NULL; /* General GMT internal parameters */
372 struct GMT_OPTION *options = NULL;
373 struct PSL_CTRL *PSL = NULL; /* General PSL internal parameters */
374 struct GMTAPI_CTRL *API = gmt_get_api_ptr (V_API); /* Cast from void to GMTAPI_CTRL pointer */
375
376 /*----------------------- Standard module initialization and parsing ----------------------*/
377
378 if (API == NULL) return (GMT_NOT_A_SESSION);
379 if (mode == GMT_MODULE_PURPOSE) return (usage (API, GMT_MODULE_PURPOSE)); /* Return the purpose of program */
380 options = GMT_Create_Options (API, mode, args); if (API->error) return (API->error); /* Set or get option list */
381
382 if ((error = gmt_report_usage (API, options, 0, usage)) != GMT_NOERROR) bailout (error); /* Give usage if requested */
383
384 /* Parse the command-line arguments */
385
386 if ((GMT = gmt_init_module (API, THIS_MODULE_LIB, THIS_MODULE_CLASSIC_NAME, THIS_MODULE_KEYS, THIS_MODULE_NEEDS, NULL, &options, &GMT_cpy)) == NULL) bailout (API->error); /* Save current state */
387 if (GMT_Parse_Common (API, THIS_MODULE_OPTIONS, options)) Return (API->error);
388 Ctrl = New_Ctrl (GMT); /* Allocate and initialize a new control structure */
389 if ((error = parse (GMT, Ctrl, options)) != 0) Return (error);
390
391 /*---------------------------- This is the grdvector main code ----------------------------*/
392
393 gmt_M_memset (&last_headpen, 1, struct GMT_PEN); /* Initilaize members to zero */
394 GMT_Report (API, GMT_MSG_INFORMATION, "Processing input grids\n");
395 d_col = d_row = 1;
396 col_0 = row_0 = 0;
397
398 if (!(strcmp (Ctrl->In.file[0], "=") || strcmp (Ctrl->In.file[1], "="))) {
399 GMT_Report (API, GMT_MSG_ERROR, "Piping of grid files not supported!\n");
400 Return (GMT_RUNTIME_ERROR);
401 }
402
403 for (k = 0; k < 2; k++) {
404 if ((Grid[k] = GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_ONLY, NULL, Ctrl->In.file[k], NULL)) == NULL) { /* Get header only */
405 Return (API->error);
406 }
407 if ((API->error = gmt_img_sanitycheck (GMT, Grid[k]->header))) { /* Used map projection on a Mercator (cartesian) grid */
408 Return (API->error);
409 }
410 gmt_grd_init (GMT, Grid[k]->header, options, true);
411 }
412
413 if (!(gmt_M_grd_same_shape (GMT, Grid[0], Grid[1]) && gmt_M_grd_same_region (GMT, Grid[0], Grid[1]) && gmt_M_grd_same_inc (GMT, Grid[0], Grid[1]))) {
414 GMT_Report (API, GMT_MSG_ERROR, "files %s and %s does not match!\n", Ctrl->In.file[0], Ctrl->In.file[1]);
415 Return (GMT_RUNTIME_ERROR);
416 }
417
418 /* Determine what wesn to pass to map_setup */
419
420 if (!GMT->common.R.active[RSET]) /* -R was not set so we use the grid domain */
421 gmt_set_R_from_grd (GMT, Grid[0]->header);
422
423 if (gmt_map_setup (GMT, GMT->common.R.wesn)) Return (GMT_PROJECTION_ERROR);
424
425 /* Determine the wesn to be used to read the grid file */
426
427 if (!gmt_grd_setregion (GMT, Grid[0]->header, wesn, BCR_BILINEAR) || !gmt_grd_setregion (GMT, Grid[1]->header, wesn, BCR_BILINEAR)) {
428 /* No grid to plot; just do empty map and return */
429 if (GMT_End_IO (API, GMT_IN, 0) != GMT_NOERROR) { /* Disables further data input */
430 Return (API->error);
431 }
432 GMT_Report (API, GMT_MSG_WARNING, "No data within specified region\n");
433 if ((PSL = gmt_plotinit (GMT, options)) == NULL) Return (GMT_RUNTIME_ERROR);
434 gmt_plane_perspective (GMT, GMT->current.proj.z_project.view_plane, GMT->current.proj.z_level);
435 gmt_set_basemap_orders (GMT, GMT_BASEMAP_FRAME_AFTER, GMT_BASEMAP_GRID_AFTER, GMT_BASEMAP_ANNOT_AFTER);
436 GMT->current.map.frame.order = GMT_BASEMAP_AFTER; /* Move to last order since only calling gmt_map_basemap once */
437 gmt_plotcanvas (GMT); /* Fill canvas if requested */
438 gmt_map_basemap (GMT);
439 gmt_plane_perspective (GMT, -1, 0.0);
440 gmt_plotend (GMT);
441 Return (GMT_NOERROR);
442 }
443
444 /* Read data */
445
446 for (k = 0; k < 2; k++) if (GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_DATA_ONLY, wesn, Ctrl->In.file[k], Grid[k]) == NULL) { /* Get data */
447 Return (API->error);
448 }
449
450 if (Ctrl->C.active) {
451 double v_data_min, v_data_max;
452 if (Ctrl->A.active) { /* Polar grid, just use min/max of radius grid */
453 v_data_min = Grid[0]->header->z_min;
454 v_data_max = Grid[0]->header->z_max;
455 }
456 else { /* Find min/max vector lengths from the components */
457 v_data_min = DBL_MAX; v_data_max = 0.0;
458 gmt_M_grd_loop (GMT, Grid[GMT_X], row, col, ij) {
459 vec_data_length = hypot (Grid[GMT_X]->data[ij], Grid[GMT_Y]->data[ij]);
460 if (vec_data_length < v_data_min) v_data_min = vec_data_length;
461 if (vec_data_length > v_data_max) v_data_max = vec_data_length;
462 }
463 }
464 if ((P = gmt_get_palette (GMT, Ctrl->C.file, GMT_CPT_OPTIONAL, v_data_min, v_data_max, Ctrl->C.dz)) == NULL) {
465 Return (API->error);
466 }
467 }
468
469 Geographic = (gmt_M_is_geographic (GMT, GMT_IN)); /* Will be overridden if c|i|p units for scaling is selected */
470
471 if (!Ctrl->S.invert) {
472 double was = Ctrl->S.factor;
473 Ctrl->S.factor = 1.0 / Ctrl->S.factor; /* Turn length into a scale */
474 GMT_Report (API, GMT_MSG_INFORMATION, "Vector scale of %g <data-unit>/%c converts to %g %c/<data-unit>.\n", was, Ctrl->S.unit, Ctrl->S.factor, Ctrl->S.unit);
475 }
476
477 switch (Ctrl->S.unit) { /* Adjust for possible unit selection */
478 /* First three choices will give straight vectors scaled from user length to plot lengths */
479 case 'c':
480 Ctrl->S.factor *= GMT->session.u2u[GMT_CM][GMT_INCH];
481 Geographic = false;
482 break;
483 case 'i':
484 Ctrl->S.factor *= GMT->session.u2u[GMT_INCH][GMT_INCH];
485 Geographic = false;
486 break;
487 case 'p':
488 Ctrl->S.factor *= GMT->session.u2u[GMT_PT][GMT_INCH];
489 Geographic = false;
490 break;
491 /* Remaining choices will give geo vectors scaled from user length to distance lengths [dmsefkMnu] */
492 case 'd':
493 Ctrl->S.factor *= GMT->current.proj.KM_PR_DEG;
494 break;
495 case 'm':
496 Ctrl->S.factor *= (GMT->current.proj.KM_PR_DEG / GMT_DEG2MIN_F);
497 break;
498 case 's':
499 Ctrl->S.factor *= (GMT->current.proj.KM_PR_DEG / GMT_DEG2SEC_F);
500 break;
501 case 'e':
502 Ctrl->S.factor *= (1 / METERS_IN_A_KM); /* 1 is "meters in a meter" ... */
503 break;
504 case 'f':
505 Ctrl->S.factor *= (METERS_IN_A_FOOT / METERS_IN_A_KM);
506 break;
507 case 'k':
508 /* Already in km */
509 break;
510 case 'M':
511 Ctrl->S.factor *= (METERS_IN_A_MILE / METERS_IN_A_KM);
512 break;
513 case 'n':
514 Ctrl->S.factor *= (METERS_IN_A_NAUTICAL_MILE / METERS_IN_A_KM);
515 break;
516 case 'u':
517 Ctrl->S.factor *= (METERS_IN_A_SURVEY_FOOT / METERS_IN_A_KM);
518 break;
519 default:
520 GMT_Report (API, GMT_MSG_ERROR, "Bad scale unit %c\n", Ctrl->S.unit);
521 Return (GMT_RUNTIME_ERROR);
522 break;
523 }
524
525 if (Geographic) { /* Now that we know this we make sure -T is disabled if given */
526 if (Ctrl->T.active) { /* This is a mistake */
527 Ctrl->T.active = false;
528 GMT_Report (API, GMT_MSG_ERROR, "-T does not apply to geographic grids - ignored\n");
529 }
530 GMT_Report (API, GMT_MSG_DEBUG, "Great-circle geo-vectors will be drawn. Scale converting user lengths to km is %g\n", Ctrl->S.factor);
531 }
532 else {
533 GMT_Report (API, GMT_MSG_DEBUG, "Cartesian straight vectors will be drawn\n");
534 GMT_Report (API, GMT_MSG_DEBUG, "Cartesian straight vectors will be drawn. Scale converting user lengths to inches is %g\n", Ctrl->S.factor);
535 }
536
537 if (Ctrl->Q.active) { /* Prepare vector parameters */
538 if (Ctrl->Q.S.symbol == PSL_VECTOR) Ctrl->Q.S.v.v_width = (float)(Ctrl->W.pen.width * GMT->session.u2u[GMT_PT][GMT_INCH]);
539 gmt_init_vector_param (GMT, &Ctrl->Q.S, true, Ctrl->W.active, &Ctrl->W.pen, Ctrl->G.active, &Ctrl->G.fill);
540 }
541 if ((PSL = gmt_plotinit (GMT, options)) == NULL) Return (GMT_RUNTIME_ERROR);
542 gmt_plane_perspective (GMT, GMT->current.proj.z_project.view_plane, GMT->current.proj.z_level);
543 gmt_set_basemap_orders (GMT, Ctrl->N.active ? GMT_BASEMAP_FRAME_BEFORE : GMT_BASEMAP_FRAME_AFTER, GMT_BASEMAP_GRID_BEFORE, GMT_BASEMAP_ANNOT_AFTER);
544 gmt_plotcanvas (GMT); /* Fill canvas if requested */
545 gmt_map_basemap (GMT);
546
547 gmt_setpen (GMT, &Ctrl->W.pen);
548 if (!Ctrl->C.active) gmt_setfill (GMT, &Ctrl->G.fill, Ctrl->W.active);
549
550 if (!Ctrl->N.active) gmt_map_clip_on (GMT, GMT->session.no_rgb, 3);
551 if (Ctrl->I.mode) { /* Gave multiplier so get actual strides */
552 GMT->common.R.inc[GMT_X] *= Grid[0]->header->inc[GMT_X];
553 GMT->common.R.inc[GMT_Y] *= Grid[0]->header->inc[GMT_Y];
554 }
555 if (GMT->common.R.inc[GMT_X] != 0.0 && GMT->common.R.inc[GMT_Y] != 0.0) { /* Coarsen the output interval. The new -Idx/dy must be integer multiples of the grid dx/dy */
556 struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (Grid[0]->header);
557 double val = GMT->common.R.inc[GMT_Y] * HH->r_inc[GMT_Y]; /* Should be ~ an integer within 1 ppm */
558 d_row = urint (val);
559 if (d_row == 0 || fabs ((d_row - val)/d_row) > GMT_CONV6_LIMIT) {
560 GMT_Report (API, GMT_MSG_ERROR, "New y grid spacing (%.12lg) is not a multiple of actual grid spacing (%.12g) [within %g]\n", GMT->common.R.inc[GMT_Y], Grid[0]->header->inc[GMT_Y], GMT_CONV6_LIMIT);
561 Return (GMT_RUNTIME_ERROR);
562 }
563 GMT->common.R.inc[GMT_Y] = d_row * Grid[0]->header->inc[GMT_Y]; /* Get exact y-increment in case of slop */
564 val = GMT->common.R.inc[GMT_X] * HH->r_inc[GMT_X];
565 d_col = urint (val);
566 if (d_col == 0 || fabs ((d_col - val)/d_col) > GMT_CONV6_LIMIT) {
567 GMT_Report (API, GMT_MSG_ERROR, "New x grid spacing (%.12g) is not a multiple of actual grid spacing (%.12g) [within %g]\n", GMT->common.R.inc[GMT_X], Grid[0]->header->inc[GMT_X], GMT_CONV6_LIMIT);
568 Return (GMT_RUNTIME_ERROR);
569 }
570 GMT->common.R.inc[GMT_X] = d_col * Grid[0]->header->inc[GMT_X]; /* Get exact x-increment in case of slop */
571
572 /* Determine starting row/col for straddled access */
573 tmp = ceil (Grid[0]->header->wesn[YHI] / GMT->common.R.inc[GMT_Y]) * GMT->common.R.inc[GMT_Y];
574 if (tmp > Grid[0]->header->wesn[YHI]) tmp -= GMT->common.R.inc[GMT_Y];
575 row_0 = urint ((Grid[0]->header->wesn[YHI] - tmp) * HH->r_inc[GMT_Y]);
576 tmp = floor (Grid[0]->header->wesn[XLO] / GMT->common.R.inc[GMT_X]) * GMT->common.R.inc[GMT_X];
577 if (tmp < Grid[0]->header->wesn[XLO]) tmp += GMT->common.R.inc[GMT_X];
578 col_0 = urint ((tmp - Grid[0]->header->wesn[XLO]) * HH->r_inc[GMT_X]);
579 }
580
581 dim[PSL_VEC_HEAD_SHAPE] = Ctrl->Q.S.v.v_shape; /* head shape, type, and status do not change inside the loop */
582 dim[PSL_VEC_STATUS] = (double)Ctrl->Q.S.v.status;
583 dim[PSL_VEC_HEAD_TYPE_BEGIN] = (double)Ctrl->Q.S.v.v_kind[0];
584 dim[PSL_VEC_HEAD_TYPE_END] = (double)Ctrl->Q.S.v.v_kind[1];
585
586 if (Ctrl->Q.S.v.status & PSL_VEC_OUTLINE2) { /* Vector head outline pen specified separately */
587 PSL_defpen (PSL, "PSL_vecheadpen", Ctrl->Q.S.v.pen.width, Ctrl->Q.S.v.pen.style, Ctrl->Q.S.v.pen.offset, Ctrl->Q.S.v.pen.rgb);
588 headpen_width = Ctrl->Q.S.v.pen.width;
589 }
590 else { /* Reset to default pen */
591 if (Ctrl->W.active) { /* Vector head outline pen default is half that of stem pen */
592 PSL_defpen (PSL, "PSL_vecheadpen", Ctrl->W.pen.width, Ctrl->W.pen.style, Ctrl->W.pen.offset, Ctrl->W.pen.rgb);
593 headpen_width = 0.5 * Ctrl->W.pen.width;
594 }
595 }
596 if (Ctrl->W.cpt_effect) { /* Should color apply to pen, fill, or both [fill] */
597 if ((Ctrl->W.pen.cptmode & 2) == 0 && !Ctrl->G.active) /* Turn off CPT fill */
598 gmt_M_rgb_copy (Ctrl->G.fill.rgb, GMT->session.no_rgb);
599 }
600
601 if (gmt_M_is_verbose (GMT, GMT_MSG_INFORMATION)) { /* Report min/max/mean scaled vector length */
602 double v_scaled_min = DBL_MAX, v_scaled_max = -DBL_MAX, v_scaled_mean = 0.0;
603 double v_data_min = DBL_MAX, v_data_max = -DBL_MAX, v_data_mean = 0.0;
604 uint64_t v_n = 0;
605 char v_unit[GMT_LEN8] = {""};
606
607 for (row = row_0; row < Grid[1]->header->n_rows; row += d_row) {
608 y = gmt_M_grd_row_to_y (GMT, row, Grid[0]->header); /* Latitude OR y OR radius */
609 for (col = col_0; col < Grid[1]->header->n_columns; col += d_col) {
610
611 ij = gmt_M_ijp (Grid[0]->header, row, col);
612 if (gmt_M_is_fnan (Grid[0]->data[ij]) || gmt_M_is_fnan (Grid[1]->data[ij])) continue; /* Cannot plot NaN-vectors */
613 x = gmt_M_grd_col_to_x (GMT, col, Grid[0]->header); /* Longitude OR x OR theta [or azimuth] */
614 if (!Ctrl->N.active) { /* Throw out vectors whose node is outside */
615 gmt_map_outside (GMT, x, y);
616 if (abs (GMT->current.map.this_x_status) > 1 || abs (GMT->current.map.this_y_status) > 1) continue;
617 }
618
619 if (Ctrl->A.active) { /* Got r,theta grids */
620 vec_data_length = Grid[0]->data[ij];
621 if (vec_data_length == 0.0) continue; /* No length = no plotting */
622 if (vec_data_length < 0.0) /* Interpret negative lengths to mean pointing in opposite direction 180-degrees off */
623 vec_data_length = -vec_data_length;
624 }
625 else { /* Cartesian component grids: Convert to polar form of radius, theta */
626 vec_data_length = hypot (Grid[GMT_X]->data[ij], Grid[GMT_Y]->data[ij]);
627 if (vec_data_length == 0.0) continue; /* No length = no plotting */
628 }
629 scaled_vec_length = (Ctrl->S.constant) ? Ctrl->S.factor : vec_data_length * Ctrl->S.factor;
630 /* scaled_vec_length is now in inches (Cartesian) or km (Geographic) */
631 if (vec_data_length < v_data_min) v_data_min = vec_data_length;
632 if (vec_data_length > v_data_max) v_data_max = vec_data_length;
633 v_data_mean += vec_data_length;
634 if (scaled_vec_length < v_scaled_min) v_scaled_min = scaled_vec_length;
635 if (scaled_vec_length > v_scaled_max) v_scaled_max = scaled_vec_length;
636 v_scaled_mean += scaled_vec_length;
637 v_n++;
638 }
639 }
640 if (v_n) { /* Compute the means */
641 v_data_mean /= v_n;
642 v_scaled_mean /= v_n;
643 }
644 if (Geographic) /* Since regardless of unit chosen, we end up with a length in km */
645 sprintf (v_unit, "km");
646 else { /* Report length in selected unit and scale results to match */
647 strcpy (v_unit, API->GMT->session.unit_name[API->GMT->current.setting.proj_length_unit]);
648 v_scaled_min *= GMT->session.u2u[GMT_INCH][GMT->current.setting.proj_length_unit];
649 v_scaled_max *= GMT->session.u2u[GMT_INCH][GMT->current.setting.proj_length_unit];
650 v_scaled_mean *= GMT->session.u2u[GMT_INCH][GMT->current.setting.proj_length_unit];
651 }
652
653 GMT_Report (API, GMT_MSG_INFORMATION, "Minimum length of data vector (user unit) : %g\n", v_data_min);
654 GMT_Report (API, GMT_MSG_INFORMATION, "Maximum length of data vector (user unit) : %g\n", v_data_max);
655 GMT_Report (API, GMT_MSG_INFORMATION, "Mean length of the data vector (user unit) : %g\n", v_data_mean);
656
657 if (!Ctrl->S.constant) { /* No point reporting the mean of n identical vectors */
658 GMT_Report (API, GMT_MSG_INFORMATION, "Minimum length of scaled vector in %4s : %g\n", v_unit, v_scaled_min);
659 GMT_Report (API, GMT_MSG_INFORMATION, "Maximum length of scaled vector in %4s : %g\n", v_unit, v_scaled_max);
660 GMT_Report (API, GMT_MSG_INFORMATION, "Mean length of the scaled vector in %4s : %g\n", v_unit, v_scaled_mean);
661 }
662 }
663
664 PSL_command (GMT->PSL, "V\n");
665 for (row = row_0; row < Grid[1]->header->n_rows; row += d_row) {
666 y = gmt_M_grd_row_to_y (GMT, row, Grid[0]->header); /* Latitude OR y OR radius */
667 for (col = col_0; col < Grid[1]->header->n_columns; col += d_col) {
668
669 ij = gmt_M_ijp (Grid[0]->header, row, col);
670 if (gmt_M_is_fnan (Grid[0]->data[ij]) || gmt_M_is_fnan (Grid[1]->data[ij])) continue; /* Cannot plot NaN-vectors */
671 x = gmt_M_grd_col_to_x (GMT, col, Grid[0]->header); /* Longitude OR x OR theta [or azimuth] */
672 if (!Ctrl->N.active) { /* Throw out vectors whose node is outside */
673 gmt_map_outside (GMT, x, y);
674 if (abs (GMT->current.map.this_x_status) > 1 || abs (GMT->current.map.this_y_status) > 1) continue;
675 }
676
677 if (Ctrl->A.active) { /* Got r,theta grids */
678 vec_data_length = Grid[0]->data[ij];
679 if (vec_data_length == 0.0) continue; /* No length = no plotting */
680 vec_azim = Grid[1]->data[ij];
681 value = vec_data_length;
682 if (vec_data_length < 0.0) { /* Interpret negative lengths to mean pointing in opposite direction 180-degrees off */
683 vec_data_length = -vec_data_length;
684 vec_azim += 180.0;
685 }
686 if (!Ctrl->Z.active) vec_azim = 90.0 - vec_azim; /* Convert theta to azimuth */
687 }
688 else { /* Cartesian component grids: Convert to polar form of radius, theta */
689 vec_data_length = hypot (Grid[GMT_X]->data[ij], Grid[GMT_Y]->data[ij]);
690 if (vec_data_length == 0.0) continue; /* No length = no plotting */
691 vec_azim = 90.0 - atan2d (Grid[GMT_Y]->data[ij], Grid[GMT_X]->data[ij]); /* Convert dy,dx to azimuth */
692 value = vec_data_length;
693 }
694
695 if (Ctrl->C.active) { /* Get color based on the vector length */
696 gmt_get_fill_from_z (GMT, P, value, &Ctrl->G.fill);
697 }
698
699 if (Ctrl->W.cpt_effect) { /* Should color apply to pen, fill, or both [fill] */
700 if (Ctrl->W.pen.cptmode & 1) { /* Change pen color via CPT */
701 gmt_M_rgb_copy (Ctrl->W.pen.rgb, Ctrl->G.fill.rgb);
702 if (!gmt_M_same_pen (Ctrl->W.pen, last_headpen)) { /* Since color may have changed */
703 PSL_defpen (PSL, "PSL_vecheadpen", Ctrl->W.pen.width, Ctrl->W.pen.style, Ctrl->W.pen.offset, Ctrl->W.pen.rgb);
704 last_headpen = Ctrl->W.pen;
705 }
706 }
707 }
708 if (Ctrl->C.active) { /* Update pen and fill color settings */
709 if (!Ctrl->Q.active) /* Must update stick pen */
710 gmt_M_rgb_copy (Ctrl->W.pen.rgb, Ctrl->G.fill.rgb);
711 gmt_setpen (GMT, &Ctrl->W.pen);
712 if (Ctrl->Q.active) gmt_setfill (GMT, &Ctrl->G.fill, Ctrl->W.active);
713 gmt_init_vector_param (GMT, &Ctrl->Q.S, true, Ctrl->W.active, &Ctrl->W.pen, true, &Ctrl->G.fill);
714 }
715
716 scaled_vec_length = (Ctrl->S.constant) ? Ctrl->S.factor : vec_data_length * Ctrl->S.factor;
717 /* scaled_vec_length is now in inches (Cartesian) or km (Geographic) */
718
719 if (Geographic) { /* Draw great-circle geo-vectors */
720 warn = gmt_geo_vector (GMT, x, y, vec_azim, scaled_vec_length, &Ctrl->W.pen, &Ctrl->Q.S);
721 n_warn[warn]++;
722 }
723 else { /* Draw straight Cartesian vectors */
724 gmt_geo_to_xy (GMT, x, y, &plot_x, &plot_y);
725 if (gmt_M_is_geographic (GMT, GMT_IN)) /* Must align azimuth with local north */
726 vec_azim = 90.0 - gmt_azim_to_angle (GMT, x, y, 0.1, vec_azim);
727 if (Ctrl->T.active) /* Deal with negative scales in x and/or y which affect the azimuths */
728 gmt_flip_azim_d (GMT, &vec_azim);
729 vec_azim = 90.0 - vec_azim; /* Transform azimuths to plot angle */
730 if (GMT->current.proj.projection_GMT == GMT_POLAR) { /* Must rotate azimuth since circular projection */
731 double x_orient;
732 x_orient = (GMT->current.proj.got_azimuths) ? -(x + GMT->current.proj.p_base_angle) : x - GMT->current.proj.p_base_angle - 90.0;
733 vec_azim += x_orient;
734 }
735 vec_azim *= D2R; /* vec_azim is now in radians */
736 sincos (vec_azim, &s, &c);
737 x2 = plot_x + scaled_vec_length * c;
738 y2 = plot_y + scaled_vec_length * s;
739
740 justify = PSL_vec_justify (Ctrl->Q.S.v.status); /* Return justification as 0-2 */
741 if (justify) { /* Justify vector at center, or tip [beginning] */
742 x_off = justify * 0.5 * (x2 - plot_x); y_off = justify * 0.5 * (y2 - plot_y);
743 plot_x -= x_off; plot_y -= y_off;
744 x2 -= x_off; y2 -= y_off;
745 }
746 n_warn[0]++;
747 if (!Ctrl->Q.active) { /* Just a vector stem: line segment */
748 PSL_plotsegment (PSL, plot_x, plot_y, x2, y2);
749 continue;
750 }
751 /* Must plot a vector head */
752 dim[PSL_VEC_XTIP] = x2;
753 dim[PSL_VEC_YTIP] = y2;
754 dim[PSL_VEC_TAIL_WIDTH] = Ctrl->Q.S.v.v_width;
755 dim[PSL_VEC_HEAD_LENGTH] = Ctrl->Q.S.v.h_length;
756 dim[PSL_VEC_HEAD_WIDTH] = Ctrl->Q.S.v.h_width;
757 dim[PSL_VEC_HEAD_PENWIDTH] = headpen_width; /* Possibly shrunk head pen width */
758 if (scaled_vec_length < Ctrl->Q.S.v.v_norm) { /* Scale arrow attributes down with length */
759 f = scaled_vec_length / Ctrl->Q.S.v.v_norm;
760 for (k = 2; k <= 4; k++) dim[k] *= f;
761 dim[PSL_VEC_HEAD_PENWIDTH] *= f;
762 }
763 if (Ctrl->Q.S.symbol == GMT_SYMBOL_VECTOR_V4) { /* Do the deprecated GMT4 vector polygon instead */
764 int v4_outline = Ctrl->W.active;
765 double *this_rgb = NULL;
766 if (Ctrl->G.active || Ctrl->C.active)
767 this_rgb = Ctrl->G.fill.rgb;
768 else
769 this_rgb = GMT->session.no_rgb;
770 if (v4_outline) gmt_setpen (GMT, &Ctrl->W.pen);
771 if (Ctrl->Q.S.v.status & PSL_VEC_BEGIN) v4_outline += 8; /* Double-headed */
772 psl_vector_v4 (PSL, plot_x, plot_y, dim, this_rgb, v4_outline);
773 }
774 else
775 PSL_plotsymbol (PSL, plot_x, plot_y, dim, PSL_VECTOR);
776 }
777 }
778 }
779 PSL_command (GMT->PSL, "U\n");
780 PSL->current.linewidth = 0.0; /* Since we changed things under clip; this will force it to be set next */
781
782 if (!Ctrl->N.active) gmt_map_clip_off (GMT);
783
784 gmt_map_basemap (GMT);
785
786 gmt_plane_perspective (GMT, -1, 0.0);
787
788 gmt_plotend (GMT);
789
790 GMT_Report (API, GMT_MSG_INFORMATION, "%d vectors plotted successfully\n", n_warn[0]);
791 if (n_warn[1]) GMT_Report (API, GMT_MSG_INFORMATION, "%d vector heads had length exceeding the vector length and were skipped. Consider the +n<norm> modifier to -Q\n", n_warn[1]);
792 if (n_warn[2]) GMT_Report (API, GMT_MSG_INFORMATION, "%d vector heads had to be scaled more than implied by +n<norm> since they were still too long. Consider changing the +n<norm> modifier to -Q\n", n_warn[2]);
793
794
795 Return (GMT_NOERROR);
796 }
797