/*----------------------------------------------------------------- * $Id: x2sys_cross.c 12822 2014-01-31 23:39:56Z remko $ * * Copyright (c) 1999-2014 by P. Wessel * See LICENSE.TXT file for copying and redistribution conditions. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by * the Free Software Foundation; version 3 or any later version. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU Lesser General Public License for more details. * * Contact info: www.soest.hawaii.edu/pwessel *--------------------------------------------------------------------*/ /* x2sys_cross will calculate crossovers generated by the * intersections of two tracks. Optionally, it will also evaluate * the interpolated datafields at the crossover locations. * * Author: Paul Wessel * Date: 15-JUN-2004 * Version: 1.0, based on the spirit of the old xsystem code, * but with a smarter algorithm based on the book * "Algorithms in C" by R. Sedgewick. * 31-MAR-2006: Changed -O to -L to avoid clash with GMT. * */ #define THIS_MODULE_NAME "x2sys_cross" #define THIS_MODULE_LIB "x2sys" #define THIS_MODULE_PURPOSE "Calculate crossovers between track data files" #include "x2sys.h" #define GMT_PROG_OPTIONS "->JRVb" /* Control structure for x2sys_cross */ #define HHI 0 #define VLO 1 #define VHI 2 struct X2SYS_CROSS_CTRL { struct X2S_CROSS_A { /* -A */ bool active; char *file; } A; struct X2S_CROSS_C { /* -C */ bool active; char *file; } C; struct X2S_CROSS_I { /* -I */ bool active; int mode; } I; struct X2S_CROSS_S { /* -S */ bool active[2]; double limit[3]; } S; struct X2S_CROSS_T { /* -T */ bool active; char *TAG; } T; struct X2S_CROSS_W { /* -W */ bool active; unsigned int width; } W; struct X2S_CROSS_Q { /* -Q */ bool active; int mode; } Q; struct X2S_CROSS_Z { /* -Z */ bool active; } Z; }; struct PAIR { /* Used with -Kkombinations.lis option */ char *id1, *id2; }; void *New_x2sys_cross_Ctrl (struct GMT_CTRL *GMT) { /* Allocate and initialize a new control structure */ struct X2SYS_CROSS_CTRL *C; C = GMT_memory (GMT, NULL, 1, struct X2SYS_CROSS_CTRL); /* Initialize values whose defaults are not 0/false/NULL */ C->S.limit[VHI] = DBL_MAX; /* Ignore crossovers on segments that implies speed higher than this */ C->W.width = 3; /* Number of points on either side in the interpolation */ return (C); } void Free_x2sys_cross_Ctrl (struct GMT_CTRL *GMT, struct X2SYS_CROSS_CTRL *C) { /* Deallocate control structure */ if (!C) return; if (C->A.file) free (C->A.file); if (C->C.file) free (C->C.file); if (C->T.TAG) free (C->T.TAG); GMT_free (GMT, C); } int GMT_x2sys_cross_usage (struct GMTAPI_CTRL *API, int level) { GMT_show_name_and_purpose (API, THIS_MODULE_LIB, THIS_MODULE_NAME, THIS_MODULE_PURPOSE); if (level == GMT_MODULE_PURPOSE) return (GMT_NOERROR); GMT_Message (API, GMT_TIME_NONE, "usage: x2sys_cross -T [-A] [-C[]] [-Il|a|c] [%s] [-Qe|i]\n", GMT_J_OPT); GMT_Message (API, GMT_TIME_NONE, "\t[%s] [-Sl|h|u] [%s] [-W] [-Z]\n", GMT_Rgeo_OPT, GMT_V_OPT); GMT_Message (API, GMT_TIME_NONE, "\t[%s]\n\n", GMT_bo_OPT); GMT_Message (API, GMT_TIME_NONE, "\tOutput is x y t1 t2 d1 d2 az1 az2 v1 v2 xval1 xmean1 xval2 xmean2 ...\n"); GMT_Message (API, GMT_TIME_NONE, "\tIf time is not selected (or present) we use record numbers as proxies i1 i2\n"); if (level == GMT_SYNOPSIS) return (EXIT_FAILURE); GMT_Message (API, GMT_TIME_NONE, "\t is one or more datafiles, or give = for a file with a list of datafiles.\n"); GMT_Message (API, GMT_TIME_NONE, "\t-T is the system tag for the data set.\n"); GMT_Message (API, GMT_TIME_NONE, "\n\tOPTIONS:\n"); GMT_Message (API, GMT_TIME_NONE, "\t-A Give list of file pairs that are ok to compare [Default is all combinations].\n"); GMT_Message (API, GMT_TIME_NONE, "\t-C Print run time for each pair. Optionally append to save them in file.\n"); GMT_Message (API, GMT_TIME_NONE, "\t-I Set the interpolation mode. Choose among:\n"); GMT_Message (API, GMT_TIME_NONE, "\t l Linear interpolation [Default].\n"); GMT_Message (API, GMT_TIME_NONE, "\t a Akima spline interpolation.\n"); GMT_Message (API, GMT_TIME_NONE, "\t c Acubic spline interpolation.\n"); GMT_Option (API, "J-"); GMT_Message (API, GMT_TIME_NONE, "\t-Q Append e for external crossovers.\n"); GMT_Message (API, GMT_TIME_NONE, "\t Append i for internal crossovers [Default is all crossovers].\n"); GMT_Option (API, "R"); GMT_Message (API, GMT_TIME_NONE, "\t-S Set limits on lower and upper speeds (units determined by -Ns):\n"); GMT_Message (API, GMT_TIME_NONE, "\t -Sl sets lower speed [Default is 0].\n"); GMT_Message (API, GMT_TIME_NONE, "\t -Sh no headings should be computed if velocity drops below this value [0].\n"); GMT_Message (API, GMT_TIME_NONE, "\t -Su sets upper speed [Default is Infinity].\n"); GMT_Option (API, "V"); GMT_Message (API, GMT_TIME_NONE, "\t-W Set maximum points on either side of xover to use in interpolation [Default is 3].\n"); GMT_Message (API, GMT_TIME_NONE, "\t-Z Return z-values for each track [Default is crossover and mean value].\n"); GMT_Option (API, "bo,."); return (EXIT_FAILURE); } int GMT_x2sys_cross_parse (struct GMT_CTRL *GMT, struct X2SYS_CROSS_CTRL *Ctrl, struct GMT_OPTION *options) { /* This parses the options provided to grdcut and sets parameters in CTRL. * Any GMT common options will override values set previously by other commands. * It also replaces any file names specified as input or output with the data ID * returned when registering these sources/destinations with the API. */ unsigned int n_errors = 0; struct GMT_OPTION *opt = NULL; struct GMTAPI_CTRL *API = GMT->parent; for (opt = options; opt; opt = opt->next) { /* Process all the options given */ switch (opt->option) { /* Common parameters */ case '<': /* Skip input files */ if (!GMT_check_filearg (GMT, '<', opt->arg, GMT_IN)) n_errors++; break; /* Processes program-specific parameters */ case 'A': /* Get list of approved filepair combinations to check */ if ((Ctrl->A.active = GMT_check_filearg (GMT, 'A', opt->arg, GMT_IN))) Ctrl->A.file = strdup (opt->arg); else n_errors++; break; case 'C': Ctrl->C.active = true; if (strlen(opt->arg)) Ctrl->C.file = strdup (opt->arg); break; case 'I': Ctrl->I.active = true; switch (opt->arg[0]) { case 'l': Ctrl->I.mode = 0; break; case 'a': Ctrl->I.mode = 1; break; case 'c': Ctrl->I.mode = 2; break; case 'n': Ctrl->I.mode = 3; break; default: n_errors++; break; } break; case 'S': /* Speed checks */ switch (opt->arg[0]) { case 'L': case 'l': /* Lower cutoff speed */ Ctrl->S.limit[VLO] = atof (&opt->arg[1]); Ctrl->S.active[VLO] = true; break; case 'U': case 'u': /* Upper cutoff speed */ Ctrl->S.limit[VHI] = atof (&opt->arg[1]); Ctrl->S.active[VLO] = true; break; case 'H': case 'h': /* Heading calculation cutoff speed */ Ctrl->S.limit[HHI] = atof (&opt->arg[1]); Ctrl->S.active[HHI] = true; break; default: GMT_Report (API, GMT_MSG_NORMAL, "Syntax error: -S\n"); n_errors++; break; } break; case 'T': Ctrl->T.active = true; Ctrl->T.TAG = strdup (opt->arg); break; case 'W': /* Get new window half-width as number of points */ Ctrl->W.active = true; Ctrl->W.width = atoi (opt->arg); break; case 'Q': /* Specify internal or external only */ Ctrl->Q.active = true; if (opt->arg[0] == 'e') Ctrl->Q.mode = 1; else if (opt->arg[0] == 'i') Ctrl->Q.mode = 2; else Ctrl->Q.mode = 3; break; case 'Z': /* Return z1, z1 rather than (z1-z1) and 0.5 * (z1 + z2) */ Ctrl->Z.active = true; break; default: /* Report bad options */ n_errors += GMT_default_error (GMT, opt->option); break; } } n_errors += GMT_check_condition (GMT, !Ctrl->T.active || !Ctrl->T.TAG, "Syntax error: -T must be used to set the TAG\n"); n_errors += GMT_check_condition (GMT, Ctrl->W.width < 1, "Syntax error: Error -W: window must be at least 1\n"); n_errors += GMT_check_condition (GMT, Ctrl->S.limit[VLO] > Ctrl->S.limit[VHI], "Syntax error: Error -S: lower speed cutoff higher than upper cutoff!\n"); n_errors += GMT_check_condition (GMT, Ctrl->Q.mode == 3, "Syntax error: Error -Q: Only one of -Qe -Qi can be specified!\n"); return (n_errors ? GMT_PARSE_ERROR : GMT_OK); } int combo_ok (char *name_1, char *name_2, struct PAIR *pair, uint64_t n_pairs) { uint64_t i; /* Return true if this particular combination is found in the list of pairs */ for (i = 0; i < n_pairs; i++) { if (!(strcmp (name_1, pair[i].id1) || strcmp (name_2, pair[i].id2))) return (true); if (!(strcmp (name_2, pair[i].id1) || strcmp (name_1, pair[i].id2))) return (true); } return (false); } #define bailout(code) {GMT_Free_Options (mode); return (code);} #define Return(code) {Free_x2sys_cross_Ctrl (GMT, Ctrl); GMT_end_module (GMT, GMT_cpy); bailout (code);} int GMT_x2sys_cross (void *V_API, int mode, void *args) { char **trk_name = NULL; /* Name of tracks */ char line[GMT_BUFSIZ] = {""}; /* buffer */ char item[GMT_BUFSIZ] = {""}; /* buffer */ char t_or_i; /* t = time, i = dummy node time */ char name1[80] = {""}, name2[80] = {""}; /* Name of two files to be examined */ char *x2sys_header = "%s %d %s %d %s"; uint64_t n_rec[2]; /* Number of data records for both files */ uint64_t window_width; /* Max number of points to use in the interpolation */ uint64_t n_tracks = 0; /* Total number of data sets to compare */ uint64_t nx; /* Number of crossovers found for this pair */ uint64_t *col_number = NULL; /* Array with the column numbers of the data fields */ unsigned int n_output; /* Number of columns on output */ uint64_t n_pairs = 0; /* Number of acceptable combinations */ uint64_t A, B, i, j, col, k, start, n_bad; /* Misc. counters and local variables */ uint64_t end, first, n_ok; uint64_t n_data_col, left[2], t_left; uint64_t n_left, right[2], t_right, n_right; uint64_t n_duplicates, n_errors; uint64_t add_chunk; int scol; int error = 0; /* nonzero for invalid arguments */ unsigned int *ok = NULL; bool xover_locations_only = false; /* true if only x,y (and possible indices) to be output */ bool internal = true; /* false if only external xovers are needed */ bool external = true; /* false if only internal xovers are needed */ bool do_project = false; /* true if we must mapproject first */ bool got_time = false; /* true if there is a time column */ bool first_header = true; /* true for very first crossover */ bool first_crossover; /* true for first crossover between two data sets */ bool same = false; /* true when the two cruises we compare have the same name */ bool has_time[2]; /* true for each cruises that actually has a time column */ bool *duplicate = NULL; /* Array, true for any cruise that is already listed */ bool cmdline_files = false; /* true if files where given directly on the command line */ bool wrap = false; /* true if data wraps so -Rg was given */ size_t n_alloc = 1; double dt; /* Time between crossover and previous node */ double dist_x[2]; /* Distance(s) along track at the crossover point */ double time_x[2]; /* Time(s) along track at the crossover point */ double deld, delt; /* Differences in dist and time across xpoint */ double speed[2]; /* speed across the xpoint ( = deld/delt) */ double **data[2] = {NULL, NULL}; /* Data matrices for the two data sets to be checked */ double *xdata[2] = {NULL, NULL}; /* Data vectors with estimated values at crossover points */ double *dist[2] = {NULL, NULL}; /* Data vectors with along-track distances */ double *time[2] = {NULL, NULL}; /* Data vectors with along-track times (or dummy node indices) */ double *t = NULL, *y = NULL; /* Interpolation y(t) arrays */ double *out = NULL; /* Output record array */ double X2SYS_NaN; /* Value to write out when result is NaN */ double xx, yy; /* Temporary projection variables */ double dist_scale; /* Scale to give selected distance units */ double vel_scale; /* Scale to give selected velocity units */ double t_scale; /* Scale to give time in seconds */ clock_t tic = 0, toc = 0; struct X2SYS_INFO *s = NULL; /* Data format information */ struct GMT_XSEGMENT *ylist_A = NULL, *ylist_B = NULL; /* y-indices sorted in increasing order */ struct GMT_XOVER XC; /* Structure with resulting crossovers */ struct X2SYS_FILE_INFO data_set[2]; /* File information */ struct X2SYS_BIX Bix; struct PAIR *pair = NULL; /* Used with -Akombinations.lis option */ FILE *fp = NULL, *fpC = NULL; struct X2SYS_CROSS_CTRL *Ctrl = NULL; struct GMT_CTRL *GMT = NULL, *GMT_cpy = NULL; struct GMT_OPTION *options = NULL; struct GMTAPI_CTRL *API = GMT_get_API_ptr (V_API); /* Cast from void to GMTAPI_CTRL pointer */ /*----------------------------------END OF VARIBLE DECLARATIONS---------------------------------------------*/ /*----------------------- Standard module initialization and parsing ----------------------*/ if (API == NULL) return (GMT_NOT_A_SESSION); if (mode == GMT_MODULE_PURPOSE) return (GMT_x2sys_cross_usage (API, GMT_MODULE_PURPOSE)); /* Return the purpose of program */ options = GMT_Create_Options (API, mode, args); if (API->error) return (API->error); /* Set or get option list */ if (!options || options->option == GMT_OPT_USAGE) bailout (GMT_x2sys_cross_usage (API, GMT_USAGE)); /* Return the usage message */ if (options->option == GMT_OPT_SYNOPSIS) bailout (GMT_x2sys_cross_usage (API, GMT_SYNOPSIS)); /* Return the synopsis */ /* Parse the command-line arguments */ GMT = GMT_begin_module (API, THIS_MODULE_LIB, THIS_MODULE_NAME, &GMT_cpy); /* Save current state */ if (GMT_Parse_Common (API, GMT_PROG_OPTIONS, options)) Return (API->error); Ctrl = New_x2sys_cross_Ctrl (GMT); /* Allocate and initialize a new control structure */ if ((error = GMT_x2sys_cross_parse (GMT, Ctrl, options))) Return (error); /*---------------------------- This is the x2sys_cross main code ----------------------------*/ x2sys_err_fail (GMT, x2sys_set_system (GMT, Ctrl->T.TAG, &s, &Bix, &GMT->current.io), Ctrl->T.TAG); if (!s->geographic) GMT->current.io.col_type[GMT_IN][GMT_X] = GMT->current.io.col_type[GMT_IN][GMT_Y] = GMT->current.io.col_type[GMT_OUT][GMT_X] = GMT->current.io.col_type[GMT_OUT][GMT_Y] = GMT_IS_UNKNOWN; if (s->x_col == -1 || s->y_col == -1) { GMT_Report (API, GMT_MSG_NORMAL, "Error: lon,lat or x,y are not among data columns!\n"); Return (EXIT_FAILURE); } if ((n_tracks = x2sys_get_tracknames (GMT, options, &trk_name, &cmdline_files)) == 0) { GMT_Report (API, GMT_MSG_NORMAL, "Error: Must give at least one data set!\n"); Return (EXIT_FAILURE); } GMT->current.setting.interpolant = Ctrl->I.mode; if (Ctrl->Q.active) { if (Ctrl->Q.mode == 1) internal = false; if (Ctrl->Q.mode == 2) external = false; } GMT_Report (API, GMT_MSG_VERBOSE, "Files found: %" PRIu64 "\n", n_tracks); duplicate = GMT_memory (GMT, NULL, n_tracks, bool); GMT_Report (API, GMT_MSG_VERBOSE, "Checking for duplicates : "); /* Make sure there are no duplicates */ for (A = n_duplicates = 0; A < n_tracks; A++) { /* Loop over all files */ if (duplicate[A]) continue; for (B = A + 1; B < n_tracks; B++) { if (duplicate[B]) continue; same = !strcmp (trk_name[A], trk_name[B]); if (same) { GMT_Report (API, GMT_MSG_NORMAL, "File %s repeated on command line - skipped\n", trk_name[A]); duplicate[B] = true; n_duplicates++; } } } GMT_Report (API, GMT_MSG_VERBOSE, "%" PRIu64 " found\n", n_duplicates); if (Ctrl->A.active) { /* Read list of acceptable trk_name combinations */ GMT_Report (API, GMT_MSG_VERBOSE, "Explicit combinations found: "); if ((fp = fopen (Ctrl->A.file, "r")) == NULL) { GMT_Report (API, GMT_MSG_NORMAL, "Error: Could not open combinations file %s!\n", Ctrl->A.file); Return (EXIT_FAILURE); } n_alloc = add_chunk = GMT_CHUNK; pair = GMT_memory (GMT, NULL, n_alloc, struct PAIR); while (fgets (line, GMT_BUFSIZ, fp)) { if (line[0] == '#' || line[0] == '\n') continue; /* Skip comments and blanks */ GMT_chop (line); /* Get rid of CR, LF stuff */ if (sscanf (line, "%s %s", name1, name2) != 2) { GMT_Report (API, GMT_MSG_NORMAL, "Error: Error decoding combinations file for pair %" PRIu64 "!\n", n_pairs); Return (EXIT_FAILURE); } pair[n_pairs].id1 = strdup (name1); pair[n_pairs].id2 = strdup (name2); n_pairs++; if (n_pairs == n_alloc) { size_t old_n_alloc = n_alloc; add_chunk *= 2; n_alloc += add_chunk; pair = GMT_memory (GMT, pair, n_alloc, struct PAIR); GMT_memset (&(pair[old_n_alloc]), n_alloc - old_n_alloc, struct PAIR); } } fclose (fp); if (!n_pairs) { GMT_Report (API, GMT_MSG_NORMAL, "Error: No combinations found in file %s!\n", Ctrl->A.file); Return (EXIT_FAILURE); } if (n_pairs < n_alloc) pair = GMT_memory (GMT, pair, n_pairs, struct PAIR); GMT_Report (API, GMT_MSG_VERBOSE, "%" PRIu64 "\n", n_pairs); } if (Ctrl->C.file) { /* Open file to store the per pair run time */ if ((fpC = fopen (Ctrl->C.file, "w")) == NULL) { GMT_Report (API, GMT_MSG_NORMAL, "Error: Could not open save times file %s!\n", Ctrl->C.file); free(Ctrl->C.file); Ctrl->C.file = NULL; } } X2SYS_NaN = GMT->session.d_NaN; if (GMT->current.setting.interpolant == 0) Ctrl->W.width = 1; window_width = 2 * Ctrl->W.width; n_data_col = x2sys_n_data_cols (GMT, s); got_time = (s->t_col >= 0); if (!got_time) Ctrl->S.active[VLO] = false; /* Cannot check speed if there is no time */ n_output = (unsigned int)(10 + 2 * n_data_col); GMT->current.io.col_type[GMT_OUT][GMT_X] = (!strcmp (s->info[s->x_col].name, "lon")) ? GMT_IS_LON : GMT_IS_FLOAT; GMT->current.io.col_type[GMT_OUT][GMT_Y] = (!strcmp (s->info[s->x_col].name, "lat")) ? GMT_IS_LAT : GMT_IS_FLOAT; GMT->current.io.col_type[GMT_OUT][GMT_Z] = GMT->current.io.col_type[GMT_OUT][3] = (got_time) ? GMT_IS_ABSTIME : GMT_IS_FLOAT; for (i = 0; i < n_data_col+2; i++) GMT->current.io.col_type[GMT_OUT][4+2*i] = GMT->current.io.col_type[GMT_OUT][5+2*i] = GMT_IS_FLOAT; if (n_data_col == 0) { xover_locations_only = true; n_output = 2; } else { /* Set the actual column numbers with data fields */ t = GMT_memory (GMT, NULL, window_width, double); y = GMT_memory (GMT, NULL, window_width, double); col_number = GMT_memory (GMT, NULL, n_data_col, uint64_t); ok = GMT_memory (GMT, NULL, n_data_col, unsigned int); for (col = k = scol = 0; col < s->n_out_columns; col++, scol++) { if (scol == s->x_col || scol == s->y_col || scol == s->t_col) continue; col_number[k++] = col; } if (s->t_col < 0) GMT_Report (API, GMT_MSG_VERBOSE, "No time column, use dummy times\n"); } out = GMT_memory (GMT, NULL, n_output, double); xdata[0] = GMT_memory (GMT, NULL, s->n_out_columns, double); xdata[1] = GMT_memory (GMT, NULL, s->n_out_columns, double); GMT_set_segmentheader (GMT, GMT_OUT, true); /* Turn on segment headers on output */ GMT->current.setting.io_header[GMT_OUT] = true; /* Turn on -ho explicitly */ if (GMT->common.R.active && GMT->current.proj.projection != GMT_NO_PROJ) { do_project = true; s->geographic = false; /* Since we then have x,y projected coordinates, not lon,lat */ s->dist_flag = 0; GMT_err_fail (GMT, GMT_map_setup (GMT, GMT->common.R.wesn), ""); } GMT_init_distaz (GMT, s->dist_flag ? GMT_MAP_DIST_UNIT : 'X', s->dist_flag, GMT_MAP_DIST); MGD77_Set_Unit (GMT, s->unit[X2SYS_DIST_SELECTION], &dist_scale, -1); /* Gets scale which multiplies meters to chosen distance unit */ MGD77_Set_Unit (GMT, s->unit[X2SYS_SPEED_SELECTION], &vel_scale, -1); /* Sets output scale for distances using in velocities */ switch (s->unit[X2SYS_SPEED_SELECTION][0]) { case 'c': vel_scale = 1.0; break; case 'e': vel_scale /= dist_scale; /* Must counteract any distance scaling to get meters. dt is in sec so we get m/s */ break; case 'f': vel_scale /= (METERS_IN_A_FOOT * dist_scale); /* Must counteract any distance scaling to get feet. dt is in sec so we get ft/s */ break; case 'k': vel_scale *= (3600.0 / dist_scale); /* Must counteract any distance scaling to get km. dt is in sec so 3600 gives km/hr */ break; case 'm': vel_scale *= (3600.0 / dist_scale); /* Must counteract any distance scaling to get miles. dt is in sec so 3600 gives miles/hr */ break; case 'n': vel_scale *= (3600.0 / dist_scale); /* Must counteract any distance scaling to get miles. dt is in sec so 3600 gives miles/hr */ break; case 'u': vel_scale /= (METERS_IN_A_SURVEY_FOOT * dist_scale); /* Must counteract any distance scaling to get survey feet. dt is in sec so we get ft/s */ break; default: /*Cartesian */ break; } t_scale = GMT->current.setting.time_system.scale; /* Convert user's TIME_UNIT to seconds */ wrap = (GMT_is_geographic (GMT, GMT_IN) && GMT->common.R.active && GMT_360_RANGE (GMT->common.R.wesn[XLO], GMT->common.R.wesn[XHI])); if ((error = GMT_set_cols (GMT, GMT_OUT, n_output)) != GMT_OK) { Return (error); } if (GMT_Init_IO (API, GMT_IS_DATASET, GMT_IS_POINT, GMT_OUT, GMT_ADD_DEFAULT, 0, options) != GMT_OK) { /* Registers default output destination, unless already set */ Return (API->error); } if (GMT_Begin_IO (API, GMT_IS_DATASET, GMT_OUT, GMT_HEADER_ON) != GMT_OK) { /* Enables data output and sets access mode */ Return (API->error); } for (A = 0; A < n_tracks; A++) { /* Loop over all files */ if (duplicate[A]) continue; if (s->x_col < 0 || s->x_col < 0) { GMT_Report (API, GMT_MSG_NORMAL, "Error: x and/or y column not found for track %s!\n", trk_name[A]); Return (EXIT_FAILURE); } x2sys_err_fail (GMT, (s->read_file) (GMT, trk_name[A], &data[0], s, &data_set[0], &GMT->current.io, &n_rec[0]), trk_name[A]); if (n_rec[0] == 0) { /* No data in track A */ x2sys_free_data (GMT, data[0], s->n_out_columns, &data_set[0]); continue; } has_time[0] = false; if (got_time) { /* Check to make sure we do in fact have time */ for (i = n_bad = 0; i < n_rec[0]; i++) n_bad += GMT_is_dnan (data[0][s->t_col][i]); if (n_bad < n_rec[0]) has_time[0] = true; } if (do_project) { /* Convert all the coordinates */ for (i = 0; i < n_rec[0]; i++) { GMT_geo_to_xy (GMT, data[0][s->x_col][i], data[0][s->y_col][i], &xx, &yy); data[0][s->x_col][i] = xx; data[0][s->y_col][i] = yy; } } if ((dist[0] = GMT_dist_array_2 (GMT, data[0][s->x_col], data[0][s->y_col], n_rec[0], dist_scale, s->dist_flag)) == NULL) GMT_err_fail (GMT, GMT_MAP_BAD_DIST_FLAG, ""); time[0] = (has_time[0]) ? data[0][s->t_col] : x2sys_dummytimes (GMT, n_rec[0]) ; GMT_init_track (GMT, data[0][s->y_col], n_rec[0], &ylist_A); for (B = A; B < n_tracks; B++) { if (duplicate[B]) continue; same = !strcmp (trk_name[A], trk_name[B]); if (same && !(A == B)) { GMT_Report (API, GMT_MSG_NORMAL, "File %s repeated on command line - skipped\n", trk_name[A]); continue; } if (!internal && same) continue; /* Only do external errors */ if (!external && !same) continue; /* Only do internal errors */ if (Ctrl->A.active && !combo_ok (trk_name[A], trk_name[B], pair, n_pairs)) continue; /* Do not want this combo */ if (Ctrl->C.active) tic = clock(); /* To report execution time from this pair */ GMT_Report (API, GMT_MSG_VERBOSE, "Processing %s - %s : ", trk_name[A], trk_name[B]); if (same) { /* Just set pointers */ data[1] = data[0]; dist[1] = dist[0]; time[1] = time[0]; has_time[1] = has_time[0]; n_rec[1] = n_rec[0]; ylist_B = ylist_A; data_set[1] = data_set[0]; } else { /* Must read a second file */ x2sys_err_fail (GMT, (s->read_file) (GMT, trk_name[B], &data[1], s, &data_set[1], &GMT->current.io, &n_rec[1]), trk_name[B]); if (n_rec[1] == 0) { /* No data in track B */ x2sys_free_data (GMT, data[1], s->n_out_columns, &data_set[1]); continue; } has_time[1] = false; if (got_time) { /* Check to make sure we do in fact have time */ for (i = n_bad = 0; i < n_rec[1]; i++) n_bad += GMT_is_dnan (data[1][s->t_col][i]); if (n_bad < n_rec[1]) has_time[1] = true; } if (do_project) { /* Convert all the coordinates */ for (i = 0; i < n_rec[0]; i++) { GMT_geo_to_xy (GMT, data[1][s->x_col][i], data[1][s->y_col][i], &xx, &yy); data[1][s->x_col][i] = xx; data[1][s->y_col][i] = yy; } } if ((dist[1] = GMT_dist_array_2 (GMT, data[1][s->x_col], data[1][s->y_col], n_rec[1], dist_scale, s->dist_flag)) == NULL) GMT_err_fail (GMT, GMT_MAP_BAD_DIST_FLAG, ""); time[1] = (has_time[1]) ? data[1][s->t_col] : x2sys_dummytimes (GMT, n_rec[1]); GMT_init_track (GMT, data[1][s->y_col], n_rec[1], &ylist_B); } /* Calculate all possible crossover locations */ nx = GMT_crossover (GMT, data[0][s->x_col], data[0][s->y_col], data_set[0].ms_rec, ylist_A, n_rec[0], data[1][s->x_col], data[1][s->y_col], data_set[1].ms_rec, ylist_B, n_rec[1], (A == B), wrap, &XC); if (nx && xover_locations_only) { /* Report crossover locations only */ sprintf (line, "%s - %s", trk_name[A], trk_name[B]); GMT_Put_Record (API, GMT_WRITE_SEGMENT_HEADER, line); for (i = 0; i < nx; i++) { out[0] = XC.x[i]; out[1] = XC.y[i]; if (s->geographic) GMT_lon_range_adjust (s->geodetic, &out[0]); GMT_Put_Record (API, GMT_WRITE_DOUBLE, out); /* Write this to output */ } GMT_x_free (GMT, &XC); } else if (nx) { /* Got crossovers, now estimate crossover values */ int64_t start_s; first_crossover = true; for (i = 0; i < nx; i++) { /* For each potential crossover */ GMT_memset (ok, n_data_col, unsigned int); n_ok = 0; for (k = 0; k < 2; k++) { /* For each of the two data sets involved */ /* Get node number to each side of crossover location */ /* --o----------o--------o------X-------o-------o----------o-- ----> time ^ ^ ^ left xover right */ left[k] = lrint (floor (XC.xnode[k][i])); right[k] = lrint (ceil (XC.xnode[k][i])); if (left[k] == right[k]) { /* Crosses exactly on a node; move left or right so interpolation will work */ if (left[k] > 0) left[k]--; /* Move back so cross occurs at right[k] */ else right[k]++; /* Move forward so cross occurs at left[k] */ } deld = dist[k][right[k]] - dist[k][left[k]]; delt = time[k][right[k]] - time[k][left[k]]; /* Check if speed is outside accepted domain */ speed[k] = (delt == 0.0) ? GMT->session.d_NaN : vel_scale * (deld / (delt * t_scale)); if (Ctrl->S.active[VLO] && !GMT_is_dnan (speed[k]) && (speed[k] < Ctrl->S.limit[VLO] || speed[k] > Ctrl->S.limit[VHI])) continue; /* Linearly estimate the crossover times and distances */ dt = XC.xnode[k][i] - left[k]; time_x[k] = time[k][left[k]]; dist_x[k] = dist[k][left[k]]; if (dt > 0.0) { time_x[k] += dt * delt; dist_x[k] += dt * deld; } for (j = 0; j < n_data_col; j++) { /* Evaluate each field at the crossover */ col = col_number[j]; start = t_right = left[k]; end = t_left = right[k]; n_left = n_right = 0; xdata[k][col] = GMT->session.d_NaN; /* In case of nuthin' */ /* First find the required points to the left of the xover */ start_s = start; while (start_s >= 0 && n_left < Ctrl->W.width) { if (!GMT_is_dnan (data[k][col][start])) { n_left++; if (t_left > left[k]) t_left = start; y[Ctrl->W.width-n_left] = data[k][col][start]; t[Ctrl->W.width-n_left] = time[k][start]; } start--; start_s--; } if (!n_left) continue; if (got_time && ((time_x[k] - time[k][t_left]) > Bix.time_gap)) continue; if ((dist_x[k] - dist[k][t_left]) > Bix.dist_gap) continue; /* Ok, that worked. Now for the right side: */ while (end < n_rec[k] && n_right < Ctrl->W.width) { if (!GMT_is_dnan (data[k][col][end])) { y[Ctrl->W.width+n_right] = data[k][col][end]; t[Ctrl->W.width+n_right] = time[k][end]; n_right++; if (t_right < right[k]) t_right = end; } end++; } if (!n_right) continue; /* See if we pass any gap criteria */ if (got_time && ((time[k][t_right] - time_x[k]) > Bix.time_gap)) continue; /* Exceeded time gap */ if ((dist[k][t_right] - dist_x[k]) > Bix.dist_gap) continue; /* Exceeded distance gap */ /* Ok, got enough data to interpolate at xover */ first = Ctrl->W.width - n_left; n_errors = GMT_intpol (GMT, &t[first], &y[first], (n_left + n_right), 1, &time_x[k], &xdata[k][col], GMT->current.setting.interpolant); if (n_errors == 0) { /* OK */ ok[j]++; n_ok++; } } } /* Only output crossover if there are any data there */ if (n_ok == 0) continue; for (j = n_ok = 0; j < n_data_col; j++) if (ok[j] == 2) n_ok++; if (n_ok == 0) continue; /* OK, got something to report */ /* Load the out array */ out[0] = XC.x[i]; /* Crossover location */ out[1] = XC.y[i]; for (k = 0; k < 2; k++) { /* Get times, distances, headings, and velocities */ /* Get time */ out[2+k] = (got_time && !has_time[k]) ? X2SYS_NaN : time_x[k]; /* Get cumulative distance at crossover */ out[k+4] = dist_x[k]; /* Estimate heading there */ j = k + 6; out[j] = (!GMT_is_dnan (speed[k]) && (!Ctrl->S.active[HHI] || speed[k] > Ctrl->S.limit[HHI])) ? (*GMT->current.map.azimuth_func) (GMT, data[k][s->x_col][right[k]], data[k][s->y_col][right[k]], data[k][s->x_col][left[k]], data[k][s->y_col][left[k]], false) : X2SYS_NaN; /* Estimate velocities there */ j = k + 8; out[j] = (has_time[k]) ? speed[k] : X2SYS_NaN; } /* Calculate crossover and mean value */ for (k = 0, j = 10; k < n_data_col; k++) { if (Ctrl->Z.active) { col = col_number[k]; out[j++] = xdata[0][col]; out[j++] = xdata[1][col]; } else { if (ok[k] == 2) { col = col_number[k]; out[j++] = xdata[0][col] - xdata[1][col]; out[j++] = 0.5 * (xdata[0][col] + xdata[1][col]); } else { out[j] = out[j+1] = X2SYS_NaN; j += 2; } } } if (first_header) { /* Write the header record */ char *cmd = NULL, *c = GMT->current.setting.io_col_separator; t_or_i = (got_time) ? 't' : 'i'; sprintf (line, "# Tag: %s", Ctrl->T.TAG); GMT_Put_Record (API, GMT_WRITE_TABLE_HEADER, line); cmd = GMT_Create_Cmd (API, options); sprintf (line, "# Command: %s %s", THIS_MODULE_NAME, cmd); /* Build command line argument string */ GMT_free (GMT, cmd); GMT_Put_Record (API, GMT_WRITE_TABLE_HEADER, line); sprintf (line, "# %s%s%s%s%c_1%s%c_2%sdist_1%sdist_2%shead_1%shead_2%svel_1%svel_2", s->info[s->out_order[s->x_col]].name, c, s->info[s->out_order[s->y_col]].name, c, t_or_i, c, t_or_i, c, c, c, c, c, c); for (j = 0; j < n_data_col; j++) { col = col_number[j]; if (Ctrl->Z.active) sprintf (item, "%s%s_1%s%s_2", c, s->info[s->out_order[col]].name, c, s->info[s->out_order[col]].name); else sprintf (item, "%s%s_X%s%s_M", c, s->info[s->out_order[col]].name, c, s->info[s->out_order[col]].name); strcat (line, item); } GMT_Put_Record (API, GMT_WRITE_TABLE_HEADER, line); first_header = false; } if (first_crossover) { char info[GMT_BUFSIZ] = {""}, start[2][GMT_LEN64], stop[2][GMT_LEN64]; for (k = 0; k < 2; k++) { if (has_time[k]) { /* Find first and last record times */ for (j = 0; j < n_rec[k] && GMT_is_dnan (time[k][j]); j++); /* Find first non-NaN time */ GMT_ascii_format_col (GMT, start[k], time[k][j], GMT_OUT, 2); for (j = n_rec[k]-1; j > 0 && GMT_is_dnan (time[k][j]); j--); /* Find last non-NaN time */ GMT_ascii_format_col (GMT, stop[k], time[k][j], GMT_OUT, 3); } else { strcpy (start[k], "NaN"); strcpy (stop[k], "NaN"); } } sprintf (info, "%s/%s/%g %s/%s/%g", start[0], stop[0], dist[0][n_rec[0]-1], start[1], stop[1], dist[1][n_rec[1]-1]); sprintf (line, x2sys_header, trk_name[A], data_set[0].year, trk_name[B], data_set[1].year, info); GMT_Put_Record (API, GMT_WRITE_SEGMENT_HEADER, line); first_crossover = false; } if (s->geographic) GMT_lon_range_adjust (s->geodetic, &out[0]); GMT_Put_Record (API, GMT_WRITE_DOUBLE, out); /* Write this to output */ } GMT_x_free (GMT, &XC); } if (!same) { /* Must free up memory for B */ x2sys_free_data (GMT, data[1], s->n_out_columns, &data_set[1]); GMT_free (GMT, dist[1]); if (!got_time) GMT_free (GMT, time[1]); GMT_free (GMT, ylist_B); } if (!Ctrl->C.active) GMT_Report (API, GMT_MSG_VERBOSE, "%" PRIu64 "\n", nx); else { toc = clock(); GMT_Report (API, GMT_MSG_VERBOSE, "%" PRIu64 "\t%.3f sec\n", nx, (double)(toc - tic)/1000); if (fpC) /* Save also the run time in file */ fprintf (fpC, "%s\t%s\t%d\t%.3f\n", trk_name[A], trk_name[B], (int)nx, (double)(toc - tic)/1000); } } /* Must free up memory for A */ x2sys_free_data (GMT, data[0], s->n_out_columns, &data_set[0]); GMT_free (GMT, dist[0]); if (!got_time) GMT_free (GMT, time[0]); GMT_free (GMT, ylist_A); } if (fpC) fclose (fpC); if (GMT_End_IO (API, GMT_OUT, 0) != GMT_OK) { /* Disables further data output */ Return (API->error); } /* Free up other arrays */ GMT_free (GMT, xdata[0]); GMT_free (GMT, xdata[1]); GMT_free (GMT, out); GMT_free (GMT, duplicate); if (n_data_col) { GMT_free (GMT, t); GMT_free (GMT, y); GMT_free (GMT, col_number); GMT_free (GMT, ok); } x2sys_free_list (GMT, trk_name, n_tracks); x2sys_end (GMT, s); Return (GMT_OK); }