tecmaps.cpp

根据GPS观测数据

      if(help) cout << "Option d appears " << dashd.getCount() << " times\n";
   }

   if(Rest.getCount() && help) {
      cout << "Remaining options:" << endl;
      values = Rest.getValue();
      for (i=0; i<values.size(); i++) cout << values[i] << endl;
   }
   if(verbose && help) {
      cout << "\nTokens on command line (" << Args.size() << ") are:" << endl;
      for(j=0; j<Args.size(); j++) cout << Args[j] << endl;
      cout << endl;
   }
   if(help) cout << "Configuration summarized in log file\n";

      // print configuration
   if(verbose || help || debug) {
      oflog << "\nInput configuration for TECMaps:\n";
      if(!InputPath.empty()) oflog << " Path for input files is "
         << InputPath << endl;
      oflog << " Input Rinex obs file names are:\n";
      for(i=0; i<Stations.size(); i++)
         oflog << "   " << Stations[i].filename << endl;
      if(!NavDir.empty()) oflog << " Path for input nav files is "
         << NavDir  << endl;
      if(NavFiles.size() > 0) {
         oflog << " Input Rinex nav file names are:\n";
         for(i=0; i<NavFiles.size(); i++)
            oflog << "   " << NavFiles[i] << endl;
      }
      if(BegTime > DayTime::BEGINNING_OF_TIME) oflog << " Begin time is "
         << BegTime.printf("%Y/%m/%d_%H:%M:%6.3f=%F/%10.3g") << endl;
      if(EndTime < DayTime::END_OF_TIME) oflog << " End   time is "
         << EndTime.printf("%Y/%m/%d_%H:%M:%6.3f=%F/%10.3g") << endl;
      oflog << " Processing:\n";
      oflog << "  Primary Title is " << Title1 << endl;
      oflog << "  Secondary Title is " << Title2 << endl;
      oflog << "  Reference " << refSite << endl;
      oflog << "  Do " << (doVTECmap ? "":"not ") << "create VTEC map\n";
      oflog << "  Do " << (doMUFmap ? "":"not ") << "create MUF map\n";
      oflog << "  Do " << (doF0F2map ? "":"not ") << "create F0F2 map\n";
      oflog << "  Grid type is "
         << (typegrid == VTECMap::UniformLatLon ? "Uniform " : "Uniform Space ")
         << typegrid << endl;
      oflog << "  Beginning latitude (deg) is " << BeginLat << endl;
      oflog << "  Beginning longitude (deg E) is " << BeginLon << endl;
      oflog << "  Number of latitude grid points is " << NumLat << endl;
      oflog << "  Number of latitude grid points is " << NumLat << endl;
      oflog << "  Grid step in latitude (deg) is " << DeltaLat << endl;
      oflog << "  Grid step in longitude (deg) is " << DeltaLon << endl;
      oflog << "  Minimum elevation (deg) is " << ElevThresh << endl;
      oflog << "  Minimum acquisition time (sec) is " << MinAcqTime << endl;
      if(BiasFile.length() > 0)
         oflog << "  Receiver biases input from file " << BiasFile << endl;
      else
         oflog << "  Do not input sat+rx biases" << endl;
      oflog << "  Decorrelation error rate (TECU/1000km) is " << DecorrelError
         << endl;
      oflog << "  Ionosphere height = " << IonoHt << " km" << endl;
      oflog << "  Base name for output files is " << BaseName << endl;
      cout << (GridOut ? "O":"Do NOT o") << "utput grid in file named " <<
         BaseName << ".LL" << endl;
      if(GridOut) cout << "Output grid in " << (GnuplotFormat ? "gnuplot" : "Matlab")
         << " format" << endl;
      if(ExSV.size() > 0) {
         oflog << "  Exclude these satellites:";
         for(i=0; i<ExSV.size(); i++) oflog << " " << ExSV[i];
         oflog << endl;
      }
      if(help) oflog << "  Help (-h) is set\n";
      if(debug) oflog << "  Debug (-d) is set\n";
      oflog << "End of input configuration for TECMaps" << endl << endl;
   }

   if(help) return 1;

   return 0;
}
catch(Exception& e) { GPSTK_RETHROW(e); }
catch(exception& e) { Exception E("std except: "+string(e.what())); GPSTK_THROW(E); }
catch(...) { Exception e("Unknown exception"); GPSTK_THROW(e); }
}

