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

Last change on this file since 8561 was 8561, checked in by mervart, 19 months ago

Analyze more than two signals

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