tecmaps.cpp

根据GPS观测数据

#pragma ident "$Id: TECMaps.cpp 378 2007-01-26 16:11:34Z 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 TECMaps.cpp
 * Program TECMaps reads a set of Rinex files containing observation types
 *    EL, AZ, and VR or SR and fits the ionospheric vertical TEC data to a model
 *    of the ionosphere. There are input options for the type of grid, the type of
 *    model, and the type of data (VTEC, MUF or F0F2) to be used. TD ...
 */

#include "StringUtils.hpp"
#include "DayTime.hpp"
#include "RinexSatID.hpp"
#include "CommandOption.hpp"
#include "CommandOptionWithTimeArg.hpp"
#include "CommandOptionParser.hpp"

#include "BCEphemerisStore.hpp"
#include "SP3EphemerisStore.hpp"
#include "WGS84Geoid.hpp"
#include "Position.hpp"

#include "VTECMap.hpp"
#include "RinexUtilities.hpp"

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

//------------------------------------------------------------------------------------
using namespace std;
using namespace gpstk;
using namespace StringUtils;

//------------------------------------------------------------------------------------
// input data
bool verbose,debug;
   // log file
string LogFile;
ofstream oflog;
string Title;
   // input path
string InputPath;
string Title1,Title2,BaseName,BiasFile;
double DecorrelError,ElevThresh,MinAcqTime;
double BeginLat,DeltaLat,BeginLon,DeltaLon;
int NumLat,NumLon;
VTECMap::FitType typefit;
VTECMap::GridType typegrid;
bool doVTECmap,doMUFmap,doF0F2map;
Station refSite;
string KnownPos;         // string holding position x,y,z or l,l,h
bool KnownLLH;           // if true, KnownPos is l,l,h
bool GridOut;            // if true, write grid to file 'basename.LL'
bool GnuplotFormat;      // if true, write grid in format for gnuplot
   // excluded satellites
vector<RinexSatID> ExSV;
   // ephemeris
string NavDir;
vector<string> NavFiles;
EphemerisStore *pEph;         // is this used?
   // obs types needed
RinexObsHeader::RinexObsType ELot,AZot,VRot,SRot,TPot;
RinexObsHeader::RinexObsType LAot,LOot;      // TEMP
   // geoid
WGS84Geoid WGS84;
   // Start and stop times
DayTime BegTime,EndTime;
   // processing
double IonoHt;
DayTime EarliestTime;
VTECMap vtecmap;
MUFMap mufmap;
F0F2Map f0f2map;
   // map of input sat+rx biases
map<string,map<RinexSatID,double> > BiasMap;
   // Data structures for all receivers
vector<Station> Stations;
RinexObsStream *instream; // array of streams, parallell to Stations

//------------------------------------------------------------------------------------
// prototypes
void ConfigureAndDefaults(void);
int GetCommandLine(int argc, char **argv) throw(Exception);
void PreProcessArgs(const char *arg, vector<string>& Args) throw(Exception);
int Initialize(void) throw(Exception);
int ProcessStations(void) throw(Exception);
void ProcessObsAndComputeMap(void) throw(Exception);
void OutputGridToFile(VTECMap& vmap, string filename) throw(Exception);
void OutputMapToFile(VTECMap& vtmap, string filename, DayTime t, int n)
   throw(Exception);
void AddStation(string& filename) throw(Exception);
int ProcessHeader(Station& S) throw(Exception);
int ReadNextObs(Station& S) throw(Exception);
int ProcessObs(Station& S, vector<ObsData>& obsvect) throw(Exception);