//------------------------------------------------------------------------------------
// Pull out --debug --verbose options.
void PreProcessArgs(const char *arg, vector<string>& Args) throw(Exception)
{
try {
   static bool found_cfg_file=false;

   if(found_cfg_file || (arg[0]=='-' && arg[1]=='f')) {
      string filename(arg);
      if(!found_cfg_file) filename.erase(0,2); else found_cfg_file = false;
      if(debug) cout << "Found a file of options: " << filename << endl;
      ifstream infile(filename.c_str());
      if(!infile) {
         cout << "Error: could not open options file " << filename << endl;
         return;
      }

      bool again_cfg_file=false;
      char c;
      string buffer,word;
      while(1) {
         getline(infile,buffer);
         stripTrailing(buffer,'\r');

         // process the buffer before checking eof or bad b/c there can be
         // a line at EOF that has no CRLF...
         while(!buffer.empty()) {
            word = firstWord(buffer);
            if(again_cfg_file) {
               word = "-f" + word;
               again_cfg_file = false;
               PreProcessArgs(word.c_str(),Args);
            }
            else if(word[0] == '#') { // skip to end of line
               buffer.clear();
            }
            else if(word == "--file" || word == "-f")
               again_cfg_file = true;
            else if(word[0] == '"') {
               word = stripFirstWord(buffer,'"');
               buffer = "dummy " + buffer;            // to be stripped later
               PreProcessArgs(word.c_str(),Args);
            }
            else
               PreProcessArgs(word.c_str(),Args);

            word = stripFirstWord(buffer);      // now remove it from buffer
         }
         if(infile.eof() || !infile.good()) break;
      }
   }
   else if((arg[0]=='-' && arg[1]=='d') || string(arg)==string("--debug"))
      debug = true;
   else if((arg[0]=='-' && arg[1]=='v') || string(arg)==string("--verbose"))
      verbose = true;
   else if(string(arg) == "--file" || string(arg) == "-f")
      found_cfg_file = true;
   else Args.push_back(arg);
}
catch(Exception& e) { GPSTK_RETHROW(e); }
catch(exception& e) { Exception E("std except: "+string(e.what())); GPSTK_THROW(E); }
catch(...) { Exception e("Unknown exception"); GPSTK_THROW(e); }
}

//------------------------------------------------------------------------------------
int Initialize(void) throw(Exception)
{
try {
   static SP3EphemerisStore SP3EphList;
   static BCEphemerisStore BCEphList;

      // open nav files and read EphemerisStore
   if(!NavDir.empty())
      for(int i=0; i<NavFiles.size(); i++)
         NavFiles[i] = NavDir + "/" + NavFiles[i];
   FillEphemerisStore(NavFiles, SP3EphList, BCEphList);
   if(SP3EphList.size()) {
      if(verbose) SP3EphList.dump(0,oflog);
   }
   else if(verbose) oflog << "SP3 Ephemeris list is empty\n";

   if(BCEphList.size()) {
      BCEphList.SearchNear();
      if(verbose) BCEphList.dump(0,oflog);
   }
   else if(verbose) oflog << "BC Ephemeris list is empty\n";

   if(SP3EphList.size()) pEph = &SP3EphList;
   else if(BCEphList.size()) pEph = &BCEphList;
   else {
      cerr << "TECMaps abort -- no ephemeris\n";
      oflog << "TECMaps abort -- no ephemeris\n";
      return -1;
   }

      // read the sat+rx biases file and build the map
   if(BiasFile.length() > 0) {
      ifstream inf(BiasFile.c_str());
      if(!inf) {
         cerr << "Error: could not open biases file " << BiasFile << endl;
         oflog << "Error: could not open biases file " << BiasFile << endl;
         return -1;
      }
      else {
         int nbiases,n=0;
         double bias;
         RinexSatID sat;
         string line,station;
         vector<string> words;
         while(!inf.eof() && inf.good()) {
            {        // parse line into words
               istringstream iss(line);
               string wd;
               words.clear();
               while(iss >> wd) { words.push_back(wd); }
            }
            while(1) {
               if(words.size() != 5) break;
               if(debug) {
                  oflog << "Biases file:";
                  for(int i=0; i<words.size(); i++) oflog << " " << words[i];
                  oflog << endl;
               }
               if(words[0] == string("IonoBias,")) break;
               if(words[1] == string("Number")) {
                  nbiases = asInt(words[0]);
                  break;
               }
               station = words[1];
               sat.fromString(words[2]);
               bias = asDouble(words[3]);
               BiasMap[station][sat] = bias;
               n++;
               break;
            }
            getline(inf,line);
         }
         inf.close();
         if(verbose) {
            oflog << "Read sat+rx biases file " << BiasFile << ":" << endl;
            oflog << " Expected " << nbiases << " biases, read " << n << "." << endl;
            oflog << " Here is the biases map:" << endl;
            map<string,map<RinexSatID,double> >::const_iterator it;
            map<RinexSatID,double>::const_iterator jt;
            for(it=BiasMap.begin(); it!=BiasMap.end(); it++) {
               for(jt=it->second.begin(); jt!=it->second.end(); jt++) {
                  oflog << "  " << it->first
                     << " " << jt->first
                     << " " << setw(12) << setprecision(6) << jt->second << endl;
               }
            }
         }
      }
   }

      // create the obs types for later use
   ELot = RinexObsHeader::convertObsType("EL");
   AZot = RinexObsHeader::convertObsType("AZ");
   LAot = RinexObsHeader::convertObsType("LA"); // TEMP
   LOot = RinexObsHeader::convertObsType("LO"); // TEMP
   SRot = RinexObsHeader::convertObsType("SR");
   VRot = RinexObsHeader::convertObsType("VR");
   TPot = RinexObsHeader::convertObsType("TP");

      // initialize the map(s)
      // set parameters
   vtecmap.Decorrelation = DecorrelError;
   vtecmap.MinElevation = ElevThresh;
   vtecmap.IonoHeight = IonoHt*1000;
   vtecmap.gridtype = typegrid;
   vtecmap.fittype = typefit;
   vtecmap.BeginLat = BeginLat;
   vtecmap.DeltaLat = DeltaLat;
   vtecmap.NumLat = NumLat;
   vtecmap.BeginLon = BeginLon;
   vtecmap.DeltaLon = DeltaLon;
   vtecmap.NumLon = NumLon;
   if(doMUFmap) mufmap.CopyInputData(vtecmap);
   if(doF0F2map) f0f2map.CopyInputData(vtecmap);

   return 0;
}
catch(Exception& e) { GPSTK_RETHROW(e); }
catch(exception& e) { Exception E("std except: "+string(e.what())); GPSTK_THROW(E); }
catch(...) { Exception e("Unknown exception"); GPSTK_THROW(e); }
}

