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: read a file of lon, lat, zvalue[, distance, azimuth]
19 * and split it into profile segments.
20 *
21 * Author: W. H. F. Smith
22 * Date: 1 JAN 2010
23 * Version: 6 API
24 */
25
26 #include "gmt_dev.h"
27
28 #define THIS_MODULE_CLASSIC_NAME "gmtsplit"
29 #define THIS_MODULE_MODERN_NAME "gmtsplit"
30 #define THIS_MODULE_LIB "core"
31 #define THIS_MODULE_PURPOSE "Split xyz[dh] data tables into individual segments"
32 #define THIS_MODULE_KEYS "<D{,>D}"
33 #define THIS_MODULE_NEEDS ""
34 #define THIS_MODULE_OPTIONS "-:>Vbdefghiqs" GMT_OPT("H")
35
36 #define SPLITXYZ_F_RES 1000 /* Number of points in filter halfwidth */
37 #define SPLITXYZ_N_OUTPUT_CHOICES 5
38
39 struct SPLITXYZ_CTRL {
40 struct SPLITXYZ_Out { /* -> */
41 bool active;
42 char *file;
43 } Out;
44 struct SPLITXYZ_A { /* -A<azimuth>/<tolerance> */
45 bool active;
46 double azimuth, tolerance;
47 } A;
48 struct SPLITXYZ_C { /* -C<course_change> */
49 bool active;
50 double value;
51 } C;
52 struct SPLITXYZ_D { /* -D<mindist> */
53 bool active;
54 double value;
55 } D;
56 struct SPLITXYZ_F { /* -F<xy_filter>/<z_filter> */
57 bool active;
58 double xy_filter, z_filter;
59 } F;
60 struct SPLITXYZ_N { /* -N<namestem> */
61 bool active;
62 char *name;
63 } N;
64 struct SPLITXYZ_Q { /* -Q[<xyzdg>] */
65 bool active;
66 bool z_selected;
67 char col[SPLITXYZ_N_OUTPUT_CHOICES]; /* Character codes for desired output in the right order */
68 } Q;
69 struct SPLITXYZ_S { /* -S */
70 bool active;
71 } S;
72 };
73
gmtsplit_filterxy_setup(struct GMT_CTRL * GMT)74 GMT_LOCAL double *gmtsplit_filterxy_setup (struct GMT_CTRL *GMT) {
75 unsigned int i;
76 double tmp, sum = 0.0, *fwork = NULL;
77
78 fwork = gmt_M_memory (GMT, NULL, SPLITXYZ_F_RES, double); /* Initialized to zeros */
79 tmp = M_PI / SPLITXYZ_F_RES;
80 for (i = 0; i < SPLITXYZ_F_RES; i++) {
81 fwork[i] = 1.0 + cos (i * tmp);
82 sum += fwork[i];
83 }
84 for (i = 1; i < SPLITXYZ_F_RES; i++) fwork[i] /= sum;
85 return (fwork);
86 }
87
gmtsplit_filter_cols(struct GMT_CTRL * GMT,double * data[],uint64_t begin,uint64_t end,unsigned int d_col,unsigned int n_cols,unsigned int cols[],double filter_width,double * fwork)88 GMT_LOCAL void gmtsplit_filter_cols (struct GMT_CTRL *GMT, double *data[], uint64_t begin, uint64_t end, unsigned int d_col, unsigned int n_cols, unsigned int cols[], double filter_width, double *fwork) {
89 uint64_t i, j, k, p, istart, istop, ndata;
90 int64_t kk;
91 bool hilow;
92 double half_width, dt, sum, **w = NULL;
93
94 if (filter_width == 0.0) return; /* No filtering */
95 hilow = (filter_width < 0.0);
96 half_width = 0.5 * fabs (filter_width);
97 dt = SPLITXYZ_F_RES / half_width;
98 ndata = end - begin;
99 w = gmt_M_memory (GMT, NULL, n_cols, double *);
100 for (k = 0; k < n_cols; k++) w[k] = gmt_M_memory (GMT, NULL, ndata, double); /* Initialized to zeros */
101 j = istart = istop = begin;
102 while (j < end) {
103 while (istart < end && data[d_col][istart] - data[d_col][j] <= -half_width) istart++;
104 while (istop < end && data[d_col][istop] - data[d_col][j] < half_width) istop++;
105 for (i = istart, sum = 0.0; i < istop; i++) {
106 kk = lrint (floor (dt * fabs (data[d_col][i] - data[d_col][j])));
107 if (kk < 0 || kk >= SPLITXYZ_F_RES) continue; /* Safety valve */
108 k = kk;
109 sum += fwork[k];
110 for (p = 0; p < n_cols; p++) w[p][j] += (data[cols[p]][i] * fwork[k]);
111 }
112 for (p = 0; p < n_cols; p++) w[p][j] /= sum;
113 j++;
114 }
115 if (hilow) {
116 for (i = begin; i < end; i++) for (p = 0; p < n_cols; p++) data[cols[p]][i] -= w[p][i];
117 }
118 else {
119 for (i = begin; i < end; i++) for (p = 0; p < n_cols; p++) data[cols[p]][i] = w[p][i];
120 }
121 for (p = 0; p < n_cols; p++) gmt_M_free (GMT, w[p]);
122 gmt_M_free (GMT, w);
123 }
124
New_Ctrl(struct GMT_CTRL * GMT)125 static void *New_Ctrl (struct GMT_CTRL *GMT) { /* Allocate and initialize a new control structure */
126 struct SPLITXYZ_CTRL *C = NULL;
127
128 C = gmt_M_memory (GMT, NULL, 1, struct SPLITXYZ_CTRL);
129
130 /* Initialize values whose defaults are not 0/false/NULL */
131 C->A.azimuth = 90.0;
132 C->A.tolerance = 360.0;
133 C->C.value = 180.0; /* Tolerate any course change */
134 return (C);
135 }
136
Free_Ctrl(struct GMT_CTRL * GMT,struct SPLITXYZ_CTRL * C)137 static void Free_Ctrl (struct GMT_CTRL *GMT, struct SPLITXYZ_CTRL *C) { /* Deallocate control structure */
138 if (!C) return;
139 gmt_M_str_free (C->Out.file);
140 gmt_M_str_free (C->N.name);
141 gmt_M_free (GMT, C);
142 }
143
usage(struct GMTAPI_CTRL * API,int level)144 static int usage (struct GMTAPI_CTRL *API, int level) {
145 /* This displays the gmtsplit synopsis and optionally full usage information */
146
147 const char *name = gmt_show_name_and_purpose (API, THIS_MODULE_LIB, THIS_MODULE_CLASSIC_NAME, THIS_MODULE_PURPOSE);
148 if (level == GMT_MODULE_PURPOSE) return (GMT_NOERROR);
149 GMT_Usage (API, 0, "usage: %s [<table>] [-A<azimuth>/<tolerance>] [-C<course_change>] [-D<minimum_distance>] "
150 "[-F<xy_filter>/<z_filter>] [-N<template>] [-Q<flags>] [-S] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s]\n",
151 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_q_OPT, GMT_s_OPT, GMT_colon_OPT, GMT_PAR_OPT);
152
153 if (level == GMT_SYNOPSIS) return (GMT_MODULE_SYNOPSIS);
154
155 GMT_Message (API, GMT_TIME_NONE, " REQUIRED ARGUMENTS:\n");
156 GMT_Usage (API, 1, "\n<table>");
157 GMT_Usage (API, -2, "One or more data files (in ASCII, binary, netCDF) with 2, 3 or 5 columns. "
158 "If no files are given, standard input is read.");
159 GMT_Message (API, GMT_TIME_NONE, "\n OPTIONAL ARGUMENTS:\n");
160 GMT_Usage (API, 1, "\n-A<azimuth>/<tolerance>");
161 GMT_Usage (API, -2, "Only write profile if mean direction is within +/- <tolerance> "
162 "degrees of <azimuth> [Default = All].");
163 GMT_Usage (API, 1, "\n-C<course_change>");
164 GMT_Usage (API, -2, "Profile ends when change of heading exceeds <course_change> [ignore course changes].");
165 GMT_Usage (API, 1, "\n-D<minimum_distance>");
166 GMT_Usage (API, -2, "Only write profile if length is at least <minimum_distance> [0].");
167 GMT_Usage (API, 1, "\n-F<xy_filter>/<z_filter>");
168 GMT_Usage (API, -2, "Filter the data: Give full widths of cosine arch filters for xy and z. "
169 "Use negative width for high-pass filter [Default is no filtering].");
170 GMT_Usage (API, 1, "\n-N<template>");
171 GMT_Usage (API, -2, "Write individual segments to separate files [Default writes one "
172 "multisegment file to standard output]. Append file name template which MUST "
173 "contain a C-style format code for a long integer (e.g., %%d) that represents "
174 "a sequential segment number across all tables (if more than one table) "
175 "[Default uses gmtsplit_segment_%%d.txt (or .bin for binary)]. "
176 "Alternatively, supply a template with two long format codes and we will "
177 "replace them with the table number and table segment numbers.");
178 GMT_Usage (API, 1, "\n-Q<flags>");
179 GMT_Usage (API, -2, "Indicate what output you want as one or more of xyzdh in any order, "
180 "where x,y,z refer to input data locations and optional z-value(s), "
181 "and d,h are the distance and heading along track "
182 "[Default is all fields, i.e., -Qxyzdh (or -Qxydh if no z-column in the input)].");
183 GMT_Usage (API, 1, "\n-S d,h is supplied: Input is 5 col x,y,z,d,h with d non-decreasing "
184 "[Default input is 3 col x,y,z only and computes d,h from the data].");
185 GMT_Option (API, "V,bi");
186 if (gmt_M_showusage (API)) GMT_Usage (API, -2, "Default input columns is set via -S.");
187 GMT_Option (API, "bo,d,e,f,g,h,i,q,s,:,.");
188
189 return (GMT_MODULE_USAGE);
190 }
191
parse(struct GMT_CTRL * GMT,struct SPLITXYZ_CTRL * Ctrl,struct GMT_OPTION * options)192 static int parse (struct GMT_CTRL *GMT, struct SPLITXYZ_CTRL *Ctrl, struct GMT_OPTION *options) {
193 /* This parses the options provided to gmtsplit and sets parameters in Ctrl.
194 * Note Ctrl has already been initialized and non-zero default values set.
195 * Any GMT common options will override values set previously by other commands.
196 * It also replaces any file names specified as input or output with the data ID
197 * returned when registering these sources/destinations with the API.
198 */
199
200 unsigned int j, n_errors = 0, n_outputs = 0, n_files = 0;
201 char txt_a[GMT_LEN256] = {""};
202 struct GMT_OPTION *opt = NULL;
203 struct GMTAPI_CTRL *API = GMT->parent;
204
205 for (opt = options; opt; opt = opt->next) { /* Process all the options given */
206
207 switch (opt->option) {
208
209 case '<': /* Skip input files */
210 if (GMT_Get_FilePath (API, GMT_IS_DATASET, GMT_IN, GMT_FILE_REMOTE, &(opt->arg))) n_errors++;;
211 break;
212 case '>': /* Got named output file */
213 if (n_files++ == 0) Ctrl->Out.file = strdup (opt->arg);
214 break;
215
216 /* Processes program-specific parameters */
217
218 case 'A':
219 n_errors += gmt_M_repeated_module_option (API, Ctrl->A.active);
220 Ctrl->A.active = true;
221 n_errors += gmt_M_check_condition (GMT, (sscanf(opt->arg, "%lf/%lf", &Ctrl->A.azimuth, &Ctrl->A.tolerance)) != 2,
222 "Option -A: Can't decipher values\n");
223 break;
224 case 'C':
225 n_errors += gmt_M_repeated_module_option (API, Ctrl->C.active);
226 Ctrl->C.active = true;
227 n_errors += gmt_M_check_condition (GMT, (sscanf(opt->arg, "%lf", &Ctrl->C.value)) != 1,
228 "Option -C: Can't decipher value\n");
229 break;
230 case 'D':
231 n_errors += gmt_M_repeated_module_option (API, Ctrl->D.active);
232 Ctrl->D.active = true;
233 n_errors += gmt_M_check_condition (GMT, (sscanf(opt->arg, "%lf", &Ctrl->D.value)) != 1,
234 "Option -D: Can't decipher value\n");
235 break;
236 case 'F':
237 n_errors += gmt_M_repeated_module_option (API, Ctrl->F.active);
238 Ctrl->F.active = true;
239 n_errors += gmt_M_check_condition (GMT, (sscanf(opt->arg, "%lf/%lf", &Ctrl->F.xy_filter, &Ctrl->F.z_filter)) != 2,
240 "Option -F: Can't decipher values\n");
241 break;
242 case 'G':
243 if (gmt_M_compat_check (GMT, 4)) {
244 GMT_Report (API, GMT_MSG_COMPAT, "-G option is deprecated; use -g instead.\n");
245 GMT->common.g.active = true;
246 if (gmt_M_is_geographic (GMT, GMT_IN))
247 sprintf (txt_a, "D%sk", opt->arg); /* Hardwired to be km */
248 else
249 sprintf (txt_a, "d%s", opt->arg); /* Cartesian */
250 n_errors += gmt_parse_g_option (GMT, txt_a);
251 }
252 else
253 n_errors += gmt_default_error (GMT, opt->option);
254 break;
255 case 'M':
256 if (gmt_M_compat_check (GMT, 4)) {
257 GMT_Report (API, GMT_MSG_COMPAT, "Option -M is deprecated; -fg was set instead, use this in the future.\n");
258 if (gmt_M_is_cartesian (GMT, GMT_IN)) gmt_parse_common_options (GMT, "f", 'f', "g"); /* Set -fg unless already set */
259 }
260 else
261 n_errors += gmt_default_error (GMT, opt->option);
262 break;
263 case 'N':
264 n_errors += gmt_M_repeated_module_option (API, Ctrl->N.active);
265 Ctrl->N.active = true;
266 if (opt->arg[0]) Ctrl->N.name = strdup (opt->arg);
267 break;
268 case 'Q':
269 n_errors += gmt_M_repeated_module_option (API, Ctrl->Q.active);
270 Ctrl->Q.active = true;
271 for (j = 0; opt->arg[j]; j++) {
272 if (j < SPLITXYZ_N_OUTPUT_CHOICES) {
273 Ctrl->Q.col[j] = opt->arg[j];
274 if (!strchr ("xyzdh", Ctrl->Q.col[j])) {
275 GMT_Report (API, GMT_MSG_ERROR, "Option -Q: Unrecognized output choice %c\n", Ctrl->Q.col[j]);
276 n_errors++;
277 }
278 if (opt->arg[j] == 'z') Ctrl->Q.z_selected = true;
279 n_outputs++;
280 }
281 else {
282 GMT_Report (API, GMT_MSG_ERROR, "Option -Q: Too many output columns selected: Choose from -Qxyzdg\n");
283 n_errors++;
284 }
285 }
286 break;
287 case 'S':
288 n_errors += gmt_M_repeated_module_option (API, Ctrl->S.active);
289 Ctrl->S.active = true;
290 break;
291 case 'Z':
292 if (gmt_M_compat_check (GMT, 5)) /* Warn and pass through */
293 GMT_Report (API, GMT_MSG_COMPAT, "-Z option is deprecated and not longer required.\n");
294 break;
295
296 default: /* Report bad options */
297 n_errors += gmt_default_error (GMT, opt->option);
298 break;
299 }
300 }
301
302 n_errors += gmt_M_check_condition (GMT, Ctrl->D.value < 0.0, "Option -D: Minimum segment distance must be positive\n");
303 n_errors += gmt_M_check_condition (GMT, Ctrl->C.value <= 0.0, "Option -C: Course change tolerance must be positive\n");
304 n_errors += gmt_M_check_condition (GMT, Ctrl->A.tolerance < 0.0, "Option -A: Azimuth tolerance must be positive\n");
305 n_errors += gmt_M_check_condition (GMT, GMT->common.b.active[GMT_OUT] && !Ctrl->N.name,
306 "Binary output requires a namestem in -N\n");
307 n_errors += gmt_check_binary_io (GMT, (Ctrl->S.active) ? 5 : 3);
308 n_errors += gmt_M_check_condition (GMT, n_files > 1, "Only one output destination can be specified\n");
309 n_errors += gmt_M_check_condition (GMT, Ctrl->N.active && Ctrl->N.name && !strstr (Ctrl->N.name, "%"),
310 "Option -N: Output template must contain %%d\n");
311
312 return (n_errors ? GMT_PARSE_ERROR : GMT_NOERROR);
313 }
314
315 #define bailout(code) {gmt_M_free_options (mode); return (code);}
316 #define Return(code) {Free_Ctrl (GMT, Ctrl); gmt_end_module (GMT, GMT_cpy); bailout (code);}
317
GMT_gmtsplit(void * V_API,int mode,void * args)318 EXTERN_MSC int GMT_gmtsplit (void *V_API, int mode, void *args) {
319 bool ok, first = true, no_z_column;
320
321 unsigned int i, j, d_col, h_col, z_cols, xy_cols[2] = {0, 1}, io_mode = 0;
322 unsigned int output_choice[SPLITXYZ_N_OUTPUT_CHOICES], n_outputs = 0, n_in;
323
324 int error = 0;
325
326 uint64_t dim[GMT_DIM_SIZE] = {1, 0, 0, 0}; /* Output dataset will have one table only */
327 uint64_t tbl, n_out = 0, k, n, row, seg, seg2 = 0, begin, end, n_total = 0, n_columns = 0, nprofiles = 0, *rec = NULL;
328
329 size_t n_alloc_seg = 0, n_alloc = 0;
330
331 double dy, dx, last_c, last_s, csum, ssum, this_c, this_s, dotprod, mean_azim, *fwork = NULL;
332
333 char header[GMT_LEN64] = {""};
334
335 struct GMT_DATASET *D[2] = {NULL, NULL};
336 struct GMT_DATATABLE *T = NULL;
337 struct GMT_DATASEGMENT *S = NULL, *S_out = NULL;
338 struct GMT_DATASEGMENT_HIDDEN *SH = NULL;
339 struct SPLITXYZ_CTRL *Ctrl = NULL;
340 struct GMT_CTRL *GMT = NULL, *GMT_cpy = NULL;
341 struct GMT_OPTION *options = NULL;
342 struct GMTAPI_CTRL *API = gmt_get_api_ptr (V_API); /* Cast from void to GMTAPI_CTRL pointer */
343
344 /*----------------------- Standard module initialization and parsing ----------------------*/
345
346 if (API == NULL) return (GMT_NOT_A_SESSION);
347 if (mode == GMT_MODULE_PURPOSE) return (usage (API, GMT_MODULE_PURPOSE)); /* Return the purpose of program */
348 options = GMT_Create_Options (API, mode, args); if (API->error) return (API->error); /* Set or get option list */
349
350 if ((error = gmt_report_usage (API, options, 0, usage)) != GMT_NOERROR) bailout (error); /* Give usage if requested */
351
352 /* Parse the command-line arguments; return if errors are encountered */
353
354 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 */
355 if (GMT_Parse_Common (API, THIS_MODULE_OPTIONS, options)) Return (API->error);
356 Ctrl = New_Ctrl (GMT); /* Allocate and initialize a new control structure */
357 if ((error = parse (GMT, Ctrl, options)) != GMT_NOERROR) Return (error);
358
359 /*---------------------------- This is the gmtsplit main code ----------------------------*/
360
361 GMT_Report (API, GMT_MSG_INFORMATION, "Processing input table data\n");
362 n_in = (Ctrl->S.active) ? 5 : 3;
363 if ((error = GMT_Set_Columns (API, GMT_IN, n_in, GMT_COL_VAR)) != GMT_NOERROR) {
364 Return (error);
365 }
366 if (GMT_Init_IO (API, GMT_IS_DATASET, GMT_IS_LINE, GMT_IN, GMT_ADD_DEFAULT, 0, options) != GMT_NOERROR) { /* Establishes data input */
367 Return (API->error);
368 }
369 if ((D[GMT_IN] = GMT_Read_Data (API, GMT_IS_DATASET, GMT_IS_FILE, 0, GMT_READ_FILEBREAK, NULL, NULL, NULL)) == NULL) {
370 Return (API->error);
371 }
372 if (D[GMT_IN]->n_columns < (n_in-1)) {
373 GMT_Report (API, GMT_MSG_ERROR, "Input data have %d column(s) but at least %u are needed\n", (int)D[GMT_IN]->n_columns, n_in);
374 Return (GMT_DIM_TOO_SMALL);
375 }
376
377 gmt_M_memset (output_choice, SPLITXYZ_N_OUTPUT_CHOICES, int);
378 no_z_column = (D[GMT_IN]->n_columns == 2);
379
380 if (no_z_column && Ctrl->S.active) {
381 GMT_Report (API, GMT_MSG_ERROR, "The -S option requires a 3rd data column\n");
382 Return (GMT_PARSE_ERROR);
383 }
384 if (no_z_column && Ctrl->F.z_filter != 0.0) {
385 GMT_Report (API, GMT_MSG_ERROR, "The -F option requires a 3rd data column\n");
386 Return (GMT_PARSE_ERROR);
387 }
388 if (Ctrl->Q.z_selected && no_z_column) {
389 GMT_Report (API, GMT_MSG_ERROR, "Option -Q: Cannot request z if unless data have a 3rd column\n");
390 Return (GMT_PARSE_ERROR);
391 }
392
393 /* Determine output choices and order */
394 for (k = n_outputs = 0; k < SPLITXYZ_N_OUTPUT_CHOICES && Ctrl->Q.col[k]; k++) {
395 switch (Ctrl->Q.col[k]) {
396 case 'x':
397 output_choice[k] = GMT_X;
398 break;
399 case 'y':
400 output_choice[k] = GMT_Y;
401 break;
402 case 'z':
403 if (no_z_column) {
404 GMT_Report (API, GMT_MSG_ERROR, "Cannot specify z when there is no z column!\n");
405 Return (-1);
406 }
407 output_choice[k] = GMT_Z;
408 break;
409 case 'd':
410 output_choice[k] = 3 - no_z_column;
411 break;
412 case 'h':
413 output_choice[k] = 4 - no_z_column;
414 break;
415 }
416 n_outputs++;
417 }
418 if (gmt_M_is_geographic (GMT, GMT_IN))
419 gmt_set_geographic (GMT, GMT_OUT);
420
421 if (n_outputs == 0) { /* No -Q given, generate default -Q setting (all) */
422 n_outputs = 5 - no_z_column;
423 for (i = 0; i < 2; i++) output_choice[i] = i;
424 if (!no_z_column) output_choice[2] = 2;
425 for (i = 3-no_z_column; i < n_outputs; i++) output_choice[i] = i;
426 }
427
428 /* Convert to radians */
429 Ctrl->A.tolerance *= D2R;
430 Ctrl->A.azimuth = D2R * (90.0 - Ctrl->A.azimuth); /* Work in Cartesian angle and radians */
431 Ctrl->C.value *= D2R;
432 if (Ctrl->F.active) fwork = gmtsplit_filterxy_setup (GMT);
433 if (!Ctrl->N.active)
434 gmt_set_segmentheader (GMT, GMT_OUT, true); /* Turn on segment headers on output */
435
436 if ((error = GMT_Set_Columns (API, GMT_OUT, n_outputs, (D[GMT_IN]->type == GMT_READ_DATA) ? GMT_COL_FIX_NO_TEXT : GMT_COL_FIX)) != GMT_NOERROR) {
437 Return (error);
438 }
439 /* Registers default output destination, unless already set */
440 if (GMT_Init_IO (API, GMT_IS_DATASET, GMT_IS_PLP, GMT_OUT, GMT_ADD_DEFAULT, 0, options) != GMT_NOERROR) {
441 Return (API->error);
442 }
443
444 if (!Ctrl->S.active) { /* Must extend table with 2 cols to hold distance and azimuth */
445 n_columns = D[GMT_IN]->n_columns + 2;
446 d_col = (unsigned int)D[GMT_IN]->n_columns;
447 h_col = d_col + 1;
448 gmt_adjust_dataset (GMT, D[GMT_IN], n_columns);
449 }
450 else { /* Comes with d and az in file */
451 d_col = no_z_column + 2;
452 h_col = no_z_column + 3;
453 }
454 z_cols = 2;
455 S_out = gmt_get_segment (GMT, n_outputs); /* We will use a single segment to hold all output records and keep track of start/stop via rec[] array */
456
457 nprofiles = 0;
458 for (tbl = 0; tbl < D[GMT_IN]->n_tables; tbl++) {
459 T = D[GMT_IN]->table[tbl]; /* Shorthand for current table */
460 if (gmt_M_is_verbose (GMT, GMT_MSG_INFORMATION)) {
461 struct GMT_DATATABLE_HIDDEN *TH = gmt_get_DT_hidden (T);
462 GMT_Report (API, GMT_MSG_INFORMATION, "Working on file %s\n", TH->file[GMT_IN]);
463 }
464
465 for (seg = 0; seg < T->n_segments; seg++) { /* For each segment in this table */
466 S = T->segment[seg];
467 if (Ctrl->S.active) S->data[h_col][0] = D2R * (90.0 - S->data[h_col][0]); /* Angles are stored as CCW angles in radians */
468 for (row = 1; row < S->n_rows; row++) {
469 if (!Ctrl->S.active) { /* Must extend table with 2 cols to hold distance and azimuth */
470 dy = S->data[GMT_Y][row] - S->data[GMT_Y][row-1];
471 if (gmt_M_x_is_lon (GMT, GMT_IN)) { /* We do flat earth approximation for distances, which assumes points are close to each other */
472 gmt_M_set_delta_lon (S->data[GMT_X][row-1], S->data[GMT_X][row], dx);
473 dy *= GMT->current.proj.DIST_KM_PR_DEG;
474 dx *= (GMT->current.proj.DIST_KM_PR_DEG * cosd (0.5 * (S->data[GMT_Y][row] + S->data[GMT_Y][row-1])));
475 }
476 else
477 dx = S->data[GMT_X][row] - S->data[GMT_X][row-1];
478 if (dy == 0.0 && dx == 0.0) {
479 S->data[d_col][row] = S->data[d_col][row-1];
480 S->data[h_col][row] = S->data[h_col][row-1];
481 }
482 else {
483 S->data[d_col][row] = S->data[d_col][row-1] + hypot (dx,dy);
484 S->data[h_col][row] = d_atan2(dy,dx); /* Angles are stored as CCW angles in radians */
485 }
486 }
487 else
488 S->data[h_col][row] = D2R * (90.0 - S->data[h_col][row]); /* Angles are stored as CCW angles in radians */
489 }
490 if (!Ctrl->S.active) S->data[h_col][0] = S->data[h_col][1];
491
492 /* Here a complete segment is ready for further processing */
493 /* Now we have read the data and can filter z, if necessary. */
494
495 if (Ctrl->F.active)
496 gmtsplit_filter_cols (GMT, S->data, 0, S->n_rows, d_col, 1, &z_cols, Ctrl->F.z_filter, fwork);
497
498 /* Now we are ready to search for segments */
499
500 begin = end = 0;
501 while (end < S->n_rows-1) {
502 sincos (S->data[h_col][begin], &last_s, &last_c);
503 csum = last_c; ssum = last_s;
504 ok = true;
505 while (ok && end < S->n_rows-1) {
506 end++;
507 sincos (S->data[h_col][end], &this_s, &this_c);
508 dotprod = this_c * last_c + this_s * last_s;
509 if (fabs (dotprod) > 1.0) dotprod = copysign (1.0, dotprod);
510 if (d_acos (dotprod) > Ctrl->C.value) { /* Fails due to too much change in azimuth */
511 ok = false;
512 continue;
513 }
514 /* Get here when this point belongs with last one */
515 csum += this_c;
516 ssum += this_s;
517 last_c = this_c;
518 last_s = this_s;
519 }
520
521 /* Get here when we have found a beginning and end */
522
523 if (ok) end++; /* Last point in input should be included in this segment */
524
525 if (end - begin - 1) { /* There are at least two points in the list. */
526 if ((S->data[d_col][end-1] - S->data[d_col][begin]) >= Ctrl->D.value) {
527 /* List is long enough. Check strike. Compute mean_azim in range [-pi/2, pi/2] */
528
529 mean_azim = d_atan2 (ssum, csum);
530 mean_azim = fabs (mean_azim - Ctrl->A.azimuth);
531 if (mean_azim <= Ctrl->A.tolerance) { /* List has acceptable strike */
532 if (Ctrl->F.active)
533 gmtsplit_filter_cols (GMT, S->data, begin, end, d_col, 2, xy_cols, Ctrl->F.xy_filter, fwork);
534 nprofiles++;
535
536 n_out = end - begin;
537 if ((n_total + n_out) >= n_alloc) { /* Allocate more memory for temporary segment */
538 n_alloc = (first) ? D[GMT_IN]->n_records : n_alloc * 2;
539 gmt_alloc_segment (GMT, S_out, n_alloc, n_outputs, 0U, first);
540 first = false;
541 }
542
543 for (row = begin; row < end; row++, k++) {
544 for (j = 0; j < n_outputs; j++) { /* Remember to convert CCW angles back to azimuths */
545 S_out->data[j][k] = (output_choice[j] == h_col) ? 90.0 - R2D * S->data[h_col][row] : S->data[output_choice[j]][row];
546 }
547 }
548 if (seg2 == n_alloc_seg) {
549 n_alloc_seg = (n_alloc_seg == 0) ? D[GMT_IN]->n_segments : n_alloc_seg * 2;
550 rec = gmt_M_memory (GMT, rec, n_alloc_seg, uint64_t);
551 }
552 rec[seg2++] = k;
553 n_total += n_out;
554 }
555 }
556 }
557 begin = end;
558 }
559 }
560 }
561
562 if (!Ctrl->S.active) /* Must remove the 2 cols we added */
563 gmt_adjust_dataset (GMT, D[GMT_IN], n_columns-2);
564
565 if (Ctrl->F.active) gmt_M_free (GMT, fwork);
566
567 /* Get here when all profiles have been found and written */
568
569 if (nprofiles > 1)
570 gmt_set_segmentheader (GMT, GMT_OUT, true); /* Turn on segment headers on output */
571
572 dim[GMT_SEG] = seg2; dim[GMT_COL] = n_outputs; /* dim[GMT_ROW] is zero so segment pointers have no data */
573 if ((D[GMT_OUT] = GMT_Create_Data (API, GMT_IS_DATASET, GMT_IS_LINE, 0, dim, NULL, NULL, 0, 0, NULL)) == NULL) {
574 gmt_M_free (GMT, rec);
575 goto end;
576 }
577 for (seg = 0; seg < seg2; seg++) { /* We fake a table by setting the data pointers to point to various points in our single S_out arrays */
578 S = D[GMT_OUT]->table[0]->segment[seg];
579 SH = gmt_get_DS_hidden (S);
580 k = (seg == 0) ? 0 : rec[seg-1]; /* Record number of first row in the new segment */
581 n = (seg == 0) ? rec[seg] : rec[seg] - rec[seg-1]; /* How many rows in the new segment */
582 if (gmt_alloc_segment (GMT, S, n, n_outputs, S->n_columns, true)) {
583 GMT_Report (API, GMT_MSG_ERROR, "Failed to allocate output memory for a new segment");
584 goto end;
585 }
586
587 for (j = 0; j < n_outputs; j++) /* Duplicate the rows */
588 gmt_M_memcpy (S->data[j], &S_out->data[j][k], n, double);
589 S->n_rows = n;
590 sprintf (header, "Profile %" PRIu64" -I%" PRIu64, seg, seg);
591 S->header = strdup (header);
592 SH->id = seg;
593 }
594 gmt_M_free (GMT, rec);
595
596 GMT_Report (API, GMT_MSG_INFORMATION, " Split %" PRIu64 " data into %" PRIu64 " segments.\n", D[GMT_IN]->n_records, nprofiles);
597 if (Ctrl->N.active) { /* Want tables or segments written to separate files */
598 int n_formats = 0;
599 for (k = 0; Ctrl->N.name[k]; k++) if (Ctrl->N.name[k] == '%') n_formats++;
600 io_mode = (n_formats == 2) ? GMT_WRITE_TABLE_SEGMENT: GMT_WRITE_SEGMENT;
601 /* The io_mode tells the i/o function to split segments into files */
602 gmt_M_str_free (Ctrl->Out.file);
603 Ctrl->Out.file = strdup (Ctrl->N.name);
604 }
605 if (GMT_Write_Data (API, GMT_IS_DATASET, GMT_IS_FILE, GMT_IS_PLP, io_mode, NULL, Ctrl->Out.file, D[GMT_OUT]) != GMT_NOERROR) {
606 goto end;
607 }
608
609 if (GMT_End_IO (API, GMT_OUT, 0) != GMT_NOERROR) { /* Disables further data output */
610 goto end;
611 }
612
613 end: /* Clean up memory and return */
614
615 gmt_free_segment (GMT, &S_out);
616
617 Return (API->error);
618 }
619
GMT_splitxyz(void * V_API,int mode,void * args)620 EXTERN_MSC int GMT_splitxyz (void *V_API, int mode, void *args) {
621 /* This was the GMT6.1.1 name */
622 struct GMTAPI_CTRL *API = gmt_get_api_ptr (V_API); /* Cast from void to GMTAPI_CTRL pointer */
623 if (gmt_M_compat_check (API->GMT, 6)) {
624 GMT_Report (API, GMT_MSG_COMPAT, "Module splitxyz is deprecated; use gmtsplit.\n");
625 return (GMT_Call_Module (API, "gmtsplit", mode, args));
626 }
627 GMT_Report (API, GMT_MSG_ERROR, "Shared GMT module not found: splitxyz\n");
628 return (GMT_NOT_A_VALID_MODULE);
629 }
630