//------------------------------------------------------------------------------------
//------------------------------------------------------------------------------------
int main(int argc, char **argv)
{
try {
   int iret;
   clock_t totaltime=clock(); // timer
   DayTime CurrEpoch;

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

      // Title description and run time
   CurrEpoch.setLocalTime();
   Title = "TECMaps, built on the GPSTK ToolKit, Ver 1.0 8/12/04, Run ";
   Title += CurrEpoch.printf("%04Y/%02m/%02d %02H:%02M:%02S\n");
   cout << Title;

      // define extended types
   iret = RegisterARLUTExtendedTypes();
   if(iret) goto quit;
   iret = RegisterExtendedRinexObsType("TP","Acquisition time","seconds", 0);
   if(iret) goto quit;

      // set configuration and default values
   ConfigureAndDefaults();

      // get command line arguments
   iret = GetCommandLine(argc, argv);
   if(iret) goto quit;

      // initialize
   iret = Initialize();
   if(iret) goto quit;

      // make the grid
   if(doVTECmap) {
      vtecmap.MakeGrid(refSite);
      if(GridOut) OutputGridToFile(vtecmap, BaseName+string(".LL"));
   }
   if(doMUFmap) {
      mufmap.MakeGrid(refSite);
      if(GridOut) OutputGridToFile(mufmap, BaseName+string(".MUF.LL"));
   }
   if(doF0F2map) {
      f0f2map.MakeGrid(refSite);
      if(GridOut) OutputGridToFile(f0f2map, BaseName+string(".F0F2.LL"));
   }

      // process the headers, filling the Stations array
   iret = ProcessStations();
   if(iret) goto quit;

      // process the all the observation data
   ProcessObsAndComputeMap();

quit:
      // compute run time
   totaltime = clock()-totaltime;
   cout << "TECMaps timing: " << fixed << setprecision(3)
      << double(totaltime)/double(CLOCKS_PER_SEC) << " seconds.\n";
   oflog << "TECMaps timing: " << fixed << setprecision(3)
      << double(totaltime)/double(CLOCKS_PER_SEC) << " seconds.\n";
   oflog.close();

   return iret;
}
catch(gpstk::Exception& e) { cerr << e; }
catch (...) { cerr << "Unknown error.  Abort." << endl; }
   return 1;
}   // end main()

//------------------------------------------------------------------------------------
// set defaults
void ConfigureAndDefaults(void)
{
   Title1=string("TECMaps main title");
   Title2=string("TECMaps sub title");
   BaseName=string("tecmap_out");
   DecorrelError = 3.0;
   BiasFile = string("");
   ElevThresh = 10.0;
   MinAcqTime = 0.0;
   BeginLat = 21.0;
   DeltaLat = 0.25;
   BeginLon = 230;
   DeltaLon = 1.0;
   NumLat = 40;
   NumLon = 40;
   typefit = VTECMap::Constant;
   typegrid = VTECMap::UniformLatLon;
   doVTECmap = true;
   doMUFmap = false;
   doF0F2map = false;
   LogFile = string("vtm.log");
   IonoHt = 350.0;         // km
   verbose = false;
   debug = false;
   KnownPos = string("");
   GridOut = false;
   GnuplotFormat = false;
}

