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

Last change on this file since 4575 was 4575, checked in by mervart, 12 years ago
File size: 15.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
44#include <qwt_plot.h>
45#include <qwt_plot_curve.h>
46
47#include "reqcanalyze.h"
48#include "bncapp.h"
49#include "bncsettings.h"
50#include "reqcedit.h"
51#include "bncutils.h"
52#include "bncpostprocess.h"
53#include "graphwin.h"
54#include "polarplot.h"
55
56using namespace std;
57
58const double SLIPTRESH = 5.0; // cycle-slip threshold (meters)
59
60// Constructor
61////////////////////////////////////////////////////////////////////////////
62t_reqcAnalyze::t_reqcAnalyze(QObject* parent) : QThread(parent) {
63
64 bncSettings settings;
65
66 _logFileName = settings.value("reqcOutLogFile").toString(); expandEnvVar(_logFileName);
67 _logFile = 0;
68 _log = 0;
69 _obsFileNames = settings.value("reqcObsFile").toString().split(",", QString::SkipEmptyParts);
70 _navFileNames = settings.value("reqcNavFile").toString().split(",", QString::SkipEmptyParts);
71
72 _currEpo = 0;
73
74 connect(this, SIGNAL(dspSkyPlot(const QString&,
75 const QByteArray&,
76 QVector<t_polarPoint*>*,
77 const QByteArray&,
78 QVector<t_polarPoint*>*,
79 const QByteArray&, double)),
80 this, SLOT(slotDspSkyPlot(const QString&,
81 const QByteArray&,
82 QVector<t_polarPoint*>*,
83 const QByteArray&,
84 QVector<t_polarPoint*>*,
85 const QByteArray&, double)));
86
87 connect(this, SIGNAL(dspAvailPlot(const QString&,
88 const QByteArray&,
89 QMap<QString, QVector<int> >*)),
90 this, SLOT(slotDspAvailPlot(const QString&,
91 const QByteArray&,
92 QMap<QString, QVector<int> >*)));
93
94}
95
96// Destructor
97////////////////////////////////////////////////////////////////////////////
98t_reqcAnalyze::~t_reqcAnalyze() {
99 for (int ii = 0; ii < _rnxObsFiles.size(); ii++) {
100 delete _rnxObsFiles[ii];
101 }
102 for (int ii = 0; ii < _ephs.size(); ii++) {
103 delete _ephs[ii];
104 }
105 delete _log; _log = 0;
106 delete _logFile; _logFile = 0;
107 bncApp* app = (bncApp*) qApp;
108 if ( app->mode() != bncApp::interactive) {
109 app->exit(0);
110 }
111}
112
113//
114////////////////////////////////////////////////////////////////////////////
115void t_reqcAnalyze::slotDspSkyPlot(const QString& fileName,
116 const QByteArray& title1,
117 QVector<t_polarPoint*>* data1,
118 const QByteArray& title2,
119 QVector<t_polarPoint*>* data2,
120 const QByteArray& scaleTitle,
121 double maxValue) {
122
123 bncApp* app = dynamic_cast<bncApp*>(qApp);
124 if (app->GUIenabled()) {
125
126 if (maxValue == 0.0) {
127 if (data1) {
128 for (int ii = 0; ii < data1->size(); ii++) {
129 double val = data1->at(ii)->_value;
130 if (maxValue < val) {
131 maxValue = val;
132 }
133 }
134 }
135 if (data2) {
136 for (int ii = 0; ii < data2->size(); ii++) {
137 double val = data2->at(ii)->_value;
138 if (maxValue < val) {
139 maxValue = val;
140 }
141 }
142 }
143 }
144
145 QwtInterval scaleInterval(0.0, maxValue);
146
147 QVector<QWidget*> plots;
148 if (data1) {
149 t_polarPlot* plot1 = new t_polarPlot(QwtText(title1), scaleInterval,
150 app->mainWindow());
151 plot1->addCurve(data1);
152 plots << plot1;
153 }
154 if (data2) {
155 t_polarPlot* plot2 = new t_polarPlot(QwtText(title2), scaleInterval,
156 app->mainWindow());
157 plot2->addCurve(data2);
158 plots << plot2;
159 }
160
161 t_graphWin* graphWin = new t_graphWin(0, fileName, plots,
162 &scaleTitle, &scaleInterval);
163
164 graphWin->show();
165
166 bncSettings settings;
167 QString dirName = settings.value("reqcPlotDir").toString();
168 if (!dirName.isEmpty()) {
169 QByteArray ext = scaleTitle.isEmpty() ? "_SNR.png" : "_MP.png";
170 graphWin->savePNG(dirName, ext);
171 }
172 }
173}
174
175//
176////////////////////////////////////////////////////////////////////////////
177void t_reqcAnalyze::run() {
178
179 // Open Log File
180 // -------------
181 _logFile = new QFile(_logFileName);
182 if (_logFile->open(QIODevice::WriteOnly | QIODevice::Text)) {
183 _log = new QTextStream();
184 _log->setDevice(_logFile);
185 }
186
187 // Initialize RINEX Observation Files
188 // ----------------------------------
189 t_reqcEdit::initRnxObsFiles(_obsFileNames, _rnxObsFiles, _log);
190
191 // Read Ephemerides
192 // ----------------
193 t_reqcEdit::readEphemerides(_navFileNames, _ephs);
194
195 // Loop over all RINEX Files
196 // -------------------------
197 for (int ii = 0; ii < _rnxObsFiles.size(); ii++) {
198 analyzeFile(_rnxObsFiles[ii]);
199 }
200
201 // Exit
202 // ----
203 emit finished();
204 deleteLater();
205}
206
207//
208////////////////////////////////////////////////////////////////////////////
209void t_reqcAnalyze::analyzeFile(t_rnxObsFile* obsFile) {
210
211 if (_log) {
212 *_log << "\nAnalyze File\n"
213 << "------------\n"
214 << obsFile->fileName().toAscii().data() << endl << endl;
215 }
216
217 _satStat.clear();
218
219 // A priori Coordinates
220 // --------------------
221 ColumnVector xyz = obsFile->xyz();
222
223 // Loop over all Epochs
224 // --------------------
225 try {
226 while ( (_currEpo = obsFile->nextEpoch()) != 0) {
227
228 // Loop over all satellites
229 // ------------------------
230 for (unsigned iObs = 0; iObs < _currEpo->rnxSat.size(); iObs++) {
231 const t_rnxObsFile::t_rnxSat& rnxSat = _currEpo->rnxSat[iObs];
232 t_obs obs;
233 t_postProcessing::setObsFromRnx(obsFile, _currEpo, rnxSat, obs);
234
235 if (obs.satSys == 'R') {
236 continue; // TODO: set channel number
237 }
238
239 QString prn = QString("%1%2").arg(obs.satSys)
240 .arg(obs.satNum, 2, 10, QChar('0'));
241
242 t_satStat& satStat = _satStat[prn];
243
244 satStat.addObs(obs);
245 }
246
247 } // while (_currEpo)
248 }
249 catch (QString str) {
250 if (_log) {
251 *_log << "Exception " << str << endl;
252 }
253 else {
254 qDebug() << str;
255 }
256 return;
257 }
258
259 // Analyze the Multipath
260 // ---------------------
261 QVector<t_polarPoint*>* dataMP1 = new QVector<t_polarPoint*>;
262 QVector<t_polarPoint*>* dataMP2 = new QVector<t_polarPoint*>;
263 QVector<t_polarPoint*>* dataSNR1 = new QVector<t_polarPoint*>;
264 QVector<t_polarPoint*>* dataSNR2 = new QVector<t_polarPoint*>;
265 QMap<QString, QVector<int> >* availL1 = new QMap<QString, QVector<int> >;
266
267
268 QMapIterator<QString, t_satStat> it(_satStat);
269 while (it.hasNext()) {
270 it.next();
271 QString prn = it.key();
272 const t_satStat& satStat = it.value();
273 QVector<int>& dataL1 = (*availL1)[prn];
274 preparePlotData(prn, satStat, xyz, obsFile->interval(),
275 dataMP1, dataMP2, dataSNR1, dataSNR2, dataL1);
276 }
277
278 emit dspSkyPlot(obsFile->fileName(), "MP1", dataMP1, "MP2", dataMP2,
279 "Meters", 2.0);
280
281 emit dspSkyPlot(obsFile->fileName(), "SNR1", dataSNR1, "SNR2", dataSNR2,
282 "", 9.0);
283
284 ///// emit dspAvailPlot(obsFile->fileName(), "Availability L1", availL1);
285
286 if (_log) {
287 _log->flush();
288 }
289}
290
291//
292////////////////////////////////////////////////////////////////////////////
293void t_reqcAnalyze::t_satStat::addObs(const t_obs& obs) {
294
295 t_anaObs* newObs = new t_anaObs(obs.GPSWeek, obs.GPSWeeks);
296 bool okFlag = false;
297
298 // Compute the Multipath
299 // ----------------------
300 double L1 = obs.measdata("L1", 3.0);
301 if (L1 != 0) {
302 newObs->_hasL1 = true;
303 }
304 double L2 = obs.measdata("L2", 3.0);
305 if (L2 != 0) {
306 newObs->_hasL2 = true;
307 }
308 if (L1 != 0.0 && L2 != 0.0) {
309 double f1 = t_CST::f1(obs.satSys, obs.slotNum);
310 double f2 = t_CST::f2(obs.satSys, obs.slotNum);
311
312 L1 = L1 * t_CST::c / f1;
313 L2 = L2 * t_CST::c / f2;
314
315 double P1 = obs.measdata("C1", 3.0);
316 if (P1 != 0.0) {
317 newObs->_MP1 = P1 - L1 - 2.0*f2*f2/(f1*f1-f2*f2) * (L1 - L2);
318 okFlag = true;
319 }
320 double P2 = obs.measdata("C2", 3.0);
321 if (P2 != 0.0) {
322 newObs->_MP2 = P2 - L2 - 2.0*f1*f1/(f1*f1-f2*f2) * (L1 - L2);
323 okFlag = true;
324 }
325 }
326
327 // Signal-to-Noise
328 // ---------------
329 double S1 = obs.measdata("S1", 3.0);
330 if (S1 != 0.0) {
331 newObs->_SNR1 = floor(S1/6);
332 if (newObs->_SNR1 > 9.0) {
333 newObs->_SNR1 = 9.0;
334 }
335 if (newObs->_SNR1 < 1.0) {
336 newObs->_SNR1 = 1.0;
337 }
338 okFlag = true;
339 }
340 else {
341 if (obs.snrL1 > 0) {
342 newObs->_SNR1 = obs.snrL1;
343 okFlag = true;
344 }
345 }
346 double S2 = obs.measdata("S2", 3.0);
347 if (S2 != 0.0) {
348 newObs->_SNR2 = floor(S2/6);
349 if (newObs->_SNR2 > 9.0) {
350 newObs->_SNR2 = 9.0;
351 }
352 if (newObs->_SNR2 < 1.0) {
353 newObs->_SNR2 = 1.0;
354 }
355 okFlag = true;
356 }
357 else {
358 if (obs.snrL2 > 0) {
359 newObs->_SNR2 = obs.snrL2;
360 okFlag = true;
361 }
362 }
363
364 // Remember the Observation
365 // ------------------------
366 if (okFlag) {
367 anaObs << newObs;
368 }
369 else {
370 delete newObs;
371 }
372}
373
374//
375////////////////////////////////////////////////////////////////////////////
376void t_reqcAnalyze::preparePlotData(const QString& prn,
377 const t_satStat& satStat,
378 const ColumnVector& xyz,
379 double obsInterval,
380 QVector<t_polarPoint*>* dataMP1,
381 QVector<t_polarPoint*>* dataMP2,
382 QVector<t_polarPoint*>* dataSNR1,
383 QVector<t_polarPoint*>* dataSNR2,
384 QVector<int>& dataL1) {
385
386 const int chunkStep = int( 30.0 / obsInterval); // chunk step (30 sec)
387 const int numEpo = int(600.0 / obsInterval); // # epochs in one chunk (10 min)
388
389 for (int chunkStart = 0; chunkStart + numEpo < satStat.anaObs.size();
390 chunkStart += chunkStep) {
391
392 bncTime firstEpoch;
393
394 // Compute Mean
395 // ------------
396 bool slipFlag = false;
397 double mean1 = 0.0;
398 double mean2 = 0.0;
399 double SNR1 = 0.0;
400 double SNR2 = 0.0;
401
402 for (int ii = 0; ii < numEpo; ii++) {
403 int iEpo = chunkStart + ii;
404 const t_anaObs* anaObs = satStat.anaObs[iEpo];
405
406 if (ii == 0) {
407 firstEpoch.set(anaObs->_GPSWeek, anaObs->_GPSWeeks);
408 }
409
410 mean1 += anaObs->_MP1;
411 mean2 += anaObs->_MP2;
412
413 if ( anaObs->_SNR1 > 0 && (SNR1 == 0 || SNR1 > anaObs->_SNR1) ) {
414 SNR1 = anaObs->_SNR1;
415 }
416 if ( anaObs->_SNR2 > 0 && (SNR2 == 0 || SNR2 > anaObs->_SNR2) ) {
417 SNR2 = anaObs->_SNR2;
418 }
419
420 // Check Slip
421 // ----------
422 if (ii > 0) {
423 double diff1 = anaObs->_MP1 - satStat.anaObs[iEpo-1]->_MP1;
424 double diff2 = anaObs->_MP2 - satStat.anaObs[iEpo-1]->_MP2;
425 if (fabs(diff1) > SLIPTRESH || fabs(diff2) > SLIPTRESH) {
426 slipFlag = true;
427 break;
428 }
429 }
430 }
431
432 dataL1 << int(firstEpoch.gpssec());
433
434 if (slipFlag) {
435 continue;
436 }
437
438 mean1 /= numEpo;
439 mean2 /= numEpo;
440
441 // Compute Standard Deviation
442 // --------------------------
443 double stddev1 = 0.0;
444 double stddev2 = 0.0;
445 for (int ii = 0; ii < numEpo; ii++) {
446 int iEpo = chunkStart + ii;
447 const t_anaObs* anaObs = satStat.anaObs[iEpo];
448 double diff1 = anaObs->_MP1 - mean1;
449 double diff2 = anaObs->_MP2 - mean2;
450 stddev1 += diff1 * diff1;
451 stddev2 += diff2 * diff2;
452 }
453 double MP1 = sqrt(stddev1 / (numEpo-1));
454 double MP2 = sqrt(stddev2 / (numEpo-1));
455
456 const t_anaObs* anaObs0 = satStat.anaObs[chunkStart];
457
458 // Compute the Azimuth and Zenith Distance
459 // ---------------------------------------
460 double az = 0.0;
461 double zen = 0.0;
462 if (xyz.size()) {
463 t_eph* eph = 0;
464 for (int ie = 0; ie < _ephs.size(); ie++) {
465 if (_ephs[ie]->prn() == prn) {
466 eph = _ephs[ie];
467 break;
468 }
469 }
470
471 if (eph) {
472 double xSat, ySat, zSat, clkSat;
473 eph->position(anaObs0->_GPSWeek, anaObs0->_GPSWeeks,
474 xSat, ySat, zSat, clkSat);
475
476 double rho, eleSat, azSat;
477 topos(xyz(1), xyz(2), xyz(3), xSat, ySat, zSat, rho, eleSat, azSat);
478
479 az = azSat * 180.0/M_PI;
480 zen = 90.0 - eleSat * 180.0/M_PI;
481 }
482 }
483
484 // Add new Point
485 // -------------
486 (*dataMP1) << (new t_polarPoint(az, zen, MP1));
487 (*dataMP2) << (new t_polarPoint(az, zen, MP2));
488 (*dataSNR1) << (new t_polarPoint(az, zen, SNR1));
489 (*dataSNR2) << (new t_polarPoint(az, zen, SNR2));
490
491 if (_log) {
492 _log->setRealNumberNotation(QTextStream::FixedNotation);
493
494 _log->setRealNumberPrecision(2);
495 *_log << "MP1 " << prn << " " << az << " " << zen << " ";
496 _log->setRealNumberPrecision(3);
497 *_log << MP1 << endl;
498
499 _log->setRealNumberPrecision(2);
500 *_log << "MP2 " << prn << " " << az << " " << zen << " ";
501 _log->setRealNumberPrecision(3);
502 *_log << MP2 << endl;
503
504 _log->flush();
505 }
506 }
507}
508
509//
510////////////////////////////////////////////////////////////////////////////
511void t_reqcAnalyze::slotDspAvailPlot(const QString& fileName,
512 const QByteArray& title,
513 QMap<QString, QVector<int> >* prnAvail){
514
515 if (dynamic_cast<bncApp*>(qApp)->GUIenabled()) {
516
517 QwtPlot* plot = new QwtPlot();
518
519 int iC = 0;
520 QMapIterator<QString, QVector<int> > it(*prnAvail);
521 while (it.hasNext()) {
522 it.next();
523 const QString& prn = it.key();
524 const QVector<int>& epochs = it.value();
525 }
526
527 QVector<QWidget*> plots;
528 plots << plot;
529 t_graphWin* graphWin = new t_graphWin(0, fileName, plots, 0, 0);
530 graphWin->show();
531
532 bncSettings settings;
533 QString dirName = settings.value("reqcPlotDir").toString();
534 if (!dirName.isEmpty()) {
535 QByteArray ext = "_AVAIL.png";
536 graphWin->savePNG(dirName, ext);
537 }
538 }
539}
Note: See TracBrowser for help on using the repository browser.