editdds.cpp

linux的gps应用

#pragma ident "$Id: EditDDs.cpp 185 2006-10-05 18:21:39Z 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
//
//============================================================================

//============================================================================
//
//This software developed by Applied Research Laboratories at the University of
//Texas at Austin, under contract to an agency or agencies within the U.S. 
//Department of Defense. The U.S. Government retains all rights to use,
//duplicate, distribute, disclose, or release this software. 
//
//Pursuant to DoD Directive 523024 
//
// DISTRIBUTION STATEMENT A: This software has been approved for public 
//                           release, distribution is unlimited.
//
//=============================================================================

/**
 * @file EditDDs.cpp
 * Edit buffered double differences for outliers, cycle slips and isolated
 * points for program DDBase.
 */

//------------------------------------------------------------------------------------
// TD EditDDs.cpp make various numbers input parameters

//------------------------------------------------------------------------------------
// includes
// system
#include <vector>

// GPSTk
#include "Matrix.hpp"
#include "Stats.hpp"
#include "RobustStats.hpp"
//#include "SRIFilter.hpp"

// DDBase
#include "DDBase.hpp"

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

//------------------------------------------------------------------------------------
static int ngood,nbad;           // number good data, number of data marked bad
static vector<int> mark;         // parallel to count and data vectors, mark bad data
static ofstream tddofs;          // output stream for OutputTDDFile

//------------------------------------------------------------------------------------
// prototypes -- this module only
int EditDDResets(const DDid& ddid, DDData& dddata);
int EditDDIsolatedPoints(const DDid& ddid, DDData& dddata);
int EditDDSlips(const DDid& ddid, DDData& dddata, int frequency);
int EditDDOutliers(const DDid& ddid, DDData& dddata, int frequency);
//void LSPolyFunc(Vector<double>& X, Vector<double>& f, Matrix<double>& P);
// prototypes -- DataOutput.cpp
int OutputRawDData(const DDid& ddid, const DDData& dddata, const vector<int>& mark);
int OutputDDData(void);

