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

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

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