ensemble.h

Go to the documentation of this file.

#ifndef _ENSEMBLE_H_
#define _ENSEMBLE_H_
/* This include file violates the normal rule of good style recommended
 * for files in OOP languages where each object type should have it's own 
 * include.  I intentionally violate this here because there are multiple
 * types of useful data seismic data "ensembles".  The code is so parallel
 * between them that I felt it preferable to emphasize this by keeping 
 * it together. At this writing the ensembles are either TimeSeries or
 * ThreeComponentSeismogram object, but the concept is general enough 
 * that a member vector container could hold a lot of other things like
 * complex valued time series, multiwavelet transformed data, or 
 * spectra.  
 */

#include <vector>
#include "coords.h"
#include "perf.h"

#include "pfstream.h"
#include "TimeWindow.h"
#include "TimeSeries.h"
#include "ThreeComponentSeismogram.h"
#include "ComplexTimeSeries.h"
#include "Hypocenter.h"
#include "SeisppKeywords.h"
namespace SEISPP {
class TimeSeriesEnsemble : public Metadata
{
public:  
        vector <TimeSeries> member;

        TimeSeriesEnsemble();
        TimeSeriesEnsemble(int ntsin, int nsampin);
        TimeSeriesEnsemble(DatabaseHandle& dbhi,
                MetadataList& data_mdl,
                        MetadataList& ensemble_mdl,
                                AttributeMap& am);
        TimeSeriesEnsemble(DatabaseHandle& db,
                TimeWindow twin,
                        const string sta_expression="none",
                                const string chan_expression="none",
                bool require_coords=true,
                        bool require_orientation=true,
                                bool require_response=false);
        TimeSeriesEnsemble(Pf_ensemble *pfe);
        TimeSeriesEnsemble(const TimeSeriesEnsemble& tseold);
        TimeSeriesEnsemble(Metadata& md,int nmembers);
        TimeSeriesEnsemble& operator=(const TimeSeriesEnsemble& tseold);
};

class ThreeComponentEnsemble : public Metadata
{
public:
        vector <ThreeComponentSeismogram> member;
        ThreeComponentEnsemble();
        ThreeComponentEnsemble(int nsta, int nsamp);
        ThreeComponentEnsemble(int nsta, int nsamp,
                                MetadataList& mdl);
        ThreeComponentEnsemble(DatabaseHandle& db,
                MetadataList& station_mdl,
                MetadataList& ensemble_mdl,
                 AttributeMap& am);
        ThreeComponentEnsemble(DatabaseHandle& dbhi,
                TimeWindow twin,
                        MetadataList& ensemble_mdl,
                                MetadataList& data_mdl,
                                        AttributeMap& am,
                                                vector<string>chanmap);
        ThreeComponentEnsemble(const ThreeComponentEnsemble& tseold);
        ThreeComponentEnsemble(Metadata& md,int nmembers);
        ThreeComponentEnsemble& operator=(const ThreeComponentEnsemble& tseold);
};
template <class Tmember>
        void remove_trace(Tmember& ensemble, int no)
{
    if(no >= ensemble.member.size()) return;
    int i;
    typedef typename std::vector<Tmember> vector_type;
    typename vector_type::iterator it;

    it=ensemble.member.begin();
    for(i=0;i<no;++i) ++it;
    ensemble.member.erase(it);
}
/* Current implementations of PeakAmplitude */

double PeakAmplitude(TimeSeries *p);
double PeakAmplitude(ThreeComponentSeismogram *p);
double PeakAmplitude(ComplexTimeSeries *p);

template <class Tensemble,class Tmember> 
        void MeasureEnsemblePeakAmplitudes(Tensemble& d,
                string scale_attribute)
{
        int i;
        double amplitude;
        for(i=0;i<d.member.size();++i)
        {
                Tmember *dp=&(d.member[i]);
                if(dp->live)
                {
                        amplitude=PeakAmplitude(dp);
                        dp->put(scale_attribute,amplitude);
                }
        }
}

template <class Tensemble,class Tmember> 
        void ScaleEnsemble(Tensemble& d,
                string scale_attribute,
                        bool invert=false)
{
        int i;
        double amplitude;
        for(i=0;i<d.member.size();++i)
        {
                Tmember *dp=&(d.member[i]);
                if(dp->live)
                {
                        try {
                                amplitude=dp->get_double(scale_attribute);
                                // invalid amplitudes are signaled as negative
                                if(amplitude>0)
                                {
                                        if(invert)
                                                ScaleMember(dp,1.0/amplitude);
                                        else
                                                ScaleMember(dp,amplitude);
                                }
                        }
                        catch (MetadataGetError mde)
                        {
                                cerr << "Warning (ScaleEnsemble):  "
                                        <<"scaling of data failed"<<endl;
                                mde.log_error();
                                cerr << "Continuing with data left unscaled"
                                        <<endl;
                        }
                }
        }
}
/* Current implementations of ScaleMember. */

void ScaleMember(TimeSeries *p,double scale);
void ScaleMember(ThreeComponentSeismogram *p,double scale);
void ScaleMember(ComplexTimeSeries *p,double scale);

/* Companion to above */

template <class Tensemble> void ScaleCalib(Tensemble d,
        string calib_attribute, 
                string scale_attribute, 
                        bool invert=false)
{
        double calib;
        double scale;
        for(int i=0;i<d.member.size();++i)
        {
                if(d.member[i].live)
                {
                        try {
                                calib=d.member[i].get_double(calib_attribute);
                        } catch (MetadataGetError mde) 
                        {
                                calib=1.0;
                        }
                                
                        try {
                                scale=d.member[i].get_double(scale_attribute);
                        } catch (MetadataGetError mde) 
                        {
                                scale=1.0;
                        }
                        if(invert) scale=1.0/scale;
                        d.member[i].put(calib_attribute,scale*calib);
                }
        }
}
template <class Te, class Ta> 
        void InitializeEnsembleAttribute(Te& d, string name,
                Ta val)
{
        try {
                int i;
                for(i=0;i<d.member.size();++i)
                        d.member[i].put(name,val);
        } catch(...){throw;};
}
template <class Tensemble> void LoadEventArrivals(Tensemble& d, 
                DatabaseHandle& dbi,
                        string phase,
                                string predarrkeyword,
                                        string atkeyword,
                                                double tpad,
                                                        double nullvalue)
{
        DatascopeHandle dbh=dynamic_cast<DatascopeHandle&> (dbi);
        DatascopeHandle dbhss(dbh);   /* working copy */
        /* First scan the ensemble for the range of predicted arrival times.
        Initialize with negative values to provide error check below */
        TimeWindow atrange(-2.0,-1.0);
        double testval;
        bool firstpass(true);
        int i;
        for(i=0;i<d.member.size();++i)
        {
                if(d.member[i].is_attribute_set(predarrkeyword))
                {
                        testval=d.member[i].get_double(predarrkeyword);
                        if(firstpass)
                        {
                                atrange.start=testval;
                                atrange.end=testval;
                                firstpass=false;
                        }
                        else
                        {
                                if(atrange.start>testval)
                                        atrange.start=testval;
                                if(atrange.end<testval)
                                        atrange.end=testval;
                        }
                }
        }
        if(atrange.start<0.0)
                throw SeisppError(string("LoadEventArrivals:  no predicted arrivals are defined")
                        + string(" in ensemble passed to this procedure\n")
                        + string("Coding error or problem in the way metadata were loaded.") );
        atrange.start -= tpad;
        atrange.end += tpad;
        char ssexpression[256];
        sprintf(ssexpression,"(time >= %lf && time<=%lf && iphase=~/%s/)",
                atrange.start,atrange.end,phase.c_str());
        dbhss.subset(string(ssexpression));
        /* We always initialize the full ensemble to null values.  This is needed
        in case an exception is thrown by one of the get routines */
        for(i=0;i<d.member.size();++i)
                        d.member[i].put(atkeyword,nullvalue);
        if(dbhss.number_tuples()>0)
        {
                try {
                        map<string,double> atimes;
                        map<string,double>::iterator atptr;
                        double t;
                        string sta;
                        dbhss.rewind();
                        for(i=0;i<dbhss.number_tuples();++i,++dbhss)
                        {
                                sta=dbhss.get_string("sta");
                                t=dbhss.get_double("arrival.time");
                                atimes[sta]=t;
                        }
                        for(i=0;i<d.member.size();++i)
                        {
                                sta=d.member[i].get_string("sta");
                                atptr=atimes.find(sta);
                                if(atptr!=atimes.end())
                                {
                                        d.member[i].put(atkeyword,atptr->second);
                                }
                        }
                }
                catch (...) 
                {
                        throw SeisppError(string("LoadEventArrivals:  ")
                        + string("Problems arrivals.  Not all arrivals were loaded") );
                }
        }
        /*WARNING:  the handle really should automatically do this, but in the 
        current implementation it does not.  This will produce a memory leak if
        we don't do this.  If changed this next line must be removed */
        //dbfree(dbhss.db);
}
template <class Tensemble> int LoadPredictedArrivalTimes(Tensemble& d,
        string phase,
                string predarr_keyword=predicted_time_key,
                        string ttmethod="tttaup",
                                string ttmodel="iasp91",
                                        bool verbose=false)
{
        int i,nmembers;
        double slat,slon,sz,stime;
        double rlat,rlon,relev;
        double phasetime;
        int nfailures(0);

        nmembers=d.member.size();
        for(i=0;i<nmembers;++i)
        {
                if(d.member[i].live)
                {
                        try{
                                slat=d.member[i].get_double("source_lat");
                                slon=d.member[i].get_double("source_lon");
                                sz=d.member[i].get_double("source_depth");
                                stime=d.member[i].get_double("source_time");
                                rlat=d.member[i].get_double("sta_lat");
                                rlon=d.member[i].get_double("sta_lon");
                                relev=d.member[i].get_double("sta_elev");
                                slat=rad(slat); 
                                slon=rad(slon);
                                rlat=rad(rlat);
                                rlon=rad(rlon);
                                Hypocenter h(slat,slon,sz,stime,
                                        ttmethod,ttmodel);
                                phasetime=h.phasetime(rlat,rlon,relev,phase);
                                phasetime+=stime;
                        }
                        catch (MetadataGetError mderr)
                        {
                                phasetime=0.0;
                                ++nfailures;
                                if(verbose)
                                {
                                        cerr << "LoadPredictedArrivalTimes:  "
                                         << "get failed on attribute name="
                                        << mderr.name
                                        <<" for ensemble member number "
                                        << i <<endl
                                        << "Arrival time set to 0.0"<<endl;
                                }
                        }
                        catch (SeisppError serr)
                        {
                                phasetime=0.0;
                                ++nfailures;
                                if(verbose)
                                {
                                        cerr << "LoadPredictedArrivalTimes:  "
                                         << "travel time calculator error "
                                         << "for member="<<i <<endl
                                         << "SeisppError message:"<<endl;
                                        serr.log_error();
                                        cerr << "Arrival time set to zero"<<endl;
                                }
                        }

                }
                else
                {
                        phasetime=0.0;
                }
                d.member[i].put(predarr_keyword,phasetime);
        }
        return(nfailures);
}

auto_ptr<TimeSeriesEnsemble> ExtractComponent(ThreeComponentEnsemble& tcs,int component);
} // End SEISPP namespace declaration
#endif

---