source: ntrip/branches/BNC_2.12/src/rinex/reqcanalyze.cpp @ 7908

Last change on this file since 7908 was 7908, checked in by stuerze, 2 years ago

possible obs added in rinex QC

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