seg2segy.c

完整的seg2 segy 转换程序

/* Copyright (c) Colorado School of Mines, 2008.*/
/* All rights reserved.                       */

/* Copyright (c) Colorado School of Mines, 2002.*/
/* All rights reserved.                       */

/*
	FROM THE PUBLISHED VERSION IN GEOPHYSICS, SEPT '90
	Original Header:
	NAME:            SEG2SEGY.C
	VERSION:         V1.0
	CREATION DATE:   6/15/90
	AUTHOR:          Brett Bennett
	PURPOSE:         Conversion of the SEG2-PC standard data set 
		 to SEG-Y 32-BIT format.
	LIMITATIONS:     In this implementation of the conversion only
	the data generated in the IBM byte ordering method can be
	converted.  i.e.  byte ordering must be LSB-MSB Data generated
	by a 68000 pro- cessor, for instance, will not convert.  A flag
	is however in the code for future addition of this capability.
	The SEG-Y file also differs slightly from Barry et al standard.
	Word 88 is defined as "last-trace-flag" rather than "Geophone
	group number of last trace within original file."  This
	simplification greatly reduces the amount of code needed for
	conversion.  Header word 92 is defined here as the
	Shot-Sequence-Number.  This is used to assign each trace group a
	unique number.

	EXTENSIVE REWRITE/HACK AND PORT: 11/95 
	AUTHOR: Ian Kay.  
	-Added code to handle running on big-endian machines.
	-Replaced all non-unix portable code and library calls.
	-Typed all the ints so they are the right types
	    (on DOS int=short, on Unix int=long)
	-removed calls to "ieeeibm" since workstations are ieee floating point
	    format as well.
	-partial addition of structures to handle the keywords and headers.
	but I eventually gave up.

	$Id: seg2segy.c,v 1.2 2003/02/06 19:05:45 pm Exp $

	Changes:
	Feb 96.  Trace header was not being written out with swap on bytes above
	89 although info is stored in 92,93,94.  Error in original, fixed.

	Oct 97.  Changed the input to command line, inspired by the seg2segy
	version in CWP/SU package.
	
	Jan 98.  
	fixed malloc problem; not enough space was malloced to hold the 
	string defpath. 
	freed the malloced space when finished. 
	closed the segykeyword file when finished.
	changed curf and lastf to int from short; this should allow digits
	in the file name prefix (but what will happen in DOS?).

	Sep 98.
	Ok.  So no one can figure out the segykeyword file mechanism, 
	so I've changed it.  Now, if the environment variable SEGKEYPATH
	is set , then that file is used for the segykeywords. If that
	isn't set then a system default is used /usr/local/lib/segykeyword
	or "segykeyw.ord" in the local directory in DOS.  Finally, neither
	of those exist, then the default values are built in.  

	NOTE that this means that you can do away with the segykeyword
	file altogether!

	$Log: seg2segy.c,v $
	Revision 1.2  2003/02/06 19:05:45  pm
	fixed sawtooth on 20bit packed data conversion case 3
	
	Revision 1.1  2003/02/06 18:44:59  pm
	working on 3rdparty seismic unix
	
	Revision 1.2  2002/11/12 16:01:43  john
	type changed on i,j,k to in
	t

	Revision 1.1  2001/04/12 17:55:11  john
	Initial revision

	Revision 1.7  1998/10/08 20:17:15  kay
	A little more cleanup...^[

	Revision 1.6  1998/10/07 18:54:40  kay
	Built in the defaults for the segykeyword file.  This file is no longer
	needed.

	Revision 1.5  1998/01/07 13:36:15  kay
	added new changes and a Log for the future

 */

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>
#ifndef __MSDOS__   
#include <unistd.h>
#define MAXSAMP 16384
#define MAXTRACES 16384
#else
#include <io.h>
#define MAXSAMP 4096
#define MAXTRACES 512
#endif

#ifndef NULL
#define NULL    ((void *)0)
#endif 

#define STRINGWIDTH 100
#define MAXKEYWORDS 100
#define MAXPARMS 10

/* global data... */
short int segyreelheader[1800];
int totalkeys;
char string1[200];
short int outhead[120];


