/*-------------------------------------------------------------------- * $Id: pswiggle.c,v 1.4.4.2 2002/02/27 17:41:10 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 *--------------------------------------------------------------------*/ /* * pswiggle reads x,y,z from GMT_stdin and plots a wiggleplot using the * specified map projection. If the distance between 2 consecutive points * exceeds the value of dist_gap, a data gap is assumed. The user may * select a preferred direction where s/he wants the positive anomalies * to be pointing towards. Positive normal vectors to the trackline * will then always be withing +- 90 degrees of this direction. * Separate colors may be specified for positive anomaly, outline of * anomaly, and trackline. Plotting of the outline and track is optional. * * Author: Paul Wessel * Date: 20-JUN-2000 * Version: 3.4.1 * */ #include "gmt.h" #define MAX_POINTS 400 double *xx, *yy, *zz, *lon, *lat, *z; double start_az, stop_az; void GMT_draw_z_scale (double x0, double y0, double length, double zscale, BOOLEAN gave_xy, char *units); void plot_wiggle (double *x, double *y, double *z, int np, double zscale, BOOLEAN fixed, double fix_az, struct GMT_FILL *fill, struct GMT_PEN *pen_o, struct GMT_PEN *pen_t, BOOLEAN paint_wiggle, BOOLEAN negative, BOOLEAN outline, BOOLEAN track); main (int argc, char **argv) { int i, j, n, k, n_files = 0, fno, n_alloc = GMT_CHUNK, n_args, n_expected_fields, n_fields; BOOLEAN error = FALSE, nofile = TRUE, outline = FALSE, track = FALSE, paint_wiggle, draw_scale = FALSE; BOOLEAN do_fill = FALSE, negative = FALSE, fix = FALSE, done, gave_xy = FALSE; double west = 0.0, east = 360.0, south = 0.0, north = 0.0, x_2, y_2, bias = 0.0, fix_az; double zscale = 0.0, dx, dy, dz, ds, x0, y0, length, dist_gap = DBL_MAX, *in, pref_az = 0.0; char line[BUFSIZ], *units, txt_a[32], txt_b[32], unit = 0; FILE *fp = NULL; struct GMT_FILL fill; struct GMT_PEN pen_o, pen_t; argc = GMT_begin (argc, argv); GMT_init_pen (&pen_o, GMT_PENWIDTH); GMT_init_pen (&pen_t, GMT_PENWIDTH); GMT_init_fill (&fill, gmtdefs.basemap_frame_rgb[0], gmtdefs.basemap_frame_rgb[1], gmtdefs.basemap_frame_rgb[2]); for (i = 1; i < argc; i++) { if (argv[i][0] == '-') { switch(argv[i][1]) { /* Common parameters */ case 'B': case 'H': case 'J': case 'K': case 'O': case 'P': case 'R': case 'U': case 'V': case 'X': case 'x': case 'Y': case 'y': case 'c': case ':': case '\0': error += GMT_get_common_args (argv[i], &west, &east, &south, &north); break; /* Supplemental parameters */ case 'b': error += GMT_io_selection (&argv[i][2]); break; case 'A': pref_az = atof (&argv[i][2]); break; case 'D': dist_gap = atof (&argv[i][2]); break; case 'E': sscanf (&argv[i][2], "%lf/%lf", &z_project.view_azimuth, &z_project.view_elevation); break; case 'G': if (GMT_getfill (&argv[i][2], &fill)) { GMT_fill_syntax ('G'); error++; } do_fill = TRUE; break; case 'I': fix_az = atof (&argv[i][2]); fix = TRUE; break; case 'M': GMT_multisegment (&argv[i][2]); break; case 'N': negative = TRUE; break; case 'C': bias = atof (&argv[i][2]); break; case 'T': if (GMT_getpen (&argv[i][2], &pen_t)) { GMT_pen_syntax ('T'); error++; } track = TRUE; break; case 'W': if (GMT_getpen (&argv[i][2], &pen_o)) { GMT_pen_syntax ('W'); error++; } outline = TRUE; break; case 'S': j = 0; if (argv[i][2] == 'x') gave_xy = TRUE, j = 1; k = sscanf (&argv[i][2+j], "%[^/]/%[^/]/%lf", txt_a, txt_b, &length); x0 = GMT_ddmmss_to_degree (txt_a); y0 = GMT_ddmmss_to_degree (txt_b); units = strrchr (argv[i], '/'); units++; if (k != 3) { fprintf (stderr, "%s: GMT SYNTAX ERROR -S option: Correct syntax\n", GMT_program); fprintf (stderr, " -S[x]//[/]\n"); error++; } draw_scale = TRUE; break; case 'Z': j = strlen (argv[i]) - 1; if (strchr ("cimpCIMP", (int)argv[i][j])) unit = argv[i][j]; zscale = atof (&argv[i][2]); break; /* Illegal options */ default: error = TRUE; GMT_default_error (argv[i][1]); break; } } else n_files++; } if (argc == 1 || GMT_quick) { fprintf (stderr,"pswiggle %s - Plot xyz-series along tracks\n\n", GMT_VERSION); fprintf (stderr, "usage: pswiggle -J -R -Z\n"); fprintf (stderr, "\t[-A] [-B] [-C
] [-D] [-E] [-G] [-H[]]\n"); fprintf (stderr, "\t[-I] [-K] [-M[]] [-N] [-O] [-P] [-S[x]///] [-T]\n"); fprintf (stderr, "\t[-U] [-V] [-W] [-X] [-Y] [-c] [-:] [-bi[s][]\n"); if (GMT_quick) exit (EXIT_FAILURE); fprintf (stderr, "\t is one or more files. If none, read standard input\n"); GMT_explain_option ('j'); GMT_explain_option ('R'); fprintf (stderr, "\n\tOPTIONS:\n"); fprintf (stderr, "\t-A set azimuth for preferred positive wiggle orientation [0.0 (north)]\n"); GMT_explain_option ('b'); fprintf (stderr, "\t-C sets center value to be removed from z before plotting. [0]\n"); fprintf (stderr, "\t-D means there is a datagap if 2 points are more than distance units apart\n"); fprintf (stderr, "\t gap is in km for map projections, user-units otherwise\n"); fprintf (stderr, "\t-E set azimuth and elevation of viewpoint for 3-D perspective [180/90]\n"); fprintf (stderr, "\t-G Specify color or pattern for positive areas. fill can be either\n"); fprintf (stderr, "\t 1) (each 0-255) for color or (0-255) for gray-shade [0]\n"); fprintf (stderr, "\t 2) p[or P]/ for predefined patterns (1-90).\n"); GMT_explain_option ('H'); fprintf (stderr, "\t-I set fixed projection azimuths for wiggles\n"); GMT_explain_option ('K'); GMT_explain_option ('M'); fprintf (stderr, "\t-N Fill negative wiggles instead [Default is positive]\n"); GMT_explain_option ('O'); GMT_explain_option ('P'); fprintf (stderr, "\t-S draws a simple vertical scale centered on /. Use -Sx to specify cartesian coordinates instead.\n"); fprintf (stderr, "\t is in z-units, append unit name for labeling\n"); fprintf (stderr, "\t-T specifies track pen attributes. [Default is no track]\n"); GMT_explain_option ('U'); GMT_explain_option ('V'); fprintf (stderr, "\t-W specifies outline pen attributes. [Default is no outline]\n"); GMT_explain_option ('X'); fprintf (stderr, "\t-Z gives the wiggle scale in data-units per %s\n", GMT_unit_names[gmtdefs.measure_unit]); GMT_explain_option ('c'); GMT_explain_option (':'); GMT_explain_option ('i'); GMT_explain_option ('n'); fprintf (stderr, "\t Default is 3 input columns\n"); GMT_explain_option ('.'); exit (EXIT_FAILURE); } if (!project_info.region_supplied) { fprintf (stderr, "%s: GMT SYNTAX ERROR: Must specify -R option\n", GMT_program); error++; } if (!(outline || do_fill)) { fprintf (stderr, "%s: GMT SYNTAX ERROR: Must specify at least one of -G, -W\n", GMT_program); error++; } if (z_project.view_azimuth > 360.0 || z_project.view_elevation <= 0.0 || z_project.view_elevation > 90.0) { fprintf (stderr, "%s: GMT SYNTAX ERROR -E option: Enter azimuth in 0-360 range, elevation in 0-90 range\n", GMT_program); error++; } if (zscale == 0.0) { fprintf (stderr, "%s: GMT SYNTAX ERROR -Z option: scale must be nonzero\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[0] && GMT_io.ncol[0] == 0) GMT_io.ncol[0] = 3; if (GMT_io.binary[0] && GMT_io.ncol[0] < 3) { fprintf (stderr, "%s: GMT SYNTAX ERROR. Binary input data (-bi) must have at least 3 columns\n", GMT_program); 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]]); } GMT_map_setup (west, east, south, north); if (MAPPING) dist_gap *= 1000.0; /* Distance in meters */ ps_plotinit (CNULL, gmtdefs.overlay, gmtdefs.page_orientation, gmtdefs.x_origin, gmtdefs.y_origin, gmtdefs.global_x_scale, gmtdefs.global_y_scale, gmtdefs.n_copies, gmtdefs.dpi, GMT_INCH, gmtdefs.paper_width, gmtdefs.page_rgb, GMT_epsinfo (argv[0])); GMT_echo_command (argc, argv); if (gmtdefs.unix_time) GMT_timestamp (argc, argv); if (project_info.three_D) ps_transrotate (-z_project.xmin, -z_project.ymin, 0.0); zscale = 1.0 / zscale; switch (unit) { /* Adjust for possible unit selection */ case 'C': case 'c': zscale *= GMT_u2u[GMT_CM][GMT_INCH]; break; case 'I': case 'i': zscale *= GMT_u2u[GMT_INCH][GMT_INCH]; break; case 'M': case 'm': zscale *= GMT_u2u[GMT_M][GMT_INCH]; break; case 'P': case 'p': zscale *= GMT_u2u[GMT_PT][GMT_INCH]; break; default: zscale *= GMT_u2u[gmtdefs.measure_unit][GMT_INCH]; break; } /* Now convert angles to radians and keep it that way */ start_az = (pref_az - 90.0) * D2R; stop_az = (pref_az + 90.0) * D2R; fix_az *= D2R; GMT_map_clip_on (GMT_no_rgb, 3); /* Allocate memory */ lon = (double *) GMT_memory (VNULL, (size_t)n_alloc, sizeof (double), GMT_program); lat = (double *) GMT_memory (VNULL, (size_t)n_alloc, sizeof (double), GMT_program); z = (double *) GMT_memory (VNULL, (size_t)n_alloc, sizeof (double), GMT_program); xx = (double *) GMT_memory (VNULL, (size_t)n_alloc, sizeof (double), GMT_program); yy = (double *) GMT_memory (VNULL, (size_t)n_alloc, sizeof (double), GMT_program); zz = (double *) GMT_memory (VNULL, (size_t)n_alloc, sizeof (double), GMT_program); if (n_files != 0) nofile = FALSE; done = FALSE; n_args = (argc > 1) ? argc : 2; for (fno = 1; !done && fno < n_args; fno++) { /* Loop over all the files */ if (!nofile && argv[fno][0] == '-') continue; if (nofile) { fp = GMT_stdin; done = TRUE; if (gmtdefs.verbose) fprintf (stderr, "%s: Reading data from standard input\n", GMT_program); #ifdef SET_IO_MODE GMT_setmode (0); #endif } else if ((fp = GMT_fopen (argv[fno], GMT_io.r_mode)) == NULL) { fprintf (stderr, "%s: Cannot open file %s\n", GMT_program, argv[fno]); continue; } if (!nofile && gmtdefs.verbose) { fprintf (stderr, "%s: Working on file %s\n", GMT_program, argv[fno]); sprintf (line, "File %s", argv[fno]); ps_comment (line); } if (gmtdefs.io_header) for (i = 0; i < gmtdefs.n_header_recs; i++) GMT_fgets (line, BUFSIZ, fp); if (GMT_io.multi_segments) { GMT_fgets (line, BUFSIZ, fp); if (gmtdefs.verbose) ps_comment (line); } n_expected_fields = (GMT_io.ncol[0]) ? GMT_io.ncol[0] : 3; n_fields = GMT_input (fp, &n_expected_fields, &in); while (! (GMT_io.status & GMT_IO_EOF)) { /* Not yet EOF */ while (GMT_io.status & GMT_IO_SEGMENT_HEADER) { if (gmtdefs.verbose) ps_comment (GMT_io.segment_header); n_fields = GMT_input (fp, &n_expected_fields, &in); } if (GMT_io.status & GMT_IO_EOF) continue; /* At EOF */ n = 0; while (! (GMT_io.status & (GMT_IO_SEGMENT_HEADER | GMT_IO_EOF))) { /* Keep going until FALSE or = 2 segment header */ 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, n); exit (EXIT_FAILURE); } lon[n] = in[0]; lat[n] = in[1]; z[n] = in[2] - bias; n++; if (n >= n_alloc) { n_alloc += GMT_CHUNK; lon = (double *) GMT_memory ((void *)lon, (size_t)n_alloc, sizeof (double), GMT_program); lat = (double *) GMT_memory ((void *)lat, (size_t)n_alloc, sizeof (double), GMT_program); z = (double *) GMT_memory ((void *)z, (size_t)n_alloc, sizeof (double), GMT_program); xx = (double *) GMT_memory ((void *)xx, (size_t)n_alloc, sizeof (double), GMT_program); yy = (double *) GMT_memory ((void *)yy, (size_t)n_alloc, sizeof (double), GMT_program); zz = (double *) GMT_memory ((void *)zz, (size_t)n_alloc, sizeof (double), GMT_program); } n_fields = GMT_input (fp, &n_expected_fields, &in); } GMT_geo_to_xy (lon[0], lat[0], &xx[0], &yy[0]); zz[0] = z[0]; paint_wiggle = (do_fill && ((negative && z[0] <= 0.0) || (!negative && z[0] >= 0.0))); j = 1; for (i = 1; i < n; i++) { /* Convert to inches/cm and get distance increments */ dx = lon[i] - lon[i-1]; dy = lat[i] - lat[i-1]; if (MAPPING) { /* Distance in meters */ dx *= cos (0.5 * (lat[i] + lat[i-1]) * D2R); ds = hypot (dx, dy) * project_info.M_PR_DEG; } else ds = hypot (dx, dy); GMT_geo_to_xy (lon[i], lat[i], &x_2, &y_2); if (j > 0 && (ds > dist_gap || GMT_is_dnan (z[i]))) { /* Data gap, plot what we have */ paint_wiggle = (do_fill && ((negative && zz[j-1] <= 0.0) || (!negative && zz[j-1] >= 0.0))); plot_wiggle (xx, yy, zz, j, zscale, fix, fix_az, &fill, &pen_o, &pen_t, paint_wiggle, negative, outline, track); j = 0; } else if (j >= MAX_POINTS) { paint_wiggle = (do_fill && ((negative && zz[j-1] <= 0.0) || (!negative && zz[j-1] >= 0.0))); plot_wiggle (xx, yy, zz, j, zscale, fix, fix_az, &fill, &pen_o, &pen_t, paint_wiggle, negative, outline, track); xx[0] = xx[j-1]; yy[0] = yy[j-1]; zz[0] = zz[j-1]; j = 1; } else if (!GMT_is_dnan (z[i-1]) && (z[i]*z[i-1] < 0.0 || z[i] == 0.0)) { /* Crossed 0, add new point and plot */ dz = z[i] - z[i-1]; xx[j] = (dz == 0.0) ? xx[j-1] : xx[j-1] + fabs (z[i-1] / dz) * (x_2 - xx[j-1]); yy[j] = (dz == 0.0) ? yy[j-1] : yy[j-1] + fabs (z[i-1] / dz) * (y_2 - yy[j-1]); zz[j] = 0.0; j++; paint_wiggle = (do_fill && ((negative && zz[j-2] <= 0.0) || (!negative && zz[j-2] >= 0.0))); plot_wiggle (xx, yy, zz, j, zscale, fix, fix_az, &fill, &pen_o, &pen_t, paint_wiggle, negative, outline, track); xx[0] = xx[j-1]; yy[0] = yy[j-1]; zz[0] = zz[j-1]; j = 1; } xx[j] = x_2; yy[j] = y_2; zz[j] = z[i]; if (!GMT_is_dnan (z[i])) j++; } if (j > 1) { paint_wiggle = (do_fill && ((negative && zz[j-1] <= 0.0) || (!negative && zz[j-1] >= 0.0))); plot_wiggle (xx, yy, zz, j, zscale, fix, fix_az, &fill, &pen_o, &pen_t, paint_wiggle, negative, outline, track); } } if (fp != GMT_stdin) GMT_fclose (fp); } GMT_map_clip_off(); if (frame_info.plot) GMT_map_basemap (); if (draw_scale) GMT_draw_z_scale (x0, y0, length, zscale, gave_xy, units); if (project_info.three_D) ps_rotatetrans (z_project.xmin, z_project.ymin, 0.0); ps_plotend (gmtdefs.last_page); GMT_free ((void *)lon); GMT_free ((void *)lat); GMT_free ((void *)z); GMT_free ((void *)xx); GMT_free ((void *)yy); GMT_free ((void *)zz); GMT_end (argc, argv); } void plot_wiggle (double *x, double *y, double *z, int np, double zscale, BOOLEAN fixed, double fix_az, struct GMT_FILL *fill, struct GMT_PEN *pen_o, struct GMT_PEN *pen_t, BOOLEAN paint_wiggle, BOOLEAN negative, BOOLEAN outline, BOOLEAN track) { double dx, dy, len, az = 0.0, s, c, x_inc, y_inc; int n, i; static char *name[2] = {"Positive", "Negative"}, comment[32]; if (fixed) { az = fix_az; sincos (az, &s, &c); } if (paint_wiggle || outline) { if (!GMT_n_alloc) GMT_get_plot_array (); GMT_x_plot[0] = x[0]; GMT_y_plot[0] = y[0]; n = 1; for (i = 0; i < np; i++) { if (!fixed && i < np-1) { dx = x[i+1] - x[i]; dy = y[i+1] - y[i]; if (!(dx == 0.0 && dy == 0.0)) az = -d_atan2 (dy, dx) - TWO_PI; /* Azimuth of normal to track */ while (az < start_az) az += TWO_PI; if (az > stop_az) az -= M_PI; } if (fabs (z[i]) > 0.0) { if (!fixed) sincos (az, &s, &c); len = zscale * z[i]; x_inc = len * s; y_inc = len * c; } else x_inc = y_inc = 0.0; GMT_x_plot[n] = x[i] + x_inc; GMT_y_plot[n] = y[i] + y_inc; n++; if (n == GMT_n_alloc) GMT_get_plot_array (); } GMT_x_plot[n] = x[np-1]; GMT_y_plot[n] = y[np-1]; n++; if (paint_wiggle) { for (i = np - 2; i >= 0; i--, n++) { /* Go back to 1st point along track */ if (n == GMT_n_alloc) GMT_get_plot_array (); GMT_x_plot[n] = x[i]; GMT_y_plot[n] = y[i]; } } if (project_info.three_D) GMT_2D_to_3D (GMT_x_plot, GMT_y_plot, n); } if (paint_wiggle) { /* First shade wiggles */ sprintf (comment, "%s wiggle", name[negative]); ps_comment (comment); GMT_fill (GMT_x_plot, GMT_y_plot, n, fill, FALSE); } if (outline) { /* Then draw wiggle outline */ sprintf (comment, "Wiggle line"); ps_comment (comment); GMT_setpen (pen_o); ps_line (&GMT_x_plot[1], &GMT_y_plot[1], np, 3, FALSE, TRUE); } if (track) { /* Finally draw track line */ if (project_info.three_D) GMT_2D_to_3D (x, y, np); sprintf (comment, "Track line"); ps_comment (comment); GMT_setpen (pen_t); ps_line (x, y, np, 3, FALSE, TRUE); } } void GMT_draw_z_scale (double x0, double y0, double length, double zscale, BOOLEAN gave_xy, char *units) { double dy, a, b, off; char txt[128]; GMT_setpen (&gmtdefs.tick_pen); if (!gave_xy) { GMT_geo_to_xy (x0, y0, &a, &b); x0 = a; y0 = b; } if (units) sprintf (txt, "%lg %s", length, units); else sprintf (txt, "%lg", length); dy = 0.5 * length * zscale; GMT_xyz_to_xy (x0 + gmtdefs.map_scale_height, y0 - dy, 0.0, &a, &b); ps_plot (a, b, 3); GMT_xyz_to_xy (x0, y0 - dy, 0.0, &a, &b); ps_plot (a, b, 2); GMT_xyz_to_xy (x0, y0 + dy, 0.0, &a, &b); ps_plot (a, b, 2); GMT_xyz_to_xy (x0 + gmtdefs.map_scale_height, y0 + dy, 0.0, &a, &b); ps_plot (a, b, 2); off = ((gmtdefs.map_scale_height > 0.0) ? gmtdefs.tick_length : 0.0) + gmtdefs.anot_offset; GMT_text3d (x0 + off, y0, project_info.z_level, gmtdefs.anot_font_size, gmtdefs.anot_font, txt, 0.0, 5, 0); }