rinsum.cpp

GPS数据预处理软件

#pragma ident "$Id: RinSum.cpp 698 2007-07-23 17:08:30Z btolman $"

//============================================================================
//
//  This file is part of GPSTk, the GPS Toolkit.
//
//  The GPSTk is free software; you can redistribute it and/or modify
//  it under the terms of the GNU Lesser General Public License as published
//  by the Free Software Foundation; either version 2.1 of the License, or
//  any later version.
//
//  The GPSTk is distributed in the hope that it will be useful,
//  but WITHOUT ANY WARRANTY; without even the implied warranty of
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//  GNU Lesser General Public License for more details.
//
//  You should have received a copy of the GNU Lesser General Public
//  License along with GPSTk; if not, write to the Free Software Foundation,
//  Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
//  
//  Copyright 2004, The University of Texas at Austin
//
//============================================================================

/**
 * @file RinSum.cpp
 * Read and summarize Rinex observation files, optionally fill header in-place.
 */

#include "MathBase.hpp"
#include "RinexObsBase.hpp"
#include "RinexObsData.hpp"
#include "RinexObsHeader.hpp"
#include "RinexObsStream.hpp"
#include "RinexNavBase.hpp"
#include "RinexNavHeader.hpp"
#include "RinexNavData.hpp"
#include "RinexNavStream.hpp"
#include "DayTime.hpp"
#include "SatID.hpp"
#include "RinexSatID.hpp"
#include "CommandOptionParser.hpp"
#include "CommandOption.hpp"
#include "CommandOptionWithTimeArg.hpp"
#include "icd_200_constants.hpp"
#include "RinexUtilities.hpp"

#include <string>
#include <vector>
#include <iostream>
#include <fstream>
#include <algorithm>
#include <time.h>

using namespace std;
using namespace gpstk;
using namespace StringUtils;

//------------------------------------------------------------------------------------
string version("2.4 1/22/07");

// data input from command line
vector<string> InputFiles;
string InputDirectory;
string OutputFile;
ostream* pout;
DayTime BegTime, EndTime;
bool ReplaceHeader=false;
bool TimeSortTable=false;
bool GPSTimeOutput=false;
bool debug=false;
bool brief=false;

//------------------------------------------------------------------------------------
// data used for computation
const int ndtmax=15;
double dt,bestdt[ndtmax];
int ndt[ndtmax]={-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1};
int nepochs,ncommentblocks;

//------------------------------------------------------------------------------------
// class used to store SAT/Obs table
class TableData {
public:
   RinexSatID sat;
   vector<int> nobs;
   double prevC1,prevP1,prevL1;
   DayTime begin,end;
   TableData(const SatID& p, const int& n)
      { sat=RinexSatID(p); nobs=vector<int>(n); prevC1=prevP1=prevL1=0; };
      // needed for find()
   inline bool operator==(const TableData& d) {return d.sat == sat;}
};
   // for sort()
class TableSATLessThan  {      
public:
   bool operator()(const TableData& d1, const TableData& d2)
      { return d1.sat < d2.sat; }
};
class TableBegLessThan  {
public:
   bool operator()(const TableData& d1, const TableData& d2)
      { return d1.begin < d2.begin; }
};

//------------------------------------------------------------------------------------
// prototypes
int GetCommandLine(int argc, char **argv) throw(Exception);
void PreProcessArgs(const char *arg, vector<string>& Args) throw(Exception);