typedef struct _SegyKeyTable
	{
	char segykeyword[STRINGWIDTH];
	int segyfunction;
	int segyheader;
	double segyparms[MAXPARMS];
	} SegyKeyTable;


void
swapshort (short *si)
{
/*      void swab(char *from, char *to, size_t nbytes); */
	short tempc;
	tempc = *si;
	swab ((char *) &tempc, (char *) si, sizeof (short));
}

void
swaplong (long *li)
{
	short sl, sh;
	long temp, temp1;
	temp = *li;
	sl = temp & 0x0000ffff;
	sh = (temp >> (8 * sizeof (short))) & 0x0000ffff;
	swapshort (&sl);
	swapshort (&sh);
	temp = (sl << (8 * sizeof (short)));
	temp1 = (sh & 0x0000ffff);
	temp = temp | temp1;

	*li = temp;
}

void
float_to_ibm(long from[], long to[], long n)
{
	register long fconv, fmant, ii, t;

	for (ii=0;ii<n;++ii) {
	fconv = from[ii];
	if (fconv) {
	fmant = (0x007fffff & fconv) | 0x00800000;
	t = (long) ((0x7f800000 & fconv) >> 23) - 126;
	while (t & 0x3) { ++t; fmant >>= 1; }
	fconv = (0x80000000 & fconv) | (((t>>2) + 64) << 24) | fmant; 
	} 
	to[ii] = fconv; 
	} 
	return; 
}

char *
strrev (char *s)
{
	char *tmpstr;
	int i, len;
	tmpstr = (char *) malloc (sizeof (char) * strlen (s));
	len = 0;
	while (s[len] != '\0')
	 len++;
	for (i = 0; i < len; i++) {
	tmpstr[i] = s[len - 1 - i];
	  tmpstr[i + 1] = '\0';
	}
	tmpstr[++i] = '\0';
	strcpy (s, tmpstr);
	free (tmpstr);
	tmpstr = (char *) NULL;
	return s;
}


void
parseSegyKeys(FILE *keyfptr, SegyKeyTable *keywords)
{
	int i, j;
	char input[STRINGWIDTH], inputbuf[STRINGWIDTH];
	char *token;
	i = 0;
	while (fgets (input, STRINGWIDTH, keyfptr)) {
	j = 0;
	if ( strlen(input) > (size_t)STRINGWIDTH) {
	fprintf (stderr, "String too long!\n (line:%d)\n", __LINE__);
	exit (-12);
	}
	/*  
	 * if left most character = "*" then 
	 * this line is a comment and should be ignored. 
	 */
	if (input[0] == '*' || input[0] == '\n')  
	continue;

	strcpy (inputbuf, input);  /* make a working copy of input */
	token = strtok (inputbuf, " ");  /* search out space. */

	/* copy keyword into array */
	strncpy (keywords[i].segykeyword, token, 1 + strlen (token));

	/* get function value. */
	token = strtok (NULL, " ");
	keywords[i].segyfunction = atoi (token);  /* convert to value. */
	token = strtok (NULL, " ");

	/* get segy header pointer */
	keywords[i].segyheader = atoi (token);

	/* now start getting any parameters on the line. */
	token = strtok (NULL, " ");
	j = 0;
	while (token != NULL) {
	/* get next token and start putting them in their array location */
	keywords[i].segyparms[j] = atof (token);
	token = strtok (NULL, " ");
	j++;
	if (j > MAXPARMS) {
	printf ("Too many parameters in %s keyword\n", keywords[i].segykeyword);
	printf ("No more than %d allowed per function\n", j - 1);
	exit (-13);
	}
	}      /* end parameter extraction while loop */
	i++;      /* inc counter to keyword number */
	}        /* end keyword string while loop */
	totalkeys = i;
	
}        /* end parseSegyKeys() */

