/*-------------------------------------------------------------------- * $Id: segy2grd.c 12822 2014-01-31 23:39:56Z remko $ * * Copyright (c) 1991-2014 by T. Henstock * 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: gmt.soest.hawaii.edu *--------------------------------------------------------------------*/ /* * segy2grd was modified from xyz2grd by T. J. Henstock, June 2002 * to read a segy file and generate a * corresponding grd-file. * * Author: Tim Henstock (then@noc.soton.ac.uk) * Date: 30-JUN-2002 * Version: 3.4.1 */ #define THIS_MODULE_NAME "segy2grd" #define THIS_MODULE_LIB "segy" #define THIS_MODULE_PURPOSE "Converting SEGY data to a GMT grid" #include "gmt_dev.h" #include "segy_io.h" #define GMT_PROG_OPTIONS "-VRr" GMT_OPT("F") #define COUNT 1 #define AVERAGE 2 #define X_ID 0 #define Y_ID 1 #define PLOT_CDP 1 #define PLOT_OFFSET 2 struct SEGY2GRD_CTRL { struct In { /* -In */ bool active; char *file; } In; struct A { /* -A */ bool active; int mode; } A; struct C { /* -C */ bool active; double value; } C; struct D { /* -D */ bool active; char *text; } D; struct G { /* -G */ bool active; char *file; } G; struct I { /* -Idx[/dy] */ bool active; double inc[2]; } I; struct L { /* -L */ bool active; int value; } L; struct M { /* -M */ bool active; unsigned int value; } M; struct N { /* -N */ bool active; double d_value; float f_value; } N; struct Q { /* -Qx|y */ bool active[2]; double value[2]; } Q; struct S { /* -S */ bool active; unsigned int mode; int value; } S; }; void *New_segy2grd_Ctrl (struct GMT_CTRL *GMT) { /* Allocate and initialize a new control structure */ struct SEGY2GRD_CTRL *C; C = GMT_memory (GMT, NULL, 1, struct SEGY2GRD_CTRL); /* Initialize values whose defaults are not 0/false/NULL */ C->A.mode = AVERAGE; C->M.value = 10000; C->N.f_value = GMT->session.f_NaN; C->N.d_value = GMT->session.d_NaN; C->Q.value[X_ID] = 1.0; return (C); } void Free_segy2grd_Ctrl (struct GMT_CTRL *GMT, struct SEGY2GRD_CTRL *C) { /* Deallocate control structure */ if (!C) return; if (C && C->In.file) free (C->In.file); if (C && C->D.text) free (C->D.text); if (C && C->G.file) free (C->G.file); GMT_free (GMT, C); } int GMT_segy2grd_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: segy2grd -G %s\n", GMT_Id_OPT); GMT_Message (API, GMT_TIME_NONE, "\t%s [-A[n|z]]\n\t[%s] [-L]\n", GMT_Rgeo_OPT, GMT_GRDEDIT); GMT_Message (API, GMT_TIME_NONE, "\t[-M] [-N][-Q] [-S
] [%s] [%s]\n\n", GMT_V_OPT, GMT_r_OPT); if (level == GMT_SYNOPSIS) return (EXIT_FAILURE); GMT_Message (API, GMT_TIME_NONE, "\tsegyfile(s) is an IEEE floating point SEGY file. Traces are all assumed to start at 0 time/depth\n"); GMT_Message (API, GMT_TIME_NONE, "\t-G to name the output grid file.\n"); GMT_Message (API, GMT_TIME_NONE, "\t-I specifies grid size(s).\n"); GMT_Option (API, "R"); GMT_Message (API, GMT_TIME_NONE, "\n\tOPTIONS:\n"); GMT_Message (API, GMT_TIME_NONE, "\t-A (or -Az): Add multiple entries at the same node.\n"); GMT_Message (API, GMT_TIME_NONE, "\t Append n (-An): Count number of multiple entries per node instead.\n"); GMT_Message (API, GMT_TIME_NONE, "\t [Default (no -A option) will compute mean values]\n"); GMT_Message (API, GMT_TIME_NONE, "\t-D to enter header information. Specify '=' to get default value\n"); GMT_Message (API, GMT_TIME_NONE, "\t-L to override number of samples\n"); GMT_Message (API, GMT_TIME_NONE, "\t-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"); GMT_Message (API, GMT_TIME_NONE, "\t-N set value for nodes without corresponding input sample [Default is NaN]\n"); GMT_Message (API, GMT_TIME_NONE, "\t-Q can be used to change two different settings:\n"); GMT_Message (API, GMT_TIME_NONE, "\t -Qx applies scalar x-scale to coordinates in trace header to match the coordinates specified in -R\n"); GMT_Message (API, GMT_TIME_NONE, "\t -Qy specifies sample interval as if incorrect in the SEGY file\n"); GMT_Message (API, GMT_TIME_NONE, "\t-S
to set variable spacing\n"); GMT_Message (API, GMT_TIME_NONE, "\t
is c for cdp, o for offset, b for 4-byte float starting at byte number\n"); GMT_Message (API, GMT_TIME_NONE, "\t\tIf -S not set, assumes even spacing of samples at dx, dy supplied with -I\n"); GMT_Option (API, "V,r,."); return (EXIT_FAILURE); } int GMT_segy2grd_parse (struct GMT_CTRL *GMT, struct SEGY2GRD_CTRL *Ctrl, struct GMT_OPTION *options) { /* This parses the options provided to segy2grd and sets parameters in Ctrl. * Note Ctrl has already been initialized and non-zero default values set. * 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, n_files = 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) { case '<': /* Input files */ if (n_files++ > 0) break; if ((Ctrl->In.active = GMT_check_filearg (GMT, '<', opt->arg, GMT_IN))) Ctrl->In.file = strdup (opt->arg); else n_errors++; break; /* Processes program-specific parameters */ case 'A': Ctrl->A.active = true; if (opt->arg[0] == 'n') Ctrl->A.mode = COUNT; else if (opt->arg[0] == '\0' || opt->arg[0] == 'z') Ctrl->A.mode = AVERAGE; else { GMT_Report (API, GMT_MSG_NORMAL, "Syntax error -A option: Select -An or -A[z]\n"); n_errors++; } break; case 'D': Ctrl->D.active = true; Ctrl->D.text = strdup (opt->arg); break; case 'G': if ((Ctrl->G.active = GMT_check_filearg (GMT, 'G', opt->arg, GMT_OUT))) Ctrl->G.file = strdup (opt->arg); else n_errors++; break; case 'I': Ctrl->I.active = true; if (GMT_getinc (GMT, opt->arg, Ctrl->I.inc)) { GMT_inc_syntax (GMT, 'I', 1); n_errors++; } break; case 'N': if (!opt->arg[0]) { GMT_Report (API, GMT_MSG_NORMAL, "Syntax error -N option: Must specify value or NaN\n"); n_errors++; } else { Ctrl->N.d_value = (opt->arg[0] == 'N' || opt->arg[0] == 'n') ? GMT->session.d_NaN : atof (opt->arg); Ctrl->N.f_value = (float)Ctrl->N.d_value; Ctrl->N.active = true; } break; case 'Q': switch (opt->arg[0]) { case 'x': /* over-rides of header info */ Ctrl->Q.active[X_ID] = true; Ctrl->Q.value[X_ID] = atof (opt->arg); break; case 'y': /* over-rides of header info */ Ctrl->Q.active[Y_ID] = true; Ctrl->Q.value[Y_ID] = atof (opt->arg); break; } break; case 'L': Ctrl->L.active = true; Ctrl->L.value = atoi (opt->arg); break; case 'M': Ctrl->M.active = true; Ctrl->M.value = atoi (opt->arg); break; /* variable spacing */ case 'S': if (Ctrl->S.active) { GMT_Report (API, GMT_MSG_NORMAL, "Syntax error -S option: Can only be set once\n"); n_errors++; } Ctrl->S.active = true; switch (opt->arg[0]) { case 'o': Ctrl->S.mode = PLOT_OFFSET; break; case 'c': Ctrl->S.mode = PLOT_CDP; break; case 'b': Ctrl->S.value = atoi (&opt->arg[1]); break; } break; default: /* Report bad options */ n_errors += GMT_default_error (GMT, opt->option); break; } } GMT_check_lattice (GMT, Ctrl->I.inc, &GMT->common.r.registration, &Ctrl->I.active); n_errors += GMT_check_condition (GMT, !GMT->common.R.active, "Syntax error: Must specify -R option\n"); n_errors += GMT_check_condition (GMT, Ctrl->I.inc[GMT_X] <= 0.0 || Ctrl->I.inc[GMT_Y] <= 0.0, "Syntax error -I option: Must specify positive increment(s)\n"); n_errors += GMT_check_condition (GMT, !Ctrl->G.active || !Ctrl->G.file, "Syntax error -G: Must specify output file\n"); n_errors += GMT_check_condition (GMT, !Ctrl->G.active || !Ctrl->G.file, "Syntax error -G: Must specify output file\n"); return (n_errors ? GMT_PARSE_ERROR : GMT_OK); } #define bailout(code) {GMT_Free_Options (mode); return (code);} #define Return(code) {Free_segy2grd_Ctrl (GMT, Ctrl); GMT_end_module (GMT, GMT_cpy); bailout (code);} int GMT_segy2grd (void *V_API, int mode, void *args) { bool read_cont = false, swap_bytes = !GMT_BIGENDIAN; int error = 0; unsigned int ii, jj, n_read = 0, n_filled = 0, n_used = 0, *flag = NULL; unsigned int n_empty = 0, n_stuffed = 0, n_bad = 0, n_confused = 0, check, ix; uint64_t ij, ij0, n_samp = 0, isamp; double idy, x0, yval; char line[GMT_BUFSIZ] = {""}; FILE *fpi = NULL; struct GMT_GRID *Grid = NULL; /* SEGY parameters */ char reelhead[3200] = {""}; float *data = NULL; SEGYHEAD *header = NULL; SEGYREEL binhead; struct SEGY2GRD_CTRL *Ctrl = NULL; struct GMT_CTRL *GMT = NULL, *GMT_cpy = NULL; /* General GMT interal parameters */ struct GMT_OPTION *options = NULL; struct GMTAPI_CTRL *API = GMT_get_API_ptr (V_API); /* Cast from void to GMTAPI_CTRL pointer */ /*----------------------- Standard module initialization and parsing ----------------------*/ if (API == NULL) return (GMT_NOT_A_SESSION); if (mode == GMT_MODULE_PURPOSE) return (GMT_segy2grd_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_segy2grd_usage (API, GMT_USAGE)); /* Return the usage message */ if (options->option == GMT_OPT_SYNOPSIS) bailout (GMT_segy2grd_usage (API, GMT_SYNOPSIS)); /* Return the synopsis */ /* Parse the command-line arguments; return if errors are encountered */ 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_segy2grd_Ctrl (GMT); /* Allocate and initialize a new control structure */ if ((error = GMT_segy2grd_parse (GMT, Ctrl, options))) Return (error); /*---------------------------- This is the segy2grd main code ----------------------------*/ read_cont = (Ctrl->S.mode != PLOT_CDP && Ctrl->S.mode != PLOT_OFFSET && !Ctrl->S.value); if ((Grid = GMT_Create_Data (API, GMT_IS_GRID, GMT_IS_SURFACE, GMT_GRID_ALL, NULL, NULL, Ctrl->I.inc, \ GMT_GRID_DEFAULT_REG, GMT_NOTSET, Ctrl->G.file)) == NULL) Return (API->error); /* Decode grd information given, if any */ if (Ctrl->D.active) GMT_decode_grd_h_info (GMT, Ctrl->D.text, Grid->header); GMT_Report (API, GMT_MSG_VERBOSE, "nx = %d ny = %d\n", Grid->header->nx, Grid->header->ny); flag = GMT_memory (GMT, NULL, Grid->header->size, unsigned int); GMT_grd_pad_off (GMT, Grid); /* Undo pad since algorithm does not expect on */ idy = 1.0 / Grid->header->inc[GMT_Y]; ij = -1; /* Will be incremented to 0 or set first time around */ /* read in reel headers from segy file */ if (Ctrl->In.active) { GMT_Report (API, GMT_MSG_VERBOSE, "Will read segy file %s\n", Ctrl->In.file); if ((fpi = fopen (Ctrl->In.file, "rb")) == NULL) { GMT_Report (API, GMT_MSG_NORMAL, "Cannot find segy file %s\n", Ctrl->In.file); Return (EXIT_FAILURE); } } else { GMT_Report (API, GMT_MSG_VERBOSE, "Will read segy file from standard input\n"); if (fpi == NULL) fpi = stdin; } 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*/ GMT_Report (API, GMT_MSG_VERBOSE, "Swapping bytes for ints in the headers\n"); binhead.num_traces = bswap16 (binhead.num_traces); binhead.nsamp = bswap16 (binhead.nsamp); binhead.dsfc = bswap16 (binhead.dsfc); binhead.sr = bswap16 (binhead.sr); } /* set parameters from the reel headers */ if (!Ctrl->M.value) Ctrl->M.value = binhead.num_traces; GMT_Report (API, GMT_MSG_VERBOSE, "Number of traces in header is %d\n", Ctrl->M.value); if (!Ctrl->L.value) { /* number of samples not overridden*/ Ctrl->L.value = binhead.nsamp; GMT_Report (API, GMT_MSG_VERBOSE, "Number of samples per trace is %d\n", Ctrl->L.value); } else if ((Ctrl->L.value != binhead.nsamp) && (binhead.nsamp)) GMT_Report (API, GMT_MSG_VERBOSE, "Warning nsampr input %d, nsampr in header %d\n", Ctrl->L.value, binhead.nsamp); if (!Ctrl->L.value) { /* no number of samples still - a problem! */ GMT_Report (API, GMT_MSG_NORMAL, "Error, number of samples per trace unknown\n"); Return (EXIT_FAILURE); } GMT_Report (API, GMT_MSG_VERBOSE, "Number of samples for reel is %d\n", Ctrl->L.value); if (binhead.dsfc != 5) GMT_Report (API, GMT_MSG_VERBOSE, "Warning: data not in IEEE format\n"); if (!Ctrl->Q.value[Y_ID]) { Ctrl->Q.value[Y_ID] = (double) binhead.sr; /* sample interval of data (microseconds) */ Ctrl->Q.value[Y_ID] /= 1000000.0; GMT_Report (API, GMT_MSG_VERBOSE,"Sample interval is %f s\n", Ctrl->Q.value[Y_ID]); } else if ((Ctrl->Q.value[Y_ID] != binhead.sr) && (binhead.sr)) /* value in header overridden by input */ GMT_Report (API, GMT_MSG_VERBOSE, "Warning s_int input %f, s_int in header %f\n", Ctrl->Q.value[Y_ID], (float)binhead.sr); if (!Ctrl->Q.value[Y_ID]) { /* still no sample interval at this point is a problem! */ GMT_Report (API, GMT_MSG_NORMAL, "Error, no sample interval in reel header\n"); exit(EXIT_FAILURE); } if (read_cont && (Ctrl->Q.value[Y_ID] != Grid->header->inc[GMT_Y])) { GMT_Report (API, GMT_MSG_VERBOSE, "Warning, grid spacing != sample interval, setting sample interval to grid spacing\n"); Ctrl->Q.value[Y_ID] = Grid->header->inc[GMT_Y]; } if (Grid->header->inc[GMT_Y] < Ctrl->Q.value[Y_ID]) GMT_Report (API, GMT_MSG_VERBOSE, "Warning, grid spacing < sample interval, expect gaps in output....\n"); /* starts reading actual data here....... */ if (read_cont) { /* old-style segy2grd */ ix = 0; for (ij = 0; ij < Grid->header->size; ij++) Grid->data[ij] = Ctrl->N.f_value; if (Grid->header->nx < Ctrl->M.value) { GMT_Report (API, GMT_MSG_VERBOSE, "Warning, number of traces in header > size of grid. Reading may be truncated\n"); Ctrl->M.value = Grid->header->nx; } while ((ix < Ctrl->M.value) && (header = get_segy_header (fpi))) { if (swap_bytes) { /* need to permanently byte-swap number of samples in the trace header */ header->num_samps = bswap32 (header->num_samps); header->sampleLength = bswap16 (header->sampleLength); } data = get_segy_data (fpi, header); /* read a trace */ /* get number of samples in _this_ trace or set to number in reel header */ if (!(n_samp = samp_rd (header))) n_samp = Ctrl->L.value; ij0 = lrint (GMT->common.R.wesn[YLO] * idy); if ((n_samp - ij0) > (uint64_t)Grid->header->ny) n_samp = Grid->header->ny + ij0; if (swap_bytes) { /* need to swap the order of the bytes in the data even though assuming IEEE format */ uint32_t tmp; for (isamp = 0; isamp < n_samp; ++isamp) { memcpy (&tmp, &data[isamp], sizeof(uint32_t)); tmp = bswap32 (tmp); memcpy (&data[isamp], &tmp, sizeof(uint32_t)); } } for (ij = ij0; ij < n_samp ; ij++) { /* n*idy is index of first sample to be included in the grid */ Grid->data[ix + Grid->header->nx*(Grid->header->ny+ij0-ij-1)] = data[ij]; } free (data); free (header); ix++; } } else { /* Get trace data and position by headers */ ix = 0; while ((ix < Ctrl->M.value) && (header = get_segy_header (fpi))) { /* read traces one by one */ if (Ctrl->S.mode == PLOT_OFFSET) { /* plot traces by offset, cdp, or input order */ int32_t tmp = header->sourceToRecDist; if (swap_bytes) { uint32_t *p = (uint32_t *)&tmp; *p = bswap32 (*p); } x0 = (double) tmp; } else if (Ctrl->S.mode == PLOT_CDP) { int32_t tmp = header->cdpEns; if (swap_bytes) { uint32_t *p = (uint32_t *)&tmp; *p = bswap32 (*p); } x0 = (double) tmp; } else if (Ctrl->S.value) { /* get value starting at Ctrl->S.value of header into a double */ uint32_t tmp; memcpy (&tmp, &header[Ctrl->S.value], sizeof (uint32_t)); x0 = (swap_bytes) ? (double) bswap32 (tmp) : (double) tmp; } else x0 = (1.0 + (double) ix); x0 *= Ctrl->Q.value[X_ID]; 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 Ctrl->S.value case overlaps a defined header in a strange way */ uint32_t *p = (uint32_t *)&header->sourceToRecDist; *p = bswap32 (*p); header->sampleLength = bswap16 (header->sampleLength); header->num_samps = bswap32 (header->num_samps); } data = 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 = Ctrl->L.value; if (swap_bytes) { /* need to swap the order of the bytes in the data even though assuming IEEE format */ uint32_t tmp; for (isamp = 0; isamp < n_samp; ++isamp) { memcpy (&tmp, &data[isamp], sizeof(uint32_t)); tmp = bswap32 (tmp); memcpy (&data[isamp], &tmp, sizeof(uint32_t)); } } if (!(x0 < GMT->common.R.wesn[XLO] || x0 > GMT->common.R.wesn[XHI])) { /* inside x-range */ /* find horizontal grid pos of this trace */ ii = (unsigned int)GMT_grd_x_to_col (GMT, x0, Grid->header); if (ii == Grid->header->nx) ii--, n_confused++; for (isamp = 0; isamp< n_samp; ++isamp) { yval = isamp*Ctrl->Q.value[Y_ID]; if (!(yval < GMT->common.R.wesn[YLO] || yval > GMT->common.R.wesn[YHI])) { /* inside y-range */ jj = (unsigned int)GMT_grd_y_to_row (GMT, yval, Grid->header); if (jj == Grid->header->ny) jj--, n_confused++; ij = GMT_IJ0 (Grid->header, jj, ii); Grid->data[ij] += data[isamp]; /* Add up incase we must average */ flag[ij]++; n_used++; } } } free (data); free (header); ix++; } for (ij = 0; ij < Grid->header->nm; ij++) { /* Check if all nodes got one value only */ if (flag[ij] == 1) { if (Ctrl->A.mode == COUNT) Grid->data[ij] = 1.0f; n_filled++; } else if (flag[ij] == 0) { n_empty++; Grid->data[ij] = Ctrl->N.f_value; } else { /* More than 1 value went to this node */ if (Ctrl->A.mode == COUNT) Grid->data[ij] = (float)flag[ij]; else if (Ctrl->A.mode == AVERAGE) Grid->data[ij] /= (float)flag[ij]; n_filled++; n_stuffed++; } } if (GMT_is_verbose (GMT, GMT_MSG_VERBOSE)) { sprintf (line, "%s\n", GMT->current.setting.format_float_out); GMT_Message (API, GMT_TIME_NONE, " n_read: %d n_used: %d n_filled: %d n_empty: %d set to ", n_read, n_used, n_filled, n_empty); (GMT_is_dnan (Ctrl->N.d_value)) ? GMT_Message (API, GMT_TIME_NONE, "NaN\n") : GMT_Message (API, GMT_TIME_NONE, line, Ctrl->N.d_value); if (n_bad) GMT_Message (API, GMT_TIME_NONE, "%d records unreadable\n", n_bad); if (n_stuffed) GMT_Message (API, GMT_TIME_NONE, "Warning - %d nodes had multiple entries that were averaged\n", n_stuffed); if (n_confused) GMT_Message (API, GMT_TIME_NONE, "Warning - %d values gave bad indices: Pixel vs gridline confusion?\n", n_confused); } } GMT_grd_pad_on (GMT, Grid, GMT->current.io.pad); /* Restore padding */ if (GMT_Set_Comment (API, GMT_IS_GRID, GMT_COMMENT_IS_OPTION | GMT_COMMENT_IS_COMMAND, options, Grid)) Return (API->error); if (GMT_Write_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_GRID_ALL, NULL, Ctrl->G.file, Grid) != GMT_OK) { Return (API->error); } GMT_free (GMT, flag); Return (GMT_OK); }