plotsourcerel.cpp

Go to the documentation of this file.

#include "SeismicStationList.h"
#include <iostream>
#include <string>
#include "NetCDFTools.h"
#include "VecMat.h"

using namespace std;
using namespace gplib;

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

    string version = "$Id: plotsourcerel.cpp 1816 2009-09-07 11:28:35Z mmoorkamp $";
    cout << endl << endl;
    cout << "Program " << version << endl;
    cout << " Print a 2D histogram of backazimuth and distance relations"
        << endl;
    cout << " Input is a list of files, output a netcdf matrix file" << endl;
    cout << endl << endl;
    if (argc > 2)
      {
        listfilename = argv[1];
        outfilename = argv[2];
      }
    else
      {
        cout << "List filename: ";
        cin >> listfilename;
        cout << "Output filename: ";
        cin >> outfilename;
      }
    Stations.ReadList(listfilename);
    const double bazbinsize = 10;
    const double gcarcbinsize = 5;
    const double minbaz = 0;
    const double maxbaz = 360;
    const double mingcarc = 30;
    const double maxgcarc = 130;
    const int nbazbins = (maxbaz - minbaz) / bazbinsize + 1;
    const int ngcarcbins = (maxgcarc - mingcarc) / gcarcbinsize + 1;
    gplib::rmat Hist(nbazbins, ngcarcbins);
    Hist *= 0.0;
    for (SeismicStationList::tseiscompvector::iterator CurrentStation =
        Stations.GetList().begin(); CurrentStation != Stations.GetList().end(); CurrentStation++)
      {
        int bazindex = 0;
        int gcarcindex = 0;
        while ((bazindex + 1) * bazbinsize < (*CurrentStation)->GetBaz())
          ++bazindex;
        while ((gcarcindex + 1) * gcarcbinsize < (*CurrentStation)->GetGcarc())
          ++gcarcindex;
        Hist(bazindex, gcarcindex) += 1.0;
      }

    NcFile mtrescdf(outfilename.c_str(), NcFile::Replace);
    NcDim* xd = mtrescdf.add_dim("Baz", nbazbins);
    NcDim* yd = mtrescdf.add_dim("Gcarc", ngcarcbins);
    NcVar* x = mtrescdf.add_var("Baz", ncFloat, xd);
    NcVar* y = mtrescdf.add_var("Gcarc", ncFloat, yd);
    NcVar* z = mtrescdf.add_var("Count", 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 i = 0; i < Hist.size1(); ++i)
      {
        xvals[i] = i * bazbinsize;
        for (size_t j = 0; j < Hist.size2(); ++j)
          {
            yvals[j] = j * gcarcbinsize;
            zvals[i * Hist.size2() + j] = Hist(i, j);
          }
      }
    x->put(xvals, xd->size());
    y->put(yvals, yd->size());
    z->put(zvals, z->edges());
    delete xvals;
    delete yvals;
    delete zvals;
  }