sac2mseed.c

C编写的格式转换程序

/***************************************************************************
 * sac2mseed.c
 *
 * Simple waveform data conversion from SAC timeseries to Mini-SEED.
 * No support is included for SAC spectral or generic X-Y data.
 *
 * Written by Chad Trabant, IRIS Data Management Center
 *
 * modified 2007.284
 ***************************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include <errno.h>
#include <ctype.h>

#include <libmseed.h>

#include "sacformat.h"

#define VERSION "1.6"
#define PACKAGE "sac2mseed"

#if defined (LWP_WIN32)
  #define strtoull _strtoui64
#endif

struct listnode {
  char *key;
  char *data;
  struct listnode *next;
};

static void packtraces (flag flush);
static int sac2group (char *sacfile, MSTraceGroup *mstg);
static int parsesac (FILE *ifp, struct SACHeader *sh, float **data, int format,
		     int verbose, char *sacfile);
static int readbinaryheader (FILE *ifp, struct SACHeader *sh, int *format,
			     int *swapflag, int verbose, char *sacfile);
static int readbinarydata (FILE *ifp, float *data, int datacnt,
			   int swapflag, int verbose, char *sacfile);
static int readalphaheader (FILE *ifp, struct SACHeader *sh);
static int readalphadata (FILE *ifp, float *data, int datacnt);
static int swapsacheader (struct SACHeader *sh);
static int writemetadata (struct SACHeader *sh);
static int parameter_proc (int argcount, char **argvec);
static char *getoptval (int argcount, char **argvec, int argopt);
static int readlistfile (char *listfile);
static void addnode (struct listnode **listroot, char *key, char *data);
static void record_handler (char *record, int reclen, void *handlerdata);
static void usage (void);

static int   verbose     = 0;
static int   packreclen  = -1;
static int   encoding    = 11;
static int   byteorder   = -1;
static int   sacformat   = 0;
static char  srateblkt   = 0;
static char *forcenet    = 0;
static char *forceloc    = 0;
static char *outputfile  = 0;
static FILE *ofp         = 0;
static char *metafile    = 0;
static FILE *mfp         = 0;
static long long int datascaling = 0;

/* A list of input files */
struct listnode *filelist = 0;

static MSTraceGroup *mstg = 0;

static int packedtraces  = 0;
static int packedsamples = 0;
static int packedrecords = 0;

int
main (int argc, char **argv)
{
  struct listnode *flp;
  
  /* Process given parameters (command line and parameter file) */
  if (parameter_proc (argc, argv) < 0)
    return -1;
  
  /* Init MSTraceGroup */
  mstg = mst_initgroup (mstg);
  
  /* Open the output file if specified */
  if ( outputfile )
    {
      if ( strcmp (outputfile, "-") == 0 )
        {
          ofp = stdout;
        }
      else if ( (ofp = fopen (outputfile, "wb")) == NULL )
        {
          fprintf (stderr, "Cannot open output file: %s (%s)\n",
                   outputfile, strerror(errno));
          return -1;
        }
    }
  
  /* Open the metadata output file if specified */
  if ( metafile )
    {
      if ( strcmp (metafile, "-") == 0 )
        {
          mfp = stdout;
        }
      else if ( (mfp = fopen (metafile, "wb")) == NULL )
        {
          fprintf (stderr, "Cannot open metadata output file: %s (%s)\n",
                   metafile, strerror(errno));
          return -1;
        }
    }
  
  /* Read input SAC files into MSTraceGroup */
  flp = filelist;
  while ( flp != 0 )
    {
      if ( verbose )
	fprintf (stderr, "Reading %s\n", flp->data);

      sac2group (flp->data, mstg);
      
      flp = flp->next;
    }
  
  fprintf (stderr, "Packed %d trace(s) of %d samples into %d records\n",
	   packedtraces, packedsamples, packedrecords);
  
  /* Make sure everything is cleaned up */
  if ( ofp )
    fclose (ofp);
  
  if ( mfp )
    fclose (mfp);
  
  return 0;
}  /* End of main() */


/***************************************************************************
 * packtraces:
 *
 * Pack all traces in a group using per-MSTrace templates.
 *
 * Returns 0 on success, and -1 on failure
 ***************************************************************************/
static void
packtraces (flag flush)
{
  MSTrace *mst;
  int trpackedsamples = 0;
  int trpackedrecords = 0;
  
  mst = mstg->traces;
  while ( mst )
    {
      
      if ( mst->numsamples <= 0 )
	{
	  mst = mst->next;
	  continue;
	}
      
      trpackedrecords = mst_pack (mst, &record_handler, 0, packreclen, encoding, byteorder,
				  &trpackedsamples, flush, verbose-2, (MSRecord *) mst->prvtptr);
      
      if ( trpackedrecords < 0 )
	{
	  fprintf (stderr, "Error packing data\n");
	}
      else
	{
	  packedrecords += trpackedrecords;
	  packedsamples += trpackedsamples;
	}
      
      mst = mst->next;
    }
}  /* End of packtraces() */