//------------------------------------------------------------------------------------
// Define, parse and evaluate command line
int GetCommandLine(int argc, char **argv) throw(Exception)
{
try {
   bool help=false;
   int i,j;
   RinexSatID sat;
   sat.setfill('0');

      // required options

   RequiredOption dashin(CommandOption::hasArgument, CommandOption::stdType,
      0,"input"," --input <file>       Input Rinex obs file name(s)"
                "                    \n(Reference site position also required)"
      );
   //dashin.setMaxCount(1);

      // optional options

   // this only so it will show up in help page...
   CommandOption dashf(CommandOption::hasArgument, CommandOption::stdType,
      'f',""," -f<file>             file containing more options");

   // reference position(s)
   CommandOption dashllh(CommandOption::hasArgument, CommandOption::stdType,0,"RxLLH",
     "Reference station position (one required):\n"
     " --RxLLH <l,l,h>      Reference site position in geodetic"
     " lat, lon (E), ht (deg,deg,m)");
   dashllh.setMaxCount(1);

   CommandOption dashxyz(CommandOption::hasArgument, CommandOption::stdType,0,"RxXYZ",
      " --RxXYZ <x,y,z>      Reference site position in ECEF coordinates (m)");
   dashxyz.setMaxCount(1);

      // require one only
   CommandOptionMutex refmutex(true);
   refmutex.addOption(&dashllh);
   refmutex.addOption(&dashxyz);

   CommandOption dashp(CommandOption::hasArgument, CommandOption::stdType,
      0,"inputdir"," --inputdir <path>    Path for input file(s)");
   dashp.setMaxCount(1);

   // ephemeris
   CommandOption dashnd(CommandOption::hasArgument, CommandOption::stdType, 0,
      "navdir","Ephemeris input:\n --navdir <dir>       Path of navigation file(s)");
   dashnd.setMaxCount(1);

   CommandOption dashn(CommandOption::hasArgument, CommandOption::stdType,
      0,"nav"," --nav <file>         Navigation (Rinex Nav OR SP3) file(s)");

   CommandOption dashl(CommandOption::hasArgument, CommandOption::stdType,
      0,"log","Output:\n --log <file>         Output log file name");
   dashl.setMaxCount(1);
   
   // time
   CommandOptionWithTimeArg dasheb(0,"BeginTime","%Y,%m,%d,%H,%M,%f",
      "Time limits:\n --BeginTime <arg>    Start time, arg is of the form "
      "YYYY,MM,DD,HH,Min,Sec");
   CommandOptionWithTimeArg dashgb(0,"BeginGPSTime","%F,%g",
      " --BeginGPSTime <arg> Start time, arg is of the form GPSweek,GPSsow");

   CommandOptionWithTimeArg dashee(0,"EndTime","%Y,%m,%d,%H,%M,%f",
      " --EndTime <arg>      End time, arg is of the form YYYY,MM,DD,HH,Min,Sec");
   CommandOptionWithTimeArg dashge(0,"EndGPSTime","%F,%g",
      " --EndGPSTime <arg>   End time, arg is of the form GPSweek,GPSsow");

   // allow ONLY one start time (use startmutex(true) if one is required)
   CommandOptionMutex startmutex(false);
   startmutex.addOption(&dasheb);
   startmutex.addOption(&dashgb);
   CommandOptionMutex stopmutex(false);
   stopmutex.addOption(&dashee);
   stopmutex.addOption(&dashge);

   CommandOptionNoArg dashVmap(0, "MUFmap","Processing:\n"
      " --noVTECmap          Do NOT create the VTEC map.");

   CommandOptionNoArg dashMUF(0, "MUFmap",
      " --MUFmap             Create MUF map as well as VTEC map.");

   CommandOptionNoArg dashF0F2(0, "F0F2map",
      " --F0F2map            Create F0F2 map as well as VTEC map.");

   CommandOption dashTitle1(CommandOption::hasArgument, CommandOption::stdType,
      0,"Title1"," --Title1 <title>     Title information");
   dashTitle1.setMaxCount(1);

   CommandOption dashTitle2(CommandOption::hasArgument, CommandOption::stdType,
      0,"Title2", " --Title2 <title>     Second title information");
   dashTitle2.setMaxCount(1);

   CommandOption dashBaseName(CommandOption::hasArgument, CommandOption::stdType,
      0,"BaseName", " --BaseName <name>    Base name for output files (a)");
   dashBaseName.setMaxCount(1);

   CommandOption dashDecor(CommandOption::hasArgument, CommandOption::stdType,
      0,"DecorrError",
      " --DecorrError <de>   Decorrelation error rate in TECU/1000km (3)");
   dashDecor.setMaxCount(1);

   CommandOption dashBiases(CommandOption::hasArgument, CommandOption::stdType,
      0,"Biases", " --Biases <file>      "
      "File containing estimated sat+rx biases (Prgm IonoBias)");
   dashBiases.setMaxCount(1);

   CommandOption dashElevThresh(CommandOption::hasArgument, CommandOption::stdType,
      0,"ElevThresh", " --ElevThresh <ele>   Minimum elevation (6 deg)");
   dashElevThresh.setMaxCount(1);

   CommandOption dashMinAcqTime(CommandOption::hasArgument, CommandOption::stdType,
      0,"MinAcqTime", " --MinAcqTime <t>     Minimum acquisition time (0 sec)");
   dashMinAcqTime.setMaxCount(1);

   CommandOptionNoArg dashFlatFit(
      0, "FlatFit", " --FlatFit            Flat fit type (default)");

   CommandOptionNoArg dashLinearFit(
      0, "LinearFit", " --LinearFit          Linear fit type");

   CommandOption dashIonoHt(CommandOption::hasArgument, CommandOption::stdType,
      0,"IonoHeight", " --IonoHeight <n>     Ionosphere height (km)");
   dashIonoHt.setMaxCount(1);

   CommandOptionNoArg dashUniSpace(0, "UniformSpacing",
      "Grid:\n --UniformSpacing     Grid uniform in space (XYZ) (default)");

   CommandOptionNoArg dashUniGrid(0, "UniformGrid",
      " --UniformGrid        Grid uniform in Lat and Lon");

   CommandOptionNoArg dashOutGrid(
      0, "OutputGrid", " --OutputGrid         Output the grid to file <basename.LL>");

   CommandOptionNoArg dashGnuOut(
      0, "GnuplotOutput", " --GnuplotOutput      Write the grid file for gnuplot"
      " (default: for Matlab)");

   CommandOption dashNumLat(CommandOption::hasArgument, CommandOption::stdType,
      0,"NumLat", " --NumLat <n>         Number of latitude grid points (40)");
   dashNumLat.setMaxCount(1);

   CommandOption dashNumLon(CommandOption::hasArgument, CommandOption::stdType,
      0,"NumLon", " --NumLon <n>         Number of longitude grid points (40)");
   dashNumLon.setMaxCount(1);

   CommandOption dashBeginLat(CommandOption::hasArgument, CommandOption::stdType,
      0,"BeginLat", " --BeginLat <lat>     Beginning latitude (21 deg)");