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

Last change on this file since 8561 was 8561, checked in by mervart, 5 years ago

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