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 * This program is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU Lesser General Public License as published by
7 * the Free Software Foundation; version 3 or any later version.
8 *
9 * This program is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 * GNU Lesser General Public License for more details.
13 *
14 * Contact info: www.generic-mapping-tools.org
15 *--------------------------------------------------------------------*/
16 /*
17 * Brief synopsis: gmtvector performs basic vector operations such as dot and
18 * cross products, sums, and conversion between Cartesian and polar/spherical
19 * coordinates systems.
20 *
21 * Author: Paul Wessel
22 * Date: 10-Aug-2010
23 * Version: 6 API
24 */
25
26 #include "gmt_dev.h"
27
28 #define THIS_MODULE_CLASSIC_NAME "gmtvector"
29 #define THIS_MODULE_MODERN_NAME "gmtvector"
30 #define THIS_MODULE_LIB "core"
31 #define THIS_MODULE_PURPOSE "Operations on Cartesian vectors in 2-D and 3-D"
32 #define THIS_MODULE_KEYS "<D{,>D}"
33 #define THIS_MODULE_NEEDS ""
34 #define THIS_MODULE_OPTIONS "-:>Vbdefghijoqs" GMT_OPT("HMm")
35
36 enum gmtvector_method { /* The available methods */
37 DO_NOTHING=0,
38 DO_AVERAGE,
39 DO_DOT2D,
40 DO_DOT3D,
41 DO_CROSS,
42 DO_SUM,
43 DO_ROT2D,
44 DO_ROT3D,
45 DO_ROTVAR2D,
46 DO_ROTVAR3D,
47 DO_DOT,
48 DO_POLE,
49 DO_TRANSLATE,
50 DO_BISECTOR};
51
52 struct GMTVECTOR_CTRL {
53 struct GMTVECTOR_Out { /* -> */
54 bool active;
55 char *file;
56 } Out;
57 struct GMTVECTOR_In { /* infile */
58 bool active;
59 unsigned int n_args;
60 char *arg;
61 } In;
62 struct GMTVECTOR_A { /* -A[m[<conf>]|<vec>] */
63 bool active;
64 unsigned int mode;
65 double conf;
66 char *arg;
67 } A;
68 struct GMTVECTOR_C { /* -C[i|o] */
69 bool active[2];
70 } C;
71 struct GMTVECTOR_E { /* -E */
72 bool active;
73 } E;
74 struct GMTVECTOR_N { /* -N */
75 bool active;
76 } N;
77 struct GMTVECTOR_S { /* -S[vec] */
78 bool active;
79 char *arg;
80 } S;
81 struct GMTVECTOR_T { /* -T[operator][<arg>] */
82 char unit;
83 bool active;
84 bool degree;
85 bool a_and_d;
86 enum gmtvector_method mode;
87 unsigned int dmode;
88 double par[3];
89 } T;
90 };
91
New_Ctrl(struct GMT_CTRL * GMT)92 static void *New_Ctrl (struct GMT_CTRL *GMT) { /* Allocate and initialize a new control structure */
93 struct GMTVECTOR_CTRL *C;
94
95 C = gmt_M_memory (GMT, NULL, 1, struct GMTVECTOR_CTRL);
96
97 /* Initialize values whose defaults are not 0/false/NULL */
98 C->A.conf = 0.95; /* 95% conf level */
99 C->T.unit = 'e'; /* Unit is meter unless specified */
100 return (C);
101 }
102
Free_Ctrl(struct GMT_CTRL * GMT,struct GMTVECTOR_CTRL * C)103 static void Free_Ctrl (struct GMT_CTRL *GMT, struct GMTVECTOR_CTRL *C) { /* Deallocate control structure */
104 if (!C) return;
105 gmt_M_str_free (C->In.arg);
106 gmt_M_str_free (C->Out.file);
107 gmt_M_str_free (C->A.arg);
108 gmt_M_str_free (C->S.arg);
109 gmt_M_free (GMT, C);
110 }
111
usage(struct GMTAPI_CTRL * API,int level)112 static int usage (struct GMTAPI_CTRL *API, int level) {
113 const char *name = gmt_show_name_and_purpose (API, THIS_MODULE_LIB, THIS_MODULE_CLASSIC_NAME, THIS_MODULE_PURPOSE);
114 if (level == GMT_MODULE_PURPOSE) return (GMT_NOERROR);
115 GMT_Usage (API, 0, "usage: %s [<table>] [-A[m][<conf>]|<vector>] [-C[i|o]] [-E] [-N] [-S<vector>] "
116 "[-Ta|b|d|D|p<azim>|r<rot>|R|s|t[<azim>/<dist>]|x] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s]\n",
117 name, GMT_V_OPT, GMT_b_OPT, GMT_d_OPT, GMT_e_OPT, GMT_f_OPT, GMT_g_OPT, GMT_h_OPT, GMT_i_OPT, GMT_o_OPT, GMT_q_OPT, GMT_s_OPT, GMT_colon_OPT, GMT_PAR_OPT);
118
119 if (level == GMT_SYNOPSIS) return (GMT_MODULE_SYNOPSIS);
120
121 GMT_Message (API, GMT_TIME_NONE, " OPTIONAL ARGUMENTS:\n");
122 GMT_Usage (API, 1, "\n<table>");
123 GMT_Usage (API, -2, "<table> is one or more data files (in ASCII, binary, netCDF) with (x,y[,z]), (r,theta) or (lon,lat) in the "
124 "first 2-3 input columns. If one item is given and it cannot be opened we will interpret it as a single x/y[/z], r/theta, or lon/lat entry. "
125 "If no file(s) is given, standard input is read.");
126 GMT_Usage (API, 1, "\n-A[m][<conf>]|<vector>");
127 GMT_Usage (API, -2, "Single primary vector, given as lon/lat, r/theta, or x/y[/z]. No tables will be read. "
128 "Alternatively, give -Am to compute a single primary vector as the mean of the input vectors. "
129 "The confidence ellipse for the mean vector is determined (95%% level); "
130 "optionally append a different confidence level in percent.");
131 GMT_Usage (API, 1, "\n-C[i|o]");
132 GMT_Usage (API, -2, "Indicate Cartesian coordinates on input/output instead of lon,lat or r/theta. "
133 "Append i or o to only affect input or output coordinates.");
134 GMT_Usage (API, 1, "\n-E Automatically convert between geodetic and geocentric coordinates [no conversion]. "
135 "Ignored unless -fg is given.");
136 GMT_Usage (API, 1, "\n-N Normalize the transformed vectors (only affects -Co output).");
137 GMT_Usage (API, 1, "\n-S<vector>");
138 GMT_Usage (API, -2, "The secondary vector (if needed by -T), given as lon/lat, r/theta, or x/y[/z].");
139 GMT_Usage (API, 1, "\n-Ta|b|d|D|p<azim>|r<angle>|R|s|t[<azim>/<dist>]|x");
140 GMT_Usage (API, -2, "Specify the desired transformation of the input data:");
141 GMT_Usage (API, 3, "a: Average of the primary and secondary vector (see -S).");
142 GMT_Usage (API, 3, "b: Find the bisector great circle pole(s) for input and secondary vector (see -S).");
143 GMT_Usage (API, 3, "d: Compute dot-product(s) with secondary vector (see -S).");
144 GMT_Usage (API, 3, "D: Same as -Td, but returns the angles in degrees between the vectors.");
145 GMT_Usage (API, 3, "p: Find pole to great circle with <azim> azimuth trend at input vector location.");
146 GMT_Usage (API, 3, "s: Find the sum of the secondary vector (see -S) and the input vector(s).");
147 GMT_Usage (API, 3, "r: Rotate the input vectors. Depending on your input (2-D or 3-D), append "
148 "<angle> or <plon/plat/angle>, respectively, to define the rotation.");
149 GMT_Usage (API, 3, "R: As r, but assumes the input vectors/angles are different rotations and repeatedly "
150 "rotate the fixed secondary vector (see -S) using the input rotations.");
151 GMT_Usage (API, 3, "t: Translate input vectors to points a distance <dist> away in the azimuth <azim>. "
152 "Append <azim>/<dist> for a fixed set of azimuth and distance for all points, "
153 "otherwise we expect to read <azim>, <dist> from the input file; append a unit [e]");
154 GMT_Usage (API, 3, "x: Compute cross-product(s) with secondary vector (see -S).");
155 GMT_Option (API, "V,bi0");
156 if (gmt_M_showusage (API)) GMT_Usage (API, -2, "Default is 2 [or 3; see -C, -fg] input columns.");
157 GMT_Option (API, "bo,d,e,f,g,h,i,o,q,s,:,.");
158
159 return (GMT_MODULE_USAGE);
160 }
161
parse(struct GMT_CTRL * GMT,struct GMTVECTOR_CTRL * Ctrl,struct GMT_OPTION * options)162 static int parse (struct GMT_CTRL *GMT, struct GMTVECTOR_CTRL *Ctrl, struct GMT_OPTION *options) {
163
164 /* This parses the options provided to grdsample and sets parameters in CTRL.
165 * Any GMT common options will override values set previously by other commands.
166 * It also replaces any file names specified as input or output with the data ID
167 * returned when registering these sources/destinations with the API.
168 */
169
170 unsigned int n_in, n_errors = 0, n_files = 0;
171 int n;
172 char txt_a[GMT_LEN64] = {""}, txt_b[GMT_LEN64] = {""}, txt_c[GMT_LEN64] = {""};
173 struct GMT_OPTION *opt = NULL;
174 struct GMTAPI_CTRL *API = GMT->parent;
175
176 for (opt = options; opt; opt = opt->next) {
177 switch (opt->option) {
178
179 case '<': /* Input files or single point */
180 Ctrl->In.active = true;
181 if (Ctrl->In.n_args++ == 0) Ctrl->In.arg = strdup (opt->arg);
182 break;
183 case '>': /* Got named output file */
184 if (n_files++ == 0) Ctrl->Out.file = strdup (opt->arg);
185 break;
186
187 /* Processes program-specific parameters */
188
189 case 'A': /* Secondary vector */
190 n_errors += gmt_M_repeated_module_option (API, Ctrl->A.active);
191 Ctrl->A.active = true;
192 if (opt->arg[0] == 'm') {
193 Ctrl->A.mode = 1;
194 if (opt->arg[1]) Ctrl->A.conf = 0.01 * atof (&opt->arg[1]);
195 }
196 else
197 Ctrl->A.arg = strdup (opt->arg);
198 break;
199 case 'C': /* Cartesian coordinates on in|out */
200 if (opt->arg[0] == 'i') {
201 n_errors += gmt_M_repeated_module_option (API, Ctrl->C.active[GMT_IN]);
202 Ctrl->C.active[GMT_IN] = true;
203 }
204 else if (opt->arg[0] == 'o') {
205 n_errors += gmt_M_repeated_module_option (API, Ctrl->C.active[GMT_OUT]);
206 Ctrl->C.active[GMT_OUT] = true;
207 }
208 else if (opt->arg[0] == '\0') {
209 n_errors += gmt_M_repeated_module_option (API, Ctrl->C.active[GMT_IN]);
210 n_errors += gmt_M_repeated_module_option (API, Ctrl->C.active[GMT_OUT]);
211 Ctrl->C.active[GMT_IN] = Ctrl->C.active[GMT_OUT] = true;
212 }
213 else {
214 GMT_Report (API, GMT_MSG_ERROR, "Bad modifier given to -C (%s)\n", opt->arg);
215 n_errors++;
216 }
217 break;
218 case 'E': /* geodetic/geocentric conversion */
219 n_errors += gmt_M_repeated_module_option (API, Ctrl->E.active);
220 Ctrl->E.active = true;
221 break;
222 case 'N': /* Normalize vectors */
223 n_errors += gmt_M_repeated_module_option (API, Ctrl->N.active);
224 Ctrl->N.active = true;
225 break;
226 case 'T': /* Selects transformation */
227 n_errors += gmt_M_repeated_module_option (API, Ctrl->T.active);
228 Ctrl->T.active = true;
229 switch (opt->arg[0]) {
230 case 'a': /* Angle between vectors */
231 Ctrl->T.mode = DO_AVERAGE;
232 break;
233 case 'b': /* Pole of bisector great circle */
234 Ctrl->T.mode = DO_BISECTOR;
235 break;
236 case 'D': /* Angle between vectors */
237 Ctrl->T.degree = true;
238 /* Intentionally fall through - to 'd' */
239 case 'd': /* dot-product of two vectors */
240 Ctrl->T.mode = DO_DOT;
241 break;
242 case 'p': /* Pole of great circle */
243 Ctrl->T.mode = DO_POLE;
244 Ctrl->T.par[0] = atof (&opt->arg[1]);
245 break;
246 case 'r': /* Rotate vectors */
247 n = sscanf (&opt->arg[1], "%[^/]/%[^/]/%s", txt_a, txt_b, txt_c);
248 if (n == 1) { /* 2-D Cartesian rotation */
249 Ctrl->T.mode = DO_ROT2D;
250 Ctrl->T.par[2] = atof (txt_a);
251 }
252 else if (n == 3) { /* 3-D spherical rotation */
253 Ctrl->T.mode = DO_ROT3D;
254 n_errors += gmt_verify_expectations (GMT, gmt_M_type (GMT, GMT_IN, GMT_X), gmt_scanf_arg (GMT, txt_a, gmt_M_type (GMT, GMT_IN, GMT_X), false, &Ctrl->T.par[0]), txt_a);
255 n_errors += gmt_verify_expectations (GMT, gmt_M_type (GMT, GMT_IN, GMT_Y), gmt_scanf_arg (GMT, txt_b, gmt_M_type (GMT, GMT_IN, GMT_Y), false, &Ctrl->T.par[1]), txt_b);
256 Ctrl->T.par[2] = atof (txt_c);
257 }
258 else {
259 GMT_Report (API, GMT_MSG_ERROR, "Bad arguments given to -Tr (%s)\n", &opt->arg[1]);
260 n_errors++;
261 }
262 break;
263 case 'R': /* Rotate secondary vector using input rotations */
264 Ctrl->T.mode = DO_ROTVAR2D; /* Change to 3D later if that is what we are doing */
265 break;
266 case 's': /* Sum of vectors */
267 Ctrl->T.mode = DO_SUM;
268 break;
269 case 't': /* Translate vectors */
270 Ctrl->T.mode = DO_TRANSLATE;
271 if (opt->arg[1]) { /* Gave azimuth/distance[<unit>] or just <unit> */
272 if (strchr (opt->arg, '/')) { /* Gave a fixed azimuth/distance[<unit>] combination */
273 if ((n = sscanf (&opt->arg[1], "%lg/%s", &Ctrl->T.par[0], txt_a)) != 2) {
274 GMT_Report (API, GMT_MSG_ERROR, "Bad arguments given to -Tr (%s)\n", &opt->arg[1]);
275 n_errors++;
276 }
277 else
278 Ctrl->T.dmode = gmt_get_distance (GMT, txt_a, &(Ctrl->T.par[1]), &(Ctrl->T.unit));
279 }
280 else if (strchr (GMT_LEN_UNITS, opt->arg[1])) { /* Gave the unit of the data in column 3 */
281 Ctrl->T.unit = opt->arg[1];
282 Ctrl->T.a_and_d = true;
283 }
284 }
285 else /* No arg means use default meters */
286 Ctrl->T.a_and_d = true;
287 break;
288 case 'x': /* Cross-product between vectors */
289 Ctrl->T.mode = DO_CROSS;
290 break;
291 default:
292 GMT_Report (API, GMT_MSG_ERROR, "Unrecognized directive given to -T (%s)\n", opt->arg);
293 n_errors++;
294 break;
295 }
296 break;
297 case 'S': /* Secondary vector */
298 n_errors += gmt_M_repeated_module_option (API, Ctrl->S.active);
299 Ctrl->S.active = true;
300 Ctrl->S.arg = strdup (opt->arg);
301 break;
302 default: /* Report bad options */
303 n_errors += gmt_default_error (GMT, opt->option);
304 break;
305 }
306 }
307
308 if ((Ctrl->T.mode == DO_ROT3D || Ctrl->T.mode == DO_POLE) && gmt_M_is_cartesian (GMT, GMT_IN)) gmt_parse_common_options (GMT, "f", 'f', "g"); /* Set -fg unless already set for 3-D rots and pole ops */
309 if (Ctrl->T.dmode >= GMT_GREATCIRCLE || GMT->common.j.active) gmt_parse_common_options (GMT, "f", 'f', "g"); /* Set -fg implicitly */
310 n_in = (Ctrl->C.active[GMT_IN] && gmt_M_is_geographic (GMT, GMT_IN)) ? 3 : 2;
311 if (Ctrl->T.a_and_d) n_in += 2; /* Must read azimuth and distance as well */
312 if (GMT->common.b.active[GMT_IN] && GMT->common.b.ncol[GMT_IN] == 0) GMT->common.b.ncol[GMT_IN] = n_in;
313 n_errors += gmt_M_check_condition (GMT, GMT->common.b.active[GMT_IN] && GMT->common.b.ncol[GMT_IN] < n_in, "Binary input data (-bi) must have at least %d columns\n", n_in);
314 n_errors += gmt_M_check_condition (GMT, Ctrl->S.active && Ctrl->S.arg && !gmt_access (GMT, Ctrl->S.arg, R_OK), "Option -S: Secondary vector cannot be a file!\n");
315 n_errors += gmt_M_check_condition (GMT, n_files > 1, "Only one output destination can be specified\n");
316 n_errors += gmt_M_check_condition (GMT, Ctrl->In.n_args && Ctrl->A.active && Ctrl->A.mode == 0, "Cannot give input files and -A<vec> at the same time\n");
317
318 return (n_errors ? GMT_PARSE_ERROR : GMT_NOERROR);
319 }
320
gmtvector_decode_vector(struct GMT_CTRL * GMT,char * arg,double coord[],int cartesian,int geocentric)321 GMT_LOCAL unsigned int gmtvector_decode_vector (struct GMT_CTRL *GMT, char *arg, double coord[], int cartesian, int geocentric) {
322 unsigned int n_out, n_errors = 0, ix, iy;
323 int n;
324 char txt_a[GMT_LEN64] = {""}, txt_b[GMT_LEN64] = {""}, txt_c[GMT_LEN64] = {""};
325
326 ix = (GMT->current.setting.io_lonlat_toggle[GMT_IN]); iy = 1 - ix;
327 n = sscanf (arg, "%[^/]/%[^/]/%s", txt_a, txt_b, txt_c);
328 assert (n >= 2);
329 n_out = n;
330 if (n == 2) { /* Got lon/lat, r/theta, or x/y */
331 if (gmt_M_is_geographic (GMT, GMT_IN)) {
332 n_errors += gmt_verify_expectations (GMT, gmt_M_type (GMT, GMT_IN, ix), gmt_scanf_arg (GMT, txt_a, gmt_M_type (GMT, GMT_IN, ix), false, &coord[ix]), txt_a);
333 n_errors += gmt_verify_expectations (GMT, gmt_M_type (GMT, GMT_IN, iy), gmt_scanf_arg (GMT, txt_b, gmt_M_type (GMT, GMT_IN, iy), false, &coord[iy]), txt_b);
334 if (geocentric) coord[GMT_Y] = gmt_lat_swap (GMT, coord[GMT_Y], GMT_LATSWAP_G2O);
335 gmt_geo_to_cart (GMT, coord[GMT_Y], coord[GMT_X], coord, true); /* get x/y/z */
336 n_out = 3;
337 }
338 else if (cartesian) { /* Cartesian x/y */
339 coord[GMT_X] = atof (txt_a);
340 coord[GMT_Y] = atof (txt_b);
341 }
342 else /* Cylindrical r/theta */
343 gmt_polar_to_cart (GMT, atof (txt_a), atof (txt_b), coord, true);
344 }
345 else if (n == 3) { /* Got x/y/z */
346 coord[GMT_X] = atof (txt_a);
347 coord[GMT_Y] = atof (txt_b);
348 coord[GMT_Z] = atof (txt_c);
349 }
350 else {
351 GMT_Report (GMT->parent, GMT_MSG_ERROR, "Bad vector argument (%s)\n", arg);
352 return (0);
353 }
354 if (n_errors) {
355 GMT_Report (GMT->parent, GMT_MSG_ERROR, "Failed to decode the geographic coordinates (%s)\n", arg);
356 return (0);
357 }
358 return (n_out);
359 }
360
gmtvector_get_bisector(struct GMT_CTRL * GMT,double A[3],double B[3],double P[3])361 GMT_LOCAL void gmtvector_get_bisector (struct GMT_CTRL *GMT, double A[3], double B[3], double P[3]) {
362 /* Given points in A and B, return the bisector pole via P */
363
364 unsigned int i;
365 double Pa[3], M[3];
366
367 /* Get mid point between A and B */
368
369 for (i = 0; i < 3; i++) M[i] = A[i] + B[i];
370 gmt_normalize3v (GMT, M);
371
372 /* Get pole for great circle through A and B */
373
374 gmt_cross3v (GMT, A, B, Pa);
375 gmt_normalize3v (GMT, Pa);
376
377 /* Then get pole for bisector of A-B through Pa */
378
379 gmt_cross3v (GMT, M, Pa, P);
380 gmt_normalize3v (GMT, P);
381 }
382
gmtvector_get_azpole(struct GMT_CTRL * GMT,double A[3],double P[3],double az)383 GMT_LOCAL void gmtvector_get_azpole (struct GMT_CTRL *GMT, double A[3], double P[3], double az) {
384 /* Given point in A and azimuth az, return the pole P to the oblique equator with given az at A */
385 double R[3][3], tmp[3], B[3] = {0.0, 0.0, 1.0}; /* set B to north pole */
386 gmt_cross3v (GMT, A, B, tmp); /* Point C is 90 degrees away from plan through A and B */
387 gmt_normalize3v (GMT, tmp); /* Get unit vector */
388 gmt_make_rot_matrix2 (GMT, A, -az, R); /* Make rotation about A of -azim degrees */
389 gmt_matrix_vect_mult (GMT, 3U, R, tmp, P);
390 }
391
gmtvector_translate_point(struct GMT_CTRL * GMT,double A[3],double B[3],double a_d[],bool geo)392 GMT_LOCAL void gmtvector_translate_point (struct GMT_CTRL *GMT, double A[3], double B[3], double a_d[], bool geo) {
393 /* Given point in A, azimuth az and distance d, return the point P away from A */
394 if (geo)
395 GMT->current.map.second_point (GMT, A[GMT_X], A[GMT_Y], a_d[0], a_d[1], &B[GMT_X], &B[GMT_Y], NULL);
396 // gmt_translate_point (GMT, A[GMT_X], A[GMT_Y], a_d[0], a_d[1], &B[GMT_X], &B[GMT_Y]);
397 else { /* Cartesian translation */
398 double s, c;
399 sincosd (90.0 - a_d[0], &s, &c);
400 B[GMT_X] = A[GMT_X] + a_d[1] * c;
401 B[GMT_Y] = A[GMT_Y] + a_d[1] * s;
402 }
403 }
404
gmtvector_mean_vector(struct GMT_CTRL * GMT,struct GMT_DATASET * D,bool cartesian,double conf,bool geocentric,double * M,double * E)405 GMT_LOCAL void gmtvector_mean_vector (struct GMT_CTRL *GMT, struct GMT_DATASET *D, bool cartesian, double conf, bool geocentric, double *M, double *E) {
406 /* Determines the mean vector M and the covariance matrix C */
407
408 unsigned int i, j, k, n_components, nrots, geo = gmt_M_is_geographic (GMT, GMT_IN);
409 uint64_t row, n, seg, tbl, p;
410 double lambda[3], V[9], work1[3], work2[3], lon, lat, lon2, lat2, scl, L, Y;
411 double *P[3], X[3], B[3], C[9];
412 struct GMT_DATASEGMENT *S = NULL;
413
414 gmt_M_memset (M, 3, double);
415 n_components = (geo || D->n_columns == 3) ? 3 : 2;
416 for (k = 0; k < n_components; k++) P[k] = gmt_M_memory (GMT, NULL, D->n_records, double);
417 for (tbl = n = 0; tbl < D->n_tables; tbl++) {
418 for (seg = 0; seg < D->table[tbl]->n_segments; seg++) {
419 S = D->table[tbl]->segment[seg];
420 for (row = 0; row < S->n_rows; row++) {
421 if (!cartesian) { /* Want to turn geographic or polar into Cartesian */
422 if (geo) {
423 lat = S->data[GMT_Y][row];
424 if (geocentric) lat = gmt_lat_swap (GMT, lat, GMT_LATSWAP_G2O); /* Get geocentric latitude */
425 gmt_geo_to_cart (GMT, lat, S->data[GMT_X][row], X, true); /* Get x/y/z */
426 }
427 else
428 gmt_polar_to_cart (GMT, S->data[GMT_X][row], S->data[GMT_Y][row], X, true);
429 for (k = 0; k < n_components; k++) P[k][n] = X[k];
430 }
431 else for (k = 0; k < n_components; k++) P[k][n] = S->data[k][row];
432 for (k = 0; k < n_components; k++) M[k] += P[k][n];
433 n++;
434 }
435 }
436 }
437 if (n > 0)
438 for (k = 0; k < n_components; k++) M[k] /= n; /* Compute mean resultant vector [not normalized] */
439 gmt_M_memcpy (X, M, 3, double);
440
441 /* Compute covariance matrix */
442
443 gmt_M_memset (C, 9, double);
444 for (j = k = 0; j < n_components; j++) for (i = 0; i < n_components; i++, k++) {
445 for (p = 0; p < n; p++) C[k] += (P[j][p] - M[j]) * (P[i][p] - M[i]);
446 C[k] /= (n - 1.0);
447 }
448 for (k = 0; k < n_components; k++) gmt_M_free (GMT, P[k]);
449
450 if (gmt_jacobi (GMT, C, n_components, n_components, lambda, V, work1, work2, &nrots)) { /* Solve eigen-system */
451 GMT_Report (GMT->parent, GMT_MSG_WARNING, "Eigenvalue routine failed to converge in 50 sweeps.\n");
452 }
453 if (n_components == 3 && geo) { /* Special case of geographic unit vectors; normalize and recover lon,lat */
454 gmt_normalize3v (GMT, M);
455 gmt_cart_to_geo (GMT, &lat, &lon, M, true);
456 if (geocentric) lat = gmt_lat_swap (GMT, lat, GMT_LATSWAP_G2O+1); /* Get geodetic latitude */
457 if (lon < 0.0) lon += 360.0;
458 }
459 else /* Cartesian, get y-component */
460 lat = X[GMT_Y];
461
462 /* Find the xy[z] point (B) of end of eigenvector 1: */
463 gmt_M_memset (B, 3, double);
464 for (k = 0; k < n_components; k++) B[k] = X[k] + sqrt (lambda[0]) * V[k];
465 L = sqrt (B[0] * B[0] + B[1] * B[1] + B[2] * B[2]); /* Length of B */
466 for (k = 0; k < n_components; k++) B[k] /= L; /* Normalize */
467 if (n_components == 3 && geo) { /* Special case of geographic unit vectors; normalize and recover lon,lat */
468 gmt_normalize3v (GMT, B);
469 gmt_cart_to_geo (GMT, &lat2, &lon2, B, true);
470 if (geocentric) lat2 = gmt_lat_swap (GMT, lat2, GMT_LATSWAP_G2O+1); /* Get geodetic latitude */
471 if (lon2 < 0.0) lon2 += 360.0;
472 gmt_M_set_delta_lon (lon, lon2, L);
473 scl = cosd (lat); /* Local flat-Earth approximation */
474 }
475 else { /* Cartesian, get y-component */
476 lat2 = B[GMT_Y];
477 scl = 1.0; /* No flat-Earth scaling here */
478 }
479 Y = lat2 - lat;
480
481 E[0] = 90.0 - atan2 (Y, L * scl) * R2D; /* Get azimuth */
482 /* Convert to 95% confidence (in km, if geographic) */
483
484 scl = sqrt (gmt_chi2crit (GMT, conf, n_components));
485 if (geo) scl *= GMT->current.proj.DIST_KM_PR_DEG * R2D;
486 E[1] = 2.0 * sqrt (lambda[0]) * scl; /* 2x since we need the major axis not semi-major */
487 E[2] = 2.0 * sqrt (lambda[1]) * scl;
488
489 GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "%g%% confidence ellipse on mean position: Major axis = %g Minor axis = %g Major axis azimuth = %g\n", 100.0 * conf, E[1], E[2], E[0]);
490 }
491
gmtvector_gmt_make_rot2d_matrix(double angle,double R[3][3])492 GMT_LOCAL void gmtvector_gmt_make_rot2d_matrix (double angle, double R[3][3]) {
493 double s, c;
494 gmt_M_memset (R, 9, double);
495 sincosd (angle, &s, &c);
496 R[0][0] = c; R[0][1] = -s;
497 R[1][0] = s; R[1][1] = c;
498 }
499
gmtvector_dist_to_meter(struct GMT_CTRL * GMT,double d_in)500 GMT_LOCAL double gmtvector_dist_to_meter (struct GMT_CTRL *GMT, double d_in) {
501 double d_out = d_in / GMT->current.map.dist[GMT_MAP_DIST].scale; /* Now in degrees or meters */
502 if (GMT->current.map.dist[GMT_MAP_DIST].arc) /* Got arc measure */
503 d_out *= GMT->current.proj.DIST_M_PR_DEG; /* Now in degrees */
504 return (d_out);
505 }
506
gmtvector_dist_to_degree(struct GMT_CTRL * GMT,double d_in)507 GMT_LOCAL double gmtvector_dist_to_degree (struct GMT_CTRL *GMT, double d_in) {
508 double d_out = d_in / GMT->current.map.dist[GMT_MAP_DIST].scale; /* Now in degrees or meters */
509 if (!GMT->current.map.dist[GMT_MAP_DIST].arc) /* Got unit distance measure */
510 d_out /= GMT->current.proj.DIST_M_PR_DEG; /* Now in degrees */
511 return (d_out);
512 }
513
514 #define bailout(code) {gmt_M_free_options (mode); return (code);}
515 #define Return(code) {Free_Ctrl (GMT, Ctrl); gmt_end_module (GMT, GMT_cpy); bailout (code);}
516
GMT_gmtvector(void * V_API,int mode,void * args)517 EXTERN_MSC int GMT_gmtvector (void *V_API, int mode, void *args) {
518 unsigned int tbl, error = 0, k, n, n_in, n_components, n_out, add_cols = 0, geo;
519 bool single = false, convert;
520
521 uint64_t row, seg;
522
523 double out[3], vector_1[3], vector_2[3], vector_3[3], R[3][3], E[3];
524
525 struct GMT_DATASET *Din = NULL, *Dout = NULL;
526 struct GMT_DATASEGMENT *Sin = NULL, *Sout = NULL;
527 struct GMT_DATASET_HIDDEN *DH = NULL;
528 struct GMTVECTOR_CTRL *Ctrl = NULL;
529 struct GMT_CTRL *GMT = NULL, *GMT_cpy = NULL;
530 struct GMT_OPTION *options = NULL;
531 struct GMTAPI_CTRL *API = gmt_get_api_ptr (V_API); /* Cast from void to GMTAPI_CTRL pointer */
532
533 /*----------------------- Standard module initialization and parsing ----------------------*/
534
535 if (API == NULL) return (GMT_NOT_A_SESSION);
536 if (mode == GMT_MODULE_PURPOSE) return (usage (API, GMT_MODULE_PURPOSE)); /* Return the purpose of program */
537 options = GMT_Create_Options (API, mode, args); if (API->error) return (API->error); /* Set or get option list */
538
539 if ((error = gmt_report_usage (API, options, 0, usage)) != GMT_NOERROR) bailout (error); /* Give usage if requested */
540
541 /* Parse the command-line arguments */
542
543 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 */
544 if (GMT_Parse_Common (API, THIS_MODULE_OPTIONS, options)) Return (API->error);
545 Ctrl = New_Ctrl (GMT); /* Allocate and initialize a new control structure */
546 if ((error = parse (GMT, Ctrl, options)) != 0) Return (error);
547 geo = gmt_M_is_geographic (GMT, GMT_IN);
548
549 /*---------------------------- This is the gmtvector main code ----------------------------*/
550
551 gmt_M_memset (vector_1, 3, double);
552 gmt_M_memset (vector_2, 3, double);
553 gmt_M_memset (vector_3, 3, double);
554 if (Ctrl->T.mode == DO_ROT3D) /* Spherical 3-D rotation */
555 gmt_make_rot_matrix (GMT, Ctrl->T.par[0], Ctrl->T.par[1], Ctrl->T.par[2], R);
556 else if (Ctrl->T.mode == DO_ROT2D) /* Cartesian 2-D rotation */
557 gmtvector_gmt_make_rot2d_matrix (Ctrl->T.par[2], R);
558 else if (!(Ctrl->T.mode == DO_NOTHING || Ctrl->T.mode == DO_POLE || Ctrl->T.mode == DO_TRANSLATE)) { /* Will need secondary vector, get that first before input file */
559 n = gmtvector_decode_vector (GMT, Ctrl->S.arg, vector_2, Ctrl->C.active[GMT_IN], Ctrl->E.active);
560 if (n == 0) Return (GMT_RUNTIME_ERROR);
561 if (Ctrl->T.mode == DO_DOT) { /* Must normalize to turn dot-product into angle */
562 if (n == 2)
563 gmt_normalize2v (GMT, vector_2);
564 else
565 gmt_normalize3v (GMT, vector_2);
566 Ctrl->C.active[GMT_OUT] = true; /* Since we just want to return the angle */
567 }
568 }
569
570 if (Ctrl->T.mode == DO_TRANSLATE && geo) { /* Initialize distance machinery */
571 int def_mode = (GMT->common.j.active) ? GMT->common.j.mode : GMT_GREATCIRCLE;
572 if (Ctrl->T.a_and_d) { /* Read az and dist from file, set units here */
573 Ctrl->T.dmode = def_mode; /* Store the setting here */
574 error = gmt_init_distaz (GMT, Ctrl->T.unit, Ctrl->T.dmode, GMT_MAP_DIST);
575 }
576 else { /* Got a fixed set, set units here */
577 if (GMT->common.j.active) Ctrl->T.dmode = def_mode; /* Override any setting in -T */
578 error = gmt_init_distaz (GMT, Ctrl->T.unit, Ctrl->T.dmode, GMT_MAP_DIST);
579 if (Ctrl->T.dmode == GMT_GREATCIRCLE)
580 Ctrl->T.par[1] = gmtvector_dist_to_degree (GMT, Ctrl->T.par[1]); /* Make sure we have degrees from whatever -Tt set */
581 else /* Get meters */
582 Ctrl->T.par[1] = gmtvector_dist_to_meter (GMT, Ctrl->T.par[1]); /* Make sure we have degrees from whatever -Tt set */
583 }
584 if (error == GMT_NOT_A_VALID_TYPE) Return (error);
585 }
586
587 /* Read input data set */
588
589 n_in = (geo && Ctrl->C.active[GMT_IN]) ? 3 : 2;
590 if (Ctrl->T.a_and_d) n_in += 2; /* Must read azimuth and distance as well */
591
592 if (Ctrl->A.active) { /* Want a single primary vector */
593 uint64_t dim[GMT_DIM_SIZE] = {1, 1, 1, 3};
594 GMT_Report (API, GMT_MSG_INFORMATION, "Processing single input vector; no files are read\n");
595 if (Ctrl->A.mode) { /* Compute the mean of all input vectors */
596 if (GMT_Init_IO (API, GMT_IS_DATASET, GMT_IS_POINT, GMT_IN, GMT_ADD_DEFAULT, 0, options) != GMT_NOERROR) { /* Registers default input sources, unless already set */
597 Return (API->error);
598 }
599 if ((Din = GMT_Read_Data (API, GMT_IS_DATASET, GMT_IS_FILE, 0, GMT_READ_NORMAL, NULL, NULL, NULL)) == NULL) {
600 Return (API->error);
601 }
602 if (Din->n_columns < 2) {
603 GMT_Report (API, GMT_MSG_ERROR, "Input data have %d column(s) but at least %u are needed\n", (int)Din->n_columns, n_in);
604 Return (GMT_DIM_TOO_SMALL);
605 }
606 n = n_out = (Ctrl->C.active[GMT_OUT] && (Din->n_columns == 3 || geo)) ? 3 : 2;
607 gmtvector_mean_vector (GMT, Din, Ctrl->C.active[GMT_IN], Ctrl->A.conf, Ctrl->E.active, vector_1, E); /* Get mean vector and confidence ellipse parameters */
608 if (GMT_Destroy_Data (API, &Din) != GMT_NOERROR) {
609 Return (API->error);
610 }
611 add_cols = 3; /* Make space for angle major minor */
612 }
613 else { /* Decode a single vector */
614 n = gmtvector_decode_vector (GMT, Ctrl->A.arg, vector_1, Ctrl->C.active[GMT_IN], Ctrl->E.active);
615 if (n == 0) Return (GMT_RUNTIME_ERROR);
616 if (Ctrl->T.mode == DO_DOT) { /* Must normalize before we turn dot-product into angle */
617 if (n == 2)
618 gmt_normalize2v (GMT, vector_1);
619 else
620 gmt_normalize3v (GMT, vector_1);
621 }
622 n_out = (Ctrl->C.active[GMT_OUT] && n == 3) ? 3 : 2;
623 }
624 if ((Din = GMT_Create_Data (API, GMT_IS_DATASET, GMT_IS_POINT, 0, dim, NULL, NULL, 0, 0, NULL)) == NULL) Return (API->error);
625 n_components = (n == 3 || geo) ? 3 : 2; /* Number of Cartesian vector components */
626 for (k = 0; k < n_components; k++) Din->table[0]->segment[0]->data[k][0] = vector_1[k];
627 Din->table[0]->segment[0]->n_rows = 1;
628 single = true;
629 }
630 else { /* Read input files or stdin */
631 GMT_Report (API, GMT_MSG_INFORMATION, "Processing input table data\n");
632 if (GMT_Init_IO (API, GMT_IS_DATASET, GMT_IS_POINT, GMT_IN, GMT_ADD_DEFAULT, 0, options) != GMT_NOERROR) { /* Registers default input sources, unless already set */
633 Return (API->error);
634 }
635 if ((Din = GMT_Read_Data (API, GMT_IS_DATASET, GMT_IS_FILE, 0, GMT_READ_NORMAL, NULL, NULL, NULL)) == NULL) {
636 Return (API->error);
637 }
638 if (Din->n_columns < 2) {
639 GMT_Report (API, GMT_MSG_ERROR, "Input data have %d column(s) but at least %u are needed\n", (int)Din->n_columns, n_in);
640 Return (GMT_DIM_TOO_SMALL);
641 }
642 n = n_out = (Ctrl->C.active[GMT_OUT] && (Din->n_columns == 3 || geo)) ? 3 : 2;
643 }
644
645 if (Ctrl->T.mode == DO_DOT) {
646 n_out = 1; /* Override prior setting since we just will report an angle in one column */
647 gmt_set_column_type (GMT, GMT_OUT, GMT_X, GMT_IS_FLOAT);
648 }
649 else if (Ctrl->T.mode == DO_ROTVAR2D) { /* 2D or 3D */
650 if (Din->n_columns == 1) n_out = 2;
651 if (Din->n_columns == 3) Ctrl->T.mode = DO_ROTVAR3D; /* OK, it is 3D */
652 }
653 else if (Ctrl->C.active[GMT_OUT] || !gmt_M_is_geographic (GMT, GMT_OUT)) /* Override types since output is Cartesian or polar coordinates, not lon/lat */
654 gmt_set_cartesian (GMT, GMT_OUT);
655 DH = gmt_get_DD_hidden (Din);
656 DH->dim[GMT_COL] = n_out + add_cols; /* State we want a different set of columns on output */
657 Dout = GMT_Duplicate_Data (API, GMT_IS_DATASET, GMT_DUPLICATE_ALLOC, Din);
658 gmt_M_memset (out, 3, double);
659
660 /* OK, with data in hand we can do some damage */
661
662 if (Ctrl->T.mode == DO_DOT) Ctrl->T.mode = (n == 3 || geo) ? DO_DOT3D : DO_DOT2D;
663 n_components = (n == 3 || geo) ? 3 : 2; /* Number of Cartesian vector components */
664 if (Ctrl->T.mode == DO_ROTVAR2D) n_components = 1; /* Override in case of 2-D Cartesian rotation angles on input */
665 convert = (!single && !Ctrl->C.active[GMT_IN] && !(Ctrl->T.mode == DO_ROTVAR2D || Ctrl->T.mode == DO_ROTVAR3D));
666 if (Ctrl->T.mode == DO_TRANSLATE) {
667 convert = false;
668 n_components = 2;
669 }
670 for (tbl = 0; tbl < Din->n_tables; tbl++) {
671 for (seg = 0; seg < Din->table[tbl]->n_segments; seg++) {
672 Sin = Din->table[tbl]->segment[seg];
673 Sout = Dout->table[tbl]->segment[seg];
674 if (Sin->header) Sout->header = strdup (Sin->header);
675 for (row = 0; row < Sin->n_rows; row++) {
676 if (convert) { /* Want to turn geographic or polar into Cartesian */
677 if (geo)
678 gmt_geo_to_cart (GMT, Sin->data[GMT_Y][row], Sin->data[GMT_X][row], vector_1, true); /* get x/y/z */
679 else
680 gmt_polar_to_cart (GMT, Sin->data[GMT_X][row], Sin->data[GMT_Y][row], vector_1, true);
681 }
682 else for (k = 0; k < n_components; k++) vector_1[k] = Sin->data[k][row];
683
684 switch (Ctrl->T.mode) {
685 case DO_AVERAGE: /* Get sum of 2-D or 3-D vectors and compute average */
686 for (k = 0; k < n_components; k++) vector_3[k] = 0.5 * (vector_1[k] + vector_2[k]);
687 break;
688 case DO_BISECTOR: /* Compute pole or bisector of vector 1 and 2 */
689 gmtvector_get_bisector (GMT, vector_1, vector_2, vector_3);
690 break;
691 case DO_DOT: /* Just here to quiet the compiler as DO_DOT has been replaced in line 552 above */
692 case DO_DOT2D: /* Get angle between 2-D vectors */
693 if (!single) gmt_normalize2v (GMT, vector_1);
694 vector_3[0] = gmt_dot2v (GMT, vector_1, vector_2);
695 if (Ctrl->T.degree) vector_3[0] = acosd (vector_3[0]);
696 break;
697 case DO_DOT3D: /* Get angle between 3-D vectors */
698 if (!single) gmt_normalize3v (GMT, vector_1);
699 vector_3[0] = gmt_dot3v (GMT, vector_1, vector_2);
700 if (Ctrl->T.degree) vector_3[0] = acosd (vector_3[0]);
701 break;
702 case DO_CROSS: /* Get cross-product of 3-D vectors */
703 gmt_cross3v (GMT, vector_1, vector_2, vector_3);
704 break;
705 case DO_SUM: /* Get sum of 2-D or 3-D vectors */
706 for (k = 0; k < n_components; k++) vector_3[k] = vector_1[k] + vector_2[k];
707 break;
708 case DO_ROT2D: /* Rotate a 2-D vector about the z_axis */
709 gmt_matrix_vect_mult (GMT, 2U, R, vector_1, vector_3);
710 break;
711 case DO_ROT3D: /* Rotate a 3-D vector about an arbitrary pole encoded in 3x3 matrix R */
712 gmt_matrix_vect_mult (GMT, 3U, R, vector_1, vector_3);
713 break;
714 case DO_ROTVAR2D: /* Rotate a 2-D vector about the z_axis */
715 gmtvector_gmt_make_rot2d_matrix (vector_1[0], R);
716 gmt_matrix_vect_mult (GMT, 2U, R, vector_2, vector_3);
717 break;
718 case DO_ROTVAR3D: /* Rotate a 3-D vector about an arbitrary pole encoded in 3x3 matrix R */
719 gmt_make_rot_matrix (GMT, vector_1[0], vector_1[1], vector_1[2], R);
720 gmt_matrix_vect_mult (GMT, 3U, R, vector_2, vector_3);
721 break;
722 case DO_POLE: /* Return pole of great circle defined by center point an azimuth */
723 gmtvector_get_azpole (GMT, vector_1, vector_3, Ctrl->T.par[0]);
724 break;
725 case DO_TRANSLATE: /* Return translated points moved a distance d in the direction of azimuth */
726 if (Ctrl->T.a_and_d) { /* Get azimuth and distance from input file */
727 Ctrl->T.par[0] = Sin->data[2][row];
728 if (geo)
729 Ctrl->T.par[1] = (Ctrl->T.dmode == GMT_GREATCIRCLE) ? gmtvector_dist_to_degree (GMT, Sin->data[3][row]) : gmtvector_dist_to_meter (GMT, Sin->data[3][row]); /* Make sure we have degrees or meters for calculations */
730 else
731 Ctrl->T.par[1] = Sin->data[3][row]; /* Pass as is */
732 }
733 gmtvector_translate_point (GMT, vector_1, vector_3, Ctrl->T.par, geo);
734 break;
735 case DO_NOTHING: /* Probably just want the effect of -C, -E, -N */
736 gmt_M_memcpy (vector_3, vector_1, n_components, double);
737 break;
738 }
739 if (Ctrl->T.mode == DO_TRANSLATE) /* Want to turn Cartesian output into something else */
740 gmt_M_memcpy (out, vector_3, 2, double);
741 else if (Ctrl->C.active[GMT_OUT]) { /* Report Cartesian output... */
742 if (Ctrl->N.active) { /* ...but normalize first */
743 if (n_out == 3)
744 gmt_normalize3v (GMT, vector_3);
745 else if (n_out == 2)
746 gmt_normalize2v (GMT, vector_3);
747 }
748 gmt_M_memcpy (out, vector_3, n_out, double);
749 }
750 else { /* Want to turn Cartesian output into something else */
751 if (gmt_M_is_geographic (GMT, GMT_OUT)) {
752 gmt_normalize3v (GMT, vector_3); /* Must always normalize before calling gmt_cart_to_geo */
753 gmt_cart_to_geo (GMT, &(out[GMT_Y]), &(out[GMT_X]), vector_3, true); /* Get lon/lat */
754 if (Ctrl->E.active) out[GMT_Y] = gmt_lat_swap (GMT, out[GMT_Y], GMT_LATSWAP_G2O + 1); /* Convert to geodetic */
755 }
756 else {
757 if (Ctrl->N.active) gmt_normalize2v (GMT, vector_3); /* Optional for cart_to_polar */
758 gmt_cart_to_polar (GMT, &(out[GMT_X]), &(out[GMT_Y]), vector_3, true); /* Get r/theta */
759 }
760 }
761 for (k = 0; k < n_out; k++) Sout->data[k][row] = out[k];
762 }
763 }
764 }
765
766 if (Ctrl->A.mode && Sout) for (k = 0; k < 3; k++) Sout->data[k+n_out][0] = E[k]; /* Place az, major, minor in the single output record */
767
768 /* Time to write out the results */
769
770 if (GMT_Init_IO (API, GMT_IS_DATASET, GMT_IS_POINT, GMT_OUT, GMT_ADD_DEFAULT, 0, options) != GMT_NOERROR) { /* Establishes data output */
771 Return (API->error);
772 }
773 if (GMT_Write_Data (API, GMT_IS_DATASET, GMT_IS_FILE, GMT_IS_POINT, GMT_WRITE_SET, NULL, Ctrl->Out.file, Dout) != GMT_NOERROR) {
774 Return (API->error);
775 }
776 if (single && GMT_Destroy_Data (API, &Din) != GMT_NOERROR) {
777 Return (API->error);
778 }
779
780 Return (GMT_NOERROR);
781 }
782