/*-------------------------------------------------------------------- * $Id: pssegy.c,v 1.3 2001/04/11 19:58:09 pwessel Exp $ * * Copyright (c) 1999-2001 by T. Henstock * See README file for copying and redistribution conditions. *--------------------------------------------------------------------*/ /* pssegy program to plot segy files in postscript with variable trace spacing option * uses the GMT pslib imagemask routines to plot a * 1-bit depth bitmap which will not obliterate stuff underneath! * * Author: Tim Henstock (then@soc.soton.ac.uk) * Date: 1-July-1996 * Version: 1.0 * * Bug fixes: 1.1, 11-6-96 remove undesirable normalization introduced by -Z option. Add -U option for reduction. * * enhancements: 1.2, 7-20-98 add option to match location of traces from file * * 1.3, 1/7/99 check number of samples trace by trace to cope with SEGY with variable trace length * * 2.0, 5/7/99 update for GMT 3.3.1 * * 2.1 10/4/2001 fix unreduce bug, modify to byte-swap if necessary, make 64-bit safe * * This program is free software and may be copied or redistributed under the terms * of the GNU public license. */ #include "gmt.h" #include "segy_io.h" /* internal function calls */ double rms( float *data, int nsamp); void wig_bmap(double x0, float data0, float data1, double y0, double y1); void shade_bmap(double x0, float data0, float data1, double y0, double y1, int negative); int paint( int ix, int iy); void plot_trace(float *data, double dy, double x0, int n_samp, int do_fill, int negative, int plot_wig, float toffset); unsigned char bmask[8]={128, 64, 32, 16, 8, 4, 2, 1}; unsigned char *bitmap; int bm_nx, bm_ny; main (int argc, char **argv) { int error = FALSE; int plot_cdp = FALSE, plot_offset = FALSE, byte_x = 0, doclip = FALSE; int normalize = FALSE, do_fill = FALSE, negative = FALSE, plot_wig = FALSE; int no_zero = FALSE, trace_file = FALSE; int swap_bytes = !WORDS_BIGENDIAN; int i, nm; int ix, iy, n_traces=10000, n_samp=0, n_sampr=0, shade[3], n_tracelist=0; int check, plot_it, polarity=1; double w, e, s, n, dx = 1.0, dy = 0.0; /* dx, dy are trace and sample interval */ double xlen, ylen, xpix, ypix; double x0, test, err_head=0.0; float bias = 0.0, scale = 1.0, clip = 0.0, deviation = 0.0; double redvel=0.0; double *tracelist; float toffset=0.0; char input[512]=""; char reelhead[3200]; float *data; SEGYHEAD *header; char *head; long *head2; SEGYREEL binhead; FILE *fpi = NULL, *fpt = NULL; input[0] = 0; w = e = s = n = 0.0; shade[0]=shade[1]=shade[2]=0; argc = GMT_begin (argc, argv); for (i = 1; i < argc; i++) { if (argv[i][0] == '-') { switch (argv[i][1]) { /* Common parameters */ case 'V': case 'R': case 'J': case 'O': case 'K': case 'P': case '\0': error += GMT_get_common_args (argv[i], &w, &e, &s, &n); break; /* parameters for wiggle mode */ case 'F': do_fill = TRUE; if (GMT_getrgb (&argv[i][2], shade)) { error++; GMT_rgb_syntax ('F'); } break; case 'I': negative = TRUE; break; case 'W': plot_wig = TRUE; break; /* trace norm., clip and bias */ case 'N': normalize = TRUE; break; case 'C': doclip = TRUE; clip = (float) atof (&argv[i][2]); break; case 'B': bias = (float) atof (&argv[i][2]); break; case 'Z': no_zero = TRUE; break; /* variable spacing */ case 'S': switch(argv[i][2]){ case 'o': plot_offset = TRUE; break; case 'c': plot_cdp = TRUE; break; case 'b': byte_x = atoi (&argv[i][3]); break; } break; /* trace scaling */ case 'D': deviation = (float) atof (&argv[i][2]); break; /* over-rides of header info */ case 'X': /* -X and -Y can be changed in gmt routines to lower case...*/ case 'x': dx = atof (&argv[i][2]); break; case 'Y': case 'y': dy = atof (&argv[i][2]); fprintf(stderr," %s: Overriding sample interval dy = %f\n", GMT_program, dy); break; case 'L': n_sampr = atoi (&argv[i][2]); break; case 'M': n_traces = atoi (&argv[i][2]); break; /* reduction velocity application */ case 'U': redvel = atof (&argv[i][2]); break; /* plot traces only at listed locations */ case 'T': strcpy (input, &argv[i][2]); if (trace_file = strlen(input)){ if ((fpt = fopen (input, "r")) == NULL){ fprintf(stderr, "%s: Cannot find trace list file %s\n", GMT_program, input); error++; } } else{ error++; fprintf(stderr, "%s: must give trace list filename\n", GMT_program); } break; case 'E': err_head = atof (&argv[i][2]); break; case 'A': swap_bytes = !swap_bytes; break; default: error = TRUE; break; } } else if ((fpi = fopen (argv[i], "rb")) == NULL) { fprintf (stderr, "%s: Cannot find segy file %s\n", GMT_program,argv[i]); exit (EXIT_FAILURE); } } if (argc == 1 || GMT_quick) { fprintf (stderr, "pssegy %s - Plot a segy file in PostScript\n\n", GMT_VERSION); fprintf (stderr, "usage: pssegy [] -J -R -D \n"); fprintf (stderr, " [-C] [-B] [-N] [-Z]\n"); fprintf (stderr, " [-F|] [-I] [-W] [-S
]\n"); fprintf (stderr, " [-X] [-Y] [-L] [-M] \n"); fprintf (stderr, " [-U] [-T] [-E] [-A] [-O] [-K] [-P]\n"); if (GMT_quick) exit (EXIT_FAILURE); fprintf (stderr, "\n\t-Jx for projection. Scale in INCH/units. Specify one:\n\t -Jx Linear projection\n\t-Jxl Log10 projection\n\t -Jxp x^power projection\n\tUse / to specify separate x/y scaling.\n\t If -JX is used then give axes lengths rather than scales\n\t regular map projections are not allowed!\n"); GMT_explain_option ('R'); fprintf (stderr, " NB units for y are s or km\n"); fprintf (stderr, " -D to give deviation in X units of plot for 1.0 on scaled trace\n"); fprintf (stderr, " IEEE SEGY file [or standard input] \n\n"); fprintf (stderr, "\n\tOPTIONS:\n"); GMT_explain_option ('V'); fprintf (stderr, " -C to clip scaled trace excursions at , applied after bias\n"); fprintf (stderr, " -B to bias scaled traces (-B-0.1 subtracts 0.1 from values)\n"); fprintf (stderr," -N to trace normalize the plot\n"); fprintf (stderr," order of operations: [normalize][bias][clip](deviation)\n"); fprintf (stderr," -Z to suppress plotting traces whose rms amplitude is 0 \n"); fprintf (stderr," -F to fill variable area with shade \n"); fprintf (stderr," -Fr/g/b to fill variable area with color\n"); fprintf (stderr," only a single color for the bitmap though!\n"); fprintf (stderr," -I to fill negative rather than positive excursions\n"); fprintf (stderr," -W to plot wiggle trace\n"); fprintf (stderr," must specify either -W or -F\n"); fprintf (stderr," -X multiply trace locations by \n"); fprintf (stderr," -Y to override sample interval\n"); fprintf (stderr," -S
to set variable spacing\n"); fprintf (stderr,"
is c for cdp or o for offset\n"); fprintf (stderr," -L to override number of samples\n"); fprintf (stderr," -M to fix number of traces. Default reads all traces.\n\t\t-M0 will read number in binary header, -Mn will attempt to read only n traces.\n"); fprintf (stderr," -U to apply reduction velocity (-ve removes reduction already present)\n"); fprintf (stderr," -T to look in filename for a list of locations to select traces\n"); fprintf (stderr," (same units as header * X, ie values printed by previous -V run)\n"); fprintf (stderr," -E slop to allow for -T. recommended in case of arithmetic errors!\n"); fprintf (stderr," -A flips the default byte-swap state (default assumes data have a bigendian byte-order)\n"); GMT_explain_option ('O'); GMT_explain_option ('K'); GMT_explain_option ('P'); exit (EXIT_FAILURE); } if (err_head < 0.0){ error++; fprintf(stderr,"%s: SYNTAX ERROR: slop cannot be negative\n", GMT_program); } if (negative && !do_fill){ /* negative with no fill */ error++; fprintf(stderr,"%s: SYNTAX ERROR: Must specify -F with -I\n", GMT_program); } if (!do_fill && !plot_wig){ /* no plotting specified */ error++; fprintf(stderr,"%s: SYNTAX ERROR: Must specify -F or -W\n", GMT_program); } if (deviation <= 0.0){ error++; fprintf(stderr,"%s: SYNTAX ERROR: Must specify a positive deviation\n",GMT_program); } if (project_info.projection >2){ fprintf(stderr,"%s: WARNING: you asked for a non-rectangular projection. \n It will probably still work, but be prepared for problems\n",GMT_program); } if (plot_cdp && plot_offset){ fprintf(stderr,"%s: SYNTAX ERROR: Can't specify more than one trace location key\n",GMT_program); error++; } if (error) exit (EXIT_FAILURE); if (fpi == NULL) fpi = stdin; if (trace_file){ /* must read in file of desired trace locations */ tracelist = (double *) GMT_memory (CNULL, GMT_CHUNK, sizeof(double), "pssegy"); n_tracelist = GMT_CHUNK; ix=0; while ((fscanf(fpt, "%lf", &test)) != EOF){ tracelist[ix] = test; ix++; if(ix == n_tracelist){ /* need more memory in array */ n_tracelist += GMT_CHUNK; tracelist = (double *) GMT_memory ((char *)tracelist, n_tracelist, sizeof(double), "pssegy"); } } n_tracelist = ix; if(gmtdefs.verbose) fprintf(stderr, "%s: read in %d trace locations\n", GMT_program, n_tracelist); } GMT_put_history (argc, argv); /* Update .gmtcommands */ /* set up map projection and PS plotting */ GMT_map_setup (w, e, s, n); 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); /* define area for plotting and size of array for bitmap */ xlen = project_info.xmax-project_info.xmin; xpix = xlen*gmtdefs.dpi; /* pixels in x direction */ /*xpix /= 8.0; bm_nx = 1 +(int) xpix;*/ bm_nx = (int) ceil (xpix/8.0); /* store 8 pixels per byte in x direction but must have whole number of bytes per scan */ ylen = project_info.ymax-project_info.ymin; ypix = ylen*gmtdefs.dpi; /* pixels in y direction */ bm_ny = (int) ypix; nm = bm_nx*bm_ny; /* read in reel headers from segy file */ if ((check = get_segy_reelhd (fpi, reelhead)) != TRUE) exit(1); if ((check = get_segy_binhd (fpi, &binhead)) != TRUE) exit(1); if(swap_bytes){ /* this is a little-endian system, and we need to byte-swap ints in the reel header - we only use a few of these*/ if (gmtdefs.verbose) fprintf(stderr, "%s: swapping bytes for ints in the headers\n",GMT_program); binhead.num_traces = GMT_swab2(binhead.num_traces); binhead.nsamp = GMT_swab2(binhead.nsamp); binhead.dsfc = GMT_swab2(binhead.dsfc); binhead.sr = GMT_swab2(binhead.sr); } /* set parameters from the reel headers */ if (!n_traces) n_traces = binhead.num_traces; if (gmtdefs.verbose) fprintf(stderr, "%s: Number of traces in header is %d\n", GMT_program, n_traces); if (!n_sampr){/* number of samples not overridden*/ n_sampr = binhead.nsamp; fprintf(stderr,"%s: Number of samples per trace is %d\n", GMT_program, n_sampr); } else if ((n_sampr != binhead.nsamp) && (binhead.nsamp)) fprintf(stderr,"%s: warning nsampr input %d, nsampr in header %d\n", GMT_program, n_sampr, binhead.nsamp); if (!n_sampr){ /* no number of samples still - a problem! */ fprintf(stderr, "%s: Error, number of samples per trace unknown\n", GMT_program); exit(-1); } if(gmtdefs.verbose) fprintf(stderr, "%s: Number of samples for reel is %d\n", GMT_program, n_sampr); if(binhead.dsfc != 5) fprintf(stderr, "pssegy: WARNING data not in IEEE format\n"); if (!dy){ dy = (double) binhead.sr; /* sample interval of data (microseconds) */ dy /= 1000000.0; fprintf(stderr,"%s: Sample interval is %f s\n", GMT_program, dy); } else if ((dy != binhead.sr) && (binhead.sr)) /* value in header overridden by input */ fprintf(stderr, "%s: Warning dy input %lf, dy in header %f\n", GMT_program, dy, (float)binhead.sr); if (!dy){ /* still no sample interval at this point is a problem! */ fprintf(stderr, "%s: Error, no sample interval in reel header\n", GMT_program); exit(-1); } bitmap = (unsigned char *) GMT_memory (NULL, nm, sizeof (unsigned char), "pssegy"); ix=0; while ((ixsourceToRecDist): header->sourceToRecDist); x0 = (double) offset; } else if (plot_cdp){ long cdpval = ((swap_bytes)? GMT_swab4(header->cdpEns): header->cdpEns); x0 = (double) cdpval; } else if (byte_x){ /* ugly code - want to get value starting at byte_x of header into a double... */ long plot_loc; head = (char *) header; head2 = (long *) &(head[byte_x]); plot_loc = ((swap_bytes)? GMT_swab4(*head2): *head2); x0 = (double) plot_loc; } else x0 = (1.0 + (double) ix); x0 *= dx; if (swap_bytes){ /* need to permanently byte-swap some things in the trace header do this after getting the location of where traces are plotted in case the general byte_x case overlaps a defined header in a strange way */ header->sourceToRecDist=GMT_swab4(header->sourceToRecDist); header->sampleLength=GMT_swab2(header->sampleLength); header->num_samps=GMT_swab4(header->num_samps); } /* now check that on list to plot if list exists */ if (n_tracelist){ plot_it = FALSE; for (i=0; i< n_tracelist; i++){ if(fabs(x0-tracelist[i])<=err_head){ plot_it = TRUE; } } } if (redvel){ toffset = (float) -(fabs((double)(header->sourceToRecDist))/redvel); if (gmtdefs.verbose) fprintf(stderr, "pssegy: time shifted by %f\n",toffset); } data = (float *) get_segy_data(fpi, header); /* read a trace */ /* get number of samples in _this_ trace (e.g. OMEGA has strange ideas about SEGY standard) or set to number in reel header */ if ( !(n_samp = samp_rd(header)) ) n_samp = n_sampr; if(swap_bytes){ /* need to swap the order of the bytes in the data even though assuming IEEE format */ int *intdata = (int *) data; for (iy=0; iy clip)) data[iy] = (float)(clip*data[iy]/fabs(data[iy])); /* apply bias and then clip */ data[iy] *= deviation; } if ((!no_zero || scale) && (plot_it || !n_tracelist)){ if (gmtdefs.verbose) fprintf(stderr, "pssegy: trace %d plotting at %f \n", ix+1, x0); plot_trace (data, dy, x0, n_samp, do_fill, negative, plot_wig, toffset); } free (data); free (header); ix++; } GMT_map_clip_on (GMT_no_rgb, 3); /* set a clip at the map boundary since the image space overlaps a little */ ps_imagemask (0.0,0.0,xlen, ylen, bitmap, 8*bm_nx, bm_ny, polarity, shade); /* have to multiply by 8 since pslib version of ps_imagemask is based on a _pixel_ count, whereas pssegy uses _byte_ count internally */ GMT_map_clip_off (); if (fpi != stdin) fclose (fpi); ps_plotend (gmtdefs.last_page); GMT_end (argc, argv); } double rms (float *data, int n_samp) {/* function to return rms amplitude of n_samp values from the array data */ int ix; double sumsq=0.0; for (ix=0; ix=0.)||(data[iy]>=0.))) || (negative && ((data[iy-1]<=0.0)||(data[iy]<=0.0))) ); if (paint_wiggle) shade_bmap (x0, data[iy-1], data[iy], y0, y1, negative); } y0=y1; } } void wig_bmap(double x0, float data0, float data1, double y0, double y1) /* apply current sample with all options to bitmap */ { double xp0, xp1, yp0, yp1, slope; int px0, px1, py0, py1, ix, iy; GMT_geo_to_xy (x0+ (double)data0, y0, &xp0, &yp0); /* returns 2 ends of line segment in plot coords */ GMT_geo_to_xy (x0+ (double)data1, y1, &xp1, &yp1); slope = (yp1-yp0)/(xp1-xp0); px0 = (int) (xp0*gmtdefs.dpi); px1 = (int) (xp1*gmtdefs.dpi); py0 = (int) (yp0*gmtdefs.dpi); py1 = (int) (yp1*gmtdefs.dpi); /* now have the pixel locations for the two samples - join with a line..... */ if (fabs(slope) <= 1.0){ /* more pixels needed in x direction */ if (px00.0)&& !negative)) { /* plot from top to zero */ y1=interp; data1=0.0; } else { y0=interp; data0=0.0; } } GMT_geo_to_xy (x0+(double)data0, y0, &xp0, &yp0); /* returns 2 ends of line segment in plot coords */ GMT_geo_to_xy (x0+(double)data1, y1, &xp1, &yp1); GMT_geo_to_xy (x0, y0, &xp00, &yp0); /* to get position of zero */ slope = (yp1-yp0)/(xp1-xp0); px0 = (int) (0.49+xp0*gmtdefs.dpi); px00 = (int) (0.49+xp00*gmtdefs.dpi); py0 = (int) (0.49+yp0*gmtdefs.dpi); py1 = (int) (0.49+yp1*gmtdefs.dpi); /* can rasterize simply by looping over values of y */ if (py0= bm_nx-1) || (iy >= bm_ny-1) || (ix < 0) || (iy < 0)) return (-1); /* outside bounds of plot array */ byte = (bm_ny-iy-1)*bm_nx + quot; /* find byte to paint - flip vertical! */ bitmap[byte] = bitmap[byte] | bmask[rem]; return (0); }