/******************************** GPR_STAK.C ********************************\ * version 1.08.09.01 * GPR_STAK reduces the size of digital ground penetrating radar data by * stacking traces. A beginning trace and end trace must be specified * and the number of traces to stack. The following computer storage * formats are recognized. * * GSSI SIR-10A version 3.x to 5.x binary files, 8-bit or 16-bit unsigned * integers with an embedded descriptive 512- or 1024-byte header. * Data files usually have the extension "DZT". * Sensors and Software pulseEKKO binary files, 16-bit signed integers * accompanied by a descriptive ASCII "HD" file. Data files must have * the extensiom "DT1". * SEG-Y formatted files, 16-bit or 32-bit signed integers or 32-bit * floating point reals. Data files usually have the extension "SGY". * * Jeff Lucius * Geophysicist * U.S. Geological Survey * Geologic Division, Branch of Geophysics * Box 25046 Denver Federal Center MS 964 * Denver, CO 80225-0046 * 303-236-1413 (office), 303-236-1212 (secretary), 303-236-1425 (fax) * email: lucius@musette.cr.usgs.gov * * Revision History: * February 15, 1996 - Completed as a modification of GPR_CNDS.C * February 16, 1996 - Added exit() if GetGprFileType() returns error. * August 11, 1997 - Added support for modified DZT formats. * August 9, 2001 - Updated time and date functions * * Supported compilers: * Intel 386/486 C Code Builder 1.1a * To compile for 80386 executable: * icc /F /xnovm /zfloatsync gpr_stak.c gpr_io.lib jl_util1.lib * To compile for 80486 executable: * icc /F /xnovm /zfloatsync /zmod486 gpr_stak.c gpr_io.lib jl_util1.lib * * /F removes the floating point emulator, reduces binary size. * /xnovm specifies no virtual memory manager, disk access is too slow. * /zfloatsync insures the FPU is operand-synchronized with the CPU. * /zmod486 generates specific 80486 instructions, speeds up program. * Jeff Lucius's GPR I/O functions (DZT, DT1, SGY) are in gpr_io.lib. * Jeff Lucius's utility functions are in jl_util1.lib. * * to remove assert() functions from code also add: * /DNDEBUG * * To run: GPR_STAK in_filename out_filename start end stack [/d /b] * \****************************************************************************/ /********************** Include required header files ***********************/ /* non-ANSI headers */ #include /* getch */ #include /* errno */ #include /* access */ #include /* _dos_getdate, _dos_gettime, struct dosdate_t, struct dostime_t */ /* ANSI-compatible headers */ #include /* isdigit */ #include /* FLT_MAX */ #include /* puts,printf,ferror,clearerr,sprintf,fopen,fclose, fseek,setvbuf */ #include /* atoi */ #include /* strcpy,strupr,strstr,strcat,strchr,strncpy */ /* application headers */ #include "gpr_io.h" /* to use gpr_io.lib functions, etc. */ #include "jl_defs.h" /* common manifest constants and macros */ #include "jl_util1.h" /* ANSI_there, etc. */ /**************************** Global variables ******************************/ int Debug = FALSE; /* if TRUE, turn debugging on */ int Batch = FALSE; /* if TRUE, supress interaction with user */ /*************************** Manifest constants ****************************/ #define GPR_STAK_VER "1.08.09.01" #define DATETIME "August 9, 2001" /********************************** main() **********************************/ void main (int argc,char *argv[]) { char in_filename[MAX_PATHLEN],out_filename[MAX_PATHLEN]; int itemp, min_args; int start, end, stack; int file_hdr_bytes, num_hdrs, trace_hdr_bytes, samples_per_trace; int num_channels, input_datatype, total_time, file_type; long num_traces; /* following in errno.h */ extern int errno; /* Verify usage */ min_args = 6; if (argc < min_args) { puts("\n"); if (ANSI_there()) printf("%c%c32;40m",0x1B,0x5B); /* green on black */ puts(""); puts("###############################################################################"); puts(" Usage: GPR_STAK in_filename out_filename start end stack [/d /b]"); puts(" Computer code written by Jeff Lucius USGS lucius@usgs.gov"); puts(" This program reduces the size of GPR files by stacking traces."); puts(" Recognized storage formats include: GSSI DZT, S&S DT1/HD, and SEG-Y."); puts("\n Required command line arguments:"); puts(" in_filename - The name of the input GPR data file."); puts(" out_filename - The name of the output GPR data file."); puts(" start - The first trace to use from the file (indexed from 0)."); puts(" If -1, then start with the first trace in the file."); puts(" end - The last trace to use from the file (indexed from 0)."); puts(" If -1, then end with the last trace in the file."); puts(" stack - The number of traces to stack to construct a new trace."); puts(" Markers are preserved for DZT files."); puts(" Trace numbers are updated for DT1 and SGY files."); puts(" Optional command line arguments (do not include brackets):"); puts(" /b - turn batch processing on (only channel 0 available for DZT files)."); puts(" /d - turn debugging on."); puts(" Examples: gpr_stak file1.dzt newfile1.dzt -1 -1 2"); puts(" gpr_stak file1.dt1 newfile1.dt1 53 870 10 /b"); printf(" Version %s - Last Revision Date: %s at %s\n",GPR_STAK_VER,DATETIME,__TIME__); puts("###############################################################################"); if (ANSI_there()) printf("%c%c37;40m",0x1B,0x5B); /* white on black */ exit(1); } /* Get command line arguments */ strcpy(in_filename,argv[1]); strcpy(out_filename,argv[2]); start = atoi(argv[3]); end = atoi(argv[4]); stack = atoi(argv[5]); for (itemp = min_args; itemp 0) { if (ANSI_there()) printf("%c%c31;40m",0x1B,0x5B); /* red on black */ printf("%s\n",GetGprFileTypeMsg[itemp]); if (ANSI_there()) printf("%c%c37;40m",0x1B,0x5B); /* white on black */ exit(1); } /* Get the data, stack, and save to disk */ printf("Stacking file %s to %s ...",in_filename,out_filename); switch (file_type) { case DZT: case MOD_DZT: { char text[400],str[160]; int max_proc_hist_size = 400; unsigned char proc_hist[400]; char name[12],antname[14], drive[_MAX_DRIVE], dir[_MAX_DIR]; char fname[_MAX_FNAME], fext[_MAX_EXT]; int i,itemp; int num_hdrs,num_samps,header_size,channel,selection,dot,icount; unsigned short rh_nchan,rg_breaks,rh_nproc; long num_traces,trace,offset,trace_offset,rh_dtype; long disk_chunk_bytes,disk_trace_bytes,block_size; float rg_values[8]; /* this limit may change in future versions */ void *v_ptr; double *d_ptr; FILE *infile,*outfile; struct DztDateStruct create_date; struct DztHdrStruct hdr; struct tm *newtime; time_t long_time; #if defined(_INTELC32_) char *inbuffer = NULL; char *outbuffer = NULL; size_t vbufsize = 64512; /* 63 KB */ #endif struct dosdate_t dosdate; struct dostime_t dostime; char *month[] = { "","January","February","March","April","May","June", "July","August","September","October", "November","December" } ; /* following in gpr_io.h */ extern const char *ReadOneDztHeaderMsg[]; /* Get current time and date */ _dos_getdate(&dosdate); _dos_gettime(&dostime); /* Get first header from the input file */ channel = 0; itemp = ReadOneDztHeader(in_filename,&num_hdrs,&num_traces,channel, &header_size,&hdr); if (itemp > 0) { if (ANSI_there()) printf("%c%c31;40m",0x1B,0x5B); /* red on black */ printf("%s\n",ReadOneDztHeaderMsg[itemp]); if (ANSI_there()) printf("%c%c37;40m",0x1B,0x5B); /* white on black */ if (strstr(strupr((char *)(ReadOneDztHeaderMsg[itemp])),"WARNING") == NULL) exit(3); } /* If more than one header ask user to select */ if (num_hdrs > 1 && !Batch) { puts("\n\nSelect SIR-10A header to use."); do { printf("\r \r"); printf("Enter a number from 1 to %d --->",num_hdrs); selection = getch(); selection -= '0'; } while ((selection < 1) || (selection > num_hdrs)); puts(""); } else selection = 1; /* Get header again if selection not first one */ channel = selection - 1; if (channel > 0) { itemp = ReadOneDztHeader(in_filename,&num_hdrs,&num_traces,channel, &header_size,&hdr); if (itemp > 0) { if (ANSI_there()) printf("%c%c31;40m",0x1B,0x5B); /* red on black */ printf("%s\n",ReadOneDztHeaderMsg[itemp]); if (ANSI_there()) printf("%c%c37;40m",0x1B,0x5B); /* white on black */ if (strstr(strupr((char *)(ReadOneDztHeaderMsg[itemp])),"WARNING") == NULL) exit(3); } } /* Check limits of arguments */ if (start <= 0) start = 0; if (end <= 0 || end <= start || end >= num_traces) end = num_traces - 1; if ((end - start) < 1) { printf("GPR_STAK - ERROR: (end - start) is less than 1.\n"); exit(4); } if (stack < 1) stack = 1; if (stack >= (end-start)) stack = end - start - 1; /* Allocate space for a block of traces (include stacked ones) */ num_samps = hdr.rh_nsamp; disk_trace_bytes = (hdr.rh_bits/8) * num_samps; disk_chunk_bytes = disk_trace_bytes * hdr.rh_nchan; block_size = disk_chunk_bytes * stack; if ((v_ptr = (void *)malloc(block_size)) == NULL) { puts("GPR_STAK - ERROR: unable to allocate storage for one chunk of DZT traces."); exit(4); } /* Allocate space for summing trace */ d_ptr = (double *)malloc(num_samps * sizeof(double)); if (d_ptr == NULL) { puts("GPR_STAK - ERROR: unable to allocate storage for summing trace."); exit(4); } /* Open input file */ if ((infile = fopen(in_filename,"rb")) == NULL) { printf("GPR_STAK - ERROR: Unable to open input file %s\n",in_filename); if (access(in_filename,4)) { if (errno == ENOENT) puts("File or path name not found."); else if (errno == EACCES) puts("Read Access denied."); else puts("Unknown access error."); } free(v_ptr); exit(5); } /* Speed up disk access by using large (conventional memory) buffer */ #if defined(_INTELC32_) realmalloc(inbuffer,vbufsize); if (inbuffer) setvbuf(infile,inbuffer,_IOFBF,vbufsize); else setvbuf(infile,NULL,_IOFBF,vbufsize); #endif /* Open output file */ if ((outfile = fopen(out_filename,"wb")) == NULL) { printf("GPR_STAK - ERROR: Unable to open output file %s\n",out_filename); if (access(out_filename,4)) { if (errno == ENOENT) puts("File or path name not found."); else if (errno == EACCES) puts("Read Access denied."); else puts("Unknown access error."); } fclose(infile); free(v_ptr); exit(5); } /* Speed up disk access by using large (conventional memory) buffer */ #if defined(_INTELC32_) realmalloc(outbuffer,vbufsize); if (outbuffer) setvbuf(outfile,outbuffer,_IOFBF,vbufsize); else setvbuf(outfile,NULL,_IOFBF,vbufsize); #endif /* Set/copy values for output header */ if (hdr.rh_dtype == 0 && hdr.rh_bits != 8) { if (hdr.rh_bits == 16) rh_dtype = 0x0002; else if (hdr.rh_bits == 32) rh_dtype = 0x0004; else rh_dtype = 0x0000; } rg_breaks = *(unsigned short *) ((char *)&hdr + hdr.rh_rgain); for (i=0; i\n", month[dosdate.month],dosdate.day,dosdate.year, dostime.hour,dostime.minute); strcat((char *)proc_hist,str); if (strlen(text) < sizeof(text)-strlen((char *)proc_hist)-2) strcat(text,(char *)proc_hist); rh_nproc = hdr.rh_nproc; if (rh_nproc > max_proc_hist_size-1) rh_nproc = max_proc_hist_size-1; for (i=0; i= stack) /* far enough from end? */ { fread(v_ptr,block_size,1,infile); } else /* get end trace */ { offset = (num_hdrs * header_size) + (end * disk_chunk_bytes) + trace_offset ; fseek(infile,offset,SEEK_SET); fread(v_ptr,disk_trace_bytes,1,infile); stack = 1; } /* check for markers and stack */ switch (rh_dtype) { case 0x0000: { unsigned char *uc_ptr; int samp; /* initialize stacking trace */ for (i=0; i 0) { if (ANSI_there()) printf("%c%c31;40m",0x1B,0x5B); /* red on black */ printf("\n%s\n",GetSsHdFileMsg[itemp]); if (ANSI_there()) printf("%c%c37;40m",0x1B,0x5B); /* white on black */ exit(3); } if (Debug) PrintSsHdInfo("stdout",hd_outfilename,&HdInfo); /* Check limits of arguments */ if (start <= 0) start = 0; if (end <= 0 || end <= start || end >= HdInfo.num_traces) end = HdInfo.num_traces - 1; if ((end - start) < 1) { printf("GPR_STAK - ERROR: (end - start) is less than 1.\n"); exit(4); } if (stack < 1) stack = 1; if (stack >= (end-start)) stack = end - start - 1; /* Allocate space for a block of traces (include skipped ones) */ num_samps = HdInfo.num_samples; data_block_bytes = num_samps * 2; /* data are signed short ints */ hdr_block_bytes = sizeof(struct SsTraceHdrStruct); trace_size = data_block_bytes + hdr_block_bytes; block_size = trace_size * stack; v_ptr = (void *)malloc(block_size); if (v_ptr == NULL) { printf("GPR_STAK - ERROR: Unable to allocate storage for DT1 data block.\n"); exit(4); } /* Allocate space for summing trace */ d_ptr = (double *)malloc(num_samps * sizeof(double)); if (d_ptr == NULL) { puts("GPR_STAK - ERROR: unable to allocate storage for summing trace."); exit(4); } /* Open input file */ if ((infile = fopen(in_filename,"rb")) == NULL) { printf("GPR_STAK - ERROR: Unable to open input file %s\n",in_filename); if (access(in_filename,4)) { if (errno == ENOENT) puts("File or path name not found."); else if (errno == EACCES) puts("Read Access denied."); else puts("Unknown access error."); } free(v_ptr); exit(5); } /* Speed up disk access by using large (conventional memory) buffer */ #if defined(_INTELC32_) realmalloc(inbuffer,vbufsize); if (inbuffer) setvbuf(infile,inbuffer,_IOFBF,vbufsize); else setvbuf(infile,NULL,_IOFBF,vbufsize); #endif /* Open output file */ if ((outfile = fopen(out_filename,"wb")) == NULL) { printf("GPR_STAK - ERROR: Unable to open output file %s\n",out_filename); if (access(out_filename,4)) { if (errno == ENOENT) puts("File or path name not found."); else if (errno == EACCES) puts("Read Access denied."); else puts("Unknown access error."); } fclose(infile); free(v_ptr); exit(5); } /* Speed up disk access by using large (conventional memory) buffer */ #if defined(_INTELC32_) realmalloc(outbuffer,vbufsize); if (outbuffer) setvbuf(outfile,outbuffer,_IOFBF,vbufsize); else setvbuf(outfile,NULL,_IOFBF,vbufsize); #endif /* Set/save values for output header */ start_pos = HdInfo.start_pos + (HdInfo.step_size * start); final_pos = HdInfo.start_pos + (HdInfo.step_size * end); step_size = HdInfo.step_size * stack; num_traces = 2 + ((end - 1 - start)/stack); traces_per_sec = 0.0; meters_per_mark = 0.0; num_gain_pts = 0; gain_pts = NULL; itemp = 0; if (HdInfo.proc_hist) /* if allocated */ itemp = strlen(HdInfo.proc_hist) + 1; itemp += proc_hist_size + 1; proc_hist = (char *)malloc((size_t)itemp); if (proc_hist) { if (HdInfo.proc_hist) /* allocated */ strcpy(proc_hist,HdInfo.proc_hist); sprintf((char *)proc_hist,"\nGPR_STAK v. 1.00 (USGS-JL) on %s %d, %d at %d:%02d ->\n", month[dosdate.month],dosdate.day,dosdate.year, dostime.hour,dostime.minute); strcat(proc_hist,str); sprintf(str," file = %s start = %d end = %d stack = %d\n", in_filename,start,end,stack); strcat(proc_hist,str); } itemp = SaveSsHdFile((int)HdInfo.day,(int)HdInfo.month,(int)HdInfo.year, num_traces,(long)HdInfo.num_samples, (long)HdInfo.time_zero_sample,(int)HdInfo.total_time_ns, start_pos,final_pos,step_size,HdInfo.pos_units, (double)HdInfo.ant_freq,(double)HdInfo.ant_sep, (double)HdInfo.pulser_voltage, (int)HdInfo.num_stacks,HdInfo.survey_mode, proc_hist,hd_outfilename,out_filename, HdInfo.job_number,HdInfo.title1,HdInfo.title2, traces_per_sec,meters_per_mark, num_gain_pts,gain_pts); free(proc_hist); /* Move input file pointer to first trace chunk to retrieve */ offset = start * trace_size; fseek(infile,offset,SEEK_SET); /* Condense file. File pointer must be at start of first trace to get. Read in desired trace plus stacked ones. */ dot = (end-start)/stack/10; if (dot < 1) dot = 1; icount = 0; trace = start; for (trace=start; trace<=end; trace+=stack,icount++) { if (!(icount % dot)) printf("."); if (end - trace >= stack) /* far enough from end? */ { fread(v_ptr,block_size,1,infile); } else /* get end trace */ { offset = end * trace_size; fseek(infile,offset,SEEK_SET); fread(v_ptr,trace_size,1,infile); stack = 1; } /* update trace header */ TraceHdr = (struct SsTraceHdrStruct *)v_ptr; TraceHdr->aux_trace_num = TraceHdr->trace_num = icount; /* initialize stacking trace */ for (i=0; i 0) { if (ANSI_there()) printf("%c%c31;40m",0x1B,0x5B); /* red on black */ printf("\n%s\n",ReadSegyReelHdrMsg[itemp]); if (ANSI_there()) printf("%c%c37;40m",0x1B,0x5B); /* white on black */ exit(3); } if (Debug) PrintSegyReelHdr(swap_bytes,num_traces,"stdout", in_filename,&ReelHdr); /* Check limits of arguments */ if (start <= 0) start = 0; if (end <= 0 || end <= start || end >= num_traces) end = num_traces - 1; if ((end - start) < 1) { printf("GPR_STAK - ERROR: (end - start) is less than 1.\n"); exit(4); } if (stack < 1) stack = 1; if (stack >= (end-start)) stack = end - start - 1; /* Allocate space for a block of traces (include skipped ones) */ num_samps = ReelHdr.hns; switch (ReelHdr.format) { case 1: case 2: case 4: data_block_bytes = num_samps * 4; break; case 3: data_block_bytes = num_samps * 2; break; default: printf("GPR_STAK - ERROR: Unvalid SEG-Y sample format code.\n"); exit(4); } hdr_block_bytes = sizeof(struct SegyTraceHdrStruct); trace_size = data_block_bytes + hdr_block_bytes; block_size = trace_size * stack; v_ptr = (void *)malloc(block_size); if (v_ptr == NULL) { printf("GPR_STAK - ERROR: Unable to allocate storage for SEG-Y data block.\n"); exit(4); } /* Allocate space for summing trace */ d_ptr = (double *)malloc(num_samps * sizeof(double)); if (d_ptr == NULL) { puts("GPR_STAK - ERROR: unable to allocate storage for summing trace."); exit(4); } /* Open input file */ if ((infile = fopen(in_filename,"rb")) == NULL) { printf("GPR_STAK - ERROR: Unable to open input file %s\n",in_filename); if (access(in_filename,4)) { if (errno == ENOENT) puts("File or path name not found."); else if (errno == EACCES) puts("Read Access denied."); else puts("Unknown access error."); } free(v_ptr); exit(5); } /* Speed up disk access by using large (conventional memory) buffer */ #if defined(_INTELC32_) realmalloc(inbuffer,vbufsize); if (inbuffer) setvbuf(infile,inbuffer,_IOFBF,vbufsize); else setvbuf(infile,NULL,_IOFBF,vbufsize); #endif /* Open output file */ if ((outfile = fopen(out_filename,"wb")) == NULL) { printf("GPR_STAK - ERROR: Unable to open output file %s\n",out_filename); if (access(out_filename,4)) { if (errno == ENOENT) puts("File or path name not found."); else if (errno == EACCES) puts("Read Access denied."); else puts("Unknown access error."); } fclose(infile); free(v_ptr); exit(5); } /* Speed up disk access by using large (conventional memory) buffer */ #if defined(_INTELC32_) realmalloc(outbuffer,vbufsize); if (outbuffer) setvbuf(outfile,outbuffer,_IOFBF,vbufsize); else setvbuf(outfile,NULL,_IOFBF,vbufsize); #endif /* Set/save values for output reel header */ ExtractSsInfoFromSegy(&GprInfo,&ReelHdr); start_pos = GprInfo.starting_position + (GprInfo.position_step_size * start); final_pos = GprInfo.starting_position + (GprInfo.position_step_size * end); step_size = GprInfo.position_step_size * stack; num_traces = 2 + ((end - 1 - start)/stack); traces_per_sec = GprInfo.traces_per_sec; meters_per_mark = GprInfo.meters_per_mark; num_gain_pts = GprInfo.num_gain_pts; gain_pts = GprInfo.gain_pts; itemp = 0; if (GprInfo.proc_hist && GprInfo.proc_hist_bytes > 0) /* if allocated */ itemp = GprInfo.proc_hist_bytes; itemp += proc_hist_size + 1; proc_hist = (char *)malloc((size_t)itemp); if (proc_hist) { if (GprInfo.proc_hist && GprInfo.proc_hist_bytes > 0) strncpy(proc_hist,GprInfo.proc_hist,GprInfo.proc_hist_bytes); sprintf((char *)proc_hist,"\nGPR_STAK v. 1.00 (USGS-JL) on %s %d, %d at %d:%02d ->\n", month[dosdate.month],dosdate.day,dosdate.year, dostime.hour,dostime.minute); strcat(proc_hist,str); sprintf(str," file=%s start=%d end=%d stack=%d\n", in_filename,start,end,stack); strcat(proc_hist,str); } SetSgyFileHeader((int)GprInfo.time_date_created.day, (int)GprInfo.time_date_created.month, (int)GprInfo.time_date_created.year+1980, num_traces,num_samps,(long)GprInfo.timezero_sample, GprInfo.in_time_window, start_pos,final_pos,step_size,GprInfo.position_units, (double)GprInfo.nominal_frequency, GprInfo.antenna_separation,GprInfo.pulser_voltage, GprInfo.number_of_stacks,GprInfo.survey_mode, proc_hist,GprInfo.text,GprInfo.job_number, GprInfo.title1,GprInfo.title2,out_filename, traces_per_sec,meters_per_mark, num_gain_pts,gain_pts,&ReelHdr); fwrite((void *)&ReelHdr,sizeof(struct SegyReelHdrStruct),1,outfile); free(proc_hist); /* Move input file pointer to first trace chunk to retrieve */ offset = sizeof(struct SegyReelHdrStruct) + (start * trace_size); fseek(infile,offset,SEEK_SET); /* Condense file. File pointer must be at start of first trace to get. Read in desired trace plus stacked ones. */ dot = (end-start)/stack/10; if (dot < 1) dot = 1; icount = 0; trace = start; for (trace=start; trace<=end; trace+=stack,icount++) { if (!(icount % dot)) printf("."); if (end - trace >= stack) /* far enough from end? */ { fread(v_ptr,block_size,1,infile); } else /* get end trace */ { offset = sizeof(struct SegyReelHdrStruct) + (end * trace_size); fseek(infile,offset,SEEK_SET); fread(v_ptr,trace_size,1,infile); stack = 1; } /* update trace header */ TraceHdr = (struct SegyTraceHdrStruct *)v_ptr; TraceHdr->tracl = TraceHdr->tracr = TraceHdr->fldr = icount; /* initialize stacking trace */ for (i=0; i