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

Last change on this file since 10945 was 10945, checked in by stuerze, 8 weeks ago

Added a Minimum Elevation parameter to BNC's RINEX Editing & QC feature
File size: 33.2 KB
Line 
1// Part of BNC, a utility for retrieving decoding and
2// converting GNSS data streams from NTRIP broadcasters.
3//
4// Copyright (C) 2007
5// German Federal Agency for Cartography and Geodesy (BKG)
6// http://www.bkg.bund.de
7// Czech Technical University Prague, Department of Geodesy
8// http://www.fsv.cvut.cz
9//
10// Email: euref-ip@bkg.bund.de
11//
12// This program is free software; you can redistribute it and/or
13// modify it under the terms of the GNU General Public License
14// as published by the Free Software Foundation, version 2.
15//
16// This program is distributed in the hope that it will be useful,
17// but WITHOUT ANY WARRANTY; without even the implied warranty of
18// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19// GNU General Public License for more details.
20//
21// You should have received a copy of the GNU General Public License
22// along with this program; if not, write to the Free Software
23// Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
24
25/* -------------------------------------------------------------------------
26 * BKG NTRIP Client
27 * -------------------------------------------------------------------------
28 *
29 * Class: t_reqcAnalyze
30 *
31 * Purpose: Analyze RINEX Files
32 *
33 * Author: L. Mervart
34 *
35 * Created: 11-Apr-2012
36 *
37 * Changes:
38 *
39 * -----------------------------------------------------------------------*/
40
41#include <iostream>
42#include <iomanip>
43#include <qwt_plot_renderer.h>
44
45#include "reqcanalyze.h"
46#include "bnccore.h"
47#include "bncsettings.h"
48#include "reqcedit.h"
49#include "bncutils.h"
50#include "graphwin.h"
51#include "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 _logStream = 0;
67 _currEpo = 0;
68 _obsFileNames = settings.value("reqcObsFile").toString().split(",", Qt::SkipEmptyParts);
69 _navFileNames = settings.value("reqcNavFile").toString().split(",", Qt::SkipEmptyParts);
70 _minEle = settings.value("reqcMinEle").toDouble();
71 _reqcPlotSignals = settings.value("reqcSkyPlotSignals").toString();
72 _defaultSignalTypes << "G:1&2&5" << "R:1&2&3" << "J:1&2" << "E:1&5" << "S:1&5" << "C:2&6" << "I:5&9";
73 if (_reqcPlotSignals.isEmpty()) {
74 _reqcPlotSignals = _defaultSignalTypes.join(" ");
75 }
76 analyzePlotSignals();
77
78 _checkEph = true;
79
80 qRegisterMetaType< QVector<t_skyPlotData> >("QVector<t_skyPlotData>");
81
82 connect(this, SIGNAL(dspSkyPlot(const QString&, QVector<t_skyPlotData>, const QByteArray&, double)),
83 this, SLOT(slotDspSkyPlot(const QString&, QVector<t_skyPlotData>, 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 _logStream; _logStream = 0;
99 delete _logFile; _logFile = 0;
100 if (BNC_CORE->mode() != t_bncCore::interactive) {
101 qApp->exit(8);
102 msleep(100); //sleep 0.1 sec
103 }
104}
105
106//
107////////////////////////////////////////////////////////////////////////////
108void t_reqcAnalyze::run() {
109
110 static const double QC_FORMAT_VERSION = 1.1;
111
112 // Open Log File
113 // -------------
114 if (!_logFileName.isEmpty()) {
115 _logFile = new QFile(_logFileName);
116 if (_logFile->open(QIODevice::WriteOnly | QIODevice::Text)) {
117 _logStream = new QTextStream();
118 _logStream->setDevice(_logFile);
119 }
120 }
121
122 if (_logStream) {
123 *_logStream << "QC Format Version : " << QString("%1").arg(QC_FORMAT_VERSION,3,'f',1) << Qt::endl << Qt::endl;
124 }
125
126 // Initialize RINEX Observation Files
127 // ----------------------------------
128 t_reqcEdit::initRnxObsFiles(_obsFileNames, _rnxObsFiles, _logStream);
129
130 // Read Ephemerides
131 // ----------------
132 t_reqcEdit::readEphemerides(_navFileNames, _ephs, _logStream, _checkEph);
133
134 // Loop over all RINEX Files
135 // -------------------------
136 for (int ii = 0; ii < _rnxObsFiles.size(); ii++) {
137 analyzeFile(_rnxObsFiles[ii]);
138 }
139
140 // Exit
141 // ----
142 emit finished();
143 deleteLater();
144}
145
146//
147////////////////////////////////////////////////////////////////////////////
148void t_reqcAnalyze::analyzePlotSignals() {
149
150 QStringList signalsOpt = _reqcPlotSignals.split(" ", Qt::SkipEmptyParts);
151
152 for (int ii = 0; ii < signalsOpt.size(); ii++) {
153 QStringList input = signalsOpt.at(ii).split(QRegExp("[:&]"), Qt::SkipEmptyParts);
154 if (input.size() > 1 && input[0].length() == 1) {
155 char system = input[0].toLatin1().constData()[0];
156 QStringList sysValid = _defaultSignalTypes.filter(QString(system));
157 if (!sysValid.isEmpty()) {
158 for (int iSig = 1; iSig < input.size(); iSig++) {
159 if (input[iSig].length() == 1 && input[iSig][0].isDigit()) {
160 _signalTypes[system].append(input[iSig][0].toLatin1());
161 }
162 }
163 }
164 }
165 }
166}
167
168//
169////////////////////////////////////////////////////////////////////////////
170void t_reqcAnalyze::analyzeFile(t_rnxObsFile* obsFile) {
171
172 _qcFile.clear();
173
174 // A priori Coordinates
175 // --------------------
176 ColumnVector xyzSta = obsFile->xyz();
177
178 // Loop over all Epochs
179 // --------------------
180 try {
181 QMap<QString, bncTime> lastObsTime;
182 bool firstEpo = true;
183 while ( (_currEpo = obsFile->nextEpoch()) != 0) {
184 if (firstEpo) {
185 firstEpo = false;
186 _qcFile._startTime = _currEpo->tt;
187 _qcFile._antennaName = obsFile->antennaName();
188 _qcFile._markerName = obsFile->markerName();
189 _qcFile._receiverType = obsFile->receiverType();
190 _qcFile._interval = obsFile->interval();
191 }
192 _qcFile._endTime = _currEpo->tt;
193
194 t_qcEpo qcEpo;
195 qcEpo._epoTime = _currEpo->tt;
196 qcEpo._PDOP = cmpDOP(xyzSta);
197
198 // Loop over all satellites
199 // ------------------------
200 for (unsigned iObs = 0; iObs < _currEpo->rnxSat.size(); iObs++) {
201 const t_rnxObsFile::t_rnxSat& rnxSat = _currEpo->rnxSat[iObs];
202 if (_navFileNames.size() && _numExpObs.find(rnxSat.prn) == _numExpObs.end()) {
203 _numExpObs[rnxSat.prn] = 0;
204 }
205 if (_signalTypes.find(rnxSat.prn.system()) == _signalTypes.end()) {
206 continue;
207 }
208 t_satObs satObs;
209 t_rnxObsFile::setObsFromRnx(obsFile, _currEpo, rnxSat, satObs);
210 t_qcSat& qcSat = qcEpo._qcSat[satObs._prn];
211 setQcObs(qcEpo._epoTime, xyzSta, satObs, lastObsTime, qcSat);
212 if (qcSat._eleSet && qcSat._eleDeg < _minEle) {
213 qcEpo._qcSat.remove(satObs._prn);
214 continue;
215 }
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 (_logStream) {
233 *_logStream << "Exception " << str << Qt::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.NormFrobenius();
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() == satObs._prn) {
332 eph = _ephs[ie];
333 break;
334 }
335 }
336 if (eph) {
337 ColumnVector xc(6);
338 ColumnVector vv(3);
339 if ( xyzSta.size() == 3 && (xyzSta[0] != 0.0 || xyzSta[1] != 0.0 || xyzSta[2] != 0.0) &&
340 eph->getCrd(epoTime, xc, vv, false) == success) {
341 double rho, eleSat, azSat;
342 topos(xyzSta(1), xyzSta(2), xyzSta(3), xc(1), xc(2), xc(3), rho, eleSat, azSat);
343 if (round(eleSat * 180.0/M_PI) >= 0.0) {
344 qcSat._eleSet = true;
345 qcSat._azDeg = azSat * 180.0/M_PI;
346 qcSat._eleDeg = eleSat * 180.0/M_PI;
347 }
348 }
349 if (satObs._prn.system() == 'R') {
350 qcSat._slotSet = true;
351 qcSat._slotNum = eph->slotNum();
352 }
353 }
354
355 // Availability and Slip Flags
356 // ---------------------------
357 for (unsigned ii = 0; ii < satObs._obs.size(); ii++) {
358 const t_frqObs* frqObs = satObs._obs[ii];
359
360 qcSat._qcFrq.push_back(t_qcFrq());
361 t_qcFrq& qcFrq = qcSat._qcFrq.back();
362
363 qcFrq._rnxType2ch = QString(frqObs->_rnxType2ch.c_str());
364 qcFrq._SNR = frqObs->_snr;
365 qcFrq._slip = frqObs->_slip;
366 qcFrq._phaseValid = frqObs->_phaseValid;
367 qcFrq._codeValid = frqObs->_codeValid;
368 // Check Gaps
369 // ----------
370 QString key = QString(satObs._prn.toString().c_str()) + qcFrq._rnxType2ch;
371 if (lastObsTime[key].valid()) {
372 double dt = epoTime - lastObsTime[key];
373 if (dt > 1.5 * _qcFile._interval) {
374 qcFrq._gap = true;
375 }
376 }
377 lastObsTime[key] = epoTime;
378
379 // Compute the Multipath Linear Combination
380 // ----------------------------------------
381 if (frqObs->_codeValid) {
382 t_frequency::type fA = t_frequency::dummy;
383 t_frequency::type fB = t_frequency::dummy;
384 char sys = satObs._prn.system();
385 if (_signalTypes.find(sys) != _signalTypes.end()) {
386 for (int iSig = 0; iSig < _signalTypes[sys].size(); iSig++) {
387 if (frqObs->_rnxType2ch[0] == _signalTypes[sys][iSig]) {
388 string frqType; frqType.push_back(sys); frqType.push_back(_signalTypes[sys][iSig]);
389 fA = t_frequency::toInt(frqType);
390 break;
391 }
392 }
393 if (fA != t_frequency::dummy) {
394 for (int iSig = 0; iSig < _signalTypes[sys].size(); iSig++) {
395 string frqType; frqType.push_back(sys); frqType.push_back(_signalTypes[sys][iSig]);
396 t_frequency::type fHlp = t_frequency::toInt(frqType);
397 if (fA != fHlp) {
398 fB = fHlp;
399 break;
400 }
401 }
402 }
403 if (fA != t_frequency::dummy && fB != t_frequency::dummy) {
404
405 double f_a = t_CST::freq(fA, qcSat._slotNum);
406 double f_b = t_CST::freq(fB, qcSat._slotNum);
407 double C_a = frqObs->_code;
408
409 bool foundA = false;
410 double L_a = 0.0;
411 bool foundB = false;
412 double L_b = 0.0;
413 for (unsigned jj = 0; jj < satObs._obs.size(); jj++) {
414 const t_frqObs* frqObsHlp = satObs._obs[jj];
415 if (frqObsHlp->_rnxType2ch[0] == t_frequency::toString(fA)[1] && frqObsHlp->_phaseValid) {
416 foundA = true;
417 L_a = frqObsHlp->_phase * t_CST::c / f_a;
418 }
419 else if (frqObsHlp->_rnxType2ch[0] == t_frequency::toString(fB)[1] && frqObsHlp->_phaseValid) {
420 foundB = true;
421 L_b = frqObsHlp->_phase * t_CST::c / f_b;
422 }
423 }
424 if (foundA && foundB && C_a) {
425 qcFrq._setMP = true;
426 qcFrq._rawMP = C_a - L_a - 2.0*f_b*f_b/(f_a*f_a-f_b*f_b) * (L_a - L_b);
427 }
428 }
429 }
430 }
431 } // satObs loop
432}
433
434//
435////////////////////////////////////////////////////////////////////////////
436void t_reqcAnalyze::analyzeMultipath() {
437
438 const double SLIPTRESH = 10.0; // cycle-slip threshold (meters)
439 const double chunkStep = 600.0; // 10 minutes
440
441 // Loop over all satellites available
442 // ----------------------------------
443 QMutableMapIterator<t_prn, t_qcSatSum> itSat(_qcFile._qcSatSum);
444 while (itSat.hasNext()) {
445 itSat.next();
446 const t_prn& prn = itSat.key();
447 t_qcSatSum& qcSatSum = itSat.value();
448
449 // Loop over all frequencies available
450 // -----------------------------------
451 QMutableMapIterator<QString, t_qcFrqSum> itFrq(qcSatSum._qcFrqSum);
452 while (itFrq.hasNext()) {
453 itFrq.next();
454 const QString& frqType = itFrq.key();
455 t_qcFrqSum& qcFrqSum = itFrq.value();
456
457 // Loop over all Chunks of Data
458 // ----------------------------
459 for (bncTime chunkStart = _qcFile._startTime;
460 chunkStart < _qcFile._endTime; chunkStart += chunkStep) {
461
462 bncTime chunkEnd = chunkStart + chunkStep;
463
464 QVector<t_qcFrq*> frqVec;
465 QVector<double> MP;
466
467 // Loop over all Epochs within one Chunk of Data
468 // ---------------------------------------------
469 for (int iEpo = 0; iEpo < _qcFile._qcEpo.size(); iEpo++) {
470 t_qcEpo& qcEpo = _qcFile._qcEpo[iEpo];
471 if (chunkStart <= qcEpo._epoTime && qcEpo._epoTime < chunkEnd) {
472 if (qcEpo._qcSat.contains(prn)) {
473 t_qcSat& qcSat = qcEpo._qcSat[prn];
474 for (int iFrq = 0; iFrq < qcSat._qcFrq.size(); iFrq++) {
475 t_qcFrq& qcFrq = qcSat._qcFrq[iFrq];
476 if (qcFrq._rnxType2ch == frqType) {
477 frqVec << &qcFrq;
478 if (qcFrq._setMP) {
479 MP << qcFrq._rawMP;
480 }
481 }
482 }
483 }
484 }
485 }
486
487 // Compute the multipath mean and standard deviation
488 // -------------------------------------------------
489 if (MP.size() > 1) {
490 double meanMP = 0.0;
491 for (int ii = 0; ii < MP.size(); ii++) {
492 meanMP += MP[ii];
493 }
494 meanMP /= MP.size();
495
496 bool slipMP = false;
497
498 double stdMP = 0.0;
499 for (int ii = 0; ii < MP.size(); ii++) {
500 double diff = MP[ii] - meanMP;
501 if (fabs(diff) > SLIPTRESH) {
502 slipMP = true;
503 break;
504 }
505 stdMP += diff * diff;
506 }
507
508 if (slipMP) {
509 stdMP = 0.0;
510 stdMP = 0.0;
511 qcFrqSum._numSlipsFound += 1;
512 }
513 else {
514 stdMP = sqrt(stdMP / (MP.size()-1));
515 qcFrqSum._numMP += 1;
516 qcFrqSum._sumMP += stdMP;
517 }
518
519 for (int ii = 0; ii < frqVec.size(); ii++) {
520 t_qcFrq* qcFrq = frqVec[ii];
521 if (slipMP) {
522 qcFrq->_slip = true;
523 }
524 else {
525 qcFrq->_stdMP = stdMP;
526 }
527 }
528 }
529 } // chunk loop
530 } // frq loop
531 } // sat loop
532}
533
534// Auxiliary Function for sorting MP data
535////////////////////////////////////////////////////////////////////////////
536bool t_reqcAnalyze::mpLessThan(const t_polarPoint* p1, const t_polarPoint* p2) {
537 return p1->_value < p2->_value;
538}
539
540//
541////////////////////////////////////////////////////////////////////////////
542void t_reqcAnalyze::preparePlotData(const t_rnxObsFile* obsFile) {
543
544 if (!BNC_CORE->GUIenabled()) {
545 return ;
546 }
547
548 QVector<t_skyPlotData> skyPlotDataMP;
549 QVector<t_skyPlotData> skyPlotDataSN;
550
551 for(QMap<char, QVector<char> >::iterator it1 = _signalTypes.begin();
552 it1 != _signalTypes.end(); it1++) {
553
554 char sys = it1.key();
555
556 for (int ii = 0; ii < it1.value().size(); ii++) {
557
558 skyPlotDataMP.append(t_skyPlotData(sys)); t_skyPlotData& dataMP = skyPlotDataMP.last();
559 skyPlotDataSN.append(t_skyPlotData(sys)); t_skyPlotData& dataSN = skyPlotDataSN.last();
560
561 dataMP._title = "Multipath\n" + QString(it1.key()) + ":" + it1.value()[ii] + " ";
562 dataSN._title = "Signal-to-Noise Ratio\n" + QString(it1.key()) + ":" + it1.value()[ii] + " ";
563
564 // Loop over all observations
565 // --------------------------
566 for (int iEpo = 0; iEpo < _qcFile._qcEpo.size(); iEpo++) {
567 t_qcEpo& qcEpo = _qcFile._qcEpo[iEpo];
568 QMapIterator<t_prn, t_qcSat> it2(qcEpo._qcSat);
569 while (it2.hasNext()) {
570 it2.next();
571 const t_qcSat& qcSat = it2.value();
572 if (qcSat._eleSet && it2.key().system() == sys) {
573 for (int iFrq = 0; iFrq < qcSat._qcFrq.size(); iFrq++) {
574 const t_qcFrq& qcFrq = qcSat._qcFrq[iFrq];
575 if (QString(it1.value()[ii]) == qcFrq._rnxType2ch.left(1)) {
576 (*dataMP._data) << (new t_polarPoint(qcSat._azDeg, 90.0 - qcSat._eleDeg, qcFrq._stdMP));
577 (*dataSN._data) << (new t_polarPoint(qcSat._azDeg, 90.0 - qcSat._eleDeg, qcFrq._SNR));
578 }
579 }
580 }
581 }
582 }
583
584 // Sort MP data (make the largest values always visible)
585 // -----------------------------------------------------
586 std::stable_sort(dataMP._data->begin(), dataMP._data->end(), mpLessThan);
587 }
588 }
589
590 // Show the plots
591 // --------------
592 QFileInfo fileInfo(obsFile->fileName());
593 QByteArray title = fileInfo.fileName().toLatin1();
594 emit dspSkyPlot(obsFile->fileName(), skyPlotDataMP, "Meters", 1.0);
595 emit dspSkyPlot(obsFile->fileName(), skyPlotDataSN, "dbHz", 54.0);
596 emit dspAvailPlot(obsFile->fileName(), title);
597}
598
599//
600////////////////////////////////////////////////////////////////////////////
601void t_reqcAnalyze::slotDspSkyPlot(const QString& fileName,
602 QVector<t_skyPlotData> skyPlotData,
603 const QByteArray& scaleTitle, double maxValue) {
604
605 if (BNC_CORE->GUIenabled()) {
606
607 if (maxValue == 0.0) {
608 QVectorIterator<t_skyPlotData> it(skyPlotData);
609 while (it.hasNext()) {
610 const t_skyPlotData& plotData = it.next();
611 const QVector<t_polarPoint*>* data = plotData._data;
612 for (int ii = 0; ii < data->size(); ii++) {
613 double val = data->at(ii)->_value;
614 if (maxValue < val) {
615 maxValue = val;
616 }
617 }
618 }
619 }
620
621 QwtInterval scaleInterval(0.0, maxValue);
622
623 QVector<QWidget*> plots;
624 QVector<int> rows;
625 QMutableVectorIterator<t_skyPlotData> it(skyPlotData);
626 int iRow = 0;
627 char sys = ' ';
628 while (it.hasNext()) {
629 t_skyPlotData& plotData = it.next();
630 QVector<t_polarPoint*>* data = plotData._data;
631 QwtText title(plotData._title);
632 QFont font = title.font(); font.setPointSize(font.pointSize()-1); title.setFont(font);
633 t_polarPlot* plot = new t_polarPlot(title, scaleInterval, BNC_CORE->mainWindow());
634 plot->addCurve(data);
635 plots << plot;
636 if (sys != ' ' && sys != plotData._sys) {
637 iRow += 1;
638 }
639 sys = plotData._sys;
640 rows << iRow;
641 }
642
643 t_graphWin* graphWin = new t_graphWin(0, fileName, plots, false,
644 &scaleTitle, &scaleInterval, &rows);
645
646 graphWin->show();
647
648 bncSettings settings;
649 QString dirName = settings.value("reqcPlotDir").toString();
650 if (!dirName.isEmpty()) {
651 QByteArray ext = (scaleTitle == "Meters") ? "_M.png" : "_S.png";
652 graphWin->savePNG(dirName, ext);
653 }
654 }
655}
656
657//
658////////////////////////////////////////////////////////////////////////////
659void t_reqcAnalyze::slotDspAvailPlot(const QString& fileName, const QByteArray& title) {
660
661 t_plotData plotData;
662 QMap<t_prn, t_plotData> plotDataMap;
663
664 for (int ii = 0; ii < _qcFile._qcEpo.size(); ii++) {
665 const t_qcEpo& qcEpo = _qcFile._qcEpo[ii];
666 double mjdX24 = qcEpo._epoTime.mjddec() * 24.0;
667
668 plotData._mjdX24 << mjdX24;
669 plotData._PDOP << qcEpo._PDOP;
670 plotData._numSat << qcEpo._qcSat.size();
671
672 QMapIterator<t_prn, t_qcSat> it(qcEpo._qcSat);
673 while (it.hasNext()) {
674 it.next();
675 const t_prn& prn = it.key();
676 const t_qcSat& qcSat = it.value();
677
678 t_plotData& data = plotDataMap[prn];
679
680 if (qcSat._eleSet) {
681 data._mjdX24 << mjdX24;
682 data._eleDeg << qcSat._eleDeg;
683 }
684
685 for (int iSig = 0; iSig < _signalTypes[prn.system()].size(); iSig++) {
686 char frqChar = _signalTypes[prn.system()][iSig];
687 QString frqType;
688 for (int iFrq = 0; iFrq < qcSat._qcFrq.size(); iFrq++) {
689 const t_qcFrq& qcFrq = qcSat._qcFrq[iFrq];
690 if (qcFrq._rnxType2ch[0] == frqChar && frqType.isEmpty()) {
691 frqType = qcFrq._rnxType2ch;
692 }
693 if (qcFrq._rnxType2ch == frqType) {
694 t_plotData::t_hlpStatus& status = data._status[frqChar];
695 if (qcFrq._slip) {
696 status._slip << mjdX24;
697 }
698 else if (qcFrq._gap) {
699 status._gap << mjdX24;
700 }
701 else {
702 status._ok << mjdX24;
703 }
704 }
705 }
706 }
707 }
708 }
709
710 if (BNC_CORE->GUIenabled()) {
711 t_availPlot* plotA = new t_availPlot(0, plotDataMap);
712 plotA->setTitle(title);
713
714 t_elePlot* plotZ = new t_elePlot(0, plotDataMap);
715
716 t_dopPlot* plotD = new t_dopPlot(0, plotData);
717
718 QVector<QWidget*> plots;
719 plots << plotA << plotZ << plotD;
720 t_graphWin* graphWin = new t_graphWin(0, fileName, plots, true);
721
722 int ww = QFontMetrics(graphWin->font()).horizontalAdvance('w');
723 graphWin->setMinimumSize(120*ww, 40*ww);
724
725 graphWin->show();
726
727 bncSettings settings;
728 QString dirName = settings.value("reqcPlotDir").toString();
729 if (!dirName.isEmpty()) {
730 QByteArray ext = "_A.png";
731 graphWin->savePNG(dirName, ext);
732 }
733 }
734}
735
736// Finish the report
737////////////////////////////////////////////////////////////////////////////
738void t_reqcAnalyze::printReport(const t_rnxObsFile* obsFile) {
739
740 if (!_logStream) {
741 return;
742 }
743
744 QFileInfo obsFi(obsFile->fileName());
745 QString obsFileName = obsFi.fileName();
746
747 // Summary
748 // -------
749 *_logStream << "Observation File : " << obsFileName << Qt::endl
750 << "RINEX Version : " << QString("%1").arg(obsFile->version(),4,'f',2) << Qt::endl
751 << "Marker Name : " << _qcFile._markerName << Qt::endl
752 << "Marker Number : " << obsFile->markerNumber() << Qt::endl
753 << "Receiver : " << _qcFile._receiverType << Qt::endl
754 << "Antenna : " << _qcFile._antennaName << Qt::endl
755 << "Position XYZ : " << QString("%1 %2 %3").arg(obsFile->xyz()(1), 14, 'f', 4)
756 .arg(obsFile->xyz()(2), 14, 'f', 4)
757 .arg(obsFile->xyz()(3), 14, 'f', 4) << Qt::endl
758 << "Antenna dH/dE/dN : " << QString("%1 %2 %3").arg(obsFile->antNEU()(3), 8, 'f', 4)
759 .arg(obsFile->antNEU()(2), 8, 'f', 4)
760 .arg(obsFile->antNEU()(1), 8, 'f', 4) << Qt::endl
761 << "Start Time : " << _qcFile._startTime.datestr().c_str() << ' '
762 << _qcFile._startTime.timestr(1,'.').c_str() << Qt::endl
763 << "End Time : " << _qcFile._endTime.datestr().c_str() << ' '
764 << _qcFile._endTime.timestr(1,'.').c_str() << Qt::endl
765 << "Interval : " << _qcFile._interval << " sec" << Qt::endl;
766
767
768 // Observation types per system
769 // -----------------------------
770 for (int iSys = 0; iSys < obsFile->numSys(); iSys++) {
771 char sys = obsFile->system(iSys);
772 if (sys != ' ') {
773 *_logStream << "Observation Types " << sys << ":";
774 for (int iType = 0; iType < obsFile->nTypes(sys); iType++) {
775 QString type = obsFile->obsType(sys, iType);
776 *_logStream << " " << type;
777 }
778 *_logStream << Qt::endl;
779 }
780 }
781
782 // Number of analysed systems
783 // --------------------------
784 QMap<QChar, QVector<const t_qcSatSum*> > systemMap;
785 QMapIterator<t_prn, t_qcSatSum> itSat(_qcFile._qcSatSum);
786 while (itSat.hasNext()) {
787 itSat.next();
788 const t_prn& prn = itSat.key();
789 const t_qcSatSum& qcSatSum = itSat.value();
790 systemMap[prn.system()].push_back(&qcSatSum);
791 }
792 *_logStream << "Analysed GNSS : " << systemMap.size() << " ";
793 QMapIterator<QChar, QVector<const t_qcSatSum*> > itSys(systemMap);
794 while (itSys.hasNext()) {
795 itSys.next();
796 *_logStream << ' ' << itSys.key();
797 }
798 *_logStream << Qt::endl;
799
800
801 // System specific summary
802 // -----------------------
803 itSys.toFront();
804 while (itSys.hasNext()) {
805 itSys.next();
806 const QChar& sys = itSys.key();
807 const QVector<const t_qcSatSum*>& qcSatVec = itSys.value();
808 int numExpectedObs = 0;
809 for(QMap<t_prn, int>::iterator it = _numExpObs.begin(); it != _numExpObs.end(); it++) {
810 if (sys == it.key().system()) {
811 numExpectedObs += it.value();
812 }
813 }
814 QString prefixSys = QString(" ") + sys + QString(": ");
815 QMap<QString, QVector<const t_qcFrqSum*> > frqMap;
816 for (int ii = 0; ii < qcSatVec.size(); ii++) {
817 const t_qcSatSum* qcSatSum = qcSatVec[ii];
818 QMapIterator<QString, t_qcFrqSum> itFrq(qcSatSum->_qcFrqSum);
819 while (itFrq.hasNext()) {
820 itFrq.next();
821 QString frqType = itFrq.key(); if (frqType.length() < 2) frqType += '?';
822 const t_qcFrqSum& qcFrqSum = itFrq.value();
823 frqMap[frqType].push_back(&qcFrqSum);
824 }
825 }
826 *_logStream << Qt::endl
827 << prefixSys << "Satellites: " << qcSatVec.size() << Qt::endl
828 << prefixSys << "Signals : " << frqMap.size() << " ";
829 QMapIterator<QString, QVector<const t_qcFrqSum*> > itFrq(frqMap);
830 while (itFrq.hasNext()) {
831 itFrq.next();
832 QString frqType = itFrq.key(); if (frqType.length() < 2) frqType += '?';
833 *_logStream << ' ' << frqType;
834 }
835 *_logStream << Qt::endl;
836 QString prefixSys2 = " " + prefixSys;
837 itFrq.toFront();
838 while (itFrq.hasNext()) {
839 itFrq.next();
840 QString frqType = itFrq.key(); if (frqType.length() < 2) frqType += '?';
841 const QVector<const t_qcFrqSum*> qcFrqVec = itFrq.value();
842 QString prefixFrq = QString(" ") + frqType + QString(": ");
843 int numObs = 0;
844 int numSlipsFlagged = 0;
845 int numSlipsFound = 0;
846 int numGaps = 0;
847 int numSNR = 0;
848 int numMP = 0;
849 double sumSNR = 0.0;
850 double sumMP = 0.0;
851 for (int ii = 0; ii < qcFrqVec.size(); ii++) {
852 const t_qcFrqSum* qcFrqSum = qcFrqVec[ii];
853 numObs += qcFrqSum->_numObs ;
854 numSlipsFlagged += qcFrqSum->_numSlipsFlagged;
855 numSlipsFound += qcFrqSum->_numSlipsFound ;
856 numGaps += qcFrqSum->_numGaps ;
857 numSNR += qcFrqSum->_numSNR;
858 numMP += qcFrqSum->_numMP;
859 sumSNR += qcFrqSum->_sumSNR;
860 sumMP += qcFrqSum->_sumMP;
861 }
862 if (numSNR > 0) {
863 sumSNR /= numSNR;
864 }
865 if (numMP > 0) {
866 sumMP /= numMP;
867 }
868
869 double ratio = (double(numObs) / double(numExpectedObs)) * 100.0;
870
871 *_logStream << Qt::endl
872 << prefixSys2 << prefixFrq << "Observations : ";
873 if(_navFileNames.isEmpty() || numExpectedObs == 0.0 || _navFileIncomplete.contains(sys.toLatin1())) {
874 *_logStream << QString("%1\n").arg(numObs, 6);
875 }
876 else {
877 *_logStream << QString("%1 (%2) %3 \%\n").arg(numObs, 6).arg(numExpectedObs, 8).arg(ratio, 8, 'f', 2);
878 }
879 *_logStream << prefixSys2 << prefixFrq << "Slips (file+found): " << QString("%1 +").arg(numSlipsFlagged, 8)
880 << QString("%1\n").arg(numSlipsFound, 8)
881 << prefixSys2 << prefixFrq << "Gaps : " << QString("%1\n").arg(numGaps, 8)
882 << prefixSys2 << prefixFrq << "Mean SNR : " << QString("%1\n").arg(sumSNR, 8, 'f', 1)
883 << prefixSys2 << prefixFrq << "Mean Multipath : " << QString("%1\n").arg(sumMP, 8, 'f', 2);
884 }
885 }
886
887 // Epoch-Specific Output
888 // ---------------------
889 bncSettings settings;
890 if (Qt::CheckState(settings.value("reqcLogSummaryOnly").toInt()) == Qt::Checked) {
891 return;
892 }
893 *_logStream << Qt::endl;
894 for (int iEpo = 0; iEpo < _qcFile._qcEpo.size(); iEpo++) {
895 const t_qcEpo& qcEpo = _qcFile._qcEpo[iEpo];
896
897 unsigned year, month, day, hour, min;
898 double sec;
899 qcEpo._epoTime.civil_date(year, month, day);
900 qcEpo._epoTime.civil_time(hour, min, sec);
901
902 QString dateStr;
903 QTextStream(&dateStr) << QString("> %1 %2 %3 %4 %5%6")
904 .arg(year, 4)
905 .arg(month, 2, 10, QChar('0'))
906 .arg(day, 2, 10, QChar('0'))
907 .arg(hour, 2, 10, QChar('0'))
908 .arg(min, 2, 10, QChar('0'))
909 .arg(sec, 11, 'f', 7);
910
911 *_logStream << dateStr << QString(" %1").arg(qcEpo._qcSat.size(), 2)
912 << QString(" %1").arg(qcEpo._PDOP, 4, 'f', 1)
913 << Qt::endl;
914
915 QMapIterator<t_prn, t_qcSat> itSat(qcEpo._qcSat);
916 while (itSat.hasNext()) {
917 itSat.next();
918 const t_prn& prn = itSat.key();
919 const t_qcSat& qcSat = itSat.value();
920
921 *_logStream << prn.toString().c_str()
922 << QString(" %1 %2").arg(qcSat._eleDeg, 6, 'f', 2).arg(qcSat._azDeg, 7, 'f', 2);
923
924 int numObsTypes = 0;
925 for (int iFrq = 0; iFrq < qcSat._qcFrq.size(); iFrq++) {
926 const t_qcFrq& qcFrq = qcSat._qcFrq[iFrq];
927 if (qcFrq._phaseValid) {
928 numObsTypes += 1;
929 }
930 if (qcFrq._codeValid) {
931 numObsTypes += 1;
932 }
933 }
934 *_logStream << QString(" %1").arg(numObsTypes, 2);
935
936 for (int iFrq = 0; iFrq < qcSat._qcFrq.size(); iFrq++) {
937 const t_qcFrq& qcFrq = qcSat._qcFrq[iFrq];
938 if (qcFrq._phaseValid) {
939 *_logStream << " L" << qcFrq._rnxType2ch << ' ';
940 if (qcFrq._slip) {
941 *_logStream << 's';
942 }
943 else {
944 *_logStream << '.';
945 }
946 if (qcFrq._gap) {
947 *_logStream << 'g';
948 }
949 else {
950 *_logStream << '.';
951 }
952 *_logStream << QString(" %1").arg(qcFrq._SNR, 4, 'f', 1);
953 }
954 if (qcFrq._codeValid) {
955 *_logStream << " C" << qcFrq._rnxType2ch << ' ';
956 if (qcFrq._gap) {
957 *_logStream << " g";
958 }
959 else {
960 *_logStream << " .";
961 }
962 *_logStream << QString(" %1").arg(qcFrq._stdMP, 3, 'f', 2);
963 }
964 }
965 *_logStream << Qt::endl;
966 }
967 }
968 _logStream->flush();
969}
970
971//
972////////////////////////////////////////////////////////////////////////////
973void t_reqcAnalyze::setExpectedObs(const bncTime& startTime, const bncTime& endTime,
974 double interval, const ColumnVector& xyzSta) {
975 for(QMap<t_prn, int>::iterator it = _numExpObs.begin(); it != _numExpObs.end(); it++) {
976 t_eph* eph = 0;
977 for (int ie = 0; ie < _ephs.size(); ie++) {
978 if (_ephs[ie]->prn() == it.key()) {
979 eph = _ephs[ie];
980 break;
981 }
982 }
983 if (eph) {
984 int numExpObs = 0;
985 bncTime epoTime;
986 for (epoTime = startTime - interval; epoTime < endTime; epoTime = epoTime + interval) {
987 ColumnVector xc(6);
988 ColumnVector vv(3);
989 if ( xyzSta.size() == 3 && (xyzSta[0] != 0.0 || xyzSta[1] != 0.0 || xyzSta[2] != 0.0) &&
990 eph->getCrd(epoTime, xc, vv, false) == success) {
991 double rho, eleSat, azSat;
992 topos(xyzSta(1), xyzSta(2), xyzSta(3), xc(1), xc(2), xc(3), rho, eleSat, azSat);
993 if (round(eleSat * 180.0/M_PI) >= _minEle) {
994 numExpObs++;
995 }
996 }
997 }
998 it.value() = numExpObs;
999 }
1000 else {
1001 if (!_navFileIncomplete.contains(it.key().system())) {
1002 _navFileIncomplete.append(it.key().system());
1003 }
1004 }
1005 }
1006}
Note: See TracBrowser for help on using the repository browser.