void
loadSegyKeyTableDefaults(SegyKeyTable * keywords)
{
char *defaults[] = {
"ACQUISITION_DATE 3 6 ",
"ACQUISITION_TIME 3 7 ",
"ALIAS_FILTER 5 0 1 71 72 ",
"CDP_NUMBER 5 1 1 12 ",
"CDP_TRACE 5 1 1 14 ",
"CHANNEL_NUMBER 5 1 1 8 ",
"CLIENT 3 1 ", 
"COMPANY 3 2 ",
"DATUM 5 1 1 28 ",
"DELAY 1 55 1000 ",
"END_OF_GROUP 1 88 1 ",
"FIXED_GAIN 1 61 1 ",
"FIXED_GAIN 1 62 1 ",
"HIGH_CUT_FILTER 5 0 1 76 78 ",
"INSTRUMENT 3 4 ",
"JOB_ID 3 9 ",
"LINE_ID 3 3 ",
"LOW_CUT_FILTER 5 0 1 75 77 ",
"NOTCH_FREQUENCY 5 0 1 73 ",
"OBSERVER 3 5 ",
"RAW_RECORD 5 1 1 6 ",
"RECEIVER_LOCATION 5 1 1 42 44 22 ",
"RECEIVER_STATION_NUMBER 1 93 1 ",
"SAMPLE_INTERVAL 1 59 1000000 ",
"SOURCE_LOCATION 5 1 1 38 40 26 ",
"SHOT_SEQUENCE_NUMBER 1 92 1 ",
"SOURCE_STATION_NUMBER 1 94 1 ",
"STACK 1 16 1 ",
"STATIC_CORRECTIONS 5 0 1 50 51 52 ",
"TRACE_SORT 2 1615 ",
"TRACE_TYPE 4 15 ",
"UNITS 3 8 ",
"AMPLITUDE_RECOVERY 0 0 ",
"BAND_REJECT_FILTER 0 0 ",
"DESCALING_FACTOR 0 0 ",
"DIGITAL_BAND_REJECT_FILTER 0 0 ",
"DIGITAL_HIGH_CUT_FILTER 0 0 ",
"DIGITAL_LOW_CUT_FILTER 0 0 ",
"GENERAL_CONSTANT 0 0 ",
"NOTE 0 0  ",
"POLARITY 0 0 ",
"PROCESSING_DATE 0 0 ",
"PROCESSING_TIME 0 0 ",
"RECEIVER 0 0  ",
"RECEIVER_GEOMETRY 0 0 ",
"RECEIVER_SPECS 0 0 ",
"SKEW 0 0  ",
"SOURCE 0 0 ",
"SOURCE_GEOMETRY 0 0 ",
NULL};
	char **defaultsptr=defaults;

	char tmpfilename[L_tmpnam];
	FILE *tmpfileptr;

	tmpnam(tmpfilename);
	tmpfileptr=fopen(tmpfilename, "w");
	/*fwrite(defaults, sizeof(char), strlen(defaults), tmpfileptr); */

	while(*defaultsptr) fprintf(tmpfileptr,"%s\n", *defaultsptr++);

	tmpfileptr=freopen(tmpfilename,"r",tmpfileptr);

	parseSegyKeys(tmpfileptr, keywords);
	fclose(tmpfileptr);
	remove(tmpfilename);
	
}        /* end loadSegyKeyTableDefaults() */


/*-----------------------------------------------------------------------
	READSEGYKEYS
	This routine reads in the valid keywords in file SEGKEYW.ORD 
	and stores the results in global arrays
----------------------------------------------------------------------*/
void
readSegyKeys (SegyKeyTable *keywords)
{
	char *envkeypath;
	char *syskeypath;
#ifdef __MSDOS__
	char *defpath = "segykeyw.ord";
#else
	char *defpath = "/usr/local/lib/segykeyword";
#endif

	FILE *keyfile;

	syskeypath = (char *) malloc ((strlen(defpath)+1) * sizeof(char));
	sprintf (syskeypath, "%s", defpath);
	envkeypath = getenv ("SEGKEYPATH");

	if (envkeypath != (char *) NULL)
	{
	if((keyfile=fopen(envkeypath,"rb")) != (FILE *)NULL)
	fprintf (stderr, "%s used for header word mapping.\n", envkeypath); 
	parseSegyKeys(keyfile, keywords);
	fclose (keyfile);
	} 
	else if ((keyfile = fopen (syskeypath, "rb")) != (FILE *)NULL)
	{
	fprintf (stderr, "Using header word mappings in %s\n",defpath);
	parseSegyKeys(keyfile, keywords);
	fclose (keyfile);
	} 
	else 
	{ 
	fprintf (stderr, "Using default header word mappings. \n");
	loadSegyKeyTableDefaults(keywords);
	}
	free(syskeypath);
	return;
}        /* end readsegyKeys() */



