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

Last change on this file since 10605 was 10605, checked in by stuerze, 6 hours ago

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