//------------------------------------------------------------------------------------
int EditDDs(void)
{
try {
   if(CI.Verbose) oflog << "BEGIN EditDDs()" << endl;

   if(!CI.OutputTDDFile.empty()) {
      tddofs.open(CI.OutputTDDFile.c_str(),ios::out);
      if(tddofs) {
         oflog << "Opened file " << CI.OutputTDDFile
               << " for triple difference and cycle slip output." << endl;
         tddofs << "# " << Title << endl;
         tddofs << "TDS  site site  sat sat freq iter cnt week sow       "
                << "dcnt     TD(m)       slip(cy) frac\n";
         tddofs << "SED  site site  sat sat freq iter cnt week sow           "
                << "DDres(m)    TDres(m)\n";
      }
      else {
         oflog << "Warning - Failed to open file " << CI.OutputTDDFile << endl;
      }
   }
   
   if(CI.Verbose) {
      oflog << " TUR site site sat sat iter  N Average  StdDev   SigYX"
         << "   Median   M-est     MAD\n";
      oflog << " SUR site site sat sat iter  N Average  StdDev   SigYX"
         << "   Median   M-est     MAD\n";
   }

   int i,j,k;
   map<DDid,DDData>::iterator it;

      // -------------------------------------------------------------------
      // delete DD buffers that are too small, or that user wants to exclude
      // also compute maxCount, the largest value of Count seen in all baselines
   maxCount = 0;
   vector<DDid> DDdelete;
   for(it = DDDataMap.begin(); it != DDDataMap.end(); it++) {

         // is it too small?
      if(it->second.count.size() < CI.MinDDSeg) {
         DDdelete.push_back(it->first);
         continue;
      }

         // prepare 'mark' vector
      mark.assign(it->second.count.size(),1);
      ngood = mark.size();
      nbad = 0;

         // remove points where bias had to be reset multiple times
      k = EditDDResets(it->first, it->second);
      if(k || ngood < CI.MinDDSeg) {
         DDdelete.push_back(it->first);
         continue;
      }

         // remove isolated points
      k = EditDDIsolatedPoints(it->first, it->second);
      if(k || ngood < CI.MinDDSeg) {
         DDdelete.push_back(it->first);
         continue;
      }

         // find and remove slips
      if(CI.Frequency != 2) {                // L1
         k = EditDDSlips(it->first, it->second,1);
         if(k || ngood < CI.MinDDSeg) {
            DDdelete.push_back(it->first);
            continue;
         }
      }
      if(CI.Frequency != 1) {                // L2
         k = EditDDSlips(it->first, it->second,2);
         if(k || ngood < CI.MinDDSeg) {
            DDdelete.push_back(it->first);
            continue;
         }
      }

         // find and remove outliers
      if(CI.Frequency != 2) {                // L1
         k = EditDDOutliers(it->first, it->second,1);
         if(k || ngood < CI.MinDDSeg) {
            DDdelete.push_back(it->first);
            continue;
         }
      }
      if(CI.Frequency != 1) {                // L2
         k = EditDDOutliers(it->first, it->second,2);
         if(k || ngood < CI.MinDDSeg) {
            DDdelete.push_back(it->first);
            continue;
         }
      }

         // output raw data with mark
      OutputRawDData(it->first, it->second, mark);

         // use vector 'mark' to delete data
      if(nbad > 0) {
         vector<double> nDDL1,nDDL2,nDDP1,nDDP2,nDDER;
         vector<int> ncount;
         for(i=0; i<it->second.count.size(); i++) {
            if(mark[i] == 1) {
               nDDL1.push_back(it->second.DDL1[i]);
               nDDL2.push_back(it->second.DDL2[i]);
               nDDP1.push_back(it->second.DDP1[i]);
               nDDP2.push_back(it->second.DDP2[i]);
               nDDER.push_back(it->second.DDER[i]);
               ncount.push_back(it->second.count[i]);
            }
         }
         it->second.DDL1 = nDDL1;
         it->second.DDL2 = nDDL2;
         it->second.DDP1 = nDDP1;
         it->second.DDP2 = nDDP2;
         it->second.DDER = nDDER;
         it->second.count = ncount;
         // ignore resets from now on...
      }

         // find the max count
      if(it->second.count[it->second.count.size()-1] > maxCount)
         maxCount = it->second.count[it->second.count.size()-1];
   }

      // close the output file
   tddofs.close();
   mark.clear();

      // now delete the ones that were marked
   for(i=0; i<DDdelete.size(); i++) {
      if(CI.Verbose) oflog << setw(2) << DDdelete[i]
         << " total = " << setw(5) << DDDataMap[DDdelete[i]].count.size()
         << ", count = " << setw(5) << DDDataMap[DDdelete[i]].count[0]
         << " - " << setw(5)
         << DDDataMap[DDdelete[i]].count[DDDataMap[DDdelete[i]].count.size()-1]
         << " -- Delete this DD dataset."
         << endl;
      DDDataMap.erase(DDdelete[i]);
   }
   DDdelete.clear();

      // -------------------------------------------------------------------
      // output DD summary
   if(CI.Screen) {
      cout << "Double differences summary:" << endl;
      for(k=1,it=DDDataMap.begin(); it!=DDDataMap.end(); it++,k++) {
         cout << " " << setw(2) << k << " " << it->first
            << " " << setw(5) << it->second.count.size()
            << " " << setw(5) << it->second.count[0]
            << " - " << setw(5) << it->second.count[it->second.count.size()-1];
         for(i=0; i<it->second.count.size()-1; i++) {
            j = it->second.count.at(i+1) - it->second.count.at(i);
            if(j > 1)
               cout << " (" << it->second.count.at(i)+1 << ":" << j-1 << ")";
         }
         cout << endl;
      }
   }
   if(CI.Verbose) {
      oflog << "Double differences summary:" << endl;
      for(k=1,it=DDDataMap.begin(); it!=DDDataMap.end(); it++,k++) {
            // output
         oflog << " " << setw(2) << k << " " << it->first
            << " " << setw(5) << it->second.count.size()
            << " " << setw(5) << it->second.count[0]
            << " - " << setw(5) << it->second.count[it->second.count.size()-1];
               // gaps - (count : number of pts)
         for(i=0; i<it->second.count.size()-1; i++) {
            j = it->second.count.at(i+1) - it->second.count.at(i);
            if(j > 1)
               oflog << " (" << it->second.count.at(i)+1 << ":" << j-1 << ")";
         }
         oflog << endl;
      }
   }

      // dump buffers to a file
   OutputDDData();

   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); }
}   // end EditDDs()

