plotrecfunc.cpp

Go to the documentation of this file.

#include "SeismicStationList.h"
#include <iostream>
#include <string>
#include <vector>
#include <functional>
#include <boost/bind.hpp>
#include <boost/function.hpp>
#include <fstream>
#include "Util.h"
#include "netcdfcpp.h"
#include "statutils.h"

using namespace std;
using namespace gplib;

void WriteToNetCDF(string filename, SeismicStationList &Stations,
    boost::function<double(SeismicDataComp *)> f)
  {
    const size_t nsamples = Stations.GetList().front()->GetData().size();
    const size_t nstations = Stations.GetList().size();

    sort(Stations.GetList().begin(), Stations.GetList().end(), boost::bind(
        less<double> (), boost::bind(f, boost::bind(&boost::shared_ptr<
            SeismicDataComp>::get, _1)), boost::bind(f, boost::bind(
            &boost::shared_ptr<SeismicDataComp>::get, _2))));
    size_t nCompatStations = 0;
    for (size_t i = 0; i < nstations; ++i)
      {
        if (nsamples == Stations.GetList().at(i)->GetData().size())
          {
            ++nCompatStations;
          }
      }

    NcFile combrescdf((filename + ".nc").c_str(), NcFile::Replace);
    NcDim* xd = combrescdf.add_dim("x", nCompatStations);
    NcDim* yd = combrescdf.add_dim("y", nsamples);
    NcVar* x = combrescdf.add_var("x", ncFloat, xd);
    NcVar* y = combrescdf.add_var("y", ncFloat, yd);
    NcVar* z = combrescdf.add_var("z", ncFloat, xd, yd);
    float *xvals = new float[xd->size()];
    float *yvals = new float[yd->size()];
    float *zvals = new float[xd->size() * yd->size()];

    for (size_t j = 0; j < nsamples; ++j)
      yvals[j] = j * Stations.GetList().front()->GetDt()
          + Stations.GetList().front()->GetB();
    size_t dataindex = 0;
    size_t stationindex = 0;
    while (dataindex < nCompatStations)
      {
        xvals[dataindex] = f(Stations.GetList().at(stationindex).get());
        if (nsamples == Stations.GetList().at(stationindex)->GetData().size())
          {
            copy(Stations.GetList().at(stationindex)->GetData().begin(),
                Stations.GetList().at(stationindex)->GetData().end(),
                &zvals[dataindex * nsamples]);
            ++dataindex;
          }
        else
          {
            cerr << "File: " << Stations.GetList().at(stationindex)->GetName()
                << " " << "has incompatible length: " << Stations.GetList().at(
                stationindex)->GetData().size() << " not " << nsamples << endl;
          }
        ++stationindex;
      }
    x->put(xvals, xd->size());
    y->put(yvals, yd->size());
    z->put(zvals, z->edges());

    delete[] xvals;
    delete[] yvals;
    delete[] zvals;
  }

int main(int argc, char* argv[])
  {
    SeismicStationList Stations;
    string listfilename, outfilename;

    string version =
        "$Id: plotrecfunc.cpp 1882 2014-06-12 08:15:20Z mmoorkamp $";
    cout << endl << endl;
    cout << "Program " << version << endl;
    cout << " Print a number of receiver functions in one plot" << endl;
    cout
        << " Input is a list of files, output an ascii file for plotting with xmgrace and a netcdf file"
        << endl;
    cout
        << " Individual receiver functions are offset in y-direction by the chosen sorting criterion"
        << endl;
    cout << endl << endl;
    if (argc > 2)
      {
        listfilename = argv[1];
        outfilename = argv[2];
      }
    else
      {
        listfilename = AskFilename("List filename: ");
        outfilename = AskFilename("Output filename: ");
      }
    boost::function<double(SeismicDataComp*)> SortFunc;
    size_t sortchoice;
    cout << "Specify Sorting Criterion: " << endl;
    cout << "  1: Distance " << endl;
    cout << "  2: Backazimuth " << endl;
    cin >> sortchoice;
    switch (sortchoice)
      {
    case 1:
      SortFunc = &SeismicDataComp::GetGcarc;
      break;
    case 2:
      SortFunc = &SeismicDataComp::GetBaz;
      break;
    default:
      cerr << "Invalid choice !";
      return 100;
      break;
      }
    Stations.ReadList(listfilename);
    vector<double> distances(Stations.GetList().size(), 0.0);
    transform(Stations.GetList().begin(), Stations.GetList().end(),
        distances.begin(), boost::bind(SortFunc, boost::bind(
            &boost::shared_ptr<SeismicDataComp>::get, _1)));
    cout << "Parameters: " << endl;
    copy(distances.begin(), distances.end(), ostream_iterator<double> (cout,
        "\n"));
    cout << endl;
    double meandist = Mean(distances.begin(), distances.end());
    cout << "Mean Parameter: " << meandist << endl;
    vector<double> maxamps;
    for (SeismicStationList::tseiscompvector::iterator CurrentStation =
        Stations.GetList().begin(); CurrentStation != Stations.GetList().end(); CurrentStation++)
      {
        double currmax = *max_element(CurrentStation->get()->GetData().begin(),
            CurrentStation->get()->GetData().end());
        if (!std::isnan(currmax))
          maxamps.push_back(currmax);
      }
    double meanamp = Mean(maxamps.begin(), maxamps.end());
    double yshift = 5.0 * meanamp / meandist;
    cout << " Mean Amp: " << meanamp << endl;
    ofstream outfile(outfilename.c_str());
    for (SeismicStationList::tseiscompvector::iterator CurrentStation =
        Stations.GetList().begin(); CurrentStation != Stations.GetList().end(); CurrentStation++)
      {
        for (size_t i = 0; i < CurrentStation->get()->GetData().size(); ++i)
          {
            outfile << i * CurrentStation->get()->GetDt()
                + CurrentStation->get()->GetB() << " "
                << CurrentStation->get()->GetData().at(i) + (meandist
                    - SortFunc(CurrentStation->get())) * yshift << endl;
          }
        outfile << endl << endl;
      }
    WriteToNetCDF(outfilename, Stations, SortFunc);
  }