SeismicDataComp.cpp

Go to the documentation of this file.

#include "SeismicDataComp.h"
#include <algorithm>
#include <fstream>
#include <cassert>
#include <boost/numeric/conversion/cast.hpp>
#include <boost/filesystem/operations.hpp>
#include <boost/cstdint.hpp>
#include "CFatalException.h"

using namespace std;

SeismicDataComp::SeismicDataComp():
  stla(0), stlo(0), stel(0), stdp(0), evla(0), evlo(0), evel(0), evdp(0),
      mag(0), dist(0), az(0), baz(0), gcarc(0), b(0), dataformat(sac)
  {
  }

SeismicDataComp::SeismicDataComp(const std::string &filename, tseismicdataformat format):
  stla(0), stlo(0), stel(0), stdp(0), evla(0), evlo(0), evel(0), evdp(0),
        mag(0), dist(0), az(0), baz(0), gcarc(0), b(0), dataformat(sac)
  {
    GetData(filename,format);
  }

SeismicDataComp::~SeismicDataComp()
  {
  }

void SeismicDataComp::CopyHeader(const SeismicDataComp& source)
  {
    this->stla = source.stla;
    this->stlo= source.stlo;
    this->stel= source.stel;
    this->stdp= source.stdp;
    this->evla= source.evla;
    this->evlo= source.evlo;
    this->evel= source.evel;
    this->evdp= source.evdp;
    this->mag= source.mag;
    this->dist= source.dist;
    this->az= source.az;
    this->baz= source.baz;
    this->gcarc= source.gcarc;
    this->b = source.b;
  }

SeismicDataComp& SeismicDataComp::operator= (const SeismicDataComp& source)
  {
    if (this == &source) return *this;
    CopyHeader(source);
    this->dataformat = source.dataformat;
    TimeSeriesComponent::operator=(source);
    return *this;
  }

void SeismicDataComp::ReadAscii(const std::string &filename)
  {
    ifstream infile(filename.c_str());
    double time, amp, oldtime;

    infile >> oldtime >> amp;
    if (infile.good())
      GetData().push_back(amp);
    else
      throw CFatalException("Cannot read from file: "+filename);
    b = oldtime;
    infile >> time >> amp;
    assert(time - oldtime > 0.0);
    GetData().push_back(amp);
    SetDt(time - oldtime);
    if (!GetData().empty() && infile.good())
      GetData().clear();
    while (infile.good())
      {
        infile >> time >> amp;
        if (infile.good())
          GetData().push_back(amp);
      }
  }

void SeismicDataComp::ReadSac(const std::string &filename)
  {
    boost::int32_t ihd[40]; // These defintions are taken from the SAC to Helm Code
    float fhd[70];
    char chd[8][24];
    ifstream infile(filename.c_str());
    infile.read((char *)(&fhd), 70*4); /*Read Sac Float Field*/
    if (infile.fail()) //if reading fails
      {
        throw CFatalException("Cannot read from file: "+filename);
      }
    else
      {

        infile.read((char *)(&ihd), 40*4); /*Read Sac Int   Field*/
        infile.read((char *)(&chd), 24*8); /*Read Sac Char. Field*/

        stla =fhd[31];
        stlo =fhd[32];
        stel =fhd[33];
        stdp =fhd[34];
        evla =fhd[35];
        evlo =fhd[36];
        evel =fhd[37];
        evdp =fhd[38];
        mag =fhd[39];
        dist =fhd[50];
        az =fhd[51];
        baz =fhd[52];
        gcarc =fhd[53];
        boost::int32_t npts=ihd[9];
        assert(fhd[0] > 0.0);
        SetDt(fhd[0]);
        b=fhd[5];
        assert(npts > 0);
        float *temp = new float[npts]; //allocate memory for data
        infile.read((char *)(temp), npts*sizeof(float)); //read data part of SAC file
        if (!GetData().empty())
          GetData().clear();
        copy(temp, temp+npts, back_inserter(GetData()));

        //Data.assign(npts,0);
        //for (int i = 0; i < npts; ++i)
        //      Data.at(i) = temp[i];
        delete []temp;
        //ShiftStart(static_cast<int>(ceil(b/dt)));
      }
  }