/***************************************************************************
 * sac2group:
 * Read a SAC file and add data samples to a MSTraceGroup.  As the SAC
 * data is read in a MSRecord struct is used as a holder for the input
 * information.
 *
 * Returns 0 on success, and -1 on failure
 ***************************************************************************/
static int
sac2group (char *sacfile, MSTraceGroup *mstg)
{
  FILE *ifp = 0;
  MSRecord *msr = 0;
  MSTrace *mst;
  struct blkt_1000_s Blkt1000;
  struct blkt_1001_s Blkt1001;
  struct blkt_100_s Blkt100;
  
  struct SACHeader sh;
  float *fdata = 0;
  int32_t *idata = 0;
  int dataidx;
  int datacnt;
  long long int scaling = datascaling;
  
  /* Open input file */
  if ( (ifp = fopen (sacfile, "rb")) == NULL )
    {
      fprintf (stderr, "Cannot open input file: %s (%s)\n",
	       sacfile, strerror(errno));
      return -1;
    }
  
  /* Parse input SAC file into a header structure and data buffer */
  if ( (datacnt = parsesac (ifp, &sh, &fdata, sacformat, verbose, sacfile)) < 0 )
    {
      fprintf (stderr, "Error parsing %s\n", sacfile);
      
      return -1;
    }

  /* Write metadata to file if requested */
  if ( mfp )
    {
      if ( verbose )
	fprintf (stderr, "[%s] Writing metadata to %s\n", sacfile, metafile);
      
      if ( writemetadata (&sh) )
	{
	  fprintf (stderr, "Error writing metadata to file '%s'\n", metafile);
	  
	  return -1;
	}
    }
  
  /* Open output file if needed */
  if ( ! ofp )
    {
      char mseedoutputfile[1024];
      int namelen;
      strncpy (mseedoutputfile, sacfile, sizeof(mseedoutputfile)-6 );
      namelen = strlen (sacfile);
      
      /* Truncate file name if .sac is at the end */
      if ( namelen > 4 )
	if ( (*(mseedoutputfile + namelen - 1) == 'c' || *(mseedoutputfile + namelen - 1) == 'C') &&
             (*(mseedoutputfile + namelen - 2) == 'a' || *(mseedoutputfile + namelen - 2) == 'A') &&
             (*(mseedoutputfile + namelen - 3) == 's' || *(mseedoutputfile + namelen - 3) == 'S') &&
             (*(mseedoutputfile + namelen - 4) == '.') )
	  
	  {
	    *(mseedoutputfile + namelen - 4) = '\0';
	  }
      
      /* Add .mseed to the file name */
      strcat (mseedoutputfile, ".mseed");
      
      if ( (ofp = fopen (mseedoutputfile, "wb")) == NULL )
        {
          fprintf (stderr, "Cannot open output file: %s (%s)\n",
                   mseedoutputfile, strerror(errno));
          return -1;
        }
    }
  
  if ( ! (msr = msr_init(msr)) )
    {
      fprintf (stderr, "Cannot initialize MSRecord strcture\n");
      return -1;
    }
  
  /* Determine autoscaling */
  if ( scaling == 0 && encoding != 4 )
    {
      float datamin, datamax;
      int fractional;
      long long int autoscale;
      
      /* Determine data sample minimum and maximum
       * Detect if scaling by 1 will result in truncation (fractional=1) */
      datamin = datamax = *fdata;
      fractional = 0;
      for ( dataidx=1; dataidx < datacnt; dataidx++ )
	{
	  if ( *(fdata+dataidx) < datamin ) datamin =  *(fdata+dataidx);
	  if ( *(fdata+dataidx) > datamax ) datamax =  *(fdata+dataidx);
	  
	  if ( ! fractional )
	    if ( *(fdata+dataidx) - (int) *(fdata+dataidx) > 0.000001 )
	      fractional = 1;
	}
      
      autoscale = 1;
      
      if ( fractional )
	{
	  for (autoscale=1; abs ((int32_t) (datamax * autoscale)) < 100000; autoscale *= 10) {}
	  
	  if ( abs ((int32_t) (datamin * autoscale)) < 10 )
	    fprintf (stderr, "WARNING Large sample value range (%g/%g), autoscaling might be a bad idea\n",
		     datamax, datamin);
	}
      
      scaling = autoscale;
    }

  /* Populate MSRecord structure with header details */
  if ( strncmp (SUNDEF, sh.knetwk, 8) ) ms_strncpclean (msr->network, sh.knetwk, 2);
  if ( strncmp (SUNDEF, sh.kstnm, 8) ) ms_strncpclean (msr->station, sh.kstnm, 5);
  if ( strncmp (SUNDEF, sh.khole, 8) ) ms_strncpclean (msr->location, sh.khole, 2);
  if ( strncmp (SUNDEF, sh.kcmpnm, 8) ) ms_strncpclean (msr->channel, sh.kcmpnm, 3);
  
  if ( forcenet )
    ms_strncpclean (msr->network, forcenet, 2);
  
  if ( forceloc )
    ms_strncpclean (msr->location, forceloc, 2);
  
  msr->starttime = ms_time2hptime (sh.nzyear, sh.nzjday, sh.nzhour, sh.nzmin, sh.nzsec, sh.nzmsec * 1000);
  
  /* Adjust for Begin ('B' SAC variable) time offset */
  msr->starttime += (double) sh.b * HPTMODULUS;

  /* Calculate sample rate from interval(period) rounding to nearest 0.000001 Hz */
  msr->samprate = (double) ((int)((1 / sh.delta) * 100000 + 0.5)) / 100000;
  
  msr->samplecnt = msr->numsamples = datacnt;
  
  /* Data sample type and sample array */
  if ( encoding == 4 )
    {
      msr->sampletype = 'f';
      msr->datasamples = fdata;
    }
  else
    {
      /* Create an array of scaled integers */
      idata = (int32_t *) malloc (datacnt * sizeof(int32_t));
      
      if ( verbose )
	fprintf (stderr, "[%s] Creating integer data scaled by: %lld\n", sacfile, scaling);
      
      for ( dataidx=0; dataidx < datacnt; dataidx++ )
	*(idata + dataidx) = (int32_t) (*(fdata + dataidx) * scaling);
      
      msr->sampletype = 'i';
      msr->datasamples = idata;
    }
  
  if ( verbose >= 1 )
    {
      fprintf (stderr, "[%s] %d samps @ %.6f Hz for N: '%s', S: '%s', L: '%s', C: '%s'\n",
	       sacfile, msr->numsamples, msr->samprate,
	       msr->network, msr->station,  msr->location, msr->channel);
    }
  
  if ( ! (mst = mst_addmsrtogroup (mstg, msr, 0, -1.0, -1.0)) )
    {
      fprintf (stderr, "[%s] Error adding samples to MSTraceGroup\n", sacfile);
    }
  
  /* Create an MSRecord template for the MSTrace by copying the current holder */
  if ( ! mst->prvtptr )
    {
      mst->prvtptr = msr_duplicate (msr, 0);
      
      if ( ! mst->prvtptr )
	{
	  fprintf (stderr, "[%s] Error duplicate MSRecord for template\n", sacfile);
	  return -1;
	}
      
      /* Add blockettes 1000 & 1001 to template */
      memset (&Blkt1000, 0, sizeof(struct blkt_1000_s));
      msr_addblockette ((MSRecord *) mst->prvtptr, (char *) &Blkt1000,
			sizeof(struct blkt_1001_s), 1000, 0);
      memset (&Blkt1001, 0, sizeof(struct blkt_1001_s));
      msr_addblockette ((MSRecord *) mst->prvtptr, (char *) &Blkt1001,
			sizeof(struct blkt_1001_s), 1001, 0);
      
      /* Add blockette 100 to template if requested */
      if ( srateblkt )
	{
	  memset (&Blkt100, 0, sizeof(struct blkt_100_s));
	  Blkt100.samprate = (float) msr->samprate;
	  msr_addblockette ((MSRecord *) mst->prvtptr, (char *) &Blkt100,
			    sizeof(struct blkt_100_s), 100, 0);
	}
    }
  
  packtraces (1);
  packedtraces += mstg->numtraces;
  
  fclose (ifp);
  
  if ( ofp  && ! outputfile )
    {
      fclose (ofp);
      ofp = 0;
    }

  if ( fdata )
    free (fdata);
  
  if ( idata )
    free (idata);
  
  msr->datasamples = 0;
  
  if ( msr )
    msr_free (&msr);
  
  return 0;
}  /* End of sac2group() */


/***************************************************************************
 * parsesac:
 *
 * Parse a SAC file, autodetecting format dialect (ALPHA,
 * binary, big or little endian).  Results will be placed in the
 * supplied SAC header struct and data (float sample array in host
 * byte order).  The data array will be allocated by this routine and
 * must be free'd by the caller.  The data array will contain the
 * number of samples indicated in the SAC header (sh->npts).
 *
 * The format argument is interpreted as:
 * 0 : Unknown, detection needed
 * 1 : ALPHA
 * 2 : Binary, byte order detection needed
 * 3 : Binary, little endian
 * 4 : Binary, big endian
 *
 * Returns number of data samples in file or -1 on failure.
 ***************************************************************************/
static int
parsesac (FILE *ifp, struct SACHeader *sh, float **data, int format,
	  int verbose, char *sacfile)