//------------------------------------------------------------------------------------
int main(int argc, char **argv)
{
try {
   int iret,i,j,k,n,ifile,nsats,nclkjumps,L1lli;
   double C1,L1,P1,clkjumpave,clkjumpvar;
   DayTime last,prev,ftime;
   vector<DayTime> clkjumpTimes;
   vector<double> clkjumpMillsecs,clkjumpUncertainty;
   vector<int> clkjumpAgree;

   BegTime = DayTime::BEGINNING_OF_TIME;
   EndTime = DayTime::END_OF_TIME;

      // Title and description
   string Title;
   Title = "RINSUM, part of the GPS ToolKit, Ver " + version + ", Run ";
   time_t timer;
   struct tm *tblock;
   timer = time(NULL);
   tblock = localtime(&timer);
   last.setYMDHMS(1900+tblock->tm_year,1+tblock->tm_mon,
               tblock->tm_mday,tblock->tm_hour,tblock->tm_min,tblock->tm_sec);
   Title += last.printf("%04Y/%02m/%02d %02H:%02M:%02S\n");
   cout << Title;

   iret=GetCommandLine(argc, argv);
   if(iret) return iret;

   iret = RegisterARLUTExtendedTypes();
   if(iret) return iret;

      // open the output file and write to it
   if(!OutputFile.empty()) {
      pout = new ofstream(OutputFile.c_str(),ios::out);
      if(pout->fail()) {
         cerr << "Could not open output file " << OutputFile << endl;
         pout = &cout;
      }
      else {
         pout->exceptions(ios::failbit);
         *pout << Title;
         cout << "Writing summary to file " << OutputFile << endl;
      }
   }
   else pout = &cout;

      // add path to input file names
   if(!InputDirectory.empty()) for(ifile=0; ifile<InputFiles.size(); ifile++) {
      InputFiles[ifile] = InputDirectory + "/" + InputFiles[ifile];
   }
      // sort the input file names on header first time
   if(InputFiles.size() > 1) sortRinexObsFiles(InputFiles);

      // now open the input files, read the headers and data
   RinexObsHeader rheader;
   RinexObsData robs;
   for(ifile=0; ifile<InputFiles.size(); ifile++) {
      string filename;
      if(!InputDirectory.empty()) filename = InputDirectory + "/";
      filename += InputFiles[ifile];
      RinexObsStream InStream(filename.c_str());
      if(!InStream) {
         *pout << "File " << filename << " could not be opened.\n";
         continue;
      }
      InStream.exceptions(ios::failbit);
      if(!isRinexObsFile(filename)) {
         *pout << "File " << filename << " is not a Rinex observation file\n";
         if(isRinexNavFile(filename))
            *pout << "This file is a Rinex navigation file - try NavMerge\n";
         continue;
      }

      prev = DayTime::BEGINNING_OF_TIME;
      ftime = DayTime::BEGINNING_OF_TIME;

      if(!brief) *pout << "+++++++++++++ RinSum summary of Rinex obs file "
         << filename << " +++++++++++++\n";
      else *pout << "\nFile name: " << filename << endl;
      
         // input header
      try {
         InStream >> rheader;
      }
      catch(gpstk::FFStreamError& e) {
         cerr << "Caught an FFStreamError while reading header: "
            << e.getText(0) << endl;
      }
      catch(gpstk::Exception& e) {
         cerr << "Caught a gpstk exception while reading header: "
            << e.getText(0) << endl;
      }

      if(!brief) {
         *pout << "Rinex header:\n";
         rheader.dump(*pout);
      }
      else *pout << "Position (XYZ,m) : " << fixed << setprecision(4)
         << rheader.antennaPosition << ".\n";

      if(!rheader.isValid()) {
         *pout << "Abort: header is invalid\n";
         if(!brief) *pout << "\n+++++++++++++ End of RinSum summary of "
            << filename << " +++++++++++++\n";
         continue;
      }

      //RinexObsStream out(argv[2], ios::out);
      //out << rheader;

      nepochs = ncommentblocks = 0;
      n = rheader.obsTypeList.size();
      vector<TableData> table;
      vector<int> totals(n);

      if(pout == &cout) *pout << "Reading the observation data..." << endl;

         // input obs
      while(InStream >> robs)
      {
         if(debug) *pout << "Epoch: " << robs.time
            << ", Flag " << robs.epochFlag
            << ", Nsat " << robs.obs.size()
            << ", clk " << fixed << robs.clockOffset << endl;

          // is this a comment?
         if(robs.epochFlag > 1) {
            ncommentblocks++;
            //*pout << "inline header info:\n";
            //robs.auxHeader.dump(*pout);
            continue;
         }

         // update first and last time seen, check time limits, count epochs
         last = robs.time;
         if(last < BegTime) continue;
         if(last > EndTime) break;
         if(ftime == DayTime::BEGINNING_OF_TIME) ftime=last;
         nepochs++;
         nsats = nclkjumps = 0;  // count sats and signs clock jumps have occurred
         clkjumpave = clkjumpvar = 0.0;

         // loop over satellites
         RinexObsData::RinexSatMap::const_iterator it;
         RinexObsData::RinexObsTypeMap::const_iterator jt;
         for(it=robs.obs.begin(); it != robs.obs.end(); ++it) {
            // update the table
            vector<TableData>::iterator ptab;
            ptab = find(table.begin(),table.end(),TableData(it->first,n));
            if(ptab == table.end()) {        // sat not found in table - create one
               table.push_back(TableData(it->first,n));
               ptab = find(table.begin(),table.end(),TableData(it->first,n));
               ptab->begin = last;
            }
            // update end time for this sat
            ptab->end = last;
            if(debug) *pout << "Sat " << setw(2) << RinexSatID(it->first);

            // loop over obs types
            C1 = P1 = L1 = 0;
            for(jt=it->second.begin(); jt!=it->second.end(); jt++) {
               // find the index for this obs type
               for(k=0; k<n; k++) if(rheader.obsTypeList[k] == jt->first) break;
               // count this obs
               if(jt->second.data != 0) {
                  ptab->nobs[k]++;      // per obs
                  totals[k]++;
               }
               // save L1 range and phase for clk jump test below
               if(jt->first==RinexObsHeader::C1) C1 = jt->second.data*1000.0/C_GPS_M;
               if(jt->first==RinexObsHeader::P1) P1 = jt->second.data*1000.0/C_GPS_M;
               if(jt->first == RinexObsHeader::L1) {
                  L1 = jt->second.data * 1000.0/C_GPS_M;
                  L1lli = jt->second.lli;