#pragma ident "$Id: example9.cpp 1907 2009-05-21 14:50:11Z architest $"
#include <iostream>
#include <iomanip>
#include <fstream>
#include "BasicFramework.hpp"
#include "RinexObsStream.hpp"
#include "SP3EphemerisStore.hpp"
#include "ProcessingList.hpp"
#include "BasicModel.hpp"
#include "TropModel.hpp"
#include "DataStructures.hpp"
#include "RequireObservables.hpp"
#include "SimpleFilter.hpp"
#include "XYZ2NEU.hpp"
#include "LICSDetector2.hpp"
#include "MWCSDetector.hpp"
#include "SolidTides.hpp"
#include "OceanLoading.hpp"
#include "PoleTides.hpp"
#include "CorrectObservables.hpp"
#include "Antenna.hpp"
#include "AntexReader.hpp"
#include "ComputeWindUp.hpp"
#include "ComputeSatPCenter.hpp"
#include "ComputeTropModel.hpp"
#include "ComputeLinear.hpp"
#include "LinearCombinations.hpp"
#include "ComputeDOP.hpp"
#include "SatArcMarker.hpp"
#include "GravitationalDelay.hpp"
#include "PhaseCodeAlignment.hpp"
#include "EclipsedSatFilter.hpp"
#include "Decimate.hpp"
#include "SolverPPP.hpp"
#include "SolverPPPFB.hpp"
#include "ConfDataReader.hpp"
using namespace std;
using namespace gpstk;
class example9 : public BasicFramework
{
public:
example9(char* arg0);
protected:
virtual void process();
virtual void spinUp();
private:
CommandOptionWithArg confFile;
ConfDataReader confReader;
void printSolution( ofstream& outfile,
const SolverLMS& solver,
const DayTime& time,
const ComputeDOP& cDOP,
bool useNEU,
int numSats,
double dryTropo,
int precision = 3 );
void printModel( ofstream& modelfile,
const gnssRinex& gData,
int precision = 4 );
};
example9::example9(char* arg0)
:
BasicFramework( arg0,
"\nThis program reads GPS receiver data from a configuration file and\n"
"process such data applying a 'Precise Point Positioning' strategy.\n\n"
"Please consult the default configuration file, 'pppconf.txt', for\n"
"further details.\n\n"
"The output file format is as follows:\n\n"
" 1) Year\n"
" 2) Day of year\n"
" 3) Seconds of day\n"
" 4) dx/dLat (m)\n"
" 5) dy/dLon (m)\n"
" 6) dz/dH (m)\n"
" 7) Zenital Tropospheric Delay - zpd (m)\n"
" 8) Covariance of dx/dLat (m*m)\n"
" 9) Covariance of dy/dLon (m*m)\n"
"10) Covariance of dz/dH (m*m)\n"
"11) Covariance of Zenital Tropospheric Delay (m*m)\n"
"12) Number of satellites\n"
"13) GDOP\n"
"14) PDOP\n"
"15) TDOP\n"
"16) HDOP\n"
"17) VDOP\n" ),
confFile( CommandOption::stdType,
'c',
"conffile",
" [-c|--conffile] Name of configuration file ('pppconf.txt' by default).",
false )
{
confFile.setMaxCount(1);
}
void example9::printSolution( ofstream& outfile,
const SolverLMS& solver,
const DayTime& time,
const ComputeDOP& cDOP,
bool useNEU,
int numSats,
double dryTropo,
int precision )
{
outfile << fixed << setprecision( precision );
outfile << time.year() << " ";
outfile << time.DOY() << " ";
outfile << time.DOYsecond() << " ";
if( useNEU )
{
outfile << solver.getSolution(TypeID::dLat) << " ";
outfile << solver.getSolution(TypeID::dLon) << " ";
outfile << solver.getSolution(TypeID::dH) << " ";
outfile << solver.getSolution(TypeID::wetMap) + 0.1 + dryTropo << " ";
outfile << solver.getVariance(TypeID::dLat) << " ";
outfile << solver.getVariance(TypeID::dLon) << " ";
outfile << solver.getVariance(TypeID::dH) << " ";
outfile << solver.getVariance(TypeID::wetMap) << " ";
}
else
{
outfile << solver.getSolution(TypeID::dx) << " ";
outfile << solver.getSolution(TypeID::dy) << " ";
outfile << solver.getSolution(TypeID::dz) << " ";
outfile << solver.getSolution(TypeID::wetMap) + 0.1 + dryTropo << " ";
outfile << solver.getVariance(TypeID::dx) << " ";
outfile << solver.getVariance(TypeID::dy) << " ";
outfile << solver.getVariance(TypeID::dz) << " ";
outfile << solver.getVariance(TypeID::wetMap) << " ";
}
outfile << numSats << " ";
outfile << cDOP.getGDOP() << " ";
outfile << cDOP.getPDOP() << " ";
outfile << cDOP.getTDOP() << " ";
outfile << cDOP.getHDOP() << " ";
outfile << cDOP.getVDOP() << " ";
outfile << endl;
return;
}
void example9::printModel( ofstream& modelfile,
const gnssRinex& gData,
int precision )
{
modelfile << fixed << setprecision( precision );
DayTime time(gData.header.epoch);
for ( satTypeValueMap::const_iterator it = gData.body.begin();
it!= gData.body.end();
it++ )
{
modelfile << time.year() << " ";
modelfile << time.DOY() << " ";
modelfile << time.DOYsecond() << " ";
modelfile << (*it).first << " ";
for( typeValueMap::const_iterator itObs = (*it).second.begin();
itObs != (*it).second.end();
itObs++ )
{
modelfile << (*itObs).first << " ";
modelfile << (*itObs).second << " ";
}
modelfile << endl;
}
}
void example9::spinUp()
{
if ( confFile.getCount() > 0 )
{
confReader.exceptions(ios::failbit);
try
{
confReader.open( confFile.getValue()[0] );
}
catch(...)
{
cerr << "Problem opening file "
<< confFile.getValue()[0]
<< endl;
cerr << "Maybe it doesn't exist or you don't have proper "
<< "read permissions." << endl;
exit (-1);
}
}
else
{
try
{
confReader.open( "pppconf.txt" );
}
catch(...)
{
cerr << "Problem opening default configuration file 'pppconf.txt'"
<< endl;
cerr << "Maybe it doesn't exist or you don't have proper read "
<< "permissions. Try providing a configuration file with "
<< "option '-c'."
<< endl;
exit (-1);
}
}
confReader.setFallback2Default(true);
}
void example9::process()
{
string station;
while ( (station = confReader.getEachSection()) != "" )
{
if( station == "DEFAULT" )
{
continue;
}
cout << "Starting processing for station: '" << station << "'." << endl;
RinexObsStream rin;
rin.exceptions(ios::failbit);
try
{
rin.open( confReader("rinexObsFile", station), std::ios::in );
}
catch(...)
{
cerr << "Problem opening file '"
<< confReader.getValue("rinexObsFile", station)
<< "'." << endl;
cerr << "Maybe it doesn't exist or you don't have "
<< "proper read permissions."
<< endl;
cerr << "Skipping receiver '" << station << "'."
<< endl;
rin.close();
continue;
}
SP3EphemerisStore SP3EphList;
SP3EphList.rejectBadPositions(true);
SP3EphList.rejectBadClocks(true);
if ( confReader.getValueAsBoolean( "checkGaps", station ) )
{
SP3EphList.enableDataGapCheck();
SP3EphList.setGapInterval(
confReader.getValueAsDouble("SP3GapInterval",station) );
}
if ( confReader.getValueAsBoolean( "checkInterval", station ) )
{
SP3EphList.enableIntervalCheck();
SP3EphList.setMaxInterval(
confReader.getValueAsDouble("maxSP3Interval",station) );
}
string sp3File;
while ( (sp3File = confReader.fetchListValue("SP3List",station) ) != "" )
{
try
{
SP3EphList.loadFile( sp3File );
}
catch (FileMissingException& e)
{
cerr << "SP3 file '" << sp3File << "' doesn't exist or you don't "
<< "have permission to read it. Skipping it." << endl;
continue;
}
}
double xn(confReader.fetchListValueAsDouble("nominalPosition",station));
double yn(confReader.fetchListValueAsDouble("nominalPosition",station));
double zn(confReader.fetchListValueAsDouble("nominalPosition",station));
Position nominalPos( xn, yn, zn );
ProcessingList pList;
RequireObservables requireObs;
requireObs.addRequiredType(TypeID::P2);
requireObs.addRequiredType(TypeID::L1);
requireObs.addRequiredType(TypeID::L2);
SimpleFilter pObsFilter;
pObsFilter.setFilteredType(TypeID::P2);
bool usingC1( confReader.getValueAsBoolean( "useC1", station ) );
if ( usingC1 )
{
requireObs.addRequiredType(TypeID::C1);
pObsFilter.addFilteredType(TypeID::C1);
}
else
{
requireObs.addRequiredType(TypeID::P1);
pObsFilter.addFilteredType(TypeID::P1);
}
pList.push_back(requireObs);
bool filterCode( confReader.getValueAsBoolean( "filterCode", station ) );
if( filterCode )
{
pList.push_back(pObsFilter);
}
LinearCombinations comb;
ComputeLinear linear1;
if ( usingC1 )
{
linear1.addLinear(comb.pdeltaCombWithC1);
linear1.addLinear(comb.mwubbenaCombWithC1);
}
else
{
linear1.addLinear(comb.pdeltaCombination);
linear1.addLinear(comb.mwubbenaCombination);
}
linear1.addLinear(comb.ldeltaCombination);
linear1.addLinear(comb.liCombination);
pList.push_back(linear1);
LICSDetector2 markCSLI2;
pList.push_back(markCSLI2);
MWCSDetector markCSMW;
pList.push_back(markCSMW);
SatArcMarker markArc;
markArc.setDeleteUnstableSats(true);
markArc.setUnstablePeriod(151.0);
pList.push_back(markArc);
Decimate decimateData(
confReader.getValueAsDouble( "decimationInterval", station ),
confReader.getValueAsDouble( "decimationTolerance", station ),
SP3EphList.getInitialTime() );
pList.push_back(decimateData);
BasicModel basic(nominalPos, SP3EphList);
basic.setMinElev(confReader.getValueAsDouble("cutOffElevation",station));
if ( !usingC1 )
{
basic.setDefaultObservable(TypeID::P1);
}
pList.push_back(basic);
EclipsedSatFilter eclipsedSV;
pList.push_back(eclipsedSV);
GravitationalDelay grDelay(nominalPos);
pList.push_back(grDelay);
double uARP(confReader.fetchListValueAsDouble( "offsetARP", station ) );
double eARP(confReader.fetchListValueAsDouble( "offsetARP", station ) );
double nARP(confReader.fetchListValueAsDouble( "offsetARP", station ) );
Triple offsetARP( uARP, eARP, nARP );
Triple offsetL1( 0.0, 0.0, 0.0 ), offsetL2( 0.0, 0.0, 0.0 );
AntexReader antexReader;
Antenna receiverAntenna;
bool useantex( confReader.getValueAsBoolean( "useAntex", station ) );
if( useantex )
{
antexReader.open( confReader.getValue( "antexFile", station ) );
receiverAntenna =
antexReader.getAntenna( confReader.getValue( "antennaModel",
station ) );
}
ComputeSatPCenter svPcenter(nominalPos);
if( useantex )
{
svPcenter.setAntexReader( antexReader );
}
pList.push_back(svPcenter);
CorrectObservables corr(SP3EphList);
corr.setNominalPosition(nominalPos);
corr.setMonument( offsetARP );
bool usepatterns(confReader.getValueAsBoolean("usePCPatterns", station ));
if( useantex && usepatterns )
{
corr.setAntenna( receiverAntenna );
corr.setUseAzimuth(confReader.getValueAsBoolean("useAzim", station));
}
else
{
offsetL1[0] = confReader.fetchListValueAsDouble("offsetL1", station);
offsetL1[1] = confReader.fetchListValueAsDouble("offsetL1", station);
offsetL1[2] = confReader.fetchListValueAsDouble("offsetL1", station);
offsetL2[0] = confReader.fetchListValueAsDouble("offsetL2", station);
offsetL2[1] = confReader.fetchListValueAsDouble("offsetL2", station);
offsetL2[2] = confReader.fetchListValueAsDouble("offsetL2", station);
corr.setL1pc( offsetL1 );
corr.setL2pc( offsetL2 );
}
pList.push_back(corr);
ComputeWindUp windup( SP3EphList,
nominalPos,
confReader.getValue( "satDataFile", station ) );
pList.push_back(windup);
NeillTropModel neillTM( nominalPos.getAltitude(),
nominalPos.getGeodeticLatitude(),
confReader.getValueAsInt("dayOfYear", station) );
double drytropo( neillTM.dry_zenith_delay() );
ComputeTropModel computeTropo(neillTM);
pList.push_back(computeTropo);
ComputeLinear linear2;
if ( usingC1 )
{
linear2.addLinear(comb.pcCombWithC1);
}
else
{
linear2.addLinear(comb.pcCombination);
}
linear2.addLinear(comb.lcCombination);
pList.push_back(linear2);
SimpleFilter pcFilter;
pcFilter.setFilteredType(TypeID::PC);
bool filterPC( confReader.getValueAsBoolean( "filterPC", station ) );
if( filterPC )
{
pList.push_back(pcFilter);
}
PhaseCodeAlignment phaseAlign;
pList.push_back(phaseAlign);
ComputeLinear linear3(comb.pcPrefit);
linear3.addLinear(comb.lcPrefit);
pList.push_back(linear3);
XYZ2NEU baseChange(nominalPos);
pList.push_back(baseChange);
ComputeDOP cDOP;
pList.push_back(cDOP);
bool isNEU( confReader.getValueAsBoolean( "USENEU", station ) );
SolverPPP pppSolver(isNEU);
SolverPPPFB fbpppSolver(isNEU);
int cycles( confReader.getValueAsInt("forwardBackwardCycles", station) );
bool isWN( confReader.getValueAsBoolean( "coordinatesAsWhiteNoise",
station ) );
WhiteNoiseModel wnM(100.0);
if ( cycles > 0 )
{
if ( isWN )
{
fbpppSolver.setCoordinatesModel(&wnM);
}
pList.push_back(fbpppSolver);
}
else
{
if ( isWN )
{
pppSolver.setCoordinatesModel(&wnM);
}
pList.push_back(pppSolver);
}
SolidTides solid;
OceanLoading ocean;
ocean.setFilename( confReader.getValue( "oceanLoadingFile", station ) );
double xp( confReader.fetchListValueAsDouble( "poleDisplacements",
station ) );
double yp( confReader.fetchListValueAsDouble( "poleDisplacements",
station ) );
PoleTides pole;
pole.setXY( xp, yp );
gnssRinex gRin;
int precision( confReader.getValueAsInt( "precision", station ) );
string outName(confReader.getValue( "outputFile", station ) );
ofstream outfile;
outfile.open( outName.c_str(), ios::out );
bool printmodel( confReader.getValueAsBoolean( "printModel", station ) );
string modelName;
ofstream modelfile;
if( printmodel )
{
modelName = confReader.getValue( "modelFile", station );
modelfile.open( modelName.c_str(), ios::out );
}
while(rin >> gRin)
{
DayTime time(gRin.header.epoch);
Triple tides( solid.getSolidTide( time, nominalPos ) +
ocean.getOceanLoading( station, time ) +
pole.getPoleTide( time, nominalPos ) );
corr.setExtraBiases(tides);
try
{
gRin >> pList;
}
catch(DecimateEpoch& d)
{
continue;
}
catch(Exception& e)
{
cerr << "Exception for receiver '" << station <<
"' at epoch: " << time << "; " << e << endl;
continue;
}
catch(...)
{
cerr << "Unknown exception for receiver '" << station <<
" at epoch: " << time << endl;
continue;
}
if ( printmodel )
{
printModel( modelfile,
gRin );
}
if ( cycles < 1 )
{
printSolution( outfile,
pppSolver,
time,
cDOP,
isNEU,
gRin.numSats(),
drytropo,
precision );
}
}
rin.close();
if ( printmodel )
{
modelfile.close();
}
SP3EphList.clear();
if ( cycles < 1 )
{
outfile.close();
cout << "Processing finished for station: '" << station
<< "'. Results in file: '" << outName << "'." << endl;
continue;
}
try
{
fbpppSolver.ReProcess(cycles);
}
catch(Exception& e)
{
cerr << "Exception at reprocessing phase: " << e << endl;
cerr << "Skipping receiver '" << station << "'." << endl;
outfile.close();
continue;
}
while( fbpppSolver.LastProcess(gRin) )
{
DayTime time(gRin.header.epoch);
printSolution( outfile,
fbpppSolver,
time,
cDOP,
isNEU,
gRin.numSats(),
drytropo,
precision );
}
outfile.close();
cout << "Processing finished for station: '" << station
<< "'. Results in file: '" << outName << "'." << endl;
}
return;
}
int main(int argc, char* argv[])
{
try
{
example9 program(argv[0]);
if ( !program.initialize(argc, argv, false) )
{
return 0;
}
if ( !program.run() )
{
return 1;
}
return 0;
}
catch(Exception& e)
{
cerr << "Problem: " << e << endl;
return 1;
}
catch(...)
{
cerr << "Unknown error." << endl;
return 1;
}
return 0;
}