las2su.c

su 的源代码库

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

/* LAS2SU: $Revision: 1.2 $ ; $Date: 2006/11/07 22:58:42 $	*/

#include "par.h"
#include "su.h"
#include "segy.h"

/*********************** self documentation **********************/
char *sdoc[] = {
" 								",
" LAS2SU - convert las2 format well log curves to su traces	",
" 								",
" las2su <stdin nskip= nlogs= >stdout [optional params] 	",
" 								",
" Required parameters:						",
" nlogs=		number of well log curves in the LAS file",
" 								",
" Optional parameters:						",
" dz=0.5		input depth sampling (ft)		",
" m=0			output (d1,f1) in feet 			",
" 			=1 output (d1,f1) in meters 		",
" ss=0			do not subsample (unless nz > 32767 ) 	",
" 			=1 pass every other sample 		",
" verbose=0		=1 to echo las header lines to screen	",
" outhdr=las_hdr.asc	name of file for las headers		",
" 								",
" Notes:							",
" 1. First log curve MUST BE depth.				",
" 2. If number of depth levels > 32767 (segy NT limit)		",
"    then input is subsampled by factor of 2 or 4 as needed	",
" 3. Logs may be isolated using tracl header word (1,2,...,nlogs) ",
"    tracl=1 is depth 						",
" 									",
NULL};

/* Credits:
 *	Chris Liner based on code by Ferhan Ahmed and a2b.c (June 2005)
 *	CWP: John Stockwell 31 Oct 2006, combining lasstrip and 
 *      CENPET: lasstrip 2006 by Werner Heigl
 */

/**************** end self doc ***********************************/

/* defined quantities */
#define LAS_MAXLINE 1000
#define LAS_HDR_DEF "las_header.asc"

/* function prototype for subroutine used internally */
int las_getnewline(char line[], int maxline);

segy tr;	/* output trace structure */ 

int
main(int argc, char **argv)
{
	int nlogs;		  /* number of log curves in las file  	*/
	float **x=NULL;		  /* binary floats		   	*/
	float len=0;		  /* line length			*/
	char string[300];	  /* one line of input ascii file	*/
	char c1[30];		  /* ascii form of log value 		*/
	int m;			  /* flag for metric (d1,f1) output     */
	int nzmax;	  	  /* max number of depth levels 	*/
	int nz,nzold;	  	  /* actual number of depth levels 	*/
	int iz,izz,icurve,i;	 /* counters 				*/
	short verbose; 	  /* if 1(yes) echo las header lines to stderr  */
	float dz;		/* depth sample rate */
	int ss;			/* subsample flag 			*/
	char line[LAS_MAXLINE];
	char *outhdr;	     /* name of file holding LAS header    */
	FILE *outhdrfp;	   /* ... its file pointer	       */


	/* Hook up getpar */
	initargs(argc, argv);
	requestdoc(1);

	/* get parameters */
	MUSTGETPARINT("nlogs",&nlogs);
	if (!getparshort("verbose", &verbose))  verbose = 0;
	if (!getparint("m", &m))  	m = 0;
	if (!getparint("ss", &ss))  	ss = 0;
	if (!getparfloat("dz", &dz))  	dz = 0.5;
	if (!getparstring("outhdr", &outhdr))  outhdr = LAS_HDR_DEF;
	outhdrfp = efopen(outhdr, "w");


	/* max depth levels: 32767 (segy limit) */
	nzmax = SU_NFLTS;

	/* alloc array to hold float log values */
	x = ealloc2float(nzmax,nlogs);

	/* zero array */
	memset((void *) x[0], 0, nlogs*nzmax*FSIZE);
		
	/* copy LAS header to outhdr */
	do {
		len = las_getnewline(line,LAS_MAXLINE);
		i = 0;
		while (line[i]!='\0') {
			fputc(line[i],outhdrfp);
			++i;
		}
		fputc('\n',outhdrfp);
	} while (strncmp("~A",line,2)!=0 );
	efclose(outhdrfp);


	/* initialize depth counter */
	nz = 0;
	iz = 0;

	/* get each line as a string */
	while(fgets(string,200,stdin) != NULL) {

		/* read first token */
    	 	strcpy(c1,strtok(string,"    \n\t"));

		/* load log value into float array */
    	 	x[0][iz] = atof(c1);

		for (icurve = 1; icurve < nlogs; ++icurve) {

			/* read next token get ascii log value */
    	 		strcpy(c1,strtok(NULL,"    \n\t"));

			/* load this log value into float array */
    	 		x[icurve][iz] = atof(c1);

  		}

		/* bump depth counter */
		iz += 1;

	}

	/* number of depth values in log */
	nz = iz;
	warn("nz=%i",nz);

	/* check that nz limit is not exceeded, or subsampling requested */
	if ( nz > nzmax || ss == 1 ) {
		nzold = nz;
		/* reset number of depth samples and sample rate */
		nz = nz/2;
		dz = 2.0 * dz;
		/* subsample */
		for (icurve=0 ; icurve < nlogs ; ++icurve){
			for (iz = 0 ; iz < nz ; ++iz){
				izz = 2*iz;
				x[icurve][iz] = x[icurve][izz];
			}
		}
		warn("New: nz=%i dz=%g ft\n",nz,dz);
	}

	/* check again (possible deep well with 0.25 ft sampling) */
	if ( nz > nzmax || ss == 1 ) {
		nzold = nz;
		/* reset number of depth samples and sample rate */
		nz = nz/2;
		dz = 2.0 * dz;
		/* subsample */
		for (icurve=0 ; icurve < nlogs ; ++icurve){
			for (iz = 0 ; iz < nz ; ++iz){
				izz = 2*iz;
				x[icurve][iz] = x[icurve][izz];
			}
		}
		warn("New: nz=%i dz=%g ft\n",nz,dz);
	}

	/* set up output trace headers */
	tr.trid = 1;			/* su time traces (trick) */
	tr.ns = nz;			/* samples per trace */
	tr.dt = 1000*dz;		/* time sample rate (trick) */
	if (m == 0) {
		tr.d1 = dz;			/* actual dz (in ft) */
		tr.f1 = x[0][0];		/* first depth value (in ft) */
	} else {
		tr.d1 = dz/3.28084;		/* actual dz (in m) */
		tr.f1 = x[0][0]/3.28084;	/* first depth value (in m) */
	}

	for (icurve=0 ; icurve < nlogs ; ++icurve){
		tr.tracl = icurve+1;

		for (iz = 0 ; iz < nz ; ++iz){
			tr.data[iz] = x[icurve][iz];
		}
			
		puttr(&tr);
	}


	return EXIT_SUCCESS;
}
int las_getnewline(char s[], int lim)
/***************************************************************************
las_getnewline: read a line from stdin into s[] and return length of line 
****************************************************************************
Input:
s[]	input string
lim	maximum length of line
Returns:
i	length of line 
****************************************************************************
Author: CENPET: Werner Heigl
****************************************************************************/
{
	int c=0,i;

	for (i=0; i<lim-1 && (c=getchar())!=EOF && c!='\n'; ++i)
		s[i] = c;
	if (c=='\n') {
		s[i] = c;
		++i;
	}
	s[i] = '\0';
	return i;
}