//------------------------------------------------------------------------------------
// Return 0 ok,
//       -1 could not allocate array of RinexObsStream
//       -2 could not open a file,
//       -3 FFStream exception,
//       -4 gpstk exception,
int ProcessStations(void) throw(Exception)
{
try {
   int iret,nfile;

   if(verbose) oflog << "\nProcess " << Stations.size()
      << " input files / Stations:\n";

      // create vector of RinexObsStream for use by the files
   instream = new RinexObsStream[Stations.size()];
   if(!instream) {
      oflog << " ERROR: Array of " << Stations.size()
         << " RinexObsStream's could not be allocated\n";
      return -1;
   }

      // loop over input file names
      // open and read header, define entry in Stations array
   for(nfile=0; nfile<Stations.size(); nfile++) {
      instream[nfile].open(Stations[nfile].filename.c_str(),ios_base::in);
      if(!instream[nfile]) {
         oflog << " Rinex file " << Stations[nfile].filename
            << " could not be opened -- abort.\n";
         Stations[nfile].nfile = -1;
         Stations[nfile].nread = -1;
         Stations[nfile].getNext = false;
         return -2;
      }
      oflog << "\nInput file #" << nfile+1 << ": "
         << Stations[nfile].filename << endl;

      instream[nfile].exceptions(ios::failbit);
      Stations[nfile].nfile = nfile;
      Stations[nfile].getNext = true;

      iret = ProcessHeader(Stations[nfile]);
      if(iret != 0) return iret;

   }  // end loop over file names -> stations
   oflog << "\nEnd of loop over stations to read headers.\n\n";

   return 0;
}
catch(Exception& e) { GPSTK_RETHROW(e); }
catch(exception& e) { Exception E("std except: "+string(e.what())); GPSTK_THROW(E); }
catch(...) { Exception e("Unknown exception"); GPSTK_THROW(e); }
}

//------------------------------------------------------------------------------------
void ProcessObsAndComputeMap(void) throw(Exception)
{
try {
   int iret,nfile,ngood,nepochs=0;
   vector<ObsData> AllObs; // ObsData vector, passed into ProcessObs and Compute

      // loop over all epochs in all files
   do {
         // read the data for the next (earliest in future) observation epoch
      EarliestTime = DayTime(DayTime::END_OF_TIME);
      for(nfile=0; nfile<Stations.size(); nfile++) {
         iret = ReadNextObs(Stations[nfile]);
         if(iret < 0) {            // set file 'inactive'
            Stations[nfile].nread = -1;
         }
         else {                    // iret is 0 or 1 - file is active
            if(Stations[nfile].robs.time < EarliestTime)
               EarliestTime = Stations[nfile].robs.time;
         }
      }

         // if no more data available, EarliestTime will never get set
      if(EarliestTime == DayTime(DayTime::END_OF_TIME)) break;

         // time limits
      if(EarliestTime < BegTime) {
         oflog << "Before begin time : "
            << EarliestTime.printf("%Y/%m/%d %H:%M:%6.3f=%F/%10.3g\n");
         continue;
      }
      if(EarliestTime > EndTime) {
         oflog << "After end time (quit) : "
            << EarliestTime.printf("%Y/%m/%d %H:%M:%6.3f=%F/%10.3g\n");
         break;
      }

         // process at EarliestTime
      oflog << "Process at time = "
         << EarliestTime.printf("%Y/%m/%d %H:%M:%6.3f=%F/%10.3g\n");
      AllObs.clear();
      for(ngood=0,nfile=0; nfile<Stations.size(); nfile++) {
            // if this data time == EarliestTime, process and set flag to read again
         if(ABS(Stations[nfile].robs.time - EarliestTime) < 1.) {
            iret = ProcessObs(Stations[nfile],AllObs);
            if(iret > 0) ngood += iret;