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

Last change on this file since 10615 was 10614, checked in by stuerze, 4 days ago

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