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

Last change on this file since 11057 was 11057, checked in by stuerze, 2 weeks ago

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