source: ntrip/branches/BNC_2.12/src/rinex/reqcanalyze.cpp@ 7920

Last change on this file since 7920 was 7920, checked in by stuerze, 8 years ago

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