//------------------------------------------------------------------------------------
// There is no provision in DDBase for resetting a bias. This would imply
// solving for different biases (separated in time) for the same DDid.
// Therefore, this routine simply deletes all but the largest unbroken segment
// separated by resets.
int EditDDResets(const DDid& ddid, DDData& dddata)
{
try {
   int i,j,ibeg,iend;

   // resets[0] will always be the initial count
   if(dddata.resets.size() <= 1) return 0;

   oflog << " Warning - DD " << ddid << " had " << dddata.resets.size()-1
      << " resets between " << dddata.count[1]
      << " and " << dddata.count[dddata.count.size()-1] << " :";
   for(i=1; i<dddata.resets.size(); i++)
      oflog << " " << dddata.count[dddata.resets[i]]
         << "[" << dddata.resets[i] << "]";
   oflog << endl;

   //for(i=1; i<dddata.resets.size(); i++) {
   //   // difference in index
   //   int di = dddata.resets[i] - dddata.resets[i-1];
   //   // difference in counts
   //   int dc = dddata.count[dddata.resets[i]] - dddata.count[dddata.resets[i-1]];
   //   j = dddata.resets[i];
   //   // mark it bad
   //   if(dc < 12 && mark[j]==1) {           // TD make 12 an input parameter
   //      mark[j] = 0;
   //      ngood--;
   //      nbad++;
   //   }
   //}

   // find the largest segment between resets
   // NB this assumes nothing yet 'marked'
   ibeg = 0;
   iend = dddata.resets[1];
   for(i=2; i<=dddata.resets.size(); i++) {
      if(i == dddata.resets.size()) // last point
         j = dddata.count.size();
      else
         j = dddata.resets[i];

      if(j - dddata.resets[i-1] > iend-ibeg) {
         ibeg = dddata.resets[i-1];
         iend = j;
      }
   }

   if(CI.Verbose) oflog << " Delete data due to reset for DD " << ddid
      << " in the range " << ibeg << " to " << iend << endl;

      // mark all points from beginning to just before the 'ibeg' reset
   for(i=0; i<ibeg; i++) if(mark[i]==1) {
      mark[i] = 0;
      ngood--;
      nbad++;
   }
   
      // mark all points from 'iend' reset to the end
   for(i=iend; i<dddata.count.size(); i++) if(mark[i]==1) {
      mark[i] = 0;
      ngood--;
      nbad++;
   }

   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); }
}

//------------------------------------------------------------------------------------
int EditDDIsolatedPoints(const DDid& ddid, DDData& dddata)
{
try {
   int i,j,gappast,gapfuture;

   // loop over all counts
   // i is current (good) point, j is the next good point
   i = 0; while(i<dddata.count.size() & mark[i]==0) i++;     // find first good pt

   gapfuture = CI.MaxGap;
   while(i < dddata.count.size()) {
      gappast = gapfuture;

      j = i+1; while(j<dddata.count.size() & mark[j]==0) j++;// find next good pt

      if(j < dddata.count.size()) gapfuture = dddata.count[j] - dddata.count[i];
      else                        gapfuture = CI.MaxGap;

      if(gappast >= CI.MaxGap && gapfuture >= CI.MaxGap) {
         if(CI.Verbose) oflog << " Mark isolated " << ddid
            << " " << dddata.count[i] << endl;
         mark[i] = 0;
         ngood--;
         nbad++;
      }

      i = j;
   }

   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); }
}

//------------------------------------------------------------------------------------
int EditDDSlips(const DDid& ddid, DDData& dddata, int frequency)
{
try {
   int i,j,k,n,m,tdcount,tddt,ii,iter;
   double slip,fslip,tol;
   vector<int> slipindex;
   vector<double> td,slipsize;
   DayTime tt;
   TwoSampleStats<double> tsstats;

      // -------------------------------------- define td tolerance for slips
   const int iter_limit = 3;       // this allows iter_limit-1 slips to be fixed