void SeismicDataComp::ReadHeaderAscii(const std::string &filename) 
  {
    ifstream infile(filename.c_str());
    double amp;
    char currentline[255];
    const int headerlines = 12;

    infile.getline(currentline, 255, ':');
    double delta;
    infile >> delta;
    SetDt(delta);
    for (int i = 0; i < headerlines; ++i)
      infile.getline(currentline, 255);
    if (!GetData().empty() && infile.good())
      GetData().clear();
    while (infile.good())
      {
        infile >> amp;
        if (infile.good())
          GetData().push_back(amp);
      }
    SetName(filename);
  }

void SeismicDataComp::ReadSKS(const std::string &filename)
  {
    ifstream infile(filename.c_str(), ios::binary);
    int intdummy;
    float floatdummy;
    int npts = 0;

    if (!infile)
      {
        cerr << "Error opening file: "<< filename << endl;
      }
    infile.read((char *)&intdummy, sizeof(int));
    infile.read((char *)&floatdummy, sizeof(float));
    infile.read((char *)&intdummy, sizeof(int));

    infile.read((char *)&intdummy, sizeof(int));
    for (int i = 0; i < 5; ++i)
      infile.read((char *)&intdummy, sizeof(int));
    infile.read((char *)&intdummy, sizeof(int));

    infile.read((char *)&intdummy, sizeof(int));
    for (int i = 0; i < 6; ++i)
      infile.read((char *)&intdummy, sizeof(int));
    infile.read((char *)&intdummy, sizeof(int));

    infile.read((char *)&intdummy, sizeof(int));
    infile.read((char *)&npts, sizeof(int));
    infile.read((char *)&intdummy, sizeof(int));
    infile.read((char *)&floatdummy, sizeof(float));
    SetDt(floatdummy);
    infile.read((char *)&intdummy, sizeof(int));

    infile.read((char *)&intdummy, sizeof(int));
    infile.read((char *)&intdummy, sizeof(int));
    infile.read((char *)&intdummy, sizeof(int));
    infile.read((char *)&intdummy, sizeof(int));

    GetData().assign(npts, 0);
    infile.read((char *)&intdummy, sizeof(int));
    infile.read((char *)&floatdummy, sizeof(float));
    for (int i = 0; i < npts; ++i)
      {
        infile.read((char *)&floatdummy, sizeof(float));
        GetData().at(i) = floatdummy;
      }
    infile.read((char *)&intdummy, sizeof(int));

  }

int SeismicDataComp::GetData(const std::string &filename,
    tseismicdataformat format)
  {
    if (!boost::filesystem::exists(filename))
      throw CFatalException("File does not exist : "+filename);
    switch (format)
      {
      case sac:
        ReadSac(filename);
        break;
      case sks:
        ReadSKS(filename);
        break;
      case head:
        ReadHeaderAscii(filename);
        break;
      case ascii:
        ReadAscii(filename);
        break;
      case unknownseis:
        cerr << "Unknown dataformat ! Cannot read "<< filename << endl;
        return -100;
        break;
      default:
        cerr << "Unknown dataformat ! Cannot read "<< filename << endl;
        return -100;
        break;
      }
    SetName(filename);
    return 0;
  }

