/*-------------------------------------------------------------------- * $Id: splitxyz.c,v 1.2.4.3 2002/02/27 17:58:55 pwessel Exp $ * * Copyright (c) 1991-2002 by P. Wessel and W. H. F. Smith * See COPYING 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 General Public License as published by * the Free Software Foundation; version 2 of the License. * * 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 General Public License for more details. * * Contact info: gmt.soest.hawaii.edu *--------------------------------------------------------------------*/ /* * * read a file of lon, lat, zvalue[, distance, azimuth] * and split it into profile segments. * * Author: W. H. F. Smith * Date: 24 July, 1991-2001 * Version: 3.4.1 */ #include "gmt.h" #define F_RES 1000 /* Number of points in filter halfwidth */ struct DATA { double d; double a; double x; double y; double z; double w; } *data; main(int argc, char **argv) { int i, ndata, n_alloc, nprofiles, begin, end, n_expected_fields, n_out, n_fields; int n_required_fields; BOOLEAN error = FALSE, dh_supplied = FALSE, ok, hilow, map_units = FALSE; BOOLEAN xy_only = FALSE; double desired_azim = 90.0, tolerance_azim = 360.0; double min_dist = 100.0, xy_filter = 0.0, z_filter = 0.0; double d_gap = 10.0, course_change = 0.0, d_to_km; double dy, dx, last_d, last_c, last_s, csum, ssum, this_c, this_s, dotprod; double mean_azim, fwork[F_RES], *in, out[3]; char buffer[BUFSIZ], basename[120], filename[120]; FILE *fp = NULL, *fpout; void filterz(struct DATA *data, int ndata, double z_filter, double *fwork, int hilow), filterxy_setup(double *fwork), filterxy(struct DATA *data, int begin, int end, double xy_filter, double *fwork, int hilow); argc = GMT_begin (argc, argv); basename[0] = '\0'; /* Later, if (basename[0]) we have to write files. */ fpout = GMT_stdout; for (i = 1; i < argc; i++) { if (argv[i][0] == '-') { switch (argv[i][1]) { /* Common parameters */ case 'H': case 'V': case ':': case '\0': error += GMT_get_common_args (argv[i], 0, 0, 0, 0); break; /* Supplemental parameters */ case 'b': error += GMT_io_selection (&argv[i][2]); break; case 'A': if ( (sscanf(&argv[i][2], "%lf/%lf", &desired_azim, &tolerance_azim)) != 2) { fprintf (stderr, "%s: GMT SYNTAX ERROR -A option: Can't dechiper values\n", GMT_program); error = TRUE; } break; case 'C': if ( (sscanf(&argv[i][2], "%lf", &course_change)) != 1) { fprintf (stderr, "%s: GMT SYNTAX ERROR -C option: Can't dechiper value\n", GMT_program); error = TRUE; } break; case 'D': if ( (sscanf(&argv[i][2], "%lf", &min_dist)) != 1) { fprintf (stderr, "%s: GMT SYNTAX ERROR -D option: Can't dechiper value\n", GMT_program); error = TRUE; } break; case 'F': if ( (sscanf(&argv[i][2], "%lf/%lf", &xy_filter, &z_filter)) != 2) { fprintf (stderr, "%s: GMT SYNTAX ERROR -F option: Can't dechiper values\n", GMT_program); error = TRUE; } break; case 'G': if ( (sscanf(&argv[i][2], "%lf", &d_gap)) != 1) { fprintf (stderr, "%s: GMT SYNTAX ERROR -G option: Can't dechiper value\n", GMT_program); error = TRUE; } break; case 'M': map_units = GMT_geographic_in = GMT_geographic_out = TRUE; break; case 'N': strcpy(basename, &argv[i][2]); break; case 'S': dh_supplied = TRUE; break; case 'Z': xy_only = TRUE; break; default: error = TRUE; GMT_default_error (argv[i][1]); break; } } else { if ( (fp = GMT_fopen(argv[i], GMT_io.r_mode)) == NULL) { fprintf(stderr,"%s: Cannot open %s\n", GMT_program, argv[i]); error = TRUE; } } } if (argc == 1 || GMT_quick) { /* Display usage */ fprintf (stderr, "splitxyz %s - Split xyz[dh] files into segments\n\n", GMT_VERSION); fprintf(stderr,"usage: splitxyz [] -C [-A/]\n"); fprintf(stderr,"\t[-D] [-F/] [-G] [-H[]] [-M]\n"); fprintf(stderr,"\t[-N] [-S] [-V] [-Z] [-:] [-bi[s][]] [-bo[s][]]\n\n"); if (GMT_quick) exit (EXIT_FAILURE); fprintf(stderr,"\tGive xyz[dh]file name or read stdin.\n"); fprintf(stderr,"\t-C Profile ends when change of heading exceeds .\n"); fprintf (stderr, "\n\tOPTIONS:\n"); fprintf(stderr,"\t-A Only write profile if mean direction is w/in +/- \n"); fprintf(stderr,"\t of . [Default = All].\n"); fprintf(stderr,"\t-D Only write profile if length is at least .\n"); fprintf(stderr,"\t [Default = 100 dist units].\n"); fprintf(stderr,"\t-F Filter the data. Give full widths of cosine arch filters for xy and z.\n"); fprintf(stderr,"\t Defaults are both widths = 0, giving no filtering.\n"); fprintf(stderr,"\t Use negative width to highpass.\n"); fprintf(stderr,"\t-G Do not let profiles have gaps exceeding . [Default = 10 dist units].\n"); GMT_explain_option ('H'); fprintf(stderr,"\t-M Map units TRUE; x,y in degrees, dist units in km. [Default dist unit = x,y unit].\n"); fprintf(stderr,"\t-N Write output to separate files named .profile#.\n"); fprintf(stderr,"\t [Default all to stdout, separated by >].\n"); fprintf(stderr,"\t-S dh is supplied. Input is 5 col x,y,z,d,h with d non-decreasing.\n"); fprintf(stderr,"\t [Default input is 3 col x,y,z only and computes d,h from the data].\n"); fprintf(stderr,"\t-Z No z-values. Input is 2 col x,y only\n"); GMT_explain_option ('V'); GMT_explain_option (':'); GMT_explain_option ('i'); GMT_explain_option ('n'); fprintf(stderr,"\t Default input columns is set given -S and -Z options.\n"); GMT_explain_option ('o'); GMT_explain_option ('.'); exit (EXIT_FAILURE); } if (course_change <= 0.0) { fprintf (stderr, "%s: GMT SYNTAX ERROR -C option: Course change tolerance must be positive\n", GMT_program); error++; } if (tolerance_azim < 0.0) { fprintf (stderr, "%s: GMT SYNTAX ERROR -A option: Azimuth tolerance must be positive\n", GMT_program); error++; } if (d_gap < 0.0) { fprintf (stderr, "%s: GMT SYNTAX ERROR -G option: Data gap distance must be positive\n", GMT_program); error++; } if (desired_azim < 0.0 || desired_azim > 360) { fprintf (stderr, "%s: GMT SYNTAX ERROR -A option: azimuth must be in 0-360 degree range\n", GMT_program); error++; } if (xy_only && dh_supplied) { fprintf (stderr, "%s: GMT SYNTAX ERROR -Z option: Cannot be used with -S option\n", GMT_program); error++; } if (xy_only && z_filter != 0.0) { fprintf (stderr, "%s: GMT SYNTAX ERROR -F option: Cannot specify z-filter while using -Z option\n", GMT_program); error++; } if (GMT_io.binary[0] && gmtdefs.io_header) { fprintf (stderr, "%s: GMT SYNTAX ERROR. Binary input data cannot have header -H\n", GMT_program); error++; } if (GMT_io.binary[1] && !basename[0]) { fprintf (stderr, "%s: GMT SYNTAX ERROR. Binary output requires a namestem in -N\n", GMT_program); error++; } n_required_fields = (dh_supplied) ? 5 : ((xy_only) ? 2 : 3); if (GMT_io.binary[0] && GMT_io.ncol[0] == 0) GMT_io.ncol[0] = n_required_fields; if (GMT_io.binary[0]) n_expected_fields = GMT_io.ncol[0]; else n_expected_fields = n_required_fields; if (GMT_io.binary[0] && n_expected_fields < n_required_fields) { fprintf (stderr, "%s: GMT SYNTAX ERROR. Binary input data must have at least %d columns\n", GMT_program, n_required_fields); error++; } if (error) exit (EXIT_FAILURE); GMT_put_history (argc, argv); /* Update .gmtcommands */ if (GMT_io.binary[0] && gmtdefs.verbose) { char *type[2] = {"double", "single"}; fprintf (stderr, "%s: Expects %d-column %s-precision binary data\n", GMT_program, GMT_io.ncol[0], type[GMT_io.single_precision[0]]); } #ifdef SET_IO_MODE GMT_setmode (1); #endif tolerance_azim *= D2R; if (desired_azim > 180.0) desired_azim -= 180.0; /* Put in Easterly strikes */ desired_azim = D2R * (90.0 - desired_azim); /* Work in cartesian angle and radians */ course_change *= D2R; d_to_km = 0.001 * 2.0 * M_PI * gmtdefs.ellipse[N_ELLIPSOIDS-1].eq_radius / 360.0; if (fp == NULL) { fp = GMT_stdin; #ifdef SET_IO_MODE GMT_setmode (0); #endif } n_alloc = GMT_CHUNK; ndata = 0; i = 0; data = (struct DATA *) GMT_memory (VNULL, (size_t)n_alloc, sizeof(struct DATA), GMT_program); if (!GMT_io.binary[0] && gmtdefs.io_header) for (i = 0; i < gmtdefs.n_header_recs; i++) GMT_fgets (buffer, BUFSIZ, fp); while ((n_fields = GMT_input (fp, &n_expected_fields, &in)) >= 0 && !(GMT_io.status & GMT_IO_EOF)) { /* Not yet EOF */ if (GMT_io.status & GMT_IO_MISMATCH) { fprintf (stderr, "%s: Mismatch between actual (%d) and expected (%d) fields near line %d\n", GMT_program, n_fields, n_expected_fields, i); exit (EXIT_FAILURE); } i++; if (ndata == n_alloc) { n_alloc += GMT_CHUNK; data = (struct DATA *) GMT_memory ((void *)data, (size_t)n_alloc, sizeof(struct DATA), GMT_program); } if (!dh_supplied) { data[ndata].x = in[0]; data[ndata].y = in[1]; if (!xy_only) data[ndata].z = in[2]; if (ndata) { dy = (data[ndata].y - data[ndata-1].y); dx = (data[ndata].x - data[ndata-1].x); if (map_units) { dy *= d_to_km; dx *= (d_to_km * cos (0.5 * (data[ndata].y + data[ndata-1].y) * D2R) ); } if (dy == 0.0 && dx == 0.0) { data[ndata].d = data[ndata-1].d; data[ndata].a = data[ndata-1].a; } else { data[ndata].d = data[ndata-1].d + hypot(dx,dy); data[ndata].a = d_atan2(dy,dx); } } else { data[ndata].d = data[ndata].a = 0.0; } ndata++; } else { data[ndata].x = in[0]; data[ndata].y = in[1]; data[ndata].z = in[2]; data[ndata].d = in[3]; data[ndata].a = D2R * (90.0 - in[4]); ndata++; } } data = (struct DATA *) GMT_memory ((void *)data, (size_t)ndata, sizeof(struct DATA), GMT_program); if (!dh_supplied) data[0].a = data[1].a; if (fp != GMT_stdin) GMT_fclose(fp); /* Now we have read the data and can filter z if necessary. */ if (z_filter < 0.0) { hilow = TRUE; z_filter = fabs(z_filter); filterz(data, ndata, z_filter, fwork, hilow); } else if (z_filter > 0.0) { hilow = FALSE; filterz(data, ndata, z_filter, fwork, hilow); } if (xy_filter < 0.0) { hilow = TRUE; xy_filter = fabs(xy_filter); filterxy_setup(fwork); } else if (xy_filter > 0.0) { hilow = FALSE; filterxy_setup(fwork); } /* Now we are ready to search for segments. */ n_out = (xy_only) ? 2 : 3; nprofiles = 0; begin = 0; end = begin; while (end < ndata-1) { last_d = data[begin].d; last_c = cos(data[begin].a); last_s = sin(data[begin].a); csum = last_c; ssum = last_s; ok = TRUE; while (ok && end < ndata-1) { end++; if (data[end].d - last_d > d_gap) { /* Fails due to too much distance gap */ ok = FALSE; continue; } this_c = cos(data[end].a); this_s = sin(data[end].a); dotprod = this_c * last_c + this_s * last_s; if (fabs(dotprod) > 1.0) dotprod = copysign(1.0, dotprod); if (d_acos(dotprod) > course_change) { /* Fails due to too much change in azimuth */ ok = FALSE; continue; } /* Get here when this point belongs with last one: */ csum += this_c; ssum += this_s; last_c = this_c; last_s = this_s; last_d = data[end].d; } /* Get here when we have found a beginning and end */ if (ok) { /* Last point in input should be included in this segment */ end++; } if (end - begin - 1) { /* There are at least two points in the list. */ if ((data[end-1].d - data[begin].d) >= min_dist) { /* List is long enough. Check strike. Compute mean_azim in range [-pi/2, pi/2]: */ mean_azim = d_asin(ssum/sqrt(csum*csum + ssum*ssum)); mean_azim = fabs(mean_azim - desired_azim); if (mean_azim <= tolerance_azim || mean_azim >= (M_PI-tolerance_azim) ) { /* List has acceptable strike. */ if (xy_filter != 0.0) filterxy(data, begin, end, xy_filter, fwork, hilow); nprofiles++; /* If (basename) we write separate files, else GMT_stdout with > marks: */ if (basename[0]) { sprintf(filename, "%s.profile%d", basename, nprofiles); if ( (fpout = GMT_fopen(filename, GMT_io.w_mode)) == NULL) { fprintf(stderr,"%s: Cannot create %s\n", GMT_program, argv[i]); exit (EXIT_FAILURE); } } else fprintf (GMT_stdout, "> Start profile # %d\n", nprofiles); for (i = begin; i < end; i++) { out[0] = data[i].x; out[1] = data[i].y; if (!xy_only) out[2] = data[i].z; GMT_output (fpout, n_out, out); } if (fpout != GMT_stdout) { GMT_fclose (fpout); if (gmtdefs.verbose) fprintf (stderr,"%s: Wrote %d points to file %s\n", GMT_program, end-begin, filename); } } } } begin = end; } /* Get here when all profiles have been found and written. */ if (gmtdefs.verbose) fprintf(stderr,"%s: Split %d data into %d files.\n", GMT_program, ndata, nprofiles); GMT_free ((void *)data); GMT_end (argc, argv); } void filterz(struct DATA *data, int ndata, double z_filter, double *fwork, int hilow) { double half_width, sum, dt, tmp; int i, j, k, istart, istop; half_width = 0.5 * z_filter; sum = 0.0; dt = F_RES/half_width; tmp = M_PI / F_RES; for (i = 0; i < F_RES; i++) { fwork[i] = 1.0 + cos(i*tmp); sum += fwork[i]; } for (i = 1; i < F_RES; i++) { fwork[i] /= sum; } j = 0; istart = 0; istop = 0; while(j < ndata) { while(istart < ndata && data[istart].d - data[j].d <= -half_width) istart++; while(istop < ndata && data[istop].d - data[j].d < half_width) istop++; sum = 0.0; data[j].w = 0.0; for(i = istart; i < istop; i++) { k = (int)floor(dt*fabs(data[i].d - data[j].d)); sum += fwork[k]; data[j].w += (data[i].z * fwork[k]); } data[j].w /= sum; j++; } if (hilow) { for (i = 0; i < ndata; i++) data[i].z = data[i].z - data[i].w; } else { for (i = 0; i < ndata; i++) data[i].z = data[i].w; } return; } void filterxy_setup(double *fwork) { int i; double tmp, sum = 0.0; tmp = M_PI / F_RES; for (i = 0; i < F_RES; i++) { fwork[i] = 1.0 + cos(i*tmp); sum += fwork[i]; } for (i = 1; i < F_RES; i++) { fwork[i] /= sum; } return; } void filterxy(struct DATA *data, int begin, int end, double xy_filter, double *fwork, int hilow) { int i, j, k, istart, istop; double half_width, dt, sum; half_width = 0.5 * fabs(xy_filter); dt = F_RES/half_width; j = begin; istart = begin; istop = begin; while(j < end) { while(istart < end && data[istart].d - data[j].d <= -half_width) istart++; while(istop < end && data[istop].d - data[j].d < half_width) istop++; sum = 0.0; data[j].w = 0.0; for(i = istart; i < istop; i++) { k = (int)floor(dt*fabs(data[i].d - data[j].d)); sum += fwork[k]; data[j].w += (data[i].x * fwork[k]); } data[j].w /= sum; j++; } if (hilow) { for (i = begin; i < end; i++) data[i].x = data[i].x - data[i].w; } else { for (i = begin; i < end; i++) data[i].x = data[i].w; } j = begin; istart = begin; istop = begin; while(j < end) { while(istart < end && data[istart].d - data[j].d <= -half_width) istart++; while(istop < end && data[istop].d - data[j].d < half_width) istop++; sum = 0.0; data[j].w = 0.0; for(i = istart; i < istop; i++) { k = (int)floor(dt*fabs(data[i].d - data[j].d)); sum += fwork[k]; data[j].w += (data[i].y * fwork[k]); } data[j].w /= sum; j++; } if (hilow) { for (i = begin; i < end; i++) data[i].y = data[i].y - data[i].w; } else { for (i = begin; i < end; i++) data[i].y = data[i].w; } return; }