source: ntrip/trunk/BNC/src/rinex/reqcanalyze.cpp @ 7942

Last change on this file since 7942 was 7942, checked in by stuerze, 3 years ago

minor changes to prevent QThread errors

File size: 36.1 KB
Line 
1// Part of BNC, a utility for retrieving decoding and
2// converting GNSS data streams from NTRIP broadcasters.
3//
4// Copyright (C) 2007
5// German Federal Agency for Cartography and Geodesy (BKG)
6// http://www.bkg.bund.de
7// Czech Technical University Prague, Department of Geodesy
8// http://www.fsv.cvut.cz
9//
10// Email: euref-ip@bkg.bund.de
11//
12// This program is free software; you can redistribute it and/or
13// modify it under the terms of the GNU General Public License
14// as published by the Free Software Foundation, version 2.
15//
16// This program is distributed in the hope that it will be useful,
17// but WITHOUT ANY WARRANTY; without even the implied warranty of
18// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19// GNU General Public License for more details.
20//
21// You should have received a copy of the GNU General Public License
22// along with this program; if not, write to the Free Software
23// Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
24
25/* -------------------------------------------------------------------------
26 * BKG NTRIP Client
27 * -------------------------------------------------------------------------
28 *
29 * Class:      t_reqcAnalyze
30 *
31 * Purpose:    Analyze RINEX Files
32 *
33 * Author:     L. Mervart
34 *
35 * Created:    11-Apr-2012
36 *
37 * Changes:
38 *
39 * -----------------------------------------------------------------------*/
40
41#include <iostream>
42#include <iomanip>
43#include <qwt_plot_renderer.h>
44
45#include "reqcanalyze.h"
46#include "bnccore.h"
47#include "bncsettings.h"
48#include "reqcedit.h"
49#include "bncutils.h"
50#include "graphwin.h"
51#include "polarplot.h"
52#include "availplot.h"
53#include "eleplot.h"
54#include "dopplot.h"
55#include "bncephuser.h"
56
57using namespace std;
58
59// Constructor
60////////////////////////////////////////////////////////////////////////////
61t_reqcAnalyze::t_reqcAnalyze(QObject* parent) : QThread(parent) {
62
63  bncSettings settings;
64
65  _logFileName     = settings.value("reqcOutLogFile").toString(); expandEnvVar(_logFileName);
66  _logFile         = 0;
67  _log             = 0;
68  _currEpo         = 0;
69  _obsFileNames    = settings.value("reqcObsFile").toString().split(",", QString::SkipEmptyParts);
70  _navFileNames    = settings.value("reqcNavFile").toString().split(",", QString::SkipEmptyParts);
71  _reqcPlotSignals = settings.value("reqcSkyPlotSignals").toString();
72  _defaultSignalTypes << "G:1&2" << "R:1&2" << "J:1&2" << "E:1&5" << "S:1&5" << "C:2&7";
73  if (_reqcPlotSignals.isEmpty()) {
74    _reqcPlotSignals = _defaultSignalTypes.join(" ");
75  }
76  analyzePlotSignals(_signalTypes);
77
78  connect(this, SIGNAL(dspSkyPlot(const QString&, const QString&, QVector<t_polarPoint*>*,
79                                  const QString&, QVector<t_polarPoint*>*,
80                                  const QByteArray&, double)),
81          this, SLOT(slotDspSkyPlot(const QString&, const QString&, QVector<t_polarPoint*>*,
82                                    const QString&, QVector<t_polarPoint*>*,
83                                    const QByteArray&, double)));
84
85  connect(this, SIGNAL(dspAvailPlot(const QString&, const QByteArray&)),
86          this, SLOT(slotDspAvailPlot(const QString&, const QByteArray&)));
87}
88
89// Destructor
90////////////////////////////////////////////////////////////////////////////
91t_reqcAnalyze::~t_reqcAnalyze() {
92  for (int ii = 0; ii < _rnxObsFiles.size(); ii++) {
93    delete _rnxObsFiles[ii];
94  }
95  for (int ii = 0; ii < _ephs.size(); ii++) {
96    delete _ephs[ii];
97  }
98  delete _log;     _log     = 0;
99  delete _logFile; _logFile = 0;
100  if (BNC_CORE->mode() != t_bncCore::interactive) {
101    qApp->exit(0);
102    msleep(100); //sleep 0.1 sec
103  }
104}
105
106//
107////////////////////////////////////////////////////////////////////////////
108void t_reqcAnalyze::run() {
109
110  static const double QC_FORMAT_VERSION = 1.1;
111
112  // Open Log File
113  // -------------
114  if (!_logFileName.isEmpty()) {
115    _logFile = new QFile(_logFileName);
116    if (_logFile->open(QIODevice::WriteOnly | QIODevice::Text)) {
117      _log = new QTextStream();
118      _log->setDevice(_logFile);
119    }
120  }
121
122  if (_log) {
123    *_log << "QC Format Version  : " << QString("%1").arg(QC_FORMAT_VERSION,3,'f',1)  << endl << endl;
124  }
125
126  // Check Ephemerides
127  // -----------------
128  checkEphemerides();
129
130  // Initialize RINEX Observation Files
131  // ----------------------------------
132  t_reqcEdit::initRnxObsFiles(_obsFileNames, _rnxObsFiles, _log);
133
134  // Read Ephemerides
135  // ----------------
136  t_reqcEdit::readEphemerides(_navFileNames, _ephs);
137
138  // Loop over all RINEX Files
139  // -------------------------
140  for (int ii = 0; ii < _rnxObsFiles.size(); ii++) {
141    analyzeFile(_rnxObsFiles[ii]);
142  }
143
144  // Exit
145  // ----
146  emit finished();
147  deleteLater();
148}
149
150//
151////////////////////////////////////////////////////////////////////////////
152void t_reqcAnalyze::analyzePlotSignals(QMap<char, QVector<QString> >& signalTypes) {
153
154  QStringList signalsOpt = _reqcPlotSignals.split(" ", QString::SkipEmptyParts);
155
156  for (int ii = 0; ii < signalsOpt.size(); ii++) {
157    QStringList input = signalsOpt.at(ii).split(QRegExp("[:&]"), QString::SkipEmptyParts);
158    if (input.size() > 1 && input[0].length() == 1) {
159      char system = input[0].toAscii().constData()[0];
160      QStringList sysValid       = _defaultSignalTypes.filter(QString(system));
161      QStringList defaultSignals = sysValid.at(0).split(QRegExp("[:&]"));
162      if (sysValid.isEmpty()) {continue;}
163      if (input[1][0].isDigit()) {
164        signalTypes[system].append(input[1]);
165      }
166      else {
167        signalTypes[system].append(defaultSignals[1]);
168      }
169      if (input.size() > 2) {
170        if (input[2][0].isDigit()) {
171          signalTypes[system].append(input[2]);
172        }
173        else {
174          signalTypes[system].append(defaultSignals[2]);
175        }
176      } else {
177        signalTypes[system].append(defaultSignals[2]);
178        if (signalTypes[system][0] == signalTypes[system][1]) {
179          signalTypes[system][0] = defaultSignals[1];
180        }
181      }
182    }
183  }
184}
185
186//
187////////////////////////////////////////////////////////////////////////////
188void t_reqcAnalyze::analyzeFile(t_rnxObsFile* obsFile) {
189
190  _qcFile.clear();
191
192  // A priori Coordinates
193  // --------------------
194  ColumnVector xyzSta = obsFile->xyz();
195
196  // Loop over all Epochs
197  // --------------------
198  try {
199    QMap<QString, bncTime> lastObsTime;
200    bool firstEpo = true;
201    while ( (_currEpo = obsFile->nextEpoch()) != 0) {
202      if (firstEpo) {
203        firstEpo = false;
204        _qcFile._startTime    = _currEpo->tt;
205        _qcFile._antennaName  = obsFile->antennaName();
206        _qcFile._markerName   = obsFile->markerName();
207        _qcFile._receiverType = obsFile->receiverType();
208        _qcFile._interval     = obsFile->interval();
209      }
210      _qcFile._endTime = _currEpo->tt;
211
212      t_qcEpo qcEpo;
213      qcEpo._epoTime = _currEpo->tt;
214      qcEpo._PDOP    = cmpDOP(xyzSta);
215
216      // Loop over all satellites
217      // ------------------------
218      for (unsigned iObs = 0; iObs < _currEpo->rnxSat.size(); iObs++) {
219        const t_rnxObsFile::t_rnxSat& rnxSat = _currEpo->rnxSat[iObs];
220        if (_navFileNames.size() &&
221            _numExpObs.find(rnxSat.prn) == _numExpObs.end()) {
222          _numExpObs[rnxSat.prn] = 0;
223        }
224        if (_signalTypes.find(rnxSat.prn.system()) == _signalTypes.end()) {
225          continue;
226        }
227        t_satObs satObs;
228        t_rnxObsFile::setObsFromRnx(obsFile, _currEpo, rnxSat, satObs);
229        t_qcSat& qcSat = qcEpo._qcSat[satObs._prn];
230        setQcObs(qcEpo._epoTime, xyzSta, satObs, lastObsTime, qcSat);
231        updateQcSat(qcSat, _qcFile._qcSatSum[satObs._prn]);
232      }
233      _qcFile._qcEpo.push_back(qcEpo);
234    }
235
236    analyzeMultipath();
237
238    if (_navFileNames.size()) {
239      setExpectedObs(_qcFile._startTime, _qcFile._endTime, _qcFile._interval, xyzSta);
240    }
241
242    preparePlotData(obsFile);
243
244    printReport(obsFile);
245  }
246  catch (QString str) {
247    if (_log) {
248      *_log << "Exception " << str << endl;
249    }
250    else {
251      qDebug() << str;
252    }
253  }
254}
255
256// Compute Dilution of Precision
257////////////////////////////////////////////////////////////////////////////
258double t_reqcAnalyze::cmpDOP(const ColumnVector& xyzSta) const {
259
260  if ( xyzSta.size() != 3 || (xyzSta[0] == 0.0 && xyzSta[1] == 0.0 && xyzSta[2] == 0.0) ) {
261    return 0.0;
262  }
263
264  unsigned nSat = _currEpo->rnxSat.size();
265
266  if (nSat < 4) {
267    return 0.0;
268  }
269
270  Matrix AA(nSat, 4);
271
272  unsigned nSatUsed = 0;
273  for (unsigned iSat = 0; iSat < nSat; iSat++) {
274
275    const t_rnxObsFile::t_rnxSat& rnxSat = _currEpo->rnxSat[iSat];
276    const t_prn& prn = rnxSat.prn;
277
278    if (_signalTypes.find(prn.system()) == _signalTypes.end()) {
279      continue;
280    }
281
282    t_eph* eph = 0;
283    for (int ie = 0; ie < _ephs.size(); ie++) {
284      if (_ephs[ie]->prn() == prn) {
285        eph = _ephs[ie];
286        break;
287      }
288    }
289    if (eph) {
290      ColumnVector xSat(4);
291      ColumnVector vv(3);
292      if (eph->getCrd(_currEpo->tt, xSat, vv, false) == success) {
293        ++nSatUsed;
294        ColumnVector dx = xSat.Rows(1,3) - xyzSta;
295        double rho = dx.norm_Frobenius();
296        AA(nSatUsed,1) = dx(1) / rho;
297        AA(nSatUsed,2) = dx(2) / rho;
298        AA(nSatUsed,3) = dx(3) / rho;
299        AA(nSatUsed,4) = 1.0;
300      }
301    }
302  }
303
304  if (nSatUsed < 4) {
305    return 0.0;
306  }
307
308  AA = AA.Rows(1, nSatUsed);
309
310  SymmetricMatrix QQ;
311  QQ << AA.t() * AA;
312  QQ = QQ.i();
313
314  return sqrt(QQ.trace());
315}
316
317//
318////////////////////////////////////////////////////////////////////////////
319void t_reqcAnalyze::updateQcSat(const t_qcSat& qcSat, t_qcSatSum& qcSatSum) {
320
321  for (int ii = 0; ii < qcSat._qcFrq.size(); ii++) {
322    const t_qcFrq& qcFrq    = qcSat._qcFrq[ii];
323    t_qcFrqSum&    qcFrqSum = qcSatSum._qcFrqSum[qcFrq._rnxType2ch];
324    qcFrqSum._numObs += 1;
325    if (qcFrq._slip) {
326      qcFrqSum._numSlipsFlagged += 1;
327    }
328    if (qcFrq._gap) {
329      qcFrqSum._numGaps += 1;
330    }
331    if (qcFrq._SNR > 0.0) {
332      qcFrqSum._numSNR += 1;
333      qcFrqSum._sumSNR += qcFrq._SNR;
334    }
335  }
336}
337
338//
339////////////////////////////////////////////////////////////////////////////
340void t_reqcAnalyze::setQcObs(const bncTime& epoTime, const ColumnVector& xyzSta,
341                             const t_satObs& satObs, QMap<QString, bncTime>& lastObsTime,
342                             t_qcSat& qcSat) {
343
344  t_eph* eph = 0;
345  for (int ie = 0; ie < _ephs.size(); ie++) {
346    if (_ephs[ie]->prn().system() == satObs._prn.system() &&
347        _ephs[ie]->prn().number() == satObs._prn.number()) {
348      eph = _ephs[ie];
349      break;
350    }
351  }
352  if (eph) {
353    ColumnVector xc(4);
354    ColumnVector vv(3);
355    if ( xyzSta.size() == 3 && (xyzSta[0] != 0.0 || xyzSta[1] != 0.0 || xyzSta[2] != 0.0) &&
356         eph->getCrd(epoTime, xc, vv, false) == success) {
357      double rho, eleSat, azSat;
358      topos(xyzSta(1), xyzSta(2), xyzSta(3), xc(1), xc(2), xc(3), rho, eleSat, azSat);
359      qcSat._eleSet = true;
360      qcSat._azDeg  = azSat  * 180.0/M_PI;
361      qcSat._eleDeg = eleSat * 180.0/M_PI;
362    }
363    if (satObs._prn.system() == 'R') {
364      qcSat._slotSet = true;
365      qcSat._slotNum = eph->slotNum();
366    }
367  }
368
369  // Availability and Slip Flags
370  // ---------------------------
371  for (unsigned ii = 0; ii < satObs._obs.size(); ii++) {
372    const t_frqObs* frqObs = satObs._obs[ii];
373
374    qcSat._qcFrq.push_back(t_qcFrq());
375    t_qcFrq& qcFrq = qcSat._qcFrq.back();
376
377    qcFrq._rnxType2ch = QString(frqObs->_rnxType2ch.c_str());
378    qcFrq._SNR        = frqObs->_snr;
379    qcFrq._slip       = frqObs->_slip;
380    qcFrq._phaseValid = frqObs->_phaseValid;
381    qcFrq._codeValid  = frqObs->_codeValid;
382    // Check Gaps
383    // ----------
384    QString key = QString(satObs._prn.toString().c_str()) + qcFrq._rnxType2ch;
385    if (lastObsTime[key].valid()) {
386      double dt = epoTime - lastObsTime[key];
387      if (dt > 1.5 * _qcFile._interval) {
388        qcFrq._gap = true;
389      }
390    }
391    lastObsTime[key] = epoTime;
392
393    // Compute the Multipath Linear Combination
394    // ----------------------------------------
395    if (frqObs->_codeValid) {
396      t_frequency::type fA = t_frequency::dummy;
397      t_frequency::type fB = t_frequency::dummy;
398      char sys             = satObs._prn.system();
399      std::string frqType1, frqType2;
400      if (_signalTypes.find(sys) != _signalTypes.end()) {
401        frqType1.push_back(sys);
402        frqType1.push_back(_signalTypes[sys][0][0].toAscii());
403        frqType2.push_back(sys);
404        frqType2.push_back(_signalTypes[sys][1][0].toAscii());
405        if      (frqObs->_rnxType2ch[0] == frqType1[1]) {
406          fA = t_frequency::toInt(frqType1);
407          fB = t_frequency::toInt(frqType2);
408        }
409        else if (frqObs->_rnxType2ch[0] == frqType2[1]) {
410          fA = t_frequency::toInt(frqType2);
411          fB = t_frequency::toInt(frqType1);
412        }
413      }
414      if (fA != t_frequency::dummy && fB != t_frequency::dummy) {
415        double f_a = t_CST::freq(fA, qcSat._slotNum);
416        double f_b = t_CST::freq(fB, qcSat._slotNum);
417        double C_a = frqObs->_code;
418
419        bool   foundA = false;
420        double L_a    = 0.0;
421        bool   foundB = false;
422        double L_b    = 0.0;
423        for (unsigned jj = 0; jj < satObs._obs.size(); jj++) {
424          const t_frqObs* frqObsHlp = satObs._obs[jj];
425          if      (frqObsHlp->_rnxType2ch[0] == t_frequency::toString(fA)[1] &&
426              frqObsHlp->_phaseValid) {
427            foundA = true;
428            L_a    = frqObsHlp->_phase * t_CST::c / f_a;
429          }
430          else if (frqObsHlp->_rnxType2ch[0] == t_frequency::toString(fB)[1] &&
431              frqObsHlp->_phaseValid) {
432            foundB = true;
433            L_b    = frqObsHlp->_phase * t_CST::c / f_b;
434          }
435        }
436        if (foundA && foundB) {
437          qcFrq._setMP = true;
438          qcFrq._rawMP = C_a - L_a - 2.0*f_b*f_b/(f_a*f_a-f_b*f_b) * (L_a - L_b);
439        }
440      }
441    }
442  } // satObs loop
443}
444
445//
446////////////////////////////////////////////////////////////////////////////
447void t_reqcAnalyze::analyzeMultipath() {
448
449  const double SLIPTRESH = 10.0;  // cycle-slip threshold (meters)
450  const double chunkStep = 600.0; // 10 minutes
451
452  // Loop over all satellites available
453  // ----------------------------------
454  QMutableMapIterator<t_prn, t_qcSatSum> itSat(_qcFile._qcSatSum);
455  while (itSat.hasNext()) {
456    itSat.next();
457    const t_prn& prn      = itSat.key();
458    t_qcSatSum&  qcSatSum = itSat.value();
459
460    // Loop over all frequencies available
461    // -----------------------------------
462    QMutableMapIterator<QString, t_qcFrqSum> itFrq(qcSatSum._qcFrqSum);
463    while (itFrq.hasNext()) {
464      itFrq.next();
465      const QString& frqType  = itFrq.key();
466      t_qcFrqSum&    qcFrqSum = itFrq.value();
467
468      // Loop over all Chunks of Data
469      // ----------------------------
470      for (bncTime chunkStart = _qcFile._startTime;
471           chunkStart < _qcFile._endTime; chunkStart += chunkStep) {
472
473        bncTime chunkEnd = chunkStart + chunkStep;
474
475        QVector<t_qcFrq*> frqVec;
476        QVector<double>   MP;
477
478        // Loop over all Epochs within one Chunk of Data
479        // ---------------------------------------------
480        for (int iEpo = 0; iEpo < _qcFile._qcEpo.size(); iEpo++) {
481          t_qcEpo& qcEpo = _qcFile._qcEpo[iEpo];
482          if (chunkStart <= qcEpo._epoTime && qcEpo._epoTime < chunkEnd) {
483            if (qcEpo._qcSat.contains(prn)) {
484              t_qcSat& qcSat = qcEpo._qcSat[prn];
485              for (int iFrq = 0; iFrq < qcSat._qcFrq.size(); iFrq++) {
486                t_qcFrq& qcFrq = qcSat._qcFrq[iFrq];
487                if (qcFrq._rnxType2ch == frqType) {
488                  frqVec << &qcFrq;
489                  if (qcFrq._setMP) {
490                    MP << qcFrq._rawMP;
491                  }
492                }
493              }
494            }
495          }
496        }
497
498        // Compute the multipath mean and standard deviation
499        // -------------------------------------------------
500        if (MP.size() > 1) {
501          double meanMP = 0.0;
502          for (int ii = 0; ii < MP.size(); ii++) {
503            meanMP += MP[ii];
504          }
505          meanMP /= MP.size();
506
507          bool slipMP = false;
508
509          double stdMP = 0.0;
510          for (int ii = 0; ii < MP.size(); ii++) {
511            double diff = MP[ii] - meanMP;
512            if (fabs(diff) > SLIPTRESH) {
513              slipMP = true;
514              break;
515            }
516            stdMP += diff * diff;
517          }
518
519          if (slipMP) {
520            stdMP = 0.0;
521            stdMP = 0.0;
522            qcFrqSum._numSlipsFound += 1;
523          }
524          else {
525            stdMP = sqrt(stdMP / (MP.size()-1));
526            qcFrqSum._numMP += 1;
527            qcFrqSum._sumMP += stdMP;
528          }
529
530          for (int ii = 0; ii < frqVec.size(); ii++) {
531            t_qcFrq* qcFrq = frqVec[ii];
532            if (slipMP) {
533              qcFrq->_slip = true;
534            }
535            else {
536              qcFrq->_stdMP = stdMP;
537            }
538          }
539        }
540      } // chunk loop
541    } // frq loop
542  } // sat loop
543}
544
545//
546////////////////////////////////////////////////////////////////////////////
547void t_reqcAnalyze::preparePlotData(const t_rnxObsFile* obsFile) {
548
549  QString mp1Title = "Multipath\n";
550  QString mp2Title = "Multipath\n";
551  QString sn1Title = "Signal-to-Noise Ratio\n";
552  QString sn2Title = "Signal-to-Noise Ratio\n";
553
554  for(QMap<char, QVector<QString> >::iterator it = _signalTypes.begin();
555      it != _signalTypes.end(); it++) {
556      mp1Title += QString(it.key()) + ":" + it.value()[0] + " ";
557      sn1Title += QString(it.key()) + ":" + it.value()[0] + " ";
558      mp2Title += QString(it.key()) + ":" + it.value()[1] + " ";
559      sn2Title += QString(it.key()) + ":" + it.value()[1] + " ";
560  }
561
562  QVector<t_polarPoint*>* dataMP1  = new QVector<t_polarPoint*>;
563  QVector<t_polarPoint*>* dataMP2  = new QVector<t_polarPoint*>;
564  QVector<t_polarPoint*>* dataSNR1 = new QVector<t_polarPoint*>;
565  QVector<t_polarPoint*>* dataSNR2 = new QVector<t_polarPoint*>;
566
567  // Loop over all observations
568  // --------------------------
569  for (int iEpo = 0; iEpo < _qcFile._qcEpo.size(); iEpo++) {
570    t_qcEpo& qcEpo = _qcFile._qcEpo[iEpo];
571    QMapIterator<t_prn, t_qcSat> it(qcEpo._qcSat);
572    while (it.hasNext()) {
573      it.next();
574      const t_prn&   prn   = it.key();
575      const t_qcSat& qcSat = it.value();
576      if (qcSat._eleSet) {
577
578        QString frqType[2];
579
580        for (int iFrq = 0; iFrq < qcSat._qcFrq.size(); iFrq++) {
581          const t_qcFrq& qcFrq = qcSat._qcFrq[iFrq];
582
583          for (int ii = 0; ii < 2; ii++) {
584            if (frqType[ii].isEmpty()) {
585              QMapIterator<char, QVector<QString> > it(_signalTypes);
586              while (it.hasNext()) {
587                it.next();
588                if (it.key() == prn.system()) {
589                  if (it.value()[ii] == qcFrq._rnxType2ch || it.value()[ii] == qcFrq._rnxType2ch.left(1)) {
590                    frqType[ii] = qcFrq._rnxType2ch;
591                    break;
592                  }
593                }
594              }
595            }
596          }
597          if      (qcFrq._rnxType2ch == frqType[0]) {
598            (*dataSNR1) << (new t_polarPoint(qcSat._azDeg, 90.0 - qcSat._eleDeg, qcFrq._SNR));
599            (*dataMP1)  << (new t_polarPoint(qcSat._azDeg, 90.0 - qcSat._eleDeg, qcFrq._stdMP));
600          }
601          else if (qcFrq._rnxType2ch == frqType[1]) {
602            (*dataSNR2) << (new t_polarPoint(qcSat._azDeg, 90.0 - qcSat._eleDeg, qcFrq._SNR));
603            (*dataMP2)  << (new t_polarPoint(qcSat._azDeg, 90.0 - qcSat._eleDeg, qcFrq._stdMP));
604          }
605        }
606      }
607    }
608  }
609
610  // Show the plots
611  // --------------
612  if (BNC_CORE->GUIenabled()) {
613    QFileInfo  fileInfo(obsFile->fileName());
614    QByteArray title = fileInfo.fileName().toAscii();
615    emit dspSkyPlot(obsFile->fileName(), mp1Title,  dataMP1,  mp2Title,  dataMP2,  "Meters",  2.0);
616    emit dspSkyPlot(obsFile->fileName(), sn1Title, dataSNR1, sn2Title, dataSNR2, "dbHz",   54.0);
617    emit dspAvailPlot(obsFile->fileName(), title);
618  }
619  else {
620    for (int ii = 0; ii < dataMP1->size(); ii++) {
621      delete dataMP1->at(ii);
622    }
623    delete dataMP1;
624    for (int ii = 0; ii < dataMP2->size(); ii++) {
625      delete dataMP2->at(ii);
626    }
627    delete dataMP2;
628    for (int ii = 0; ii < dataSNR1->size(); ii++) {
629      delete dataSNR1->at(ii);
630    }
631    delete dataSNR1;
632    for (int ii = 0; ii < dataSNR2->size(); ii++) {
633      delete dataSNR2->at(ii);
634    }
635    delete dataSNR2;
636  }
637}
638
639//
640////////////////////////////////////////////////////////////////////////////
641void t_reqcAnalyze::slotDspSkyPlot(const QString& fileName, const QString& title1,
642                                   QVector<t_polarPoint*>* data1, const QString& title2,
643                                   QVector<t_polarPoint*>* data2, const QByteArray& scaleTitle,
644                                   double maxValue) {
645
646  if (BNC_CORE->GUIenabled()) {
647
648    if (maxValue == 0.0) {
649      if (data1) {
650        for (int ii = 0; ii < data1->size(); ii++) {
651          double val = data1->at(ii)->_value;
652          if (maxValue < val) {
653            maxValue = val;
654          }
655        }
656      }
657      if (data2) {
658        for (int ii = 0; ii < data2->size(); ii++) {
659          double val = data2->at(ii)->_value;
660          if (maxValue < val) {
661            maxValue = val;
662          }
663        }
664      }
665    }
666
667    QwtInterval scaleInterval(0.0, maxValue);
668
669    QVector<QWidget*> plots;
670    if (data1) {
671      QwtText title(title1);
672      QFont font = title.font(); font.setPointSize(font.pointSize()-1); title.setFont(font);
673      t_polarPlot* plot1 = new t_polarPlot(title, scaleInterval, BNC_CORE->mainWindow());
674      plot1->addCurve(data1);
675      plots << plot1;
676    }
677    if (data2) {
678      QwtText title(title2);
679      QFont font = title.font(); font.setPointSize(font.pointSize()-1); title.setFont(font);
680      t_polarPlot* plot2 = new t_polarPlot(title, scaleInterval, BNC_CORE->mainWindow());
681      plot2->addCurve(data2);
682      plots << plot2;
683    }
684
685    t_graphWin* graphWin = new t_graphWin(0, fileName, plots,
686                                          &scaleTitle, &scaleInterval);
687
688    graphWin->show();
689
690    bncSettings settings;
691    QString dirName = settings.value("reqcPlotDir").toString();
692    if (!dirName.isEmpty()) {
693      QByteArray ext = (scaleTitle == "Meters") ? "_M.png" : "_S.png";
694      graphWin->savePNG(dirName, ext);
695    }
696  }
697}
698
699//
700////////////////////////////////////////////////////////////////////////////
701void t_reqcAnalyze::slotDspAvailPlot(const QString& fileName, const QByteArray& title) {
702
703  t_plotData              plotData;
704  QMap<t_prn, t_plotData> plotDataMap;
705
706  for (int ii = 0; ii < _qcFile._qcEpo.size(); ii++) {
707    const t_qcEpo& qcEpo = _qcFile._qcEpo[ii];
708    double mjdX24 = qcEpo._epoTime.mjddec() * 24.0;
709
710    plotData._mjdX24 << mjdX24;
711    plotData._PDOP   << qcEpo._PDOP;
712    plotData._numSat << qcEpo._qcSat.size();
713
714    QMapIterator<t_prn, t_qcSat> it(qcEpo._qcSat);
715    while (it.hasNext()) {
716      it.next();
717      const t_prn&   prn   = it.key();
718      const t_qcSat& qcSat = it.value();
719
720      t_plotData&    data  = plotDataMap[prn];
721
722      if (qcSat._eleSet) {
723        data._mjdX24 << mjdX24;
724        data._eleDeg << qcSat._eleDeg;
725      }
726
727      char frqChar1 = _signalTypes[prn.system()][0][0].toAscii();
728      char frqChar2 = _signalTypes[prn.system()][1][0].toAscii();
729
730      QString frqType1;
731      QString frqType2;
732      for (int iFrq = 0; iFrq < qcSat._qcFrq.size(); iFrq++) {
733        const t_qcFrq& qcFrq = qcSat._qcFrq[iFrq];
734        if (qcFrq._rnxType2ch[0] == frqChar1 && frqType1.isEmpty()) {
735          frqType1 = qcFrq._rnxType2ch;
736        }
737        if (qcFrq._rnxType2ch[0] == frqChar2 && frqType2.isEmpty()) {
738          frqType2 = qcFrq._rnxType2ch;
739        }
740        if      (qcFrq._rnxType2ch == frqType1) {
741          if      (qcFrq._slip) {
742            data._L1slip << mjdX24;
743          }
744          else if (qcFrq._gap) {
745            data._L1gap << mjdX24;
746          }
747          else {
748            data._L1ok << mjdX24;
749          }
750        }
751        else if (qcFrq._rnxType2ch == frqType2) {
752          if      (qcFrq._slip) {
753            data._L2slip << mjdX24;
754          }
755          else if (qcFrq._gap) {
756            data._L2gap << mjdX24;
757          }
758          else {
759            data._L2ok << mjdX24;
760          }
761        }
762      }
763    }
764  }
765
766  if (BNC_CORE->GUIenabled()) {
767    t_availPlot* plotA = new t_availPlot(0, plotDataMap);
768    plotA->setTitle(title);
769
770    t_elePlot* plotZ = new t_elePlot(0, plotDataMap);
771
772    t_dopPlot* plotD = new t_dopPlot(0, plotData);
773
774    QVector<QWidget*> plots;
775    plots << plotA << plotZ << plotD;
776    t_graphWin* graphWin = new t_graphWin(0, fileName, plots, 0, 0);
777
778    int ww = QFontMetrics(graphWin->font()).width('w');
779    graphWin->setMinimumSize(120*ww, 40*ww);
780
781    graphWin->show();
782
783    bncSettings settings;
784    QString dirName = settings.value("reqcPlotDir").toString();
785    if (!dirName.isEmpty()) {
786      QByteArray ext = "_A.png";
787      graphWin->savePNG(dirName, ext);
788    }
789  }
790}
791
792// Finish the report
793////////////////////////////////////////////////////////////////////////////
794void t_reqcAnalyze::printReport(const t_rnxObsFile* obsFile) {
795
796  if (!_log) {
797    return;
798  }
799
800  QFileInfo obsFi(obsFile->fileName());
801  QString obsFileName = obsFi.fileName();
802
803  // Summary
804  // -------
805  *_log << "Observation File   : " << obsFileName                                   << endl
806        << "RINEX Version      : " << QString("%1").arg(obsFile->version(),4,'f',2) << endl
807        << "Marker Name        : " << _qcFile._markerName                           << endl
808        << "Marker Number      : " << obsFile->markerNumber()                       << endl
809        << "Receiver           : " << _qcFile._receiverType                         << endl
810        << "Antenna            : " << _qcFile._antennaName                          << endl
811        << "Position XYZ       : " << QString("%1 %2 %3").arg(obsFile->xyz()(1), 14, 'f', 4)
812                                                        .arg(obsFile->xyz()(2), 14, 'f', 4)
813                                                        .arg(obsFile->xyz()(3), 14, 'f', 4) << endl
814        << "Antenna dH/dE/dN   : " << QString("%1 %2 %3").arg(obsFile->antNEU()(3), 8, 'f', 4)
815                                                        .arg(obsFile->antNEU()(2), 8, 'f', 4)
816                                                        .arg(obsFile->antNEU()(1), 8, 'f', 4) << endl
817        << "Start Time         : " << _qcFile._startTime.datestr().c_str()         << ' '
818                                   << _qcFile._startTime.timestr(1,'.').c_str()    << endl
819        << "End Time           : " << _qcFile._endTime.datestr().c_str()           << ' '
820                                   << _qcFile._endTime.timestr(1,'.').c_str()      << endl
821        << "Interval           : " << _qcFile._interval                            << endl;
822
823  // Number of systems
824  // -----------------
825  QMap<QChar, QVector<const t_qcSatSum*> > systemMap;
826  QMapIterator<t_prn, t_qcSatSum> itSat(_qcFile._qcSatSum);
827  while (itSat.hasNext()) {
828    itSat.next();
829    const t_prn&      prn      = itSat.key();
830    const t_qcSatSum& qcSatSum = itSat.value();
831    systemMap[prn.system()].push_back(&qcSatSum);
832  }
833  *_log << "Navigation Systems : " << systemMap.size() << "   ";
834
835  QMapIterator<QChar, QVector<const t_qcSatSum*> > itSys(systemMap);
836  while (itSys.hasNext()) {
837    itSys.next();
838    *_log << ' ' << itSys.key();
839  }
840  *_log << endl;
841
842  // Observation types per system
843  // -----------------------------
844  for (int iSys = 0; iSys < obsFile->numSys(); iSys++) {
845    char sys = obsFile->system(iSys);
846    if (sys != ' ') {
847      *_log << "Observation Types " << sys << ":";
848      for (int iType = 0; iType < obsFile->nTypes(sys); iType++) {
849        QString type = obsFile->obsType(sys, iType);
850        *_log << " " << type;
851      }
852      *_log << endl;
853    }
854  }
855
856  // System specific summary
857  // -----------------------
858  itSys.toFront();
859  while (itSys.hasNext()) {
860    itSys.next();
861    const QChar&                      sys      = itSys.key();
862    const QVector<const t_qcSatSum*>& qcSatVec = itSys.value();
863    int numExpectedObs = 0;
864    for(QMap<t_prn, int>::iterator it = _numExpObs.begin();
865        it != _numExpObs.end(); it++) {
866      if (sys == it.key().system()) {
867        numExpectedObs += it.value();
868      }
869    }
870    QString prefixSys = QString("  ") + sys + QString(": ");
871    QMap<QString, QVector<const t_qcFrqSum*> > frqMap;
872    for (int ii = 0; ii < qcSatVec.size(); ii++) {
873      const t_qcSatSum* qcSatSum = qcSatVec[ii];
874      QMapIterator<QString, t_qcFrqSum> itFrq(qcSatSum->_qcFrqSum);
875      while (itFrq.hasNext()) {
876        itFrq.next();
877        QString           frqType  = itFrq.key(); if (frqType.length() < 2) frqType += '?';
878        const t_qcFrqSum& qcFrqSum = itFrq.value();
879        frqMap[frqType].push_back(&qcFrqSum);
880      }
881    }
882    *_log << endl
883          << prefixSys << "Satellites: " << qcSatVec.size() << endl
884          << prefixSys << "Signals   : " << frqMap.size() << "   ";
885    QMapIterator<QString, QVector<const t_qcFrqSum*> > itFrq(frqMap);
886    while (itFrq.hasNext()) {
887      itFrq.next();
888      QString frqType = itFrq.key(); if (frqType.length() < 2) frqType += '?';
889      *_log << ' ' << frqType;
890    }
891    *_log << endl;
892    QString prefixSys2 = "    " + prefixSys;
893    itFrq.toFront();
894    while (itFrq.hasNext()) {
895      itFrq.next();
896      QString                          frqType  = itFrq.key(); if (frqType.length() < 2) frqType += '?';
897      const QVector<const t_qcFrqSum*> qcFrqVec = itFrq.value();
898      QString prefixFrq = QString("  ") + frqType + QString(": ");
899      int    numObs          = 0;
900      int    numSlipsFlagged = 0;
901      int    numSlipsFound   = 0;
902      int    numGaps         = 0;
903      int    numSNR          = 0;
904      int    numMP           = 0;
905      double sumSNR          = 0.0;
906      double sumMP           = 0.0;
907      for (int ii = 0; ii < qcFrqVec.size(); ii++) {
908        const t_qcFrqSum* qcFrqSum = qcFrqVec[ii];
909        numObs          += qcFrqSum->_numObs         ;
910        numSlipsFlagged += qcFrqSum->_numSlipsFlagged;
911        numSlipsFound   += qcFrqSum->_numSlipsFound  ;
912        numGaps         += qcFrqSum->_numGaps        ;
913        numSNR          += qcFrqSum->_numSNR;
914        numMP           += qcFrqSum->_numMP;
915        sumSNR          += qcFrqSum->_sumSNR;
916        sumMP           += qcFrqSum->_sumMP;
917      }
918      if (numSNR > 0) {
919        sumSNR /= numSNR;
920      }
921      if (numMP > 0) {
922        sumMP /= numMP;
923      }
924
925      double ratio = (double(numObs) / double(numExpectedObs)) * 100.0;
926
927      *_log << endl
928            << prefixSys2 << prefixFrq << "Observations      : ";
929      if(_navFileNames.isEmpty() || _navFileIncomplete.contains(sys.toLatin1())) {
930        *_log << QString("%1\n").arg(numObs,           6);
931      }
932      else {
933        *_log << QString("%1 (%2) %3 \%\n").arg(numObs,           6).arg(numExpectedObs,           8).arg(ratio, 8, 'f', 2);
934      }
935      *_log << prefixSys2 << prefixFrq << "Slips (file+found): " << QString("%1 +").arg(numSlipsFlagged,  8)
936                                                                 << QString("%1\n").arg(numSlipsFound,    8)
937            << prefixSys2 << prefixFrq << "Gaps              : " << QString("%1\n").arg(numGaps,          8)
938            << prefixSys2 << prefixFrq << "Mean SNR          : " << QString("%1\n").arg(sumSNR,   8, 'f', 1)
939            << prefixSys2 << prefixFrq << "Mean Multipath    : " << QString("%1\n").arg(sumMP,    8, 'f', 2);
940    }
941  }
942
943  // Epoch-Specific Output
944  // ---------------------
945  bncSettings settings;
946  if (Qt::CheckState(settings.value("reqcLogSummaryOnly").toInt()) == Qt::Checked) {
947    return;
948  }
949  *_log << endl;
950  for (int iEpo = 0; iEpo < _qcFile._qcEpo.size(); iEpo++) {
951    const t_qcEpo& qcEpo = _qcFile._qcEpo[iEpo];
952
953    unsigned year, month, day, hour, min;
954    double sec;
955    qcEpo._epoTime.civil_date(year, month, day);
956    qcEpo._epoTime.civil_time(hour, min, sec);
957
958    QString dateStr;
959    QTextStream(&dateStr) << QString("> %1 %2 %3 %4 %5%6")
960      .arg(year,  4)
961      .arg(month, 2, 10, QChar('0'))
962      .arg(day,   2, 10, QChar('0'))
963      .arg(hour,  2, 10, QChar('0'))
964      .arg(min,   2, 10, QChar('0'))
965      .arg(sec,  11, 'f', 7);
966
967    *_log << dateStr << QString(" %1").arg(qcEpo._qcSat.size(), 2)
968          << QString(" %1").arg(qcEpo._PDOP, 4, 'f', 1)
969          << endl;
970
971    QMapIterator<t_prn, t_qcSat> itSat(qcEpo._qcSat);
972    while (itSat.hasNext()) {
973      itSat.next();
974      const t_prn&   prn   = itSat.key();
975      const t_qcSat& qcSat = itSat.value();
976
977      *_log << prn.toString().c_str()
978            << QString(" %1 %2").arg(qcSat._eleDeg, 6, 'f', 2).arg(qcSat._azDeg, 7, 'f', 2);
979
980      int numObsTypes = 0;
981      for (int iFrq = 0; iFrq < qcSat._qcFrq.size(); iFrq++) {
982        const t_qcFrq& qcFrq = qcSat._qcFrq[iFrq];
983        if (qcFrq._phaseValid) {
984          numObsTypes += 1;
985        }
986        if (qcFrq._codeValid) {
987          numObsTypes += 1;
988        }
989      }
990      *_log << QString("  %1").arg(numObsTypes, 2);
991
992      for (int iFrq = 0; iFrq < qcSat._qcFrq.size(); iFrq++) {
993        const t_qcFrq& qcFrq = qcSat._qcFrq[iFrq];
994        if (qcFrq._phaseValid) {
995          *_log << "  L" << qcFrq._rnxType2ch << ' ';
996          if (qcFrq._slip) {
997            *_log << 's';
998          }
999          else {
1000            *_log << '.';
1001          }
1002          if (qcFrq._gap) {
1003            *_log << 'g';
1004          }
1005          else {
1006            *_log << '.';
1007          }
1008          *_log << QString(" %1").arg(qcFrq._SNR,   4, 'f', 1);
1009        }
1010        if (qcFrq._codeValid) {
1011          *_log << "  C" << qcFrq._rnxType2ch << ' ';
1012          if (qcFrq._gap) {
1013            *_log << " g";
1014          }
1015          else {
1016            *_log << " .";
1017          }
1018          *_log << QString(" %1").arg(qcFrq._stdMP, 3, 'f', 2);
1019        }
1020      }
1021      *_log << endl;
1022    }
1023  }
1024  _log->flush();
1025}
1026
1027//
1028////////////////////////////////////////////////////////////////////////////
1029void t_reqcAnalyze::checkEphemerides() {
1030
1031  QString navFileName;
1032  QStringListIterator namIt(_navFileNames);
1033  bool firstName = true;
1034  while (namIt.hasNext()) {
1035    QFileInfo navFi(namIt.next());
1036    if (firstName) {
1037      firstName = false;
1038      navFileName += navFi.fileName();
1039    }
1040    else {
1041      navFileName += ", " + navFi.fileName();
1042    }
1043  }
1044  if (_log) {
1045    *_log << "Navigation File(s) : " << navFileName                                   << endl;
1046  }
1047  QStringListIterator it(_navFileNames);
1048  while (it.hasNext()) {
1049    const QString& fileName = it.next();
1050    unsigned numOK  = 0;
1051    unsigned numBad = 0;
1052    bncEphUser ephUser(false);
1053    t_rnxNavFile rnxNavFile(fileName, t_rnxNavFile::input);
1054    for (unsigned ii = 0; ii < rnxNavFile.ephs().size(); ii++) {
1055      t_eph* eph = rnxNavFile.ephs()[ii];
1056      ephUser.putNewEph(eph, true);
1057      if (eph->checkState() == t_eph::bad) {
1058        ++numBad;
1059      }
1060      else {
1061        ++numOK;
1062      }
1063    }
1064    if (_log) {
1065      *_log << "Ephemeris          : " << numOK << " OK   " << numBad << " BAD" << endl;
1066    }
1067    if (numBad > 0) {
1068      for (unsigned ii = 0; ii < rnxNavFile.ephs().size(); ii++) {
1069        t_eph* eph = rnxNavFile.ephs()[ii];
1070        if (eph->checkState() == t_eph::bad) {
1071          QFileInfo navFi(fileName);
1072          if (_log) {
1073            *_log << "  Bad Ephemeris   : " << navFi.fileName() << ' '
1074                  << eph->toString(3.0).left(24) << endl;
1075          }
1076        }
1077      }
1078    }
1079  }
1080  if (_log) {
1081    *_log << endl;
1082  }
1083}
1084
1085void t_reqcAnalyze::setExpectedObs(const bncTime& startTime, const bncTime& endTime,
1086                                   double interval, const ColumnVector& xyzSta) {
1087
1088  for(QMap<t_prn, int>::iterator it = _numExpObs.begin();
1089      it != _numExpObs.end(); it++) {
1090    t_eph* eph = 0;
1091    for (int ie = 0; ie < _ephs.size(); ie++) {
1092      if (_ephs[ie]->prn().system() == it.key().system() &&
1093          _ephs[ie]->prn().number() == it.key().number()) {
1094        eph = _ephs[ie];
1095        break;
1096      }
1097    }
1098    if (eph) {
1099      int numExpObs = 0;
1100      bncTime epoTime;
1101      for (epoTime = startTime - interval; epoTime < endTime;
1102           epoTime = epoTime + interval) {
1103        ColumnVector xc(4);
1104        ColumnVector vv(3);
1105        if ( xyzSta.size() == 3 && (xyzSta[0] != 0.0 || xyzSta[1] != 0.0 || xyzSta[2] != 0.0) &&
1106             eph->getCrd(epoTime, xc, vv, false) == success) {
1107          double rho, eleSat, azSat;
1108          topos(xyzSta(1), xyzSta(2), xyzSta(3), xc(1), xc(2), xc(3), rho, eleSat, azSat);
1109          if ((eleSat * 180.0/M_PI) > 0.0) {
1110            numExpObs++;
1111          }
1112        }
1113      }
1114      it.value() = numExpObs;
1115    }
1116    else {
1117      if (!_navFileIncomplete.contains(it.key().system())) {
1118        _navFileIncomplete.append(it.key().system());
1119      }
1120    }
1121  }
1122}
Note: See TracBrowser for help on using the repository browser.