int SeismicDataComp::WriteAsSac(const std::string &filename) const
  {
    int npts;
    int ihd[40]; // These definitions are taken from the SAC to Helm Code
    float fhd[70];
    char chd[8][24];
    ofstream outfile;
    float *temp= NULL;
    int i;

    outfile.open(filename.c_str());

    npts = GetData().size();

    temp=new float[npts]; //allocate memory for data

    copy(GetData().begin(), GetData().end(), temp);
    for (i=0; i<40; i++)
      ihd[i]=-12345;
    for (i=0; i<70; i++)
      fhd[i]=-12345.00;
    for (i=0; i<8; i++)
      sprintf(chd[i], "-12345  -12345  -12345  ");
    /*Set essential sac parameters*/
    ihd[35]=1; /*Sets file to evenly spaced*/
    ihd[15]=1; /*Sets file type to Timeseries*/
    ihd[6]=6; /*Variable Name Internal */

    fhd[31] = stla;
    fhd[32] = stlo;
    fhd[33] = stel;
    fhd[34] = stdp;
    fhd[35] = evla;
    fhd[36] = evlo;
    fhd[37] = evel;
    fhd[38] = evdp;
    fhd[39] = mag;
    fhd[50] = dist;
    fhd[51] = az;
    fhd[52] =baz;
    fhd[53] = gcarc;

    ihd[9] = npts;
    fhd[0] = GetDt();
    fhd[5] = b;

    outfile.write((char *)(&fhd), 70*4); /*Read Sac Float Field*/
    outfile.write((char *)(&ihd), 40*4); /*Read Sac Int   Field*/
    outfile.write((char *)(&chd), 24*8); /*Read Sac Char. Field*/
    outfile.write((char *)(temp), npts*sizeof(float)); //write data part of SAC file
    outfile.close();
    delete []temp;
    return (0);
  }

int SeismicDataComp::WriteAsHeaderAscii(const std::string &filename) const
  {
    return -100;
  }

int SeismicDataComp::WriteAsSKS(const std::string &filename) const
  {
    ofstream outfile(filename.c_str(), ios::binary);
    int intdummy = 0;
    float floatdummy = 0;
    const int npts = GetData().size();

    outfile.write((char *)&intdummy, sizeof(int));
    outfile.write((char *)&floatdummy, sizeof(float));
    outfile.write((char *)&intdummy, sizeof(int));

    outfile.write((char *)&intdummy, sizeof(int));
    for (int i = 0; i < 5; ++i)
      outfile.write((char *)&intdummy, sizeof(int));
    outfile.write((char *)&intdummy, sizeof(int));

    outfile.write((char *)&intdummy, sizeof(int));
    for (int i = 0; i < 6; ++i)
      outfile.write((char *)&intdummy, sizeof(int));
    outfile.write((char *)&intdummy, sizeof(int));

    outfile.write((char *)&intdummy, sizeof(int));
    outfile.write((char *)&npts, sizeof(int));
    outfile.write((char *)&intdummy, sizeof(int));
    floatdummy = GetDt();
    outfile.write((char *)&floatdummy, sizeof(float));
    outfile.write((char *)&intdummy, sizeof(int));

    outfile.write((char *)&intdummy, sizeof(int));
    outfile.write((char *)&intdummy, sizeof(int));
    outfile.write((char *)&intdummy, sizeof(int));
    outfile.write((char *)&intdummy, sizeof(int));

    outfile.write((char *)&intdummy, sizeof(int));
    outfile.write((char *)&floatdummy, sizeof(float));

    float *temp=new float[npts]; //allocate memory for data
    copy(GetData().begin(), GetData().end(), temp);
    outfile.write((char *)(temp), npts*sizeof(float));
    outfile.write((char *)&intdummy, sizeof(int));
    delete []temp;
    return 0;
  }

int SeismicDataComp::WriteAsAscii(const std::string &filename) const
  {
    ofstream outfile(filename.c_str());

    for (unsigned int i = 0; i < GetData().size(); ++i)
      {
        outfile << b + i * GetDt() << " " << GetData().at(i) << endl;
      }
    return 0;
  }

int SeismicDataComp::WriteBack() const
  {
    switch (dataformat)
      {
      case sac:
        WriteAsSac(GetName());
        break;
      case sks:
        WriteAsSKS(GetName());
        break;
      case head:
        WriteAsHeaderAscii(GetName());
        break;
      case ascii:
        WriteAsAscii(GetName());
        break;
      case unknownseis:
        cerr << "Unknown dataformat ! Cannot write " << GetName() << endl;
        return -100;
        break;
      default:
        cerr << "Unknown dataformat ! Cannot write !"<< GetName() << endl;
        return -100;
        break;
      }
    return 0;
  }

---