process_data.c

解吸SEED格式的源代码

/*===========================================================================*/
/* SEED reader     |             process_data              |    subprocedure */
/*===========================================================================*/
/*
	Name:	process_data
	Purpose:recover seismic data from SEED format
	Usage:	void process_data ();
			process_data ();
	Input:	none (seismic data are taken from globally-available structure)
	Output:	none (a subprocedure writes seismogram files)
	Externals:input.data - the data portion of a SEED logical record
		current_station - pointer to current station info in tables
		current_channel - pointer to current channel info in tables
		data_hdr - SEED fixed-length data header
		last_time - ending time of previous record
		last_sample_rate - sample rate of previous record
		last_nsamples - number of samples in previous record
		seismic_data - the current seismogram
		seismic_data_ptr - a pointer into "seismic_data"
		decode - data decoder key (1st 5 letters of compression name)
		numitems - number of seismograms wanted
		item - list of starting record numbers of desired seismograms
		start_record - starting logical record number of this seismogram
	Warnings:station name not found in tables
			channel name not found in tables
			data decomression type not found in tables
			unknown data compression scheme encountered
	Errors:	none
	Called b:main program
	Calls to:output_data - writes out a seismic data file
			decode_asro - decode ASRO data
			decode_cdsn - decode CDSN data
			decode_dwwssn - decode DWWSSN data
			decode_sro - decode SRO data
			decode_steim - decode Steim-compressed data
			decode_32bit - decode UCSD 32-bit integer data
			timeadd - add two times (yyyy,ddd,hh:mm:ss.ffff)
			timecmp - compare two times (yyyy,ddd,hh:mm:ss.ffff)
			convert_seedhdr - construct minimum data header
	Algorithm:recover the station, channel, and start time from the current
		data block; figure out the ending time and sample rate.
		If the current data block is station-, channel-, and time-
		continuous with the previous data block, append the data
		to the current seismogram and update the minimum data header
		appropriately.  If the current data block is not continuous
		with the previous data block, write out the previous seismogram
		and start a new one using this block's data.  Checks are only
		performed if the current data block is part of a desired 
		seismogram.
	Notes:	none
	Problems:can't derive an uncompressor from the abbreviation dictionary
		description, but they may be faulty anyway (esp. CDSN)
	References:Halbert et al, 1988; see main routine
	Language:C, hopefully ANSI standard
	Author:	Dennis O'Neill
	Revisions:07/15/88  Dennis O'Neill  Initial preliminary release 0.9
		11/11/88  Dennis O'Neill  removed extra argument from a printf
		11/11/88  Dennis O'Neill  added check for possible divide-by-0
		                          condition in time-cntity check
		11/11/88  Dennis O'Neill  corrected neglected "-" sign in
		                          handling sample rate and multiplier
		11/11/88  Dennis O'Neill  installed smart byte-swapping
		11/15/88  Dennis O'Neill  changed handling of long seisgms
		11/16/88  Dennis O'Neill  moved "output", "output_overflow"
		                          vbl to outside of proc body
		11/21/88  Dennis O'Neill  Production release 1.0
		01/24/89  Dennis O'Neill  fixed some fprintfs
		02/06/89  Dennis O'Neill  added decode_32bit call
		02/10/89  Dennis O'Neill  fixed the statement
		data_hdr->number_time_corrections += (long int) swap_4byte;
		data_hdr->number_time_corrections += (long int) temp_4byte;
		05/01/92  Allen Nnace     fixed the statement
		data_hdr->number_time_corrections += (long int) temp_4byte; to
		data_hdr->number_time_corrections = (long int) temp_4byte;
		05/01/92  Allen Nance     display invalid data compression type
		05/05/92  Allen Nnace     change continuity check from
		(timetol (newtime, this_time, data_hdr->nsamples) == 0) to
		(timetol (newtime, this_time, last_nsamples) == 0)
		01/17/95  CL - changed decoding name from 5 bytes MAX to whatever

*/

/* #include <sunmath.h> */

#include "rdseed.h"

#define TIME_PRECISION 10000		/* number of divisions in a */
					/* second in SEED times */

