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

Last change on this file since 8483 was 8483, checked in by stuerze, 11 months ago

SSR parameter clock rate, clock drift and URA are added within RTNET format

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" << "I:5&9";
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].toLatin1().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(7);
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].toLatin1());
403        frqType2.push_back(sys);
404        frqType2.push_back(_signalTypes[sys][1][0].toLatin1());
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().toLatin1();
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].toLatin1();
728      char frqChar2 = _signalTypes[prn.system()][1][0].toLatin1();
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 << " sec"                  << 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, false);
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(7);
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.