seg2segy.c

完整的seg2 segy 转换程序

	outhead[keywords[i].segyheader - 1] = 1;  
	/* token points to input char. */
	token = strtok (NULL, " "); 
	if (0 == strcmp ("SEISMIC_DATA", token))
		outhead[keywords[i].segyheader - 1] = 1;
	else if (0 == strcmp ("DEAD", token))
		outhead[keywords[i].segyheader - 1] = 2;
	else if (0 == strcmp ("TEST_DATA", token))
		outhead[keywords[i].segyheader - 1] = 3;
	else if (0 == strcmp ("UPHOLE", token))
		outhead[keywords[i].segyheader - 1] = 5;
	/* RADAR_DATA not defined in SEG-Y assume it to be normal seismic */
	else if (0 == strcmp ("RADAR_DATA", token))
		outhead[keywords[i].segyheader - 1] = 1;
	break;
	}      /* end case 4 */
	case 5:
	{
	/* this case calls for the input data to be inserted into the 
	 * sepcified header location of parms 2-10. The normal 
	 * 'header' location indicates the type of data, 0=int, 
	 * 1=long int, 2=floating point. Parm 1 is the multiplier. */
	token = strtok (string1, " ");
	/* do the integer case first */
	if (keywords[i].segyheader == 0) {
		int paramcount = 1;
		int headindex;
		token = strtok (NULL, " ");    	/* token points to input char. */
		while (token != NULL && paramcount < 10) {
		headindex = keywords[i].segyparms[paramcount] - 1;
		outhead[headindex] = atof (token) * keywords[i].segyparms[0];
		paramcount++;
		token = strtok (NULL, " ");  /* token points next input char. */
		}
	}
		/* do the long integer case next */
	else if (keywords[i].segyheader == 1) {
		int paramcount = 1;
		int headindex;
		long *outpoint;
		token = strtok (NULL, " ");    /* token points to input char. */
		while (token != NULL && paramcount < 10) {
		/* set outpoint to beginning of WORD of interest */
		/* need to subtract 2 from the location value. */
		headindex = keywords[i].segyparms[paramcount] - 2;
		outpoint = (long *) &outhead[headindex];
		/* outpoint now points to the area of interest and is 
		 * typed correctly. Compute the value, cast it and send it in. */
		outpoint[0] = (long) (atof (token) * keywords[i].segyparms[0]);
		paramcount++;
		token = strtok (NULL, " ");  /* token points next input char. */
		}
	}
	/* finally do floating point case */
	else if (keywords[i].segyheader == 2) {
		int paramcount = 1;
		int headindex;
		float *outpoint;
		token = strtok (NULL, " ");    /* token points to input char. */
		while ((token != NULL) && (paramcount < 10))
		{
		/* set outpoint to point to beginning of WORD of interest */
		/* need to subtract 2 from the location value. */
		headindex = keywords[i].segyparms[paramcount] - 2;
		outpoint = (float *) &outhead[headindex];
		/* outpoint now points to the area of interest and is typed 
		 * correctly. Compute the value, cast it and put it in 
		 * outhead (using outpoint). */
		outpoint[0] = (float) (atof (token) * keywords[i].segyparms[0]);
		paramcount++;
		/*   ieee2ibm(outpoint,0); *//* convert to ibm format if necessary */
		token = strtok (NULL, " ");  /* token points next input char. */
		}
	}
	break;
	}      /* end case 5 */
	case 6: 
	{
	short day, year;
	token=strtok(string1," "); 
	token=strtok(NULL,"/"); 
	day=atoi(token);
	token=strtok(NULL,"/");
	if (0==strcmp("FEB",token) || 0==strcmp("02",token)) day+=31;
	else if(0 == strcmp("MAR",token) || 0== strcmp("03",token)) day+=59;
	else if(0 == strcmp("APR",token) || 0== strcmp("04",token)) day+=90;
	else if(0 == strcmp("MAY",token) || 0== strcmp("05",token)) day+=120;
	else if(0 == strcmp("JUN",token) || 0== strcmp("06",token)) day+=151;
	else if(0 == strcmp("JUL",token) || 0== strcmp("07",token)) day+=181;
	else if(0 == strcmp("AUG",token) || 0== strcmp("08",token)) day+=212;
	else if(0 == strcmp("SEP",token) || 0== strcmp("09",token)) day+=243;
	else if(0 == strcmp("OCT",token) || 0== strcmp("10",token)) day+=273;
	else if(0 == strcmp("NOV",token) || 0== strcmp("11",token)) day+=304;
	else if(0 == strcmp("DEC",token) || 0== strcmp("12",token)) day+=334;
		token=strtok(NULL," ");/* token points to input char. */
	year=atoi(token);
	if(!year%4 && day>59) day+=1;  /* Yikes.  This may break! */
	outhead[keywords[i].segyheader - 2] = year;
	outhead[keywords[i].segyheader - 1] = day;
	break; 
	} 
	default:    /* case where function not found.. should never happen */
	{
	printf ("Function %d not defined.\n", keywords[i].segyfunction);
	break;
	}
	}            /* end switch */
	break;        /* don't go through rest of for() loop, go to next string */
	}              /* end if */
		      	/* loop through and see if it can be found again */
	}                /* end of i loop */
	if (!matchfound)      /* did a match occur */
	 	printf ("No match found for %s\n", string1);
}             /* end of keysegy */
 