extern struct stn_list *stn_listhead;

struct type50 *get_station_rec();
struct type52 *get_channel_rec();

struct data_blk_1000 *scan_for_blk_1000();

void do_time_correction();
void dump_seismic_buffer();
char determine_orient_code();
char *code_to_english();
void swap_fsdh();

extern struct mini_data_hdr mini_data_hdr;

extern int ignore_net_codes;	/* set in "main" -rdseed.c */
extern int read_summary_flag;	/* set in rdseed.c */
extern int strip_flag;		/* set in rdseed.c */
extern int mini_flag;		/* set in rdseed.c */



/* define these here so they keep their values between calls */
int output;			/* output-wanted flag */
int output_overflow;		/* for long seismograms */


struct data_blk_1000 *p;

static char network_code[3];

static int chn_log;

/* ------------------------------------------------------------------------ */
struct data_blk_1000 *scan_for_blk_1000();
struct data_blk_2000 *scan_for_blk_2000();
/* ------------------------------------------------------------------------ */


void process_data ()

{

	int i;									
	char *input_data_ptr;			/* ptr to input data */
	struct input_data_hdr *input_data_hdr;	/* fixed data header */
	struct type52 *old_current_channel, 
				  *new_current_channel, 
				  *p_channel;	/* pointer save locations */

	struct type50 *old_current_station, 
				  *new_current_station, 
				  *p_station;		/* pointer save locations */

	int continuous;		/* stn,chnl,time cts flag */
	double duration;			/* time duration of record */
	struct time newtime;		/* new ending time of s'gm */
	int sample_rate_multiplier;	/* temp used in calc smpl rt */
	short int temp_2byte;		/* temp for byte-swapping */
	int temp_4byte;			/* temp for byte-swapping */
	unsigned short int temp_u2byte; /* temp for byte-swapping */
	struct input_time temptime;	/* temp for byte-swapping */
	int next_index;			/* next index into s'gm array */

	/* "this" => current data blk */
	char this_station[5+1];		/* name of this station */
	char this_location[2+1];	/* location code this station */
	char this_channel[3+1];		/* name of this channel */
	struct time this_time;		/* start time of this block */
	int this_nsamples;		/* number of samples this blk */
	double this_sample_rate;		/* temp used to store smpl rt */

	struct type30 *this_type30;	/* index into SEED tables */

	chn_log = FALSE;


/*                 +=======================================+                 */
/*=================|   read the fixed part of data header  |=================*/
/*                 +=======================================+                 */

	/* point to beginning of data */
	input_data_ptr = lrecord_ptr + 8;

	/* recover the fixed data header from the input record */
	input_data_hdr = (struct input_data_hdr *) input_data_ptr;

	/* check for D, Q, R and compare to q_flag */
	if (q_flag != 'E')
		switch (*(lrecord_ptr + 6))
		{
			case 'D' :
				if (q_flag != 'D')
					return;	
				break;
			case 'R' :
				if (q_flag != 'R')
					return;
				break;
			case 'Q' :
				if (q_flag != 'Q')
					return;
				break;

			default:
				fprintf(stderr, "Bad Data Header Indicator byte. Byte=%c\n", *(lrecord_ptr + 7));

		}

/*                 +=======================================+                 */
/*=================| obtain particulars for this data blk  |=================*/
/*                 +=======================================+                 */

	/* recover station name */
	for (i = 0; i < 5; i++) 
		this_station[i] = input_data_hdr->station[i];
	this_station[5] = '\0';

	for (i = 4; this_station[i] == ' '; i--) 
	{
		this_station[i] = '\0';
		if (i == 0)
			break;
	}

	/* recover location code */
	for (i = 0; i < 2; i++) 
		this_location[i] = input_data_hdr->location[i];
	this_location[2] = '\0';

	for (i = 1; this_location[i] == ' '; i--)
	{
		this_location[i] = '\0';
		if (i == 0) 
                        break; 

	}


	/* recover channel name */
	for (i = 0; i < 3; i++) 
		this_channel[i] = input_data_hdr->channel[i];
	this_channel[3] = '\0';


	for (i = 2; this_channel[i] == ' '; i--) 
	{
		this_channel[i] = '\0';
		if (i == 0) 
                        break; 

	}

	/* scan for a LOG/AT/SOH record */
	if (strcasecmp(this_channel, "LOG") == 0)
		chn_log = 1;

	/* recover the network code */
	if (type10.version >= 2.3)
	{
		char *p;

		memcpy(network_code, input_data_hdr->network, 2);
		network_code[3] = 0;
			
	}
	else 
		strcpy(network_code, "");

/*                 +=======================================+                 */
/*=================| find info about this station/channel  |=================*/
/*                 +=======================================+                 */
	old_current_channel = current_channel; /* old current channel pointer */
	old_current_station = current_station; /* old current channel pointer */

	if (chn_log)
	{
		current_station = get_station_rec(this_station, 
						network_code, NULL);

                if (current_station == NULL)
                {
                        fprintf(stderr, "WARNING (process_data):  ");
                        fprintf (stderr, 
	"station %s/%s not found in station for network %2.2s.\n", 
				this_station, 
				type10.version >= 2.3 ? network_code : "N/A");
 
                        fprintf (stderr, "\tData Record Skipped.\n");
 
                        /* toggle = FALSE; */
 
                        current_station = old_current_station;
                        current_channel = old_current_channel;
 
                        return;
 
                }

	}
	else
	{
		get_stn_chn_rec(this_station, 
				this_channel, 
				network_code, 
				this_location,
				NULL);

		if (current_station == NULL || current_channel == NULL)
		{
			fprintf(stderr, "WARNING (process_data):  ");

			fprintf (stderr, "station/channel %s/%s not found in station/channel tables for network %2.2s, location code:%2.2s.\n", this_station, this_channel, type10.version >= 2.3 ? network_code : "N/A", this_location);

 
        		fprintf (stderr, "\tData Record Skipped.\n");
 
                	/* toggle = FALSE; */

                	current_station = old_current_station;
			current_channel = old_current_channel;
 
                	return;

		}

	}

   	/* check if byteswapping will be needed to recover data */
	/* hopefully, rdseed will not be around by 2010 */

	if (input_data_hdr->time.year < 1950 || 
			input_data_hdr->time.year > 2010)
		swap_fsdh(&input_data_hdr, &input_data_ptr);

	/* check for sanity */
	if (input_data_hdr->time.year < 1950 || 
		input_data_hdr->time.year > 2010)
	{
		fprintf(stderr, "ERROR - process_data(): Unknown word order for station %s, channel %s, network %s, location %s\n",
 				this_station, 
				this_channel, 	
 				network_code, 
				this_location);

		fprintf(stderr, "Skipping data record.\n");
		return;

	}

	/* recover this block's sample rate */
	this_sample_rate = input_data_hdr->sample_rate;

	sample_rate_multiplier = input_data_hdr->sample_rate_multiplier;

	if (this_sample_rate < 0) 
		this_sample_rate = 1 / (-this_sample_rate);

	if (sample_rate_multiplier > 0) 
		this_sample_rate = this_sample_rate * sample_rate_multiplier;
	else 
	if (sample_rate_multiplier < 0) 
		this_sample_rate = this_sample_rate / (-sample_rate_multiplier);

	parse_type100(input_data_ptr, &this_sample_rate);

	/* check for sanity */
	if (this_sample_rate < .01)
		return;

	if (this_sample_rate > 500)
		return;

	/* get start time of this data block */
	this_time.year = input_data_hdr->time.year;

	this_time.day = input_data_hdr->time.day;

	this_time.hour    = input_data_hdr->time.hour;
	this_time.minute  = input_data_hdr->time.minute;
	this_time.second  = input_data_hdr->time.second;

	this_time.fracsec = input_data_hdr->time.fracsec;

	/* Do time correction */
	if ((input_data_hdr->activity_flags & 0x02) == 0)
	{
		this_time.fracsec += 
			input_data_hdr->number_time_corrections;

		do_time_correction(&this_time);
	}
/* printf("time=%d,%d,%d:%d:%d.%d\n", this_time.year, this_time.day, this_time.hour, this_time.minute, this_time.second, this_time.fracsec);
*/

	atc_correct(&this_time, this_station,
                                this_channel,
                                network_code,
                                this_location,
				this_sample_rate);

	current_channel = NULL;

        if (chn_log)
        {
                current_station = get_station_rec(this_station,
                                                network_code, &this_time);

                if (current_station == NULL)
                {
                
			fprintf(stderr, "WARNING (process_data):  ");
                	fprintf(stderr, "station %s response effective time not found in station table.\n", this_station);

                	fprintf(stderr, "\tTrying again ignoring effective times.\n");

                	current_station = get_station_rec(this_station,
                                        		  network_code,
                                        		  NULL);

                	if (current_station == NULL)
                	{
                        	fprintf(stderr, "ERROR - unable to locate any station record for station:%s, network:%s\n\nSkipping data record!\n",
                                                this_station,
                                                type10.version >= 2.3 ?
                                                        network_code:
                                                        "N/A");
 
                        	current_station = old_current_station;
 
                        	return;
                	}	/* current_station == NULL */
		}

        }
        else
        {   
                get_stn_chn_rec(this_station,
                                this_channel,
                                network_code, 
				this_location,
				&this_time);

                if (current_station == NULL || current_channel == NULL)
                {
                	fprintf (stderr, "WARNING (process_data):  ");
                	fprintf (stderr, "station %s response effective time not found in station table.\n", this_station);
 
                	fprintf (stderr, "\tTrying again ignoring effective times.\n");
 
                	get_stn_chn_rec(this_station,
					  this_channel,
                               		  network_code,
					  this_location,
                               		  NULL);
 
                	if (current_station == NULL || current_channel == NULL)
                	{
                        	fprintf(stderr, "ERROR - unable to locate any station/channel record for station:%s, channel: %s, network:%s\n\nSkipping data record!\n", 
				   this_station,
				   this_channel,
                                   type10.version >= 2.3 ?
                                   	network_code:
                                       		"N/A");
 
                        	current_station = old_current_station;
				current_channel = old_current_channel;
 
                        	return;
 
                	}
			else
				fprintf(stderr, "Found stn/chn record.\n");

                }
	} 

	if (input_data_hdr->number_blockettes > 0)
        {
		if (scan_for_blk_2000(input_data_ptr +
                                        input_data_hdr->bofb - 8,
					input_data_ptr - 8))
		{
			process_opaque(input_data_ptr);
			return;
		}

	}

	p = 0;

	if (input_data_hdr->number_blockettes > 0)
	{
		p = scan_for_blk_1000(input_data_ptr + 
					/* points to blockettes */
					input_data_hdr->bofb - 8,	
					/* points to lrec base 00000D... */
					input_data_ptr - 8);  
	}

	if (p)
		byteswap = find_wordorder(p);
	else
		byteswap = find_wordorder (NULL);


	/* scan for ATC channels for this stn/chn */

	new_current_station = current_station;	/* new current channel */
	current_station = old_current_station;	/* restore current channel pointer */
	new_current_channel = current_channel;	/* new current channel */
	current_channel = old_current_channel;	/* restore current channel pointer */

	
/*                 +=======================================+                 */
/*=================| Is block within selected time spans   |=================*/
/*                 +=======================================+                 */

	/* if it is a LOG record, check the start time, fake an end time
	 * of infinity
	 */
	if (chn_log)
	{
		struct time fake_time;

		memset((char *)&fake_time, 0, sizeof(fake_time));

		if (!chk_time(this_time, fake_time))
			return;

		if (data_hdr->nsamples > 0)
			dump_seismic_buffer();

		if (Seed_flag)
        	{
                	struct data_hdr d_hdr; 
         
			memset((char *)&d_hdr, 0, sizeof(d_hdr));
 
                	d_hdr.station = this_station; 
                	d_hdr.channel = this_channel; 
			d_hdr.location = this_location;

                	strncpy(d_hdr.network, 
				 network_code, strlen(network_code));
 
                	d_hdr.time = this_time;

			d_hdr.bod = 48;
 
                	output_seed_data_file(input_data_ptr); 
 
                	update_type74(&d_hdr);