Changeset 6271 in ntrip


Ignore:
Timestamp:
Oct 31, 2014, 5:45:35 PM (9 years ago)
Author:
mervart
Message:
 
File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/BNC/src/rinex/reqcanalyze.cpp

    r6269 r6271  
    5656using namespace std;
    5757
    58 const double SLIPTRESH = 10.0; // cycle-slip threshold (meters)
    59 
    6058// Constructor
    6159////////////////////////////////////////////////////////////////////////////
    … …  
    6765  _logFile      = 0;
    6866  _log          = 0;
     67  _currEpo      = 0;
    6968  _obsFileNames = settings.value("reqcObsFile").toString().split(",", QString::SkipEmptyParts);
    7069  _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&, const QByteArray&)),
    88           this, SLOT(slotDspAvailPlot(const QString&, const QByteArray&)));
    8970}
    9071
    … …  
    10283  if (BNC_CORE->mode() != t_bncCore::interactive) {
    10384    qApp->exit(0);
    104   }
    105 }
    106 
    107 //
    108 ////////////////////////////////////////////////////////////////////////////
    109 void t_reqcAnalyze::slotDspSkyPlot(const QString& fileName,
    110                                    const QByteArray& title1,
    111                                    QVector<t_polarPoint*>* data1,
    112                                    const QByteArray& title2,
    113                                    QVector<t_polarPoint*>* data2,
    114                                    const QByteArray& scaleTitle,
    115                                    double maxValue) {
    116 
    117   if (BNC_CORE->GUIenabled()) {
    118 
    119     if (maxValue == 0.0) {
    120       if (data1) {
    121         for (int ii = 0; ii < data1->size(); ii++) {
    122           double val = data1->at(ii)->_value;
    123           if (maxValue < val) {
    124             maxValue = val;
    125           }
    126         }
    127       }
    128       if (data2) {
    129         for (int ii = 0; ii < data2->size(); ii++) {
    130           double val = data2->at(ii)->_value;
    131           if (maxValue < val) {
    132             maxValue = val;
    133           }
    134         }
    135       }
    136     }
    137 
    138     QwtInterval scaleInterval(0.0, maxValue);
    139 
    140     QVector<QWidget*> plots;
    141     if (data1) {
    142       t_polarPlot* plot1 = new t_polarPlot(QwtText(title1), scaleInterval,
    143                                           BNC_CORE->mainWindow());
    144       plot1->addCurve(data1);
    145       plots << plot1;
    146     }
    147     if (data2) {
    148       t_polarPlot* plot2 = new t_polarPlot(QwtText(title2), scaleInterval,
    149                                            BNC_CORE->mainWindow());
    150       plot2->addCurve(data2);
    151       plots << plot2;
    152     }
    153 
    154     t_graphWin* graphWin = new t_graphWin(0, fileName, plots,
    155                                           &scaleTitle, &scaleInterval);
    156 
    157     graphWin->show();
    158 
    159     bncSettings settings;
    160     QString dirName = settings.value("reqcPlotDir").toString();
    161     if (!dirName.isEmpty()) {
    162       QByteArray ext = (scaleTitle == "Meters") ? "_M.png" : "_S.png";
    163       graphWin->savePNG(dirName, ext);
    164     }
    16585  }
    16686}
    … …  
    202122void t_reqcAnalyze::analyzeFile(t_rnxObsFile* obsFile) {
    203123
    204   _mutex.lock();
     124  QMutexLocker lock(&_mutex);
    205125
    206126  if (_log) {
    … …  
    210130  }
    211131
    212   _allObsMap.clear();
    213   _availDataMap.clear();
    214   _obsStat.reset();
     132  _qcFile.clear();
    215133
    216134  // A priori Coordinates
    … …  
    221139  // --------------------
    222140  try {
    223     unsigned iEpo = 0;
     141    bool firstEpo = true;
    224142    while ( (_currEpo = obsFile->nextEpoch()) != 0) {
    225 
    226       if (iEpo == 0) {
    227         _obsStat._startTime    = _currEpo->tt;
    228         _obsStat._antennaName  = obsFile->antennaName();
    229         _obsStat._markerName   = obsFile->markerName();
    230         _obsStat._receiverType = obsFile->receiverType();
    231         _obsStat._interval     = obsFile->interval();
    232       }
    233       _obsStat._endTime = _currEpo->tt;
     143      if (firstEpo) {
     144        firstEpo = false;
     145        _qcFile._startTime    = _currEpo->tt;
     146        _qcFile._antennaName  = obsFile->antennaName();
     147        _qcFile._markerName   = obsFile->markerName();
     148        _qcFile._receiverType = obsFile->receiverType();
     149        _qcFile._interval     = obsFile->interval();
     150      }
     151      _qcFile._endTime = _currEpo->tt;
     152
     153      t_qcEpo qcEpo;
     154      qcEpo._epoTime = _currEpo->tt;
     155      qcEpo._PDOP    = cmpDOP(xyzSta);
    234156
    235157      // Loop over all satellites
    … …  
    237159      for (unsigned iObs = 0; iObs < _currEpo->rnxSat.size(); iObs++) {
    238160        const t_rnxObsFile::t_rnxSat& rnxSat = _currEpo->rnxSat[iObs];
    239         t_satObs obs;
    240         t_rnxObsFile::setObsFromRnx(obsFile, _currEpo, rnxSat, obs);
    241 
    242         const t_prn& prn = obs._prn;
    243 
    244         t_ephGlo* ephGlo  = 0;
    245         int       slotNum = 0;
    246         if (prn.system() == 'R') {
    247           for (int ie = 0; ie < _ephs.size(); ie++) {
    248             if (_ephs[ie]->prn() == prn) {
    249               ephGlo = dynamic_cast<t_ephGlo*>(_ephs[ie]);
    250               break;
    251             }
    252           }
    253           if (ephGlo) {
    254             slotNum = ephGlo->slotNum();
    255           }
    256         }
    257 
    258         t_irc irc = _allObsMap[prn].addObs(obs, slotNum);
    259 
    260         if (irc == success) {
    261           t_oneObs* newObs = _allObsMap[prn]._oneObsVec.last();
    262           if (ephGlo) {
    263             newObs->_slotSet = true;
    264           }
    265           if (newObs->_hasL1 && newObs->_hasL2) {
    266             _obsStat._prnStat[prn]._numObs += 1;
    267           }
    268           if (newObs->_slipL1 && newObs->_slipL2) {
    269             _obsStat._prnStat[prn]._numSlipsFlagged += 1;
    270           }
    271         }
    272       }
    273 
    274       prepareObsStat(iEpo, obsFile->interval(), xyzSta);
    275       iEpo++;
    276 
    277     } // while (_currEpo)
     161        t_satObs satObs;
     162        t_rnxObsFile::setObsFromRnx(obsFile, _currEpo, rnxSat, satObs);
     163        t_qcObs qcObs;
     164        if (setQcObs(satObs, qcObs) == success) {
     165          qcEpo._qcObs[satObs._prn] = qcObs;
     166          updateQcSat(qcObs, _qcFile._qcSat[satObs._prn]);
     167        }
     168      }
     169      _qcFile._qcEpo.push_back(qcEpo);
     170    }
     171
     172    preparePlotData(obsFile);
     173
     174    printReport();
    278175  }
    279176  catch (QString str) {
    … …  
    284181      qDebug() << str;
    285182    }
    286     _mutex.unlock();
    287     return;
    288   }
    289 
    290   // Analyze the Multipath
    291   // ---------------------
    292   QVector<t_polarPoint*>* dataMP1  = new QVector<t_polarPoint*>;
    293   QVector<t_polarPoint*>* dataMP2  = new QVector<t_polarPoint*>;
    294   QVector<t_polarPoint*>* dataSNR1 = new QVector<t_polarPoint*>;
    295   QVector<t_polarPoint*>* dataSNR2 = new QVector<t_polarPoint*>;
    296 
    297   QMutableMapIterator<t_prn, t_allObs> it(_allObsMap);
    298   while (it.hasNext()) {
    299     it.next();
    300     const t_prn& prn = it.key();
    301     preparePlotData(prn, xyzSta, obsFile->interval(),
    302                     dataMP1, dataMP2, dataSNR1, dataSNR2);
    303   }
    304 
    305   printReport(dataMP1, dataMP2, dataSNR1, dataSNR2);
    306 
    307   // Show the plots
    308   // --------------
    309   if (BNC_CORE->GUIenabled()) {
    310     QFileInfo  fileInfo(obsFile->fileName());
    311     QByteArray title = fileInfo.fileName().toAscii();
    312     emit dspSkyPlot(obsFile->fileName(), "MP1", dataMP1, "MP2", dataMP2,
    313                     "Meters", 2.0);
    314     double mean = 0.0;
    315     for (int ii = 0; ii < dataSNR1->size(); ii++) {
    316       const t_polarPoint* point = dataSNR1->at(ii);
    317       mean += point->_value;
    318     }
    319     if (dataSNR1->size() > 0) {
    320       mean /= dataSNR1->size();
    321     }
    322     double max = (mean > 9.0) ? 54.0 : 9.0;
    323     QByteArray str = (mean > 9.0) ? "dbHz" : "";
    324     emit dspSkyPlot(obsFile->fileName(), "SNR1", dataSNR1, "SNR2", dataSNR2,
    325                 str, max);
    326     emit dspAvailPlot(obsFile->fileName(), title);
    327   }
    328   else {
    329     for (int ii = 0; ii < dataMP1->size(); ii++) {
    330       delete dataMP1->at(ii);
    331     }
    332     delete dataMP1;
    333     for (int ii = 0; ii < dataMP2->size(); ii++) {
    334       delete dataMP2->at(ii);
    335     }
    336     delete dataMP2;
    337     for (int ii = 0; ii < dataSNR1->size(); ii++) {
    338       delete dataSNR1->at(ii);
    339     }
    340     delete dataSNR1;
    341     for (int ii = 0; ii < dataSNR2->size(); ii++) {
    342       delete dataSNR2->at(ii);
    343     }
    344     delete dataSNR2;
    345     _mutex.unlock();
    346   }
    347 }
    348 
    349 //
    350 ////////////////////////////////////////////////////////////////////////////
    351 t_irc t_reqcAnalyze::t_allObs::addObs(const t_satObs& obs, int slotNum) {
    352 
    353   t_oneObs* newObs = new t_oneObs(obs._time.gpsw(), obs._time.gpssec());
    354   bool      okFlag = false;
    355 
    356   // Availability and Slip Flags
    357   // ---------------------------
    358   double L1 = 0.0;
    359   double L2 = 0.0;
    360   double P1 = 0.0;
    361   double P2 = 0.0;
    362 
    363   for (unsigned iFrq = 0; iFrq < obs._obs.size(); iFrq++) {
    364     const t_frqObs* frqObs = obs._obs[iFrq];
    365     if      (frqObs->_rnxType2ch[0] == '1') {
    366       if (frqObs->_phaseValid) {
    367         L1              = frqObs->_phase;
    368         newObs->_hasL1  = true;
    369         newObs->_slipL1 = frqObs->_slip;
    370       }
    371       if (frqObs->_codeValid) {
    372         P1 = frqObs->_code;
    373       }
    374       if (frqObs->_snrValid) {
    375         newObs->_SNR1 = frqObs->_snr;
    376       }
    377     }
    378     else if ( (obs._prn.system() != 'E' && frqObs->_rnxType2ch[0] == '2') ||
    379               (obs._prn.system() == 'E' && frqObs->_rnxType2ch[0] == '5') ) {
    380       if (frqObs->_phaseValid) {
    381         L2             = frqObs->_phase;
    382         newObs->_hasL2 = true;
    383         newObs->_slipL2 = frqObs->_slip;
    384       }
    385       if (frqObs->_codeValid) {
    386         P2 = frqObs->_code;
    387       }
    388       if (frqObs->_snrValid) {
    389         newObs->_SNR2 = frqObs->_snr;
    390       }
    391     }
    392   }
    393 
    394   // Compute the Multipath
    395   // ----------------------
    396   if (L1 != 0.0 && L2 != 0.0) {
    397     double f1 = 0.0;
    398     double f2 = 0.0;
    399     if      (obs._prn.system() == 'G') {
    400       f1 = t_CST::freq(t_frequency::G1, 0);
    401       f2 = t_CST::freq(t_frequency::G2, 0);
    402     }
    403     else if (obs._prn.system() == 'R') {
    404       f1 = t_CST::freq(t_frequency::R1, slotNum);
    405       f2 = t_CST::freq(t_frequency::R2, slotNum);
    406     }
    407     else if (obs._prn.system() == 'E') {
    408       f1 = t_CST::freq(t_frequency::E1, 0);
    409       f2 = t_CST::freq(t_frequency::E5, 0);
    410     }
    411 
    412     L1 = L1 * t_CST::c / f1;
    413     L2 = L2 * t_CST::c / f2;
    414 
    415     if (P1 != 0.0) {
    416       newObs->_MP1 = P1 - L1 - 2.0*f2*f2/(f1*f1-f2*f2) * (L1 - L2);
    417       okFlag = true;
    418     }
    419     if (P2 != 0.0) {
    420       newObs->_MP2 = P2 - L2 - 2.0*f1*f1/(f1*f1-f2*f2) * (L1 - L2);
    421       okFlag = true;
    422     }
    423   }
    424 
    425   // Remember the Observation
    426   // ------------------------
    427   if (okFlag) {
    428     _oneObsVec << newObs;
    429     return success;
    430   }
    431   else {
    432     delete newObs;
    433     return failure;
    434   }
    435 }
    436 
    437 //
    438 ////////////////////////////////////////////////////////////////////////////
    439 void t_reqcAnalyze::prepareObsStat(unsigned iEpo, double obsInterval,
    440                                    const ColumnVector& xyzSta) {
    441   const int sampl = int(30.0 / obsInterval);
    442   if (iEpo % sampl == 0) {
    443     double mjdX24 = _currEpo->tt.mjddec() * 24.0;
    444     if (iEpo != 0) {
    445       _obsStat._mjdX24 << mjdX24;
    446       _obsStat._numSat << _obsStat._numSat.last();
    447       _obsStat._PDOP   << _obsStat._PDOP.last();
    448     }
    449     _obsStat._mjdX24 << mjdX24;
    450     _obsStat._numSat << _currEpo->rnxSat.size();
    451     _obsStat._PDOP   << cmpDOP(xyzSta);
    452   }
    453 }
    454 
    455 //
    456 ////////////////////////////////////////////////////////////////////////////
    457 void t_reqcAnalyze::preparePlotData(const t_prn& prn,
    458                                     const ColumnVector& xyzSta,
    459                                     double obsInterval,
    460                                     QVector<t_polarPoint*>* dataMP1,
    461                                     QVector<t_polarPoint*>* dataMP2,
    462                                     QVector<t_polarPoint*>* dataSNR1,
    463                                     QVector<t_polarPoint*>* dataSNR2) {
    464 
    465   const int chunkStep = int( 30.0 / obsInterval); // chunk step (30 sec)
    466   const int numEpo    = int(600.0 / obsInterval); // # epochs in one chunk (10 min)
    467 
    468   t_allObs& allObs = _allObsMap[prn];
    469 
    470   bncSettings settings;
    471   QString reqSkyPlotSystems = settings.value("reqcSkyPlotSystems").toString();
    472   bool plotGPS = false;
    473   bool plotGlo = false;
    474   bool plotGal = false;
    475   if      (reqSkyPlotSystems == "GPS") {
    476     plotGPS = true;
    477   }
    478   else if (reqSkyPlotSystems == "GLONASS") {
    479     plotGlo = true;
    480   }
    481   else if (reqSkyPlotSystems == "Galileo") {
    482     plotGal = true;
    483   }
    484   else {
    485     plotGPS = true;
    486     plotGlo = true;
    487     plotGal = true;
    488   }
    489 
    490   // Loop over all Chunks of Data
    491   // ----------------------------
    492   bool slipFound = false;
    493   for (int chunkStart = 0; chunkStart + numEpo < allObs._oneObsVec.size();
    494        chunkStart += chunkStep) {
    495 
    496     if (chunkStart * chunkStep == numEpo) {
    497       slipFound = false;
    498     }
    499 
    500     // Chunk-Specific Variables
    501     // ------------------------
    502     bncTime currTime;
    503     bncTime prevTime;
    504     bncTime chunkStartTime;
    505     double  mjdX24  = 0.0;
    506     bool    availL1 = false;
    507     bool    availL2 = false;
    508     bool    gapL1   = false;
    509     bool    gapL2   = false;
    510     bool    slipL1  = false;
    511     bool    slipL2  = false;
    512     double  meanMP1 = 0.0;
    513     double  meanMP2 = 0.0;
    514     double  minSNR1 = 0.0;
    515     double  minSNR2 = 0.0;
    516     double  aziDeg  = 0.0;
    517     double  zenDeg  = 0.0;
    518     bool    zenFlag = false;
    519 
    520     // Loop over all Epochs within one Chunk of Data
    521     // ---------------------------------------------
    522     bool slotSet = false;
    523     for (int ii = 0; ii < numEpo; ii++) {
    524       int iEpo = chunkStart + ii;
    525       const t_oneObs* oneObs = allObs._oneObsVec[iEpo];
    526       if (oneObs->_slotSet) {
    527         slotSet = true;
    528       }
    529 
    530       currTime.set(oneObs->_GPSWeek, oneObs->_GPSWeeks);
    531 
    532       // Compute the Azimuth and Zenith Distance
    533       // ---------------------------------------
    534       if (ii == 0) {
    535         chunkStartTime = currTime;
    536         mjdX24 = chunkStartTime.mjddec() * 24.0;
    537 
    538         if (xyzSta.size()) {
    539           t_eph* eph = 0;
    540           for (int ie = 0; ie < _ephs.size(); ie++) {
    541             if (_ephs[ie]->prn() == prn) {
    542               eph = _ephs[ie];
    543               break;
    544             }
    545           }
    546 
    547           if (eph) {
    548             ColumnVector xc(4);
    549             ColumnVector vv(3);
    550             if (eph->getCrd(bncTime(oneObs->_GPSWeek, oneObs->_GPSWeeks), xc, vv, false) == success) {
    551               double rho, eleSat, azSat;
    552               topos(xyzSta(1), xyzSta(2), xyzSta(3), xc(1), xc(2), xc(3), rho, eleSat, azSat);
    553               aziDeg = azSat * 180.0/M_PI;
    554               zenDeg = 90.0 - eleSat * 180.0/M_PI;
    555               zenFlag = true;
    556             }
    557           }
    558         }
    559       }
    560 
    561       // Check Interval
    562       // --------------
    563       if (prevTime.valid()) {
    564         double dt = currTime - prevTime;
    565         if (dt > 1.5 * obsInterval) {
    566           gapL1 = true;
    567           gapL2 = true;
    568         }
    569       }
    570       prevTime = currTime;
    571 
    572       // Check L1 and L2 availability
    573       // ----------------------------
    574       if (oneObs->_hasL1) {
    575         availL1 = true;
    576       }
    577       else {
    578         gapL1 = true;
    579       }
    580       if (oneObs->_hasL2) {
    581         availL2 = true;
    582       }
    583       else {
    584         gapL2 = true;
    585       }
    586 
    587       // Check Minimal Signal-to-Noise Ratio
    588       // -----------------------------------
    589       if ( oneObs->_SNR1 > 0 && (minSNR1 == 0 || minSNR1 > oneObs->_SNR1) ) {
    590         minSNR1 = oneObs->_SNR1;
    591       }
    592       if ( oneObs->_SNR2 > 0 && (minSNR2 == 0 || minSNR2 > oneObs->_SNR2) ) {
    593         minSNR2 = oneObs->_SNR2;
    594       }
    595 
    596       // Check Slip Flags
    597       // ----------------
    598       if (oneObs->_slipL1) {
    599         slipL1 = true;
    600       }
    601       if (oneObs->_slipL2) {
    602         slipL2 = true;
    603       }
    604 
    605       meanMP1 += oneObs->_MP1;
    606       meanMP2 += oneObs->_MP2;
    607     }
    608 
    609     // Compute the Multipath
    610     // ---------------------
    611     if ( (prn.system() == 'G' && plotGPS           ) ||
    612          (prn.system() == 'R' && plotGlo && slotSet) ||
    613          (prn.system() == 'E' && plotGal           ) ) {
    614       bool slipMP = false;
    615       meanMP1 /= numEpo;
    616       meanMP2 /= numEpo;
    617       double MP1 = 0.0;
    618       double MP2 = 0.0;
    619       for (int ii = 0; ii < numEpo; ii++) {
    620         int iEpo = chunkStart + ii;
    621         const t_oneObs* oneObs = allObs._oneObsVec[iEpo];
    622         double diff1 = oneObs->_MP1 - meanMP1;
    623         double diff2 = oneObs->_MP2 - meanMP2;
    624 
    625         // Check Slip Threshold
    626         // --------------------
    627         if (fabs(diff1) > SLIPTRESH || fabs(diff2) > SLIPTRESH) {
    628           slipMP = true;
    629           break;
    630         }
    631 
    632         MP1 += diff1 * diff1;
    633         MP2 += diff2 * diff2;
    634       }
    635       if (slipMP) {
    636         slipL1 = true;
    637         slipL2 = true;
    638         if (!slipFound) {
    639           slipFound = true;
    640           _obsStat._prnStat[prn]._numSlipsFound += 1;
    641         }
    642       }
    643       else {
    644         MP1 = sqrt(MP1 / (numEpo-1));
    645         MP2 = sqrt(MP2 / (numEpo-1));
    646         (*dataMP1)  << (new t_polarPoint(aziDeg, zenDeg, MP1));
    647         (*dataMP2)  << (new t_polarPoint(aziDeg, zenDeg, MP2));
    648       }
    649     }
    650 
    651     // Availability Plot Data
    652     // ----------------------
    653     if (availL1) {
    654       if      (slipL1) {
    655         _availDataMap[prn]._L1slip << mjdX24;
    656       }
    657       else if (gapL1) {
    658         _availDataMap[prn]._L1gap << mjdX24;
    659       }
    660       else {
    661         _availDataMap[prn]._L1ok << mjdX24;
    662       }
    663     }
    664     if (availL2) {
    665       if      (slipL2) {
    666         _availDataMap[prn]._L2slip << mjdX24;
    667       }
    668       else if (gapL2) {
    669         _availDataMap[prn]._L2gap << mjdX24;
    670       }
    671       else {
    672         _availDataMap[prn]._L2ok << mjdX24;
    673       }
    674     }
    675     if (zenFlag) {
    676       _availDataMap[prn]._eleTim << mjdX24;
    677       _availDataMap[prn]._eleDeg << 90.0 - zenDeg;
    678     }
    679 
    680     // Signal-to-Noise Ratio Plot Data
    681     // -------------------------------
    682     if ( (prn.system() == 'G' && plotGPS) ||
    683          (prn.system() == 'R' && plotGlo) ||
    684          (prn.system() == 'E' && plotGal) ) {
    685       (*dataSNR1) << (new t_polarPoint(aziDeg, zenDeg, minSNR1));
    686       (*dataSNR2) << (new t_polarPoint(aziDeg, zenDeg, minSNR2));
    687     }
    688   }
    689 }
    690 
    691 //
    692 ////////////////////////////////////////////////////////////////////////////
    693 void t_reqcAnalyze::slotDspAvailPlot(const QString& fileName,
    694                                      const QByteArray& title) {
    695 
    696   if (BNC_CORE->GUIenabled()) {
    697     t_availPlot* plotA = new t_availPlot(0, &_availDataMap);
    698     plotA->setTitle(title);
    699 
    700     t_elePlot* plotZ = new t_elePlot(0, &_availDataMap);
    701 
    702     t_dopPlot* plotD = new t_dopPlot(0, &_obsStat);
    703 
    704     QVector<QWidget*> plots;
    705     plots << plotA << plotZ << plotD;
    706     t_graphWin* graphWin = new t_graphWin(0, fileName, plots, 0, 0);
    707 
    708     int ww = QFontMetrics(graphWin->font()).width('w');
    709     graphWin->setMinimumSize(120*ww, 40*ww);
    710 
    711     graphWin->show();
    712 
    713     bncSettings settings;
    714     QString dirName = settings.value("reqcPlotDir").toString();
    715     if (!dirName.isEmpty()) {
    716       QByteArray ext = "_A.png";
    717       graphWin->savePNG(dirName, ext);
    718     }
    719   }
    720   _mutex.unlock();
     183  }
    721184}
    722185
    … …  
    778241}
    779242
     243//
     244////////////////////////////////////////////////////////////////////////////
     245void t_reqcAnalyze::updateQcSat(const t_qcObs& qcObs, t_qcSat& qcSat) {
     246  if (qcObs._hasL1 && qcObs._hasL2) {
     247    qcSat._numObs += 1;
     248  }
     249  if (qcObs._slipL1 && qcObs._slipL2) {
     250    qcSat._numSlipsFlagged += 1;
     251  }
     252}
     253
     254//
     255////////////////////////////////////////////////////////////////////////////
     256t_irc t_reqcAnalyze::setQcObs(const t_satObs& satObs, t_qcObs& qcObs) {
     257
     258  if (satObs._prn.system() == 'R') {
     259    t_ephGlo* ephGlo  = 0;
     260    for (int ie = 0; ie < _ephs.size(); ie++) {
     261      if (_ephs[ie]->prn() == satObs._prn) {
     262        ephGlo = dynamic_cast<t_ephGlo*>(_ephs[ie]);
     263        break;
     264      }
     265    }
     266    if (ephGlo) {
     267      qcObs._slotSet = true;
     268      qcObs._slotNum = ephGlo->slotNum();
     269    }
     270  }
     271
     272  bool okFlag = false;
     273
     274  // Availability and Slip Flags
     275  // ---------------------------
     276  double L1 = 0.0;
     277  double L2 = 0.0;
     278  double P1 = 0.0;
     279  double P2 = 0.0;
     280
     281  for (unsigned iFrq = 0; iFrq < satObs._obs.size(); iFrq++) {
     282    const t_frqObs* frqObs = satObs._obs[iFrq];
     283    if      (frqObs->_rnxType2ch[0] == '1') {
     284      if (frqObs->_phaseValid) {
     285        L1              = frqObs->_phase;
     286        qcObs._hasL1  = true;
     287        qcObs._slipL1 = frqObs->_slip;
     288      }
     289      if (frqObs->_codeValid) {
     290        P1 = frqObs->_code;
     291      }
     292      if (frqObs->_snrValid) {
     293        qcObs._SNR1 = frqObs->_snr;
     294      }
     295    }
     296    else if ( (satObs._prn.system() != 'E' && frqObs->_rnxType2ch[0] == '2') ||
     297              (satObs._prn.system() == 'E' && frqObs->_rnxType2ch[0] == '5') ) {
     298      if (frqObs->_phaseValid) {
     299        L2             = frqObs->_phase;
     300        qcObs._hasL2 = true;
     301        qcObs._slipL2 = frqObs->_slip;
     302      }
     303      if (frqObs->_codeValid) {
     304        P2 = frqObs->_code;
     305      }
     306      if (frqObs->_snrValid) {
     307        qcObs._SNR2 = frqObs->_snr;
     308      }
     309    }
     310  }
     311
     312  // Compute the Multipath
     313  // ----------------------
     314  if (L1 != 0.0 && L2 != 0.0) {
     315    double f1 = 0.0;
     316    double f2 = 0.0;
     317    if      (satObs._prn.system() == 'G') {
     318      f1 = t_CST::freq(t_frequency::G1, 0);
     319      f2 = t_CST::freq(t_frequency::G2, 0);
     320    }
     321    else if (satObs._prn.system() == 'R') {
     322      f1 = t_CST::freq(t_frequency::R1, qcObs._slotNum);
     323      f2 = t_CST::freq(t_frequency::R2, qcObs._slotNum);
     324    }
     325    else if (satObs._prn.system() == 'E') {
     326      f1 = t_CST::freq(t_frequency::E1, 0);
     327      f2 = t_CST::freq(t_frequency::E5, 0);
     328    }
     329
     330    L1 = L1 * t_CST::c / f1;
     331    L2 = L2 * t_CST::c / f2;
     332
     333    if (P1 != 0.0) {
     334      qcObs._MP1 = P1 - L1 - 2.0*f2*f2/(f1*f1-f2*f2) * (L1 - L2);
     335      okFlag = true;
     336    }
     337    if (P2 != 0.0) {
     338      qcObs._MP2 = P2 - L2 - 2.0*f1*f1/(f1*f1-f2*f2) * (L1 - L2);
     339      okFlag = true;
     340    }
     341  }
     342
     343  if (okFlag) {
     344    return success;
     345  }
     346  else {
     347    return failure;
     348  }
     349}
     350
     351//
     352////////////////////////////////////////////////////////////////////////////
     353void t_reqcAnalyze::preparePlotData(const t_rnxObsFile* obsFile) {
     354
     355  ColumnVector xyzSta = obsFile->xyz();
     356
     357  QVector<t_polarPoint*>* dataMP1  = new QVector<t_polarPoint*>;
     358  QVector<t_polarPoint*>* dataMP2  = new QVector<t_polarPoint*>;
     359  QVector<t_polarPoint*>* dataSNR1 = new QVector<t_polarPoint*>;
     360  QVector<t_polarPoint*>* dataSNR2 = new QVector<t_polarPoint*>;
     361
     362  const double SLIPTRESH = 10.0;  // cycle-slip threshold (meters)
     363  const double chunkStep = 600.0; // 10 minutes
     364
     365  bncSettings settings;
     366  QString reqSkyPlotSystems = settings.value("reqcSkyPlotSystems").toString();
     367  bool plotGPS = false;
     368  bool plotGlo = false;
     369  bool plotGal = false;
     370  if      (reqSkyPlotSystems == "GPS") {
     371    plotGPS = true;
     372  }
     373  else if (reqSkyPlotSystems == "GLONASS") {
     374    plotGlo = true;
     375  }
     376  else if (reqSkyPlotSystems == "Galileo") {
     377    plotGal = true;
     378  }
     379  else {
     380    plotGPS = true;
     381    plotGlo = true;
     382    plotGal = true;
     383  }
     384
     385  // Loop over all satellites available
     386  // ----------------------------------
     387  QMutableMapIterator<t_prn, t_qcSat> it(_qcFile._qcSat);
     388  while (it.hasNext()) {
     389    it.next();
     390    const t_prn& prn   = it.key();
     391    t_qcSat&     qcSat = it.value();
     392
     393    // Loop over all Chunks of Data
     394    // ----------------------------
     395    for (bncTime chunkStart = _qcFile._startTime;
     396         chunkStart < _qcFile._endTime; chunkStart += chunkStep) {
     397
     398      // Chunk (sampled) Epoch
     399      // ---------------------
     400      _qcFile._qcEpoSampled.push_back(t_qcEpo());
     401      t_qcEpo& qcEpoSampled = _qcFile._qcEpoSampled.back();
     402      t_qcObs& qcObsSampled = qcEpoSampled._qcObs[prn];
     403
     404      QVector<double> MP1;
     405      QVector<double> MP2;
     406   
     407      // Loop over all Epochs within one Chunk of Data
     408      // ---------------------------------------------
     409      bncTime prevTime;
     410      for (int iEpo = 0; iEpo < _qcFile._qcEpo.size(); iEpo++) {
     411        const t_qcEpo& qcEpo = _qcFile._qcEpo[iEpo];
     412        if (qcEpo._epoTime < chunkStart) {
     413          continue;
     414        }
     415        if (!qcEpo._qcObs.contains(prn)) {
     416          continue;
     417        }
     418     
     419        const t_qcObs& qcObs = qcEpo._qcObs[prn];
     420     
     421        // Compute the Azimuth and Zenith Distance
     422        // ---------------------------------------
     423        if (chunkStart == qcEpo._epoTime) {
     424          qcEpoSampled._epoTime = qcEpo._epoTime;
     425          qcEpoSampled._PDOP    = qcEpo._PDOP;
     426
     427          if (qcObs._slotSet) {
     428            qcObsSampled._slotSet = true;
     429            qcObsSampled._slotNum = qcObs._slotNum;
     430          }
     431          if (xyzSta.size()) {
     432            t_eph* eph = 0;
     433            for (int ie = 0; ie < _ephs.size(); ie++) {
     434              if (_ephs[ie]->prn() == prn) {
     435                eph = _ephs[ie];
     436                break;
     437              }
     438            }
     439            if (eph) {
     440              ColumnVector xc(4);
     441              ColumnVector vv(3);
     442              if (eph->getCrd(qcEpo._epoTime, xc, vv, false) == success) {
     443                double rho, eleSat, azSat;
     444                topos(xyzSta(1), xyzSta(2), xyzSta(3), xc(1), xc(2), xc(3), rho, eleSat, azSat);
     445                qcObsSampled._azDeg  = azSat * 180.0/M_PI;
     446                qcObsSampled._eleDeg = eleSat * 180.0/M_PI;
     447              }
     448            }
     449          }
     450        }
     451     
     452        // Check Interval
     453        // --------------
     454        if (prevTime.valid()) {
     455          double dt = qcEpo._epoTime - prevTime;
     456          if (dt > 1.5 * _qcFile._interval) {
     457            qcObsSampled._gapL1 = true;
     458            qcObsSampled._gapL2 = true;
     459          }
     460        }
     461        prevTime = qcEpo._epoTime;
     462     
     463        // Check L1 and L2 availability
     464        // ----------------------------
     465        if (qcObs._hasL1) {
     466          qcObsSampled._hasL1 = true;
     467        }
     468        if (qcObs._hasL2) {
     469          qcObsSampled._hasL2 = true;
     470        }
     471     
     472        // Check Minimal Signal-to-Noise Ratio
     473        // -----------------------------------
     474        if ( qcObs._SNR1 > 0 && (qcObsSampled._SNR1 == 0 || qcObsSampled._SNR1 > qcObs._SNR1) ) {
     475          qcObsSampled._SNR1 = qcObs._SNR1;
     476        }
     477        if ( qcObs._SNR2 > 0 && (qcObsSampled._SNR2 == 0 || qcObsSampled._SNR2 > qcObs._SNR2) ) {
     478          qcObsSampled._SNR2 = qcObs._SNR2;
     479        }
     480     
     481        // Check Slip Flags
     482        // ----------------
     483        if (qcObs._slipL1) {
     484          qcObsSampled._slipL1 = true;
     485        }
     486        if (qcObs._slipL2) {
     487          qcObsSampled._slipL2 = true;
     488        }
     489     
     490        MP1 << qcObs._MP1;
     491        MP2 << qcObs._MP2;
     492      }
     493
     494      // Compute the Multipath
     495      // ---------------------
     496      if ( MP1.size() > 0 && MP2.size() > 0 &&
     497           ( (prn.system() == 'G' && plotGPS                         ) ||
     498             (prn.system() == 'R' && plotGlo && qcObsSampled._slotSet) ||
     499             (prn.system() == 'E' && plotGal                         ) ) ) {
     500
     501        double meanMP1 = 0.0;
     502        double meanMP2 = 0.0;
     503        for (int ii = 0; ii < MP1.size(); ii++) {
     504          meanMP1 += MP1[ii];
     505          meanMP2 += MP2[ii];
     506        }
     507        meanMP1 /= MP1.size();
     508        meanMP2 /= MP2.size();
     509
     510        bool slipMP = false;
     511        double stdMP1 = 0.0;
     512        double stdMP2 = 0.0;
     513        for (int ii = 0; ii < MP1.size(); ii++) {
     514          double diff1 = MP1[ii] - meanMP1;
     515          double diff2 = MP2[ii] - meanMP2;
     516          if (fabs(diff1) > SLIPTRESH || fabs(diff2) > SLIPTRESH) {
     517            slipMP = true;
     518            break;
     519          }
     520          stdMP1 += diff1 * diff1;
     521          stdMP1 += diff2 * diff2;
     522        }
     523
     524        if (slipMP) {
     525          qcObsSampled._slipL1 = true;
     526          qcObsSampled._slipL2 = true;
     527          qcSat._numSlipsFound += 1;
     528        }
     529        else {
     530          stdMP1 = sqrt(stdMP1 / MP1.size());
     531          stdMP2 = sqrt(stdMP2 / MP2.size());
     532          (*dataMP1)  << (new t_polarPoint(qcObsSampled._azDeg, 90.0 - qcObsSampled._eleDeg, stdMP1));
     533          (*dataMP2)  << (new t_polarPoint(qcObsSampled._azDeg, 90.0 - qcObsSampled._eleDeg, stdMP2));
     534        }
     535      }
     536   
     537      // Signal-to-Noise Ratio Plot Data
     538      // -------------------------------
     539      if ( (prn.system() == 'G' && plotGPS) ||
     540           (prn.system() == 'R' && plotGlo) ||
     541           (prn.system() == 'E' && plotGal) ) {
     542        (*dataSNR1) << (new t_polarPoint(qcObsSampled._azDeg, 90.0 - qcObsSampled._eleDeg, qcObsSampled._SNR1));
     543        (*dataSNR2) << (new t_polarPoint(qcObsSampled._azDeg, 90.0 - qcObsSampled._eleDeg, qcObsSampled._SNR2));
     544      }
     545    }
     546  }
     547
     548  // Show the plots
     549  // --------------
     550  if (BNC_CORE->GUIenabled()) {
     551    QFileInfo  fileInfo(obsFile->fileName());
     552    QByteArray title = fileInfo.fileName().toAscii();
     553    dspSkyPlot(obsFile->fileName(), "MP1", dataMP1, "MP2", dataMP2, "Meters", 2.0);
     554    dspSkyPlot(obsFile->fileName(), "SNR1", dataSNR1, "SNR2", dataSNR2, "dbHz", 54.0);
     555    dspAvailPlot(obsFile->fileName(), title);
     556  }
     557  else {
     558    for (int ii = 0; ii < dataMP1->size(); ii++) {
     559      delete dataMP1->at(ii);
     560    }
     561    delete dataMP1;
     562    for (int ii = 0; ii < dataMP2->size(); ii++) {
     563      delete dataMP2->at(ii);
     564    }
     565    delete dataMP2;
     566    for (int ii = 0; ii < dataSNR1->size(); ii++) {
     567      delete dataSNR1->at(ii);
     568    }
     569    delete dataSNR1;
     570    for (int ii = 0; ii < dataSNR2->size(); ii++) {
     571      delete dataSNR2->at(ii);
     572    }
     573    delete dataSNR2;
     574  }
     575}
     576
     577//
     578////////////////////////////////////////////////////////////////////////////
     579void t_reqcAnalyze::dspSkyPlot(const QString& fileName, const QByteArray& title1,
     580                               QVector<t_polarPoint*>* data1, const QByteArray& title2,
     581                               QVector<t_polarPoint*>* data2, const QByteArray& scaleTitle,
     582                               double maxValue) {
     583
     584  if (BNC_CORE->GUIenabled()) {
     585
     586    if (maxValue == 0.0) {
     587      if (data1) {
     588        for (int ii = 0; ii < data1->size(); ii++) {
     589          double val = data1->at(ii)->_value;
     590          if (maxValue < val) {
     591            maxValue = val;
     592          }
     593        }
     594      }
     595      if (data2) {
     596        for (int ii = 0; ii < data2->size(); ii++) {
     597          double val = data2->at(ii)->_value;
     598          if (maxValue < val) {
     599            maxValue = val;
     600          }
     601        }
     602      }
     603    }
     604
     605    QwtInterval scaleInterval(0.0, maxValue);
     606
     607    QVector<QWidget*> plots;
     608    if (data1) {
     609      t_polarPlot* plot1 = new t_polarPlot(QwtText(title1), scaleInterval,
     610                                          BNC_CORE->mainWindow());
     611      plot1->addCurve(data1);
     612      plots << plot1;
     613    }
     614    if (data2) {
     615      t_polarPlot* plot2 = new t_polarPlot(QwtText(title2), scaleInterval,
     616                                           BNC_CORE->mainWindow());
     617      plot2->addCurve(data2);
     618      plots << plot2;
     619    }
     620
     621    t_graphWin* graphWin = new t_graphWin(0, fileName, plots,
     622                                          &scaleTitle, &scaleInterval);
     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::dspAvailPlot(const QString& fileName, const QByteArray& title) {
     638
     639  if (BNC_CORE->GUIenabled()) {
     640    t_availPlot* plotA = new t_availPlot(0, _qcFile._qcEpoSampled);
     641    plotA->setTitle(title);
     642
     643    t_elePlot* plotZ = new t_elePlot(0, _qcFile._qcEpoSampled);
     644
     645    t_dopPlot* plotD = new t_dopPlot(0, _qcFile._qcEpoSampled);
     646
     647    QVector<QWidget*> plots;
     648    plots << plotA << plotZ << plotD;
     649    t_graphWin* graphWin = new t_graphWin(0, fileName, plots, 0, 0);
     650
     651    int ww = QFontMetrics(graphWin->font()).width('w');
     652    graphWin->setMinimumSize(120*ww, 40*ww);
     653
     654    graphWin->show();
     655
     656    bncSettings settings;
     657    QString dirName = settings.value("reqcPlotDir").toString();
     658    if (!dirName.isEmpty()) {
     659      QByteArray ext = "_A.png";
     660      graphWin->savePNG(dirName, ext);
     661    }
     662  }
     663}
     664
    780665// Finish the report
    781666////////////////////////////////////////////////////////////////////////////
    782 void t_reqcAnalyze::printReport(QVector<t_polarPoint*>* dataMP1,
    783                                 QVector<t_polarPoint*>* dataMP2,
    784                                 QVector<t_polarPoint*>* dataSNR1,
    785                                 QVector<t_polarPoint*>* dataSNR2) {
     667void t_reqcAnalyze::printReport() {
    786668
    787669  if (!_log) {
    … …  
    789671  }
    790672
    791   *_log << "Marker name:     " << _obsStat._markerName   << endl
    792         << "Receiver:        " << _obsStat._receiverType << endl
    793         << "Antenna:         " << _obsStat._antennaName  << endl
    794         << "Start time:      " << _obsStat._startTime.datestr().c_str() << ' '
    795                                << _obsStat._startTime.timestr().c_str() << endl
    796         << "End time:        " << _obsStat._endTime.datestr().c_str() << ' '
    797                                << _obsStat._endTime.timestr().c_str() << endl
    798         << "Interval:        " << _obsStat._interval << endl
    799         << "# Sat.:          " << _obsStat._prnStat.size() << endl;
     673  *_log << "Marker name:     " << _qcFile._markerName   << endl
     674        << "Receiver:        " << _qcFile._receiverType << endl
     675        << "Antenna:         " << _qcFile._antennaName  << endl
     676        << "Start time:      " << _qcFile._startTime.datestr().c_str() << ' '
     677                               << _qcFile._startTime.timestr().c_str() << endl
     678        << "End time:        " << _qcFile._endTime.datestr().c_str() << ' '
     679                               << _qcFile._endTime.timestr().c_str() << endl
     680        << "Interval:        " << _qcFile._interval << endl
     681        << "# Sat.:          " << _qcFile._qcSat.size() << endl;
    800682
    801683  int numObs          = 0;
    802684  int numSlipsFlagged = 0;
    803685  int numSlipsFound   = 0;
    804   QMapIterator<t_prn, t_prnStat> it(_obsStat._prnStat);
     686  QMapIterator<t_prn, t_qcSat> it(_qcFile._qcSat);
    805687  while (it.hasNext()) {
    806688    it.next();
    807     const t_prnStat& prnStat = it.value();
    808     numObs          += prnStat._numObs;
    809     numSlipsFlagged += prnStat._numSlipsFlagged;
    810     numSlipsFound   += prnStat._numSlipsFound;
     689    const t_qcSat& qcSat = it.value();
     690    numObs          += qcSat._numObs;
     691    numSlipsFlagged += qcSat._numSlipsFlagged;
     692    numSlipsFound   += qcSat._numSlipsFound;
    811693  }
    812694  *_log << "# Obs.:          " << numObs          << endl
    … …  
    814696        << "# Slips (found): " << numSlipsFound   << endl;
    815697
    816   for (int kk = 1; kk <= 4; kk++) {
    817     QVector<t_polarPoint*>* data = 0;
    818     QString text;
    819     if      (kk == 1) {
    820       data = dataMP1;
    821       text = "Mean MP1:        ";
    822     }
    823     else if (kk == 2) {
    824       data = dataMP2;
    825       text = "Mean MP2:        ";
    826     }
    827     else if (kk == 3) {
    828       data = dataSNR1;
    829       text = "Mean SNR1:       ";
    830     }
    831     else if (kk == 4) {
    832       data = dataSNR2;
    833       text = "Mean SNR2:       ";
    834     }
    835     double mean = 0.0;
    836     for (int ii = 0; ii < data->size(); ii++) {
    837       const t_polarPoint* point = data->at(ii);
    838       mean += point->_value;
    839     }
    840     if (data->size() > 0) {
    841       mean /= data->size();
    842     }
    843     *_log << text << mean << endl;
    844   }
    845 
    846698  _log->flush();
    847699}
Note: See TracChangeset for help on using the changeset viewer.