int
main (int argc, char *argv[])
{

#define NFILECHAR 32
	char  prefix[NFILECHAR], seg2file[NFILECHAR], 
	segyfile[NFILECHAR], suffix[NFILECHAR], str[NFILECHAR];
	char *digits = "1234567890";
	char stringtermcount, stringterm1, stringterm2;
	char linetermcount, lineterm1, lineterm2;
	char reserved[19];
	unsigned char datatype;
	int i, j, k; 
	short int ssn; 
	short int reversed = 0;
	int curf, lastf;
	short int blockid, revnum, pointerbytecount, numtrace, stringlength;
	short int first = 1;
	size_t l,ln;
	unsigned short int blockleng;
	long tracepointers[MAXTRACES], outbuf[MAXSAMP];
	long *outheadl;
	unsigned long numsamples, datalength;
	FILE *f1, *f2;
	SegyKeyTable keywords[MAXKEYWORDS];

	double dinbuf[MAXSAMP];
	float *finbuf = (float *)dinbuf;
	long *linbuf = (long *)dinbuf;
	short int *iinbuf = (short int *)dinbuf;
	unsigned char *cinbuf = (unsigned char *)dinbuf;

	outheadl = (long *) &outhead[0];
	for (i = 0; i < 1800; i++) segyreelheader[i] = 0;

	if (argc < 3 || argc > 4) { 
	printf("Usage: seg2segy first-seg2file number-of-files [shot-number]\n");
	exit(-1);
	}
	sprintf(seg2file,"%s", argv[1]);
	if (strchr(seg2file,'.')==NULL)
	strcat(seg2file, ".dat");
	l = strcspn(seg2file,".");
	strncpy(segyfile, seg2file, l); 
	segyfile[l]='\0';
	strcpy(suffix,seg2file+l);
	l=strcspn(segyfile,digits);
	if (l==strlen(segyfile) 
	|| strspn(segyfile+l,digits)!=strlen(segyfile+l)){ 
	printf("file name seg2 %s invalid\n", seg2file);
	exit(-2);
	}
	strncpy(prefix, segyfile, l); prefix[l]='\0';
	curf=atoi(segyfile+l);
	ln=strlen(segyfile+l);
	strcat(segyfile, ".sgy");
	lastf=curf+atoi(argv[2])-1;
	if (argc==4) ssn=atoi(argv[3]); else ssn=1;

	f2 = fopen(segyfile, "wb");
	if (f2 == (FILE *) NULL) {
	fprintf (stderr, "**OUTPUT FILE OPEN FAILURE**\n **ABORTING**\n");
	  exit (-3);
	}

	readSegyKeys (keywords);

	/* start the big loop ... */
	for(; curf<=lastf; curf++) {
	strcpy(seg2file,prefix);
	sprintf(str,"%d",curf);
	l=strlen(str);
	while(l<ln) { strcat(seg2file,"0"); l++; }
	strcat(seg2file,str);
	strcat(seg2file,suffix);
	if ((f1=fopen(seg2file,"rb")) == (FILE *)NULL) {
	fprintf (stderr, "\n***ERROR OPENING FILE %s***\n", seg2file);
	fprintf (stderr, "Skipping to next file number.\n");
	continue;    /* go to end of loop, try next file */
	}
	fread (&blockid, 2, 1, f1);
	if (blockid == 0x553a) 
	reversed = 1;
   
	if (blockid != 0x3a55) {
	if (!reversed) {
	fprintf (stderr, "Not SEG-2 data can not continue\n");
	exit (-4);
	}
	}

	fread (&revnum, 2, 1, f1);
	fread (&pointerbytecount, 2, 1, f1);
	fread (&numtrace, 2, 1, f1);
	if (reversed) {
	swapshort (&revnum);
	swapshort (&pointerbytecount);
	swapshort (&numtrace);
	}

	printf ("File %s, Data Format Revision: %d, Number of traces: %d\n", seg2file, revnum, numtrace);
	fread (&stringtermcount, 1, 1, f1);
	fread (&stringterm1, 1, 1, f1);
	fread (&stringterm2, 1, 1, f1);
	fread (&linetermcount, 1, 1, f1);
	fread (&lineterm1, 1, 1, f1);
	fread (&lineterm2, 1, 1, f1);
	fread (reserved, 1, 18, f1);  /* reserved block, not used */

	if (numtrace > (pointerbytecount / 4)) {
	fprintf (stderr, "Number of traces greater than number of trace pointers\n");
	fprintf (stderr, "Number of pointers = %d\n", pointerbytecount / 4);
	fprintf (stderr, "Due to this inconsistency processing must stop\n");
	exit (-5);
	}
	fread (tracepointers, 4, numtrace, f1);
	if (reversed) {
	for (i = 0; i < numtrace; i++)
	swaplong (&tracepointers[i]);
	}

	/* now read file descriptor block.  */
	fread (&stringlength, 2, 1, f1);
	if (reversed)
	swapshort (&stringlength);

	while (0 != stringlength) {
	  fread (string1, 1, stringlength - 2, f1);
	  keycheck (keywords);
	  fread (&stringlength, 2, 1, f1);
	  if (reversed) swapshort (&stringlength);
	}
	for (j = 0; j < numtrace; j++) {
	  for (i = 0; i < 120; i++)
	    outhead[i] = 0;
	  printf ("trace-%d-\r", j + 1); fflush(stdout);
	  fseek (f1, tracepointers[j], SEEK_SET);
	  fread (&blockid, 2, 1, f1);
	  if (reversed)
	    swapshort (&blockid);
	  if (blockid == 0x2244) {
	    /* reversed=1; should already know this */
	    fprintf (stderr, "Opps, I've blown it here.... (line:%d)\n", __LINE__);
	    exit (-6);
	  }
	  if (blockid != 0x4422) {
	    fprintf (stderr, "Not a SEG-2 trace header.  Can not process %x (line %d)\n", blockid, __LINE__);
	    exit (-7);
	  }
	  fread (&blockleng, 2, 1, f1);
	  if (reversed)
	    swapshort ((short *) &blockleng);
	  fread (&datalength, 4, 1, f1);
	  if (reversed)
	    swaplong ((long *) &datalength);
	  fread (&numsamples, 4, 1, f1);
	  if (reversed)
	   swaplong ((long *) &numsamples);
	if (numsamples >= MAXSAMP){ 
	fprintf(stderr, "Your data contains more samples than I can handle\n");
	exit(-8);
	}
	  fread (&datatype, 1, 1, f1);
	  if (datatype > 5 || datatype < 1) {
	    fprintf (stderr,"Data type %d not available/valid\n", (int) datatype);
	    break;
	  }
	  outhead[57] = numsamples;
	  fread (reserved, 1, 19, f1);
	  fread (&stringlength, 2, 1, f1);
	  if (reversed)
	    swapshort (&stringlength);
	  while (0 != stringlength) {
	    fread (string1, 1, stringlength - 2, f1);
	    keycheck (keywords);
	    fread (&stringlength, 2, 1, f1);
	    if (reversed)
	      swapshort (&stringlength);
	  }
	
	  fseek (f1, blockleng + tracepointers[j], SEEK_SET);
	  switch ((int) datatype)
	  {
	case 1:
		fread (iinbuf, 2, (int) numsamples, f1);
		if (reversed) {
		for (i = 0; i < (int) numsamples; i++)
		swapshort (&iinbuf[i]);
		}
		for (k = 0; k < numsamples; k++)
		outbuf[k] = iinbuf[k];
		break;
	case 2:
		fread (linbuf, 4, (int) numsamples, f1);
		if (reversed) {
		for (i = 0; i < (int) numsamples; i++)
		swaplong (&linbuf[i]);
		}
		for (k = 0; k < numsamples; k++)
		  outbuf[k] = linbuf[k];
		break;
	case 3:
	{
		unsigned long totalbytes, subpointer;
		unsigned int expo;
		long longdat;
		short int sdata; /*modified version by PM */
		totalbytes = (numsamples * 5) / 2;
		fread (cinbuf, 1, (size_t) totalbytes, f1);
		/*  Original Code SU-36 3rd Party
		for (k = 0; k < (numsamples);) {
		subpointer = (k / 4) * 5;
		expo = (unsigned) iinbuf[subpointer++];
		for (i = 0; i < 4; i++) {
		if (0x8000 & iinbuf[subpointer]) 
		longdat = 0xffff8000;
		else 
		longdat = 0;
		longdat = longdat | ((long) iinbuf[subpointer++] << (0x000f & expo));
		expo >>= 4;
		outbuf[k++] = longdat;
		}
		}
		*/
	/*modified version by P.Michaels <[e-mail]> */
	/*fixes sawtooth conversion error on large negative values */
		for (k = 0; k < (numsamples);) {
		subpointer = (k / 4) * 5;
		expo = (unsigned) iinbuf[subpointer++];
		for (i = 0; i < 4; i++) {
	                sdata=iinbuf[subpointer++];
                     		 if (sdata>0) 
                         	  	{       
                                	longdat=(long) sdata;
                                	longdat=longdat << (0x000f & expo);
                        	   	}       
                	        	else    
                         	{       
                               	 	longdat=(long) -sdata; 
                                	longdat=longdat <<  (0x000f & expo);
                              	  	longdat = -longdat; 
                         	 	} /* endif */ 
                		expo >>= 4;
		outbuf[k++] = longdat;
                	} /* next i */
        	} /* next k */
		/* end of modifications */
		if (reversed) {
		for (i = 0; i < numsamples; i++)
		swaplong (&outbuf[i]);
		}
	}
	break;
	case 4:
	fread (finbuf, 4, (int) numsamples, f1);
		if (reversed) {
		long *buf = (long *) dinbuf;
		for (i = 0; i < numsamples; i++)
		swaplong (&buf[i]);
		}
		for (k = 0; k < numsamples; k++)
		outbuf[k] = finbuf[k];
		break;
	case 5:
	fread (dinbuf, 8, (int) numsamples, f1);
		if (reversed) {
		long *buf = (long *) dinbuf;
		for (i = 0; i < numsamples * 2; i++)
		swaplong (&buf[i]);
		}
		for (k = 0; k < numsamples; k++)
		outbuf[k] = dinbuf[k];
		break;
	}      /* end switch */
	/* set vertical stack traces=1 */
/* 
	if(outhead[15]==0) outhead[15]=1; 
	if(outheadf[59]==0.0) gain=1.; else gain=outheadf[59];
	for(i=0;i<numsamples;i++)outbuf[i]*=gain/outhead[15];
*/
	/* assign the original field record number as the current file number */
	if (outheadl[2] == 0) outheadl[2] = (long) curf;
	/* set trace type  as  sesmic data */
	if (outhead[14] == 0) outhead[14] = 1;

	/* set last trace flag (modified segy) */
	if (j == numtrace - 1 && outhead[87] == 0) {
	outhead[87] = 1;
	ssn = ssn + 1;
	}

	/* from rec-station-number and source-station-number (93 and 94) */
	/* distance from source to receiver */
	/* outheadl[9]=(long)(abs(outhead[93]-outhead[92])); */
	/* set group for trace one and roll switch position  */
	outhead[85] = outhead[86] = 
	(int) (1 + labs ((long) outhead[92] - outheadl[3]));

	/* special case, execute on first pass only... */
	if (first == 1) {
	first = 0;
	segyreelheader[1606] = numtrace;
	segyreelheader[1608] = outhead[58];
	segyreelheader[1609] = outhead[58];
	segyreelheader[1610] = numsamples;
	segyreelheader[1611] = numsamples;
	segyreelheader[1612] = 2;

	if (!reversed) /* swap only if we are on a little endian  machine */
	{  
		for (k = 1600; k < 1606; k += 2)
		swaplong ((long *)&segyreelheader[k]);
	  for (k = 1606; k < 1630; k++)
	    swapshort((short *)&segyreelheader[k]);
	}
	fwrite (segyreelheader, 1, 3600, f2);  /*  create the segy headers */
	}

	if (!reversed) { /* swap only if we are on a little endian  machine */
	/* first swap longs */
	for (k = 0; k < 7; k++)
		swaplong((long *)&outheadl[k]);
	for (k = 9; k < 17; k++)
	  swaplong((long *)&outheadl[k]);
	for (k = 18; k < 22; k++)
	  swaplong((long *)&outheadl[k]);
	/* now swap the shorts */ 
	for (k = 14; k < 18; k++)
	  swapshort((short *)&outhead[k]);
	for (k = 34; k < 36; k++)
	  swapshort((short *)&outhead[k]);
	/* for(k=44;k<90;k++)  *//* error: should have gone beyond 95 word */
	for (k = 44; k < 95; k++)
	  swapshort((short *)&outhead[k]);
	}
	if (120 != (k = fwrite (outhead, 2, 120, f2))) {      /* write header */
	fprintf (stderr,"\nWrite failure during header write\n");
	exit (-9);
	}
	if (!reversed) {
	for (k = 0; k < numsamples; k++)
		swaplong ((long *)&outbuf[k]);
	}
	if ((int) numsamples != (k = fwrite (outbuf, 4, (int) numsamples, f2))) {
	fprintf (stderr,"Write failure during trace write\n");
	exit (-10);
	}

	}      /* end trace loop */
	fclose (f1);
	outhead[87] = 0;    /* reset last trace flag. */
	}        /* end kk loop */
	fclose(f2);
	return 0;
}        /* end main */