365 lines
17 KiB
C++
365 lines
17 KiB
C++
/********************************************************
|
|
* ██████╗ ██████╗████████╗██╗
|
|
* ██╔════╝ ██╔════╝╚══██╔══╝██║
|
|
* ██║ ███╗██║ ██║ ██║
|
|
* ██║ ██║██║ ██║ ██║
|
|
* ╚██████╔╝╚██████╗ ██║ ███████╗
|
|
* ╚═════╝ ╚═════╝ ╚═╝ ╚══════╝
|
|
* Geophysical Computational Tools & Library (GCTL)
|
|
*
|
|
* Copyright (c) 2022 Yi Zhang (yizhang-geo@zju.edu.cn)
|
|
*
|
|
* GCTL is distributed under a dual licensing scheme. You can redistribute
|
|
* it and/or modify it under the terms of the GNU Lesser General Public
|
|
* License as published by the Free Software Foundation, either version 2
|
|
* of the License, or (at your option) any later version. You should have
|
|
* received a copy of the GNU Lesser General Public License along with this
|
|
* program. If not, see <http://www.gnu.org/licenses/>.
|
|
*
|
|
* If the terms and conditions of the LGPL v.2. would prevent you from using
|
|
* the GCTL, please consider the option to obtain a commercial license for a
|
|
* fee. These licenses are offered by the GCTL's original author. As a rule,
|
|
* licenses are provided "as-is", unlimited in time for a one time fee. Please
|
|
* send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget
|
|
* to include some description of your company and the realm of its activities.
|
|
* Also add information on how to contact you by electronic and paper mail.
|
|
******************************************************/
|
|
|
|
#include "sac.h"
|
|
|
|
#ifdef GCTL_SEISMIC_MATHGL
|
|
|
|
gctl::SAC_PLOT::SAC_PLOT()
|
|
{
|
|
sig_ptr = nullptr; hd_ptr = nullptr;
|
|
}
|
|
|
|
gctl::SAC_PLOT::~SAC_PLOT(){}
|
|
|
|
void gctl::SAC_PLOT::link_source(SIG_UNIT *sig, SAC_HD *hd, std::string filename)
|
|
{
|
|
sig_ptr = sig; hd_ptr = hd; file = filename;
|
|
return;
|
|
}
|
|
|
|
int gctl::SAC_PLOT::Draw(mglGraph *gr)
|
|
{
|
|
gctl::array<float> x;
|
|
x.sequence((float) 0.0, hd_ptr->delta, sig_ptr->val.size());
|
|
|
|
mglData x_d, y_d;
|
|
x_d.Set(x.get(), x.size());
|
|
y_d.Set(sig_ptr->val.get(), sig_ptr->val.size());
|
|
|
|
float mean = sig_ptr->val.mean();
|
|
std::string s_nk = hd_ptr->knetwk; s_nk.erase(s_nk.find_last_not_of(" ") + 1);
|
|
std::string s_st = hd_ptr->kstnm; s_st.erase(s_st.find_last_not_of(" ") + 1);
|
|
std::string s_cp = hd_ptr->kcmpnm; s_cp.erase(s_cp.find_last_not_of(" ") + 1);
|
|
std::string leng = s_nk + "." + s_st + "." + s_cp;
|
|
|
|
gr->SetSize(800, 400);
|
|
gr->SetRanges(x.front(), x.back(), hd_ptr->xminimum - 0.5*(mean - hd_ptr->xminimum), hd_ptr->xmaximum + 0.5*(hd_ptr->xmaximum - mean));
|
|
gr->AddLegend(leng.c_str(), "h");
|
|
gr->Title(sig_ptr->t0.time_str().c_str());
|
|
gr->Box();
|
|
gr->Axis("x");
|
|
gr->Axis("y");
|
|
gr->Grid("xy", "h:");
|
|
gr->Legend(2);
|
|
gr->Label('x', "Time (s)", 0.0);
|
|
gr->Plot(x_d, y_d, "B");
|
|
|
|
if (file != "null")
|
|
{
|
|
std::string tp = file.substr(file.find_last_of('.') + 1);
|
|
if (tp == "eps") gr->WriteEPS(file.c_str());
|
|
if (tp == "png") gr->WritePNG(file.c_str());
|
|
if (tp == "jpg") gr->WriteJPEG(file.c_str());
|
|
}
|
|
return 0;
|
|
}
|
|
|
|
#endif // GCTL_SEISMIC_MATHGL
|
|
|
|
gctl::SAC::SAC(){}
|
|
|
|
gctl::SAC::~SAC(){}
|
|
|
|
gctl::SAC::SAC(std::string filename, bool headonly) {read(filename, headonly);}
|
|
|
|
void gctl::SAC::read(std::string filename, bool headonly)
|
|
{
|
|
std::string suf_str = filename.substr(filename.find_last_of('.') + 1);
|
|
if (suf_str != "SAC" && suf_str != "SAC_ASC")
|
|
{
|
|
throw std::runtime_error("Invalid file suffix.");
|
|
}
|
|
|
|
std::ifstream fs;
|
|
if (suf_str == "SAC_ASC")
|
|
{
|
|
gctl::open_infile(fs, filename);
|
|
|
|
fs >> head.delta >> head.depmin >> head.depmax >> head.unused1 >> head.odelta; //0 - 4
|
|
fs >> head.b >> head.e >> head.o >> head.a >> head.internal; //5 - 9
|
|
fs >> head.t0 >> head.t1 >> head.t2 >> head.t3 >> head.t4; //10 - 14
|
|
fs >> head.t5 >> head.t6 >> head.t7 >> head.t8 >> head.t9; //15 - 19
|
|
fs >> head.f >> head.resp0 >> head.resp1 >> head.resp2 >> head.resp3; //20 - 24
|
|
fs >> head.resp4 >> head.resp5 >> head.resp6 >> head.resp7 >> head.resp8; //25 - 29
|
|
fs >> head.resp9 >> head.stla >> head.stlo >> head.stel >> head.stdp; //30 - 34
|
|
fs >> head.evla >> head.evlo >> head.evel >> head.evdp >> head.mag; //35 - 39
|
|
fs >> head.user0 >> head.user1 >> head.user2 >> head.user3 >> head.user4; //40 - 44
|
|
fs >> head.user5 >> head.user6 >> head.user7 >> head.user8 >> head.user9; //45 - 49
|
|
fs >> head.dist >> head.az >> head.baz >> head.gcarc >> head.sb;
|
|
fs >> head.sdelta >> head.depmen >> head.cmpaz >> head.cmpinc >> head.xminimum;
|
|
fs >> head.xmaximum >> head.yminimum >> head.ymaximum >> head.adjtm >> head.unused2;
|
|
fs >> head.unused3 >> head.unused4 >> head.unused5 >> head.unused6 >> head.unused7;
|
|
fs >> head.nzyear >> head.nzjday >> head.nzhour >> head.nzmin >> head.nzsec;
|
|
fs >> head.nzmsec >> head.nvhdr >> head.norid >> head.nevid >> head.npts;
|
|
fs >> head.nsnpts >> head.nwfid >> head.nxsize >> head.nysize >> head.unused8;
|
|
fs >> head.iftype >> head.idep >> head.iztype >> head.unused9 >> head.iinst;
|
|
fs >> head.istreg >> head.ievreg >> head.ievtyp >> head.iqual >> head.isynth;
|
|
fs >> head.imagtyp >> head.imagsrc >> head.ibody >> head.unused10 >> head.unused11;
|
|
fs >> head.unused12 >> head.unused13 >> head.unused14 >> head.unused15 >> head.unused16;
|
|
fs >> head.leven >> head.lpspol >> head.lovrok >> head.lcalda >> head.unused17;
|
|
fs >> head.kstnm >> head.kevnm;
|
|
fs >> head.khole >> head.ko >> head.ka;
|
|
fs >> head.kt0 >> head.kt1 >> head.kt2;
|
|
fs >> head.kt3 >> head.kt4 >> head.kt5;
|
|
fs >> head.kt6 >> head.kt7 >> head.kt8;
|
|
fs >> head.kt9 >> head.kf >> head.kuser0;
|
|
fs >> head.kuser1 >> head.kuser2 >> head.kcmpnm;
|
|
fs >> head.knetwk >> head.kdatrd >> head.kinst;
|
|
|
|
if (!headonly)
|
|
{
|
|
signal.val.resize(head.npts);
|
|
for (int i = 0; i < head.npts; i++)
|
|
{
|
|
fs >> signal.val[i];
|
|
}
|
|
}
|
|
fs.close();
|
|
}
|
|
else
|
|
{
|
|
gctl::open_infile(fs, filename, "", std::ios::binary);
|
|
|
|
fs.read((char*)&head, sizeof(SAC_HD));
|
|
head.kstnm[7] = '\0'; head.kevnm[15] = '\0';
|
|
head.khole[7] = '\0'; head.ko[7] = '\0'; head.ka[7] = '\0';
|
|
head.kt0[7] = '\0'; head.kt1[7] = '\0'; head.kt2[7] = '\0';
|
|
head.kt3[7] = '\0'; head.kt4[7] = '\0'; head.kt5[7] = '\0';
|
|
head.kt6[7] = '\0'; head.kt7[7] = '\0'; head.kt8[7] = '\0';
|
|
head.kt9[7] = '\0'; head.kf[7] = '\0'; head.kuser0[7] = '\0';
|
|
head.kuser1[7] = '\0'; head.kuser2[7] = '\0'; head.kcmpnm[7] = '\0';
|
|
head.knetwk[7] = '\0'; head.kdatrd[7] = '\0'; head.kinst[7] = '\0';
|
|
|
|
if (!headonly)
|
|
{
|
|
signal.val.resize(head.npts);
|
|
array<float> tmp_s(head.npts);
|
|
|
|
fs.read((char*)tmp_s.get(), head.npts*sizeof(float));
|
|
for (size_t i = 0; i < head.npts; i++)
|
|
{
|
|
signal.val[i] = (double) tmp_s[i];
|
|
}
|
|
}
|
|
fs.close();
|
|
}
|
|
|
|
signal.delta = head.delta;
|
|
signal.t0.set_time(head.nzyear, head.nzjday, head.nzhour, head.nzmin, head.nzsec, head.nzmsec);
|
|
signal.te = signal.t0;
|
|
signal.te.add_duration(head.e);
|
|
sync();
|
|
return;
|
|
}
|
|
|
|
void gctl::SAC::save(std::string filename) const
|
|
{
|
|
std::string suf_str = filename.substr(filename.find_last_of('.') + 1);
|
|
if (suf_str != "SAC" && suf_str != "SAC_ASC")
|
|
{
|
|
throw std::runtime_error("[GCTL_SEISMIC] Invalid file suffix.");
|
|
}
|
|
|
|
if (head.npts != signal.val.size())
|
|
{
|
|
throw std::runtime_error("[GCTL_SEISMIC] Record length and the head section don't match.");
|
|
}
|
|
|
|
std::ofstream fs;
|
|
if (suf_str == "SAC_ASC")
|
|
{
|
|
gctl::open_outfile(fs, filename);
|
|
|
|
fs << std::setprecision(7) << std::setw(15) << head.delta << std::setw(15) << head.depmin << std::setw(15) << head.depmax << std::setw(15) << head.unused1 << std::setw(15) << head.odelta << "\n"; //0 - 4
|
|
fs << std::setprecision(7) << std::setw(15) << head.b << std::setw(15) << head.e << std::setw(15) << head.o << std::setw(15) << head.a << std::setw(15) << head.internal << "\n"; //5 - 9
|
|
fs << std::setprecision(7) << std::setw(15) << head.t0 << std::setw(15) << head.t1 << std::setw(15) << head.t2 << std::setw(15) << head.t3 << std::setw(15) << head.t4 << "\n"; //10 - 14
|
|
fs << std::setprecision(7) << std::setw(15) << head.t5 << std::setw(15) << head.t6 << std::setw(15) << head.t7 << std::setw(15) << head.t8 << std::setw(15) << head.t9 << "\n"; //15 - 19
|
|
fs << std::setprecision(7) << std::setw(15) << head.f << std::setw(15) << head.resp0 << std::setw(15) << head.resp1 << std::setw(15) << head.resp2 << std::setw(15) << head.resp3 << "\n"; //20 - 24
|
|
fs << std::setprecision(7) << std::setw(15) << head.resp4 << std::setw(15) << head.resp5 << std::setw(15) << head.resp6 << std::setw(15) << head.resp7 << std::setw(15) << head.resp8 << "\n"; //25 - 29
|
|
fs << std::setprecision(7) << std::setw(15) << head.resp9 << std::setw(15) << head.stla << std::setw(15) << head.stlo << std::setw(15) << head.stel << std::setw(15) << head.stdp << "\n"; //30 - 34
|
|
fs << std::setprecision(7) << std::setw(15) << head.evla << std::setw(15) << head.evlo << std::setw(15) << head.evel << std::setw(15) << head.evdp << std::setw(15) << head.mag << "\n"; //35 - 39
|
|
fs << std::setprecision(7) << std::setw(15) << head.user0 << std::setw(15) << head.user1 << std::setw(15) << head.user2 << std::setw(15) << head.user3 << std::setw(15) << head.user4 << "\n"; //40 - 44
|
|
fs << std::setprecision(7) << std::setw(15) << head.user5 << std::setw(15) << head.user6 << std::setw(15) << head.user7 << std::setw(15) << head.user8 << std::setw(15) << head.user9 << "\n"; //45 - 49
|
|
fs << std::setprecision(7) << std::setw(15) << head.dist << std::setw(15) << head.az << std::setw(15) << head.baz << std::setw(15) << head.gcarc << std::setw(15) << head.sb << "\n";
|
|
fs << std::setprecision(7) << std::setw(15) << head.sdelta << std::setw(15) << head.depmen << std::setw(15) << head.cmpaz << std::setw(15) << head.cmpinc << std::setw(15) << head.xminimum << "\n";
|
|
fs << std::setprecision(7) << std::setw(15) << head.xmaximum << std::setw(15) << head.yminimum << std::setw(15) << head.ymaximum << std::setw(15) << head.adjtm << std::setw(15) << head.unused2 << "\n";
|
|
fs << std::setprecision(7) << std::setw(15) << head.unused3 << std::setw(15) << head.unused4 << std::setw(15) << head.unused5 << std::setw(15) << head.unused6 << std::setw(15) << head.unused7 << "\n";
|
|
fs << std::setprecision(7) << std::setw(10) << head.nzyear << std::setw(10) << head.nzjday << std::setw(10) << head.nzhour << std::setw(10) << head.nzmin << std::setw(10) << head.nzsec << "\n";
|
|
fs << std::setprecision(7) << std::setw(10) << head.nzmsec << std::setw(10) << head.nvhdr << std::setw(10) << head.norid << std::setw(10) << head.nevid << std::setw(10) << head.npts << "\n";
|
|
fs << std::setprecision(7) << std::setw(10) << head.nsnpts << std::setw(10) << head.nwfid << std::setw(10) << head.nxsize << std::setw(10) << head.nysize << std::setw(10) << head.unused8 << "\n";
|
|
fs << std::setprecision(7) << std::setw(10) << head.iftype << std::setw(10) << head.idep << std::setw(10) << head.iztype << std::setw(10) << head.unused9 << std::setw(10) << head.iinst << "\n";
|
|
fs << std::setprecision(7) << std::setw(10) << head.istreg << std::setw(10) << head.ievreg << std::setw(10) << head.ievtyp << std::setw(10) << head.iqual << std::setw(10) << head.isynth << "\n";
|
|
fs << std::setprecision(7) << std::setw(10) << head.imagtyp << std::setw(10) << head.imagsrc << std::setw(10) << head.ibody << std::setw(10) << head.unused10 << std::setw(10) << head.unused11 << "\n";
|
|
fs << std::setprecision(7) << std::setw(10) << head.unused12 << std::setw(10) << head.unused13 << std::setw(10) << head.unused14 << std::setw(10) << head.unused15 << std::setw(10) << head.unused16 << "\n";
|
|
fs << std::setprecision(7) << std::setw(10) << head.leven << std::setw(10) << head.lpspol << std::setw(10) << head.lovrok << std::setw(10) << head.lcalda << std::setw(10) << head.unused17 << "\n";
|
|
fs << head.kstnm << " " << head.kevnm << " \n";
|
|
fs << head.khole << " " << head.ko << " " << head.ka << " \n";
|
|
fs << head.kt0 << " " << head.kt1 << " " << head.kt2 << " \n";
|
|
fs << head.kt3 << " " << head.kt4 << " " << head.kt5 << " \n";
|
|
fs << head.kt6 << " " << head.kt7 << " " << head.kt8 << " \n";
|
|
fs << head.kt9 << " " << head.kf << " " << head.kuser0 << " \n";
|
|
fs << head.kuser1 << " " << head.kuser2 << " " << head.kcmpnm << " \n";
|
|
fs << head.knetwk << " " << head.kdatrd << " " << head.kinst << " \n";
|
|
|
|
int i = 0;
|
|
while (i < head.npts)
|
|
{
|
|
fs << std::setprecision(7) << std::setw(15) << signal.val[i];
|
|
if ((i+1)%5 == 0) fs << "\n";
|
|
i++;
|
|
}
|
|
fs.close();
|
|
}
|
|
else
|
|
{
|
|
gctl::open_outfile(fs, filename, "", std::ios::binary);
|
|
|
|
fs.write((char*)&head, sizeof(SAC_HD));
|
|
|
|
array<float> tmp_s(head.npts);
|
|
for (size_t i = 0; i < head.npts; i++)
|
|
{
|
|
tmp_s[i] = (float) signal.val[i];
|
|
}
|
|
|
|
fs.write((char*)tmp_s.get(), head.npts*sizeof(float));
|
|
fs.close();
|
|
}
|
|
return;
|
|
}
|
|
|
|
void gctl::SAC::info() const
|
|
{
|
|
std::clog << "SAC Info\n-------------\n";
|
|
std::clog << "Station: " << head.kstnm << "\nNetwork: " << head.knetwk << "\nCompnent: " << head.kcmpnm << "\n";
|
|
std::clog << "Start-Time: " << signal.t0.time_str() << "\nEnd-Time: " << signal.te.time_str() << "\n";
|
|
std::clog << "Delta: " << head.delta << " s\nDuration: " << (double) (head.npts-1)*head.delta << " s\nSingal-Length: " << head.npts << "\n";
|
|
std::clog << "-------------\n";
|
|
return;
|
|
}
|
|
|
|
void gctl::SAC::sync()
|
|
{
|
|
head.nzyear = signal.t0.year;
|
|
head.nzjday = signal.t0.julday;
|
|
head.nzhour = signal.t0.hour;
|
|
head.nzmin = signal.t0.mint;
|
|
head.nzsec = signal.t0.sec;
|
|
head.nzmsec = signal.t0.msec;
|
|
head.npts = signal.val.size();
|
|
head.delta = signal.delta;
|
|
head.e = (head.npts - 1)*head.delta;
|
|
|
|
double xmin = signal.val[0];
|
|
double xmax = signal.val[0];
|
|
for (size_t i = 1; i < signal.val.size(); i++)
|
|
{
|
|
xmin = std::min(xmin, signal.val[i]);
|
|
xmax = std::max(xmax, signal.val[i]);
|
|
}
|
|
|
|
head.xmaximum = xmax;
|
|
head.xminimum = xmin;
|
|
return;
|
|
}
|
|
|
|
void gctl::SAC::cut_section(SAC &ret, double st, double et) const
|
|
{
|
|
signal.cut_section(ret.signal, st, et);
|
|
ret.head = head;
|
|
ret.sync();
|
|
return;
|
|
}
|
|
|
|
void gctl::SAC::cut_section(SAC &ret, const UTC_TIME &st, double len) const
|
|
{
|
|
double stt = st.diff_sec(signal.t0);
|
|
double ett = stt + len;
|
|
return cut_section(ret, stt, ett);
|
|
}
|
|
|
|
void gctl::SAC::stack_sac(const gctl::array<SAC*> &list)
|
|
{
|
|
int ls = list.size();
|
|
if (list.empty())
|
|
{
|
|
throw std::runtime_error("[GCTL_SEISMIC] Invalid stacking list for gctl::SAC::stack_sac(...)");
|
|
}
|
|
|
|
array<SIG_UNIT*> sigs(list.size());
|
|
for (size_t i = 0; i < ls; i++)
|
|
{
|
|
sigs[i] = &list[i]->signal;
|
|
}
|
|
|
|
signal.stack_sig_units(sigs);
|
|
|
|
head = list[0]->head;
|
|
strcpy(head.kstnm, "STACK ");
|
|
sync();
|
|
return;
|
|
}
|
|
|
|
void gctl::SAC::down_sampling(int factor)
|
|
{
|
|
signal.down_sampling(factor);
|
|
sync();
|
|
return;
|
|
}
|
|
|
|
void gctl::SAC::remove_mean()
|
|
{
|
|
signal.remove_mean();
|
|
sync();
|
|
return;
|
|
}
|
|
|
|
void gctl::SAC::remove_outliers(double factor)
|
|
{
|
|
signal.remove_outliers(factor);
|
|
sync();
|
|
return;
|
|
}
|
|
|
|
void gctl::SAC::normalize(double norm)
|
|
{
|
|
signal.normalize(norm);
|
|
sync();
|
|
return;
|
|
}
|
|
|
|
#ifdef GCTL_SEISMIC_MATHGL
|
|
|
|
int gctl::SAC::plot(std::string filename)
|
|
{
|
|
plt.link_source(&signal, &head, filename);
|
|
mglFLTK gr(&plt, "SAC plot");
|
|
return gr.Run();
|
|
}
|
|
|
|
#endif // GCTL_SEISMIC_MATHGL
|