/*------------------------------------------------------------
 * KEYCHECK
 * This routine compares list of keywords with the global string 1 and inserts
 * values found in the input string into specified header locations (see
 * segykeyw.ord).  KEYCHECK checks EVERY keyword for a match.  Repeated 
 * keywords are valid.
 * GLOBALS USED:
 * totalkeys int
 *------------------------------------------------------------*/
void
keycheck (SegyKeyTable *keywords)
{
	int i;
	int matchfound;
	char string2[STRINGWIDTH];
	char *token;

	/* start a loop for total keys and compare the input with the list */
	/* make a copy of string1. strtok destroys it's input during token search*/
	/* therefore need to keep a copy for each keyword checked */
	strcpy (string2, string1);

	/* clear the match found flag */
	matchfound = 0;
	for (i = 0; i < totalkeys; i++) 
	{
	/* 
	 * restore string1.  It will have been destroyed if an 
	 * earlier match has occured. 
	 */
	strcpy (string1, string2);  

	if (0 == strncmp (string1, keywords[i].segykeyword, 
		strlen (keywords[i].segykeyword) ) )
	{
	/* have found a match! Set the matchfound flag */
	matchfound = 1;
	/* look up the function and implement it. */
	switch (keywords[i].segyfunction)
	{
	case 0:
		break;    /* null case nothing to do */
	case 1:
	{
	/* function 1 keywords have a single parameter 
	 * in the input;  assumed to be a number.  */
	/* find the parameter. */
	/* find first token which will be a keyword */
	token = strtok (string1, " ");  
	/* should be a number. */
	token = strtok (NULL, " ");  	
	/* parameter found. pointed to by token. */
	/* function 1 calls for the value on the input line to be mulitplied */
	/* by segykeyword parm 1. and then inserted into header. */
	outhead[keywords[i].segyheader - 1] = 
		atof (token) * keywords[i].segyparms[0];
	break;
	}
	case 2:
	{
	/* function 2 is a special function.  It has to deal with the special 
	 * case of trace sorting. In this case the parameter on the input 
	 * line is not a number  but a char string.  
	 * Notice that spaces in the keywords (i.e. AS ACQUIRED) cause 
	 * the parsing to fail.  Words must be 'spaceless' */
	token = strtok (string1, " ");  	/* pointing to keyword. */
	token = strtok (NULL, " ");  	/* token points to input char. */
	if (0 == strcmp ("AS_ACQUIRED", token))
		segyreelheader[keywords[i].segyheader - 1] = 1;
	else if (0 == strcmp ("CDP_GATHER", token))
		segyreelheader[keywords[i].segyheader - 1] = 2;
	else if (0 == strcmp ("CDP_STACK", token))
		segyreelheader[keywords[i].segyheader - 1] = 4;
	else if (0 == strcmp ("COMMON_OFFSET", token))
		segyreelheader[keywords[i].segyheader - 1] = 3;
	else if (0 == strcmp ("COMMON_RECEIVER", token))
		segyreelheader[keywords[i].segyheader - 1] = 1;
	else if (0 == strcmp ("COMMON_SOURCE", token))
		segyreelheader[keywords[i].segyheader - 1] = 1;
	else if (0 == strcmp ("METERS",token))
		segyreelheader[keywords[i].segyheader -1] = 1;
	else if (0 == strcmp ("FEET",token))
		segyreelheader[keywords[i].segyheader -1] = 2;

	break;
	}      /* end case 2 */
	case 3:
	{
	/* 
	 * this case requires the text string found on the input
	 * line to be copied to the SEGY-REEL header at the line
	 * indicated in the segyheader index. to compute this
	 * location it is (line#-1)*80
	 */
	strncpy ((char *) &segyreelheader[80 * (keywords[i].segyheader - 1)], 
		string1, 80);
	break;
	}      /* end case 3 */
	case 4:
	{
	/* this case, like 2, calls for special string parsing. */
	/* for TRACE_TYPE, see code for allowable inputs. */
	/* pointing to keyword. */
	token = strtok (string1, " "); 
	/* assume its seismic */