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

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

rinex qc: satellite visibility is considered now if expected obs is computed

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