Main Page   Class Hierarchy   Compound List   File List   Compound Members   File Members  

seiutils.c

Go to the documentation of this file.
00001 /***************************************************************************
00002  * seiutils.c - routines for writing seisan waveform files
00003  *
00004  * NB - the general parts of this file are also used by the GSE
00005  *      put away routines
00006  *
00007  * Seisan file format synopsis:
00008  * 
00009  * Record | Size | Description
00010  * -------+------+--------------------------------------------------
00011  * 1      | 80   | Network name, number of stations, event time
00012  * 2      | 80   | Not used
00013  * 3      | 80   | 1st, 2nd, 3rd channel details
00014  * 4-N    | 80   | Other channel details, 3 per record - there is
00015  *        |      | a minimum of 12 channel headers, so N has
00016  *        |      | a minimum value of 6 (but no maximum value)
00017  * N+1    | 1040 | Channel 1 header - time and response data
00018  * N+2    | X    | Channel 1 data - entire data is written
00019  *        |      | as a single record
00020  * ...    |      | Channel 2,3,... header/data follow
00021  *
00022  * Sizes are in bytes. Each record is written with its length as
00023  * a 4 byte integer immediately before and after the data:
00024  *   <length (4 bytes)> <record (x bytes)> <length (4 bytes)>
00025  *
00026  * To work correctly, you must call the routines in this file in the
00027  * following order:
00028  *
00029  *      open_seisan_file()
00030  *      foreach channel
00031  *              add_seisan_channel()
00032  *      foreach channel
00033  *              start_seisan_channel()
00034  *              add_seisan_channel_data()
00035  *              end_seisan_channel()
00036  *      close_seisan_file()
00037  *
00038  * S. Flower, Oct 2001
00039  *
00040  ***************************************************************************/
00041 
00042 #include <stdio.h>
00043 #include <time.h>
00044 #include <stdlib.h>
00045 #include <string.h>
00046 #include <errno.h>
00047 #include <math.h>
00048 #include <float.h>
00049 
00050 /* operating system specific includes */
00051 #if defined (_WINNT)
00052 #include <direct.h>
00053 #include <io.h>
00054 #elif defined  (_SOLARIS)
00055 #include <unistd.h>
00056 #else
00057 #error "_WINNT or _SOLARIS must be set before compiling"
00058 #endif
00059 
00060 #include <sys/types.h>
00061 #include <sys/stat.h>
00062 
00063 #include "earthworm.h"
00064 #include "trace_buf.h"
00065 #include "swap.h"
00066 #include "ws_clientII.h"
00067 #include "site.h"
00068 #include "time_ew.h"
00069 
00070 #include "seihead.h"
00071 
00072 /* private global variables */
00073 static char output_dir [MAX_DIR_LEN];   /* directory to use for output file */
00074 static char tmp_dirname [MAX_DIR_LEN];  /* name of the temporary directory for this request */
00075 static char network_name [50];          /* name of the network for this request */
00076 static double start_time;               /* start time of requested data */
00077 static double duration;                 /* duration of requested data */
00078 static int n_channels;                  /* number of channels in this request */
00079 static struct Sei_channel_details *channels;            /* channel details for this request */
00080 static struct Sei_channel_details *current_channel;     /* channel currently being processed */
00081 static FILE *channel_fp;                /* currently open data file */
00082 static int output_native_numbers;       /* if TRUE write numbers without byte swapping */
00083 static int tmp_dir_count = 0;           /* count of temporary directories */
00084 
00085 /* private forward declarations */
00086 static void write_seisan (FILE *fp, long length, void *data);
00087 
00088 
00089 /***************************************************************************
00090  * open_seisan_file
00091  *
00092  * Description: open a Seisan file and write the start of the header
00093  *
00094  * Input parameters: od - the directory for seisan files
00095  *                   nn - the network name
00096  *                   st - the start time for the data
00097  *                   dur - length of data in seconds
00098  *                   onn - if TRUE write numbers without byte swapping
00099  * Output parameters: channel details 
00100  * Returns: TRUE on success, FALSE on error
00101  *
00102  * Comments:
00103  *
00104  ***************************************************************************/
00105 int open_seisan_file (char *od, char *nn, double st, double dur, int onn)
00106 
00107 {
00108 
00109   /* create a temporary directory for the parts of the seisan file */
00110   sprintf (tmp_dirname, "%s%c%d_%d_sei.tmp",
00111            od, DIR_DELIM, getpid (), tmp_dir_count ++);
00112   if (CreateDir (tmp_dirname) != EW_SUCCESS) return FALSE;
00113 
00114   /* reset global variables */
00115   strcpy (output_dir, od);
00116   strcpy (network_name, nn);
00117   start_time = st;
00118   duration = dur;
00119   n_channels = 0;
00120   channels = current_channel = 0;
00121   output_native_numbers = onn;
00122 
00123   return TRUE;
00124 
00125 }
00126 
00127 /***************************************************************************
00128  * add_seisan_channel
00129  *
00130  * Description: add a channel to the current seisan file
00131  *
00132  * Input parameters: chan_name, type - the channel name and type
00133  *                                     (e.g. EKR1, SZ)
00134  * Output parameters: none
00135  * Returns: TRUE if channel added OK, FALSE otherwise
00136  *
00137  * Comments: all channels must be added using this routine before any
00138  * channel data is written - channels that are not added will be
00139  * ignored
00140  *
00141  ***************************************************************************/
00142 int add_seisan_channel (char *chan_name, char *chan_type)
00143 
00144 {
00145 
00146   struct Sei_channel_details *channel_ptr;
00147 
00148 
00149   /* allocate memory for this channel */
00150   channel_ptr = realloc (channels, sizeof (struct Sei_channel_details) * (n_channels +1));
00151   if (! channel_ptr) return FALSE;
00152   channels = channel_ptr;
00153   n_channels ++;
00154 
00155   /* fill out the new channel details - fill these in so that if a channel is
00156    * completely missing it will appear as a 1hz trace of the same duration as
00157    * the complete request */
00158   channel_ptr = channels + n_channels -1;
00159   strcpy (channel_ptr->chan_name, chan_name);
00160   strcpy (channel_ptr->chan_type, chan_type);
00161   channel_ptr->start_time = start_time;
00162   channel_ptr->n_samples = (long) duration;
00163   channel_ptr->n_written = 0;
00164   channel_ptr->sample_rate = 1.0;
00165   channel_ptr->channel_count = n_channels -1;
00166   strcpy (channel_ptr->filename, "");
00167 
00168   return TRUE;
00169 
00170 }
00171 
00172 /***************************************************************************
00173  * start_seisan_channel
00174  *
00175  * Description: start writing data for the next channel
00176  *
00177  * Input parameters: chan_name, type - the channel name and type
00178  *                                     (e.g. EKR1, SZ)
00179  *                   start_time - start time of data for this channel
00180  *                   sample_rate - the sample rate for the data
00181  *                   n_samples - the number of samples
00182  * Output parameters: none
00183  * Returns: TRUE if write succeeded, FALSE otherwise
00184  *
00185  * Comments: call this before add_seisan_channel_data() to start writing
00186  *           the next channel
00187  *
00188  ***************************************************************************/
00189 int start_seisan_channel (char *chan_name, char *chan_type,
00190                           double start_time, double sample_rate,
00191                                                   long n_samples)
00192 
00193 {
00194 
00195   int count;
00196   long value;
00197   struct Sei_channel_details *channel_ptr;
00198 
00199 
00200   /* find the correpsonding channel details - if they don't exist
00201    * then return without error - this will cause the channel to
00202    * be ignored */
00203   current_channel = 0;
00204   for (count=0, channel_ptr=channels; count<n_channels; count++, channel_ptr++)
00205   {
00206     if ((! strcmp (channel_ptr->chan_name, chan_name)) &&
00207         (! strcmp (channel_ptr->chan_type, chan_type)))
00208     {
00209       current_channel = channel_ptr;
00210       break;
00211     }
00212   }
00213   if (! current_channel) return TRUE;
00214 
00215   /* fill in the channel details */
00216   current_channel->start_time = start_time;
00217   current_channel->sample_rate = sample_rate;
00218   current_channel->n_samples = n_samples;
00219   sprintf (current_channel->filename, "%s%c%d.tmp", tmp_dirname, DIR_DELIM, current_channel->channel_count);
00220 
00221   /* open a file for this channel */
00222   channel_fp = fopen (current_channel->filename, "wb");
00223   if (! channel_fp)
00224   {
00225     strcpy (channel_ptr->filename, "");
00226     return FALSE;
00227   }
00228 
00229   /* write the length of the data at the start of the block */
00230   value = current_channel->n_samples * sizeof (long);
00231   if (! output_native_numbers) SwapLong (&value);
00232   fwrite (&value, sizeof (long), 1, channel_fp);
00233 
00234   return TRUE;
00235 
00236 }
00237 
00238 /***************************************************************************
00239  * add_seisan_channel_data
00240  * 
00241  * Input parameters: data_len - the number of samples to add
00242  *                   data - the data as an array
00243  * Output parameters: none
00244  * Returns: none
00245  *
00246  * Comments:
00247  *
00248  ***************************************************************************/
00249 void add_seisan_channel_data (long data_len, long *data)
00250 
00251 {
00252 
00253   int count;
00254   long value;
00255 
00256 
00257   /* output the data */
00258   for (count=0; count<data_len; count++)
00259   {
00260     if (current_channel->n_written >= current_channel->n_samples) break;
00261 
00262     value = *(data + count);
00263     if (! output_native_numbers) SwapLong (&value);
00264     fwrite (&value, 4, 1, channel_fp);
00265 
00266     current_channel->n_written ++;
00267   }
00268 
00269 }
00270 
00271 /***************************************************************************
00272  * end_seisan_channel
00273  *
00274  * Description: call this to complete writing data to a channel
00275  *
00276  * Input parameters: none
00277  * Output parameters: none
00278  * Returns: TRUE if channel completed OK, FALSE otherwise
00279  *
00280  * Comments:
00281  *
00282  ***************************************************************************/
00283 int end_seisan_channel (void)
00284 
00285 {
00286 
00287   int ret_val;
00288   long value;
00289 
00290 
00291   /* add any gap data to the end of the file */
00292   value = SEISAN_MISSING_DATA_FLAG;
00293   if (! output_native_numbers) SwapLong (&value);
00294   while (current_channel->n_written < current_channel->n_samples)
00295   {
00296     fwrite (&value, 4, 1, channel_fp);
00297     current_channel->n_written ++;
00298   }
00299   
00300   /* write the length of the data at the end of the block */
00301   value = current_channel->n_samples * sizeof (long);
00302   if (! output_native_numbers) SwapLong (&value);
00303   fwrite (&value, sizeof (long), 1, channel_fp);
00304 
00305   /* check the file, close it and return */
00306   if (ferror (channel_fp)) ret_val = FALSE;
00307   else ret_val = TRUE;
00308   fclose (channel_fp);
00309   return ret_val;
00310 
00311 }
00312 
00313 /***************************************************************************
00314  * close_seisan_file
00315  *
00316  * Description: check and close a Seisan event file
00317  *
00318  * Input parameters: none
00319  * Output parameters: none
00320  * Returns: TRUE/FALSE
00321  *
00322  * Comments:
00323  *
00324  ***************************************************************************/
00325 int close_seisan_file (void)
00326 
00327 {
00328 
00329   int channel_count, ret_val, length, end, channel_pos, century_code, count;
00330   long value;
00331   char padded_chan_name [50], padded_network_name [50], sei_filename [MAX_DIR_LEN];
00332   struct Sei_channel_details *channel_ptr;
00333   FILE *sei_fp, *tmp_fp;
00334   time_t ltime;
00335   struct tm conv_time;
00336 
00337   static char string [1500];
00338 
00339 
00340   ret_val = TRUE;
00341 
00342   /* work out the seisan style network name */
00343   length = strlen (network_name);
00344   for (count=0; count<5; count++)
00345   {
00346     padded_network_name [count] = '_';
00347     if (count < length)
00348     {
00349       if (isprint (*(network_name + count)))
00350         padded_network_name [count] = *(network_name + count);
00351     }
00352   }
00353   padded_network_name [5] = '\0';
00354 
00355   /* open the main seisan file */
00356   ltime = (time_t) start_time;
00357   gmtime_ew (&ltime, &conv_time);
00358   sprintf (sei_filename, "%s%c%04d-%02d-%02d-%02d%02d-%02dS.%s_%03d",
00359            output_dir, DIR_DELIM, conv_time.tm_year + 1900, conv_time.tm_mon +1,
00360            conv_time.tm_mday, conv_time.tm_hour, conv_time.tm_min,
00361            conv_time.tm_sec, padded_network_name, n_channels);
00362   sei_fp = fopen (sei_filename, "wb");
00363   if (! sei_fp) return FALSE;
00364 
00365   /* write the first two lines to the event file */
00366   ltime = (time_t) start_time;
00367   gmtime_ew (&ltime, &conv_time);
00368   sprintf (string, " %-29.29s%3d%03d %03d %02d %02d %02d %02d %6.3lf %9.3f           ",
00369                    network_name, n_channels,
00370            conv_time.tm_year, conv_time.tm_yday +1, conv_time.tm_mon +1,
00371            conv_time.tm_mday, conv_time.tm_hour, conv_time.tm_min,
00372            (double) conv_time.tm_sec + (double) ltime - start_time, duration);
00373   write_seisan (sei_fp, 80, string);
00374   memset (string, ' ', 80);
00375   write_seisan (sei_fp, 80, string);
00376 
00377   /* write the station details, minimum 12 records, 3 records per line */
00378   if (n_channels > 30) end = n_channels;
00379   else end = 30;
00380   if (end % 3) end += 3 - (end % 3);
00381   for (channel_count=0, channel_ptr = channels; channel_count<end;
00382        channel_count++, channel_ptr ++)
00383   {
00384     channel_pos = channel_count %3;
00385     if (channel_pos == 0) memset (string, ' ', 80);
00386 
00387     if (channel_count < n_channels)
00388     {
00389       strcpy (padded_chan_name, channel_ptr->chan_name);
00390       while (strlen (padded_chan_name) < 5) strcat (padded_chan_name, " ");
00391       sprintf (&string[channel_pos * 26], " %-4.4s%-4.4s%c%7.2f %8.2f",
00392                padded_chan_name, channel_ptr->chan_type, *(padded_chan_name + 4),
00393                channel_ptr->start_time - start_time, channel_ptr->n_samples / channel_ptr->sample_rate);
00394           string [(channel_pos +1) * 26] = ' ';
00395     }
00396 
00397     if (channel_pos == 2) write_seisan (sei_fp, 80, string);
00398   }
00399 
00400   /* for each channel ... */
00401   for (channel_count=0, channel_ptr = channels; channel_count<n_channels;
00402        channel_count++, channel_ptr ++)
00403   {
00404     /* write the channel header information */
00405     ltime = (time_t) channel_ptr->start_time;
00406     gmtime_ew (&ltime, &conv_time);
00407     century_code = conv_time.tm_year / 100;
00408     sprintf (string, "%-5.5s%-4.4s%d%02d %03d %02d %02d %02d %02d %6.3lf %7.2lf %6ld                          4",
00409              channel_ptr->chan_name, channel_ptr->chan_type, century_code, conv_time.tm_year % 100, 
00410              conv_time.tm_yday +1, conv_time.tm_mon +1,
00411              conv_time.tm_mday, conv_time.tm_hour, conv_time.tm_min,
00412              ((double) conv_time.tm_sec) + (channel_ptr->start_time - (double) ltime),
00413              channel_ptr->sample_rate, channel_ptr->n_samples);
00414     length = strlen (string);
00415     memset (&string[length], ' ', 1040 - length);
00416     write_seisan (sei_fp, 1040, string);
00417 
00418     /* was the channel written OK by an earlier call to 
00419      * add_seisan_channel_data() ?? */
00420     if (channel_ptr->n_written == channel_ptr->n_samples)
00421     {
00422       /* yes - copy the data from the temporary file */
00423       tmp_fp = fopen (channel_ptr->filename, "rb");
00424       if (! tmp_fp) ret_val = FALSE;
00425       else
00426       {
00427         while ((length = fread (string, 1, sizeof (string), tmp_fp)) > 0)
00428           fwrite (string, 1, length, sei_fp);
00429         if (ferror (tmp_fp)) ret_val = FALSE;
00430                 fclose (tmp_fp);
00431       }
00432     }
00433     else
00434     {
00435       /* no - invent missing data for the channel - we create a 1Hz data set
00436        * containing missing data, starting at the event start time */
00437       value = channel_ptr->n_samples * sizeof (long);
00438       if (! output_native_numbers) SwapLong (&value);
00439       fwrite (&value, 4, 1, sei_fp);
00440       value = SEISAN_MISSING_DATA_FLAG;
00441       if (! output_native_numbers) SwapLong (&value);
00442       for (count=0; count<channel_ptr->n_samples; count++)
00443         fwrite (&value, 4, 1, sei_fp);
00444       value = channel_ptr->n_samples * sizeof (long);
00445       if (! output_native_numbers) SwapLong (&value);
00446       fwrite (&value, 4, 1, sei_fp);
00447     }
00448 
00449     /* delete the temporary file for this channel */
00450     if (channel_ptr->filename [0]) remove (channel_ptr->filename);
00451   }
00452 
00453   /* check and close the seisan file */
00454   if (ferror (sei_fp)) ret_val = FALSE;
00455   fclose (sei_fp);
00456 
00457   /* free the channel details */
00458   free (channels);
00459   channels = 0;
00460   n_channels = 0;
00461 
00462   /* remove the temporary directory */
00463   rmdir (tmp_dirname);
00464 
00465   /* if there was an error remove the seisan output file */
00466   if (! ret_val) remove (sei_filename);
00467 
00468   return ret_val;
00469 
00470 }
00471 
00472 /********************************************************************
00473  * write_seisan
00474  *
00475  * Description: write a FORTRAN unformatted record (!)
00476  *
00477  * Input parameters: fp - file to write to
00478  *                   length - the length of the record in bytes
00479  *                   data - the data to write
00480  * Output parameters: none
00481  * Returns: none
00482  *
00483  * Comments:
00484  *
00485  ********************************************************************/
00486 static void write_seisan (FILE *fp, long length, void *data)
00487 
00488 {
00489 
00490   long value;
00491 
00492 
00493   value = length;
00494   if (! output_native_numbers) SwapLong (&value);
00495   fwrite (&value, 4, 1, fp);
00496   fwrite (data, length, 1, fp);
00497   fwrite (&value, 4, 1, fp);
00498 
00499 }
00500 
00501 
00502 /* --------------------------------------------------------------------
00503  * ------ General routines used by non-seisan putaway routines --------
00504  * -------------------------------------------------------------------- */
00505 
00506 /************************************************************************
00507  * pa_find_data
00508  *
00509  * Description: find data in the earthworm data message
00510  *
00511  * Input parameters: trace_req - the trace request structure
00512  *                                       start_time - start time for the data
00513  * Output parameters: data - the data found
00514  * Returns: one of the following codes:
00515  * FD_FOUND_REQUESTED   data was found starting at the requested time
00516  * FD_FOUND_GAP                 there is no data for the requested time, but
00517  *                                              there is data for a later time
00518  * FD_NO_MORE_DATA              the start time is after all data in the
00519  *                                              request structure
00520  * FD_BAD_DATATYPE              found an unreconisable data type code
00521  * FD_CHANGED_SRATE             sample rate changes between snippets
00522  *
00523  * Comments: This is how the data structures work (I think - if this comment
00524  * is wrong, then so is the code):
00525  *
00526  * The TRACE_REQ structure is in ws_clientii.h
00527  *   trace_req->pBuf points to the data buffer
00528  *   trace_req->bufLen is the length of the buffer
00529  *
00530  * The data buffer starts with a TRACE_HEADER structure (trace_buf.h)
00531  *   trace_header->starttime, endtime, samprate hold the obvious things
00532  *   trace_header->datatype[0] holds 's' or 'i' for integer data
00533  *                             or    'f' or 't' for real data
00534  *   trace_header->datatype[1] holds the number of bytes/sample (coded in ASCII)
00535  * Immediately after the TRACE_HEADER structure is the data as an array.
00536  * Another TRACE_HEADER and data block may then follow - trace_req->bufLen
00537  * can be used to determine how many of these header/data pairs there are.
00538  *
00539  ************************************************************************/
00540 int pa_find_data (TRACE_REQ *trace_req, double sample_time,
00541                                   struct Found_data *data)
00542 
00543 {
00544 
00545   int bytes_per_sample, data_type_code;
00546   long n_samples, offset;
00547   double start_time, end_time, sample_rate, first_sample_rate;
00548   char *msg_ptr;
00549   TRACE_HEADER *trace_hdr;
00550 
00551 
00552   /* walk the header/data array, looking for the given time */
00553   first_sample_rate = -1.0;
00554   msg_ptr = trace_req->pBuf;
00555   if (! msg_ptr) return FD_NO_MORE_DATA;
00556   while (msg_ptr < (trace_req->pBuf + trace_req->actLen))
00557   {
00558         /* extract details for this trace */
00559     trace_hdr = (TRACE_HEADER *) msg_ptr;
00560 
00561         /* swap bytes on incoming data if necessary (cjb 2/20/2002) */
00562         if (WaveMsgMakeLocal(trace_hdr) < 0)
00563         {
00564                 logit("e", "pa_find_next: unknown trace data type: %s\n",
00565                         trace_hdr->datatype);
00566                 return( EW_FAILURE );
00567         }
00568 
00569     bytes_per_sample  = atoi (&trace_hdr->datatype[1]);
00570     n_samples = trace_hdr->nsamp;
00571     start_time = trace_hdr->starttime;
00572     end_time = trace_hdr->endtime;
00573     sample_rate = trace_hdr->samprate;
00574 
00575     if (fabs (sample_rate) < 0.000001) sample_rate = 100.0;
00576 
00577     /* check the sample rate */
00578         if (first_sample_rate < 0.0) first_sample_rate = sample_rate;
00579         else if (sample_rate != first_sample_rate) return FD_CHANGED_SRATE;
00580 
00581     /* check the data type codes */
00582     switch (trace_hdr->datatype [0])
00583         {
00584         case 's':
00585     case 'i':
00586           switch (bytes_per_sample)
00587           {
00588           case 2: data_type_code = FD_SHORT_INT; break;
00589           case 4: data_type_code = FD_LONG_INT; break;
00590           default: return FD_BAD_DATATYPE;
00591           }
00592           break;
00593         case 't':
00594         case 'f':
00595           switch (bytes_per_sample)
00596           {
00597           case 4: data_type_code = FD_FLOAT; break;
00598           case 8: data_type_code = FD_DOUBLE; break;
00599           default: return FD_BAD_DATATYPE;
00600           }
00601           break;
00602         default:
00603           return FD_BAD_DATATYPE;
00604         }
00605 
00606     /* work out if this is the trace that we want */
00607         if ((sample_time >= start_time) && (sample_time <= end_time))
00608         {
00609           /* found the data we want - work out what to tell the caller */
00610           offset = (long) (sample_rate * (sample_time - start_time));
00611           data->n_samples = n_samples - offset;
00612           data->data = msg_ptr + sizeof (TRACE_HEADER) + (offset * bytes_per_sample);
00613           data->data_type_code = data_type_code;
00614           data->sample_rate = sample_rate;
00615           data->trace_hdr = trace_hdr;
00616           return FD_FOUND_REQUESTED;
00617         }
00618 
00619     /* is there a gap between the given time and this trace */
00620         if (sample_time < start_time)
00621         {
00622           /* calculate the gap size */
00623           data->n_samples = (long) (sample_rate * (start_time - sample_time));
00624           data->sample_rate = sample_rate;
00625           return FD_FOUND_GAP;
00626         }
00627 
00628         /* set up msg_ptr for the next header/data array */
00629     msg_ptr += sizeof (TRACE_HEADER) + (n_samples * bytes_per_sample);
00630   }
00631 
00632   /* if you get here, then the data was not found */
00633   return FD_NO_MORE_DATA;
00634 
00635 }
00636 
00637 
00638 

Generated on Tue May 6 09:16:09 2003 for Earthworm Libs by doxygen1.3-rc3