1 //==============================================================================
2 //
3 // This file is part of GPSTk, the GPS Toolkit.
4 //
5 // The GPSTk is free software; you can redistribute it and/or modify
6 // it under the terms of the GNU Lesser General Public License as published
7 // by the Free Software Foundation; either version 3.0 of the License, or
8 // any later version.
9 //
10 // The GPSTk is distributed in the hope that it will be useful,
11 // but WITHOUT ANY WARRANTY; without even the implied warranty of
12 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 // GNU Lesser General Public License for more details.
14 //
15 // You should have received a copy of the GNU Lesser General Public
16 // License along with GPSTk; if not, write to the Free Software Foundation,
17 // Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110, USA
18 //
19 // This software was developed by Applied Research Laboratories at the
20 // University of Texas at Austin.
21 // Copyright 2004-2020, The Board of Regents of The University of Texas System
22 //
23 //==============================================================================
24
25 //==============================================================================
26 //
27 // This software was developed by Applied Research Laboratories at the
28 // University of Texas at Austin, under contract to an agency or agencies
29 // within the U.S. Department of Defense. The U.S. Government retains all
30 // rights to use, duplicate, distribute, disclose, or release this software.
31 //
32 // Pursuant to DoD Directive 523024
33 //
34 // DISTRIBUTION STATEMENT A: This software has been approved for public
35 // release, distribution is unlimited.
36 //
37 //==============================================================================
38
39 /// @file RinSum.cpp
40 /// Read Rinex observation files (version 2 or 3) and output a summary of the content.
41
42 // system
43 #include <string>
44 #include <vector>
45 #include <map>
46 #include <iostream>
47 #include <fstream>
48 #include <algorithm>
49
50 // GPSTK
51 #include "Exception.hpp"
52 #include "StringUtils.hpp"
53 #include "GNSSconstants.hpp"
54
55 #include "singleton.hpp"
56 #include "expandtilde.hpp"
57 #include "stl_helpers.hpp"
58 #include "logstream.hpp"
59 #include "CommandLine.hpp"
60
61 #include "CommonTime.hpp"
62 #include "CivilTime.hpp"
63 #include "Epoch.hpp"
64 #include "TimeString.hpp"
65
66 #include "RinexSatID.hpp"
67 #include "RinexObsID.hpp"
68 #include "Rinex3ObsStream.hpp"
69 #include "Rinex3ObsHeader.hpp"
70 #include "Rinex3ObsData.hpp"
71 #include "RinexUtilities.hpp"
72
73 #include "msecHandler.hpp"
74
75 //-----------------------------------------------------------------------------
76 using namespace std;
77 using namespace gpstk;
78 using namespace StringUtils;
79
80 //-----------------------------------------------------------------------------
81 string Version(string("4.1 8/26/15"));
82
83 //-----------------------------------------------------------------------------
84 // Object for command line input and global data
85 class Configuration : public Singleton<Configuration>
86 {
87
88 public:
89
90 // Default and only constructor
Configuration()91 Configuration() throw() { SetDefaults(); }
92
93 // Create, parse and process command line options and user input
94 int ProcessUserInput(int argc, char **argv) throw();
95
96 // Design the command line
97 string BuildCommandLine(void) throw();
98
99 // Open the output file, and parse the strings used on the command line
100 // return -4 if log file could not be opened
101 int ExtraProcessing(string& errors, string& extras) throw();
102
103 private:
104
105 // Define default values
SetDefaults(void)106 void SetDefaults(void) throw()
107 {
108 defaultstartStr = string("[Beginning of dataset]");
109 defaultstopStr = string("[End of dataset]");
110 beginTime = CommonTime::BEGINNING_OF_TIME;
111 endTime = CommonTime::END_OF_TIME;
112 beginTime.setTimeSystem(TimeSystem::Any);
113 endTime.setTimeSystem(TimeSystem::Any);
114 userfmt = gpsfmt;
115 help = verbose = brief = nohead = notab = gpstime = sorttime = vistab
116 = dogaps = doms = ycode = quiet = false;
117 debug = -1;
118 dt = -1.0;
119 doCurrRversion = false;
120 Rversion = 0.0;
121 vres = 0;
122 } // end Configuration::SetDefaults()
123
124 public:
125
126 // member data
127 CommandLine opts; // command line options and syntax page
128 static const string PrgmName; // program name
129 string Title; // id line printed to screen and log
130
131 // start command line input
132 bool help, verbose, brief, nohead, notab, gpstime, sorttime, dogaps, doms,
133 vistab, ycode, quiet;
134 int debug, vres;
135 double dt;
136 double Rversion; // RINEX version of output (default=header.version)
137 bool doCurrRversion; // set Rversion to Rinex3ObsBase::currentVersion
138 string cfgfile, userfmt;
139
140 vector<string> InputObsFiles; // RINEX obs file names
141 string Obspath; // paths
142 string LogFile; // output log file (optional)
143
144 // times derived from --start and --stop
145 string defaultstartStr,startStr;
146 string defaultstopStr,stopStr;
147 CommonTime beginTime,endTime;
148 vector<RinexSatID> exSats;
149 vector<RinexSatID> onlySats;
150
151 // end of command line input
152
153 vector<int> gapcount; // for counting gaps
154 string msg;
155 static const string calfmt,gpsfmt,longfmt;
156 ofstream logstrm;
157
158 // for milliseconds
159 msecHandler msh;
160
161 }; // end class Configuration
162
163 //-----------------------------------------------------------------------------
164 // const members of Configuration
165 const string Configuration::PrgmName = string("RinSum");
166 const string Configuration::calfmt = string("%04Y/%02m/%02d %02H:%02M:%02S");
167 const string Configuration::gpsfmt = string("%4F %w %10.3g %P");
168 const string Configuration::longfmt = calfmt + " = " + gpsfmt;
169
170 //-----------------------------------------------------------------------------
171 // struct used to store SAT/Obs table
172 struct TableData
173 {
174 RinexSatID sat;
175 vector<int> nobs; // number of data for each obs type
176 vector<int> gapcount; //
177 double prevC1, prevP1, prevL1;
178 CommonTime begin, end;
179
180 // c'tor
TableDataTableData181 TableData(const SatID& p, const int& n)
182 {
183 sat = RinexSatID(p);
184 nobs = vector<int>(n);
185 prevC1 = 0;
186 prevP1 = 0;
187 prevL1 = 0;
188 }
189 // needed for find()
operator ==TableData190 inline bool operator==(const TableData& d) { return d.sat == sat; }
191 };
192
193 // needed for sort()
194 class TableSATLessThan
195 {
196 public:
operator ()(const TableData & d1,const TableData & d2)197 bool operator()(const TableData& d1, const TableData& d2)
198 { return d1.sat < d2.sat; }
199 };
200
201 class TableBegLessThan
202 {
203 public:
operator ()(const TableData & d1,const TableData & d2)204 bool operator()(const TableData& d1, const TableData& d2)
205 { return d1.begin < d2.begin; }
206 };
207
208 //-----------------------------------------------------------------------------
209 // prototypes
210 /**
211 * @throw Exception */
212 int Initialize(string& errors);
213 /**
214 * @throw Exception */
215 int ProcessFiles(void);
216
217 //-----------------------------------------------------------------------------
218 //-----------------------------------------------------------------------------
main(int argc,char ** argv)219 int main(int argc, char **argv)
220 {
221 // get the (single) instance of the configuration
222 Configuration& C(Configuration::Instance());
223
224 try
225 {
226 int iret;
227 clock_t totaltime(clock());
228 Epoch wallclkbeg;
229 wallclkbeg.setLocalTime();
230
231 // build title = first line of output
232 C.Title = "# " + C.PrgmName + ", part of the GPS Toolkit, Ver " +
233 Version + ", Run " + printTime(wallclkbeg,C.calfmt);
234
235 for(;;)
236 {
237 // get information from the command line
238 // iret -2 -3 -4
239 if((iret = C.ProcessUserInput(argc, argv)) != 0)
240 break;
241
242 // read stores, check input etc
243 string errs;
244 if((iret = Initialize(errs)) != 0)
245 {
246 LOG(ERROR) << "------- Input is not valid: ----------\n" << errs
247 << "\n------- end errors -----------";
248 break;
249 }
250 if(!errs.empty())
251 LOG(INFO) << errs; // Warnings are here too
252
253 ProcessFiles();
254 iret = 0; // successful completion.
255
256 break; // mandatory
257 }
258
259 if(iret >= 0 && !C.brief && !C.quiet)
260 {
261 // print elapsed time
262 totaltime = clock()-totaltime;
263 Epoch wallclkend;
264 wallclkend.setLocalTime();
265 ostringstream oss;
266 oss << C.PrgmName << " timing: processing " << fixed
267 << setprecision(3) << double(totaltime)/double(CLOCKS_PER_SEC)
268 << " sec, wallclock: " << setprecision(0)
269 << (wallclkend-wallclkbeg) << " sec.";
270 LOG(INFO) << oss.str();
271 }
272
273 return iret;
274 }
275 catch(FFStreamError& e)
276 {
277 cerr << "FFStreamError: " << e.what();
278 }
279 catch(Exception& e)
280 {
281 cerr << "Exception: " << e.what();
282 }
283 catch (...)
284 {
285 cerr << "Unknown exception. Abort." << endl;
286 }
287 return 1;
288
289 } // end main()
290
291 //-----------------------------------------------------------------------------
292 // return -5 if input is not valid
Initialize(string & errors)293 int Initialize(string& errors)
294 {
295 try
296 {
297 bool isValid(true);
298 Configuration& C(Configuration::Instance());
299 errors = string("");
300
301 // add path to filenames, and expand tilde (~)
302 include_path(C.Obspath, C.InputObsFiles);
303
304 // sort input obs files on time
305 if(C.InputObsFiles.size() > 1)
306 {
307 C.msg = sortRinex3ObsFiles(C.InputObsFiles);
308 if(!C.msg.empty())
309 LOG(ERROR) << C.msg;
310 }
311
312 // initialize millisecond handler
313 if(C.doms)
314 C.msh.setDT(C.dt);
315
316 // -------- save errors and output
317 //errors = oss.str();
318 //stripTrailing(errors,'\n');
319 //replaceAll(errors,"\n","\n# ");
320
321 if(!isValid)
322 return -5;
323 return 0;
324 }
325 catch(Exception& e)
326 {
327 GPSTK_RETHROW(e);
328 }
329 } // end Initialize()
330
331 //-----------------------------------------------------------------------------
ProcessUserInput(int argc,char ** argv)332 int Configuration::ProcessUserInput(int argc, char **argv) throw()
333 {
334 string PrgmDesc,cmdlineUsage, cmdlineErrors, cmdlineExtras;
335 vector<string> cmdlineUnrecognized;
336
337 // build the command line
338 opts.DefineUsageString(PrgmName + " [options]");
339 PrgmDesc = BuildCommandLine();
340
341 // let CommandLine parse options; write all errors, etc to the passed strings
342 int iret = opts.ProcessCommandLine(argc, argv, PrgmDesc,
343 cmdlineUsage, cmdlineErrors, cmdlineUnrecognized);
344
345 // handle return values
346 if(iret == -2)
347 return iret; // bad alloc
348 if(iret == -3)
349 return iret; // invalid command line
350
351 // define verbose and debug
352 verbose = (LOGlevel >= VERBOSE);
353 debug = (LOGlevel - DEBUG);
354
355 // help: print syntax page and quit
356 if(opts.hasHelp())
357 {
358 LOG(INFO) << Title;
359 LOG(INFO) << cmdlineUsage;
360 return 1;
361 }
362
363 // extra parsing (perhaps add to cmdlineErrors, cmdlineExtras)
364 iret = ExtraProcessing(cmdlineErrors, cmdlineExtras);
365 if(iret == -4)
366 return iret; // log file could not be opened
367
368 // pull out file name without --obs
369 if(cmdlineUnrecognized.size() > 0)
370 {
371 for(int i=cmdlineUnrecognized.size()-1; i >= 0; i--)
372 {
373 try
374 {
375 string filename(cmdlineUnrecognized[i]);
376 ifstream ifstrm(filename.c_str());
377 if(ifstrm.is_open())
378 {
379 LOG(DEBUG) << "# Deduce filename >" << filename << "<";
380 InputObsFiles.push_back(cmdlineUnrecognized[i]);
381 cmdlineUnrecognized.erase(cmdlineUnrecognized.begin()+i);
382 ifstrm.close();
383 continue;
384 }
385 else
386 {
387 include_path(Obspath, filename);
388 ifstrm.open(filename.c_str());
389 if(ifstrm.is_open())
390 {
391 LOG(DEBUG) << "# Deduce filename >" << filename << "<";
392 InputObsFiles.push_back(cmdlineUnrecognized[i]);
393 cmdlineUnrecognized.erase(cmdlineUnrecognized.begin()+i);
394 ifstrm.close();
395 continue;
396 }
397 }
398 }
399 catch(Exception& e) { ; }
400 }
401 }
402
403 // output warning / error messages
404 if(cmdlineUnrecognized.size() > 0)
405 {
406 LOG(WARNING) << "Warning - unrecognized arguments:";
407 for(size_t j=0; j<cmdlineUnrecognized.size(); j++)
408 LOG(WARNING) << " " << cmdlineUnrecognized[j];
409 LOG(WARNING) << "End of unrecognized arguments";
410 }
411
412 // fatal errors
413 if(!cmdlineErrors.empty())
414 {
415 stripTrailing(cmdlineErrors,'\n');
416 replaceAll(cmdlineErrors,"\n","\n ");
417 LOG(ERROR) << "Errors found on command line:\n " << cmdlineErrors
418 << "\nEnd of command line errors.";
419 return 1;
420 }
421
422 // success: dump configuration summary
423 if(verbose)
424 {
425 ostringstream oss;
426 oss << "------ Summary of " << PrgmName
427 << " command line configuration ------\n";
428 opts.DumpConfiguration(oss);
429 //if(!cmdlineExtras.empty()) oss << "# Extra Processing:\n" << cmdlineExtras;
430 oss << "------ End configuration summary ------";
431 LOG(VERBOSE) << oss.str();
432 }
433 if(!cmdlineExtras.empty())
434 {
435 stripTrailing(cmdlineExtras,'\n');
436 LOG(INFO) << cmdlineExtras;
437 }
438
439 return 0;
440
441 } // end Configuration::CommandLine()
442
443 //-----------------------------------------------------------------------------
BuildCommandLine(void)444 string Configuration::BuildCommandLine(void) throw()
445 {
446 // Program description will appear at the top of the syntax page
447 string PrgmDesc = " Program " + PrgmName + " reads one or more RINEX (v.2+) "
448 + "observation files and prints a summary of content."
449 + "\n Note verbose dumps all auxiliary headers."
450 + "\n Options:";
451 opts.DefineUsageString("RinSum <file> [options]");
452
453 // options to appear on the syntax page, and to be accepted on command line
454 //opts.Add(char, opt, arg, repeat?, required?, &target, pre-desc, desc);
455 // NB cfgfile is a dummy, but it must exist when cmdline is processed.
456 opts.Add('f', "file", "fn", true, false, &cfgfile,
457 "# Input via configuration file:",
458 "Name of file with more options [#->EOL = comment]");
459
460 opts.Add(0, "obs", "file", true, false, &InputObsFiles,
461 "# Required input [file names may appear w/o --obs if unambiguous]",
462 "Input RINEX observation file name");
463 opts.Add(0, "obspath", "p", false, false, &Obspath,
464 "# Paths of input files (optional):",
465 "Path of input RINEX observation file(s)");
466
467 startStr = defaultstartStr;
468 stopStr = defaultstopStr;
469 opts.Add(0, "start", "t[:f]", false, false, &startStr,
470 "# Editing (t,f are strings: time t; format f "
471 "defaults to wk,sow OR yr,mon,day,h,m,s",
472 "Start processing data at this epoch");
473 opts.Add(0, "stop", "t[:f]", false, false, &stopStr, "",
474 "Stop processing data at this epoch");
475 opts.Add(0, "exSat", "sat", true, false, &exSats, "",
476 "Exclude satellite (or system) <sat> e.g. G24,R");
477 opts.Add(0, "onlySat", "sat", true, false, &onlySats, "",
478 "Include ONLY satellites (or systems) <sat> e.g. G,R");
479
480
481 string rvstr("# Rinex version (current, latest version is "
482 + asString(Rinex3ObsBase::currentVersion,2) + "):");
483 opts.Add(0, "currentRinex", "", false, false, &doCurrRversion, rvstr,
484 "Output in current, not header, RINEX version");
485 opts.Add(0, "RinexVer", "V", false, false, &Rversion, "",
486 "Output in RINEX version V (default is header.version)");
487
488 opts.Add(0, "timefmt", "fmt", false, false, &userfmt, "# Output:",
489 "Format for time tags (see GPSTK::Epoch::printf) in output");
490 opts.Add('b', "brief", "", false, false, &brief, "",
491 "Produce a brief output");
492 opts.Add(0, "nohead", "", false, false, &nohead, "",
493 "Omit header from output");
494 opts.Add(0, "notable", "", false, false, ¬ab, "",
495 "Omit sat/obs table from output");
496 opts.Add(0, "dt", "sec", false, false, &dt, "",
497 "Nominal time step of data (sec); required only for gaps and millisec");
498 opts.Add(0, "milli", "", false, false, &doms, "",
499 "Find millisecond clock adjusts; req's --dt");
500 opts.Add(0, "gaps", "", false, false, &dogaps, "",
501 "Print a table of gaps in the data; req's --dt");
502 opts.Add(0, "vis", "n", false, false, &vres, "",
503 "Print graphical visibility, resolution <n> [n~20 @ 30s; req's --gaps]");
504 opts.Add(0, "vtab", "", false, false, &vistab, "",
505 "Print tabular visibility [req's --gaps and --vis]");
506
507 opts.Add(0, "ycode", "", false, false, &ycode, "# Other:",
508 "Assume v2.11 P mean Y");
509 opts.Add('q', "quiet", "", false, false, &quiet, "",
510 "Make output a little quieter");
511
512 // CommandLine adds automatically
513 //opts.Add(0, "verbose", "", false, false, &verbose, "# Help:",
514 // "Print extra output information");
515 //opts.Add(0, "debug", "", false, false, &debug, "",
516 // "Print debug output at level 0 [debug<n> for level n=1-7]");
517 //opts.Add(0, "help", "", false, false, &help, "",
518 // "Print this syntax page, and quit");
519
520 // deprecated (old,new) e.g.
521 //opts.Add_deprecated("--SP3","--eph");
522
523 return PrgmDesc;
524
525 } // end Configuration::BuildCommandLine()
526
527 //-----------------------------------------------------------------------------
ExtraProcessing(string & errors,string & extras)528 int Configuration::ExtraProcessing(string& errors, string& extras) throw()
529 {
530 int n;
531 size_t i;
532 vector<string> fld;
533 ostringstream oss,ossx; // oss for Errors, ossx for Warnings and info
534
535 // start and stop times
536 for(i=0; i<2; i++)
537 {
538 static const string fmtGPS("%F,%g"),fmtCAL("%Y,%m,%d,%H,%M,%S");
539 msg = (i==0 ? startStr : stopStr);
540 if(msg == (i==0 ? defaultstartStr : defaultstopStr))
541 continue;
542
543 bool ok(true);
544 bool hasfmt(msg.find('%') != string::npos);
545 n = numWords(msg,',');
546 if(hasfmt)
547 {
548 fld = split(msg,':');
549 if(fld.size() != 2)
550 {
551 ok = false;
552 }
553 else
554 {
555 try
556 {
557 Epoch ep;
558 stripLeading(fld[0]," \t");
559 stripLeading(fld[1]," \t");
560 ep.scanf(fld[0],fld[1]);
561 (i==0 ? beginTime : endTime) = static_cast<CommonTime>(ep);
562 }
563 catch(Exception& e)
564 {
565 ok = false;
566 LOG(INFO) << "excep " << e.what();
567 }
568 }
569 }
570 else if(n == 2 || n == 6)
571 {
572 // try the defaults
573 try
574 {
575 Epoch ep;
576 ep.scanf(msg,(n==2 ? fmtGPS : fmtCAL));
577 (i==0 ? beginTime : endTime) = static_cast<CommonTime>(ep);
578 }
579 catch(Exception& e)
580 {
581 ok = false;
582 LOG(INFO) << "excep " << e.what();
583 }
584 }
585
586 if(ok)
587 {
588 msg = printTime((i==0 ? beginTime : endTime),fmtGPS+" = "+fmtCAL);
589 if(msg.find("Error") != string::npos) ok = false;
590 }
591
592 if(!ok)
593 oss << "Error : invalid time or format in --" << (i==0 ? "start" : "stop")
594 << " " << (i==0 ? startStr : stopStr) << endl;
595 else
596 {
597 (i==0 ? beginTime : endTime).setTimeSystem(TimeSystem::Any);
598 ossx << (i==0 ? "Begin time --begin" : "End time --end") << " is "
599 << printTime((i==0 ? beginTime : endTime), fmtGPS+" = "+fmtCAL+"\n");
600 }
601 }
602
603 // open the log file (so warnings, configuration summary, etc can go there) -----
604 if(!LogFile.empty())
605 {
606 logstrm.open(LogFile.c_str(), ios::out);
607 if(!logstrm.is_open())
608 {
609 LOG(ERROR) << "Error : Failed to open log file " << LogFile;
610 return -4;
611 }
612 LOG(INFO) << "Output redirected to log file " << LogFile;
613 pLOGstrm = &logstrm;
614 }
615 if (!quiet)
616 LOG(INFO) << Title;
617
618 // check consistency of exSat and onlySat; note you CAN have --only R --ex R10,R07
619 if(exSats.size() > 0 && onlySats.size() > 0)
620 {
621 for(i=0; i<onlySats.size(); i++)
622 {
623 RinexSatID sat(onlySats[i]);
624 RinexSatID sys(-1,sat.system);
625 if((find(exSats.begin(), exSats.end(), sat) != exSats.end()) ||
626 (find(exSats.begin(), exSats.end(), sys) != exSats.end()))
627 oss << "Error : satellite " << asString(sat)
628 << " found in both --exSat and --onlySat" << endl;
629 }
630 }
631
632 // gaps and vis options
633 if(vres < 0)
634 {
635 ossx << "Warning - Option --vis, must have n positive\n";
636 vres = 0;
637 }
638 if(dt < 0.0 && dt != -1.0)
639 {
640 ossx << "Warning - Option --dt, must have dt positive\n";
641 dt = -1.0;
642 }
643 // milli requires dt
644 if(doms && dt == -1.0)
645 {
646 ossx << "Warning - Option --milli requires --dt be given\n";
647 doms = false;
648 }
649 // gaps requires dt
650 if(dogaps && dt == -1.0)
651 {
652 ossx << "Warning - Option --gaps requires --dt be given\n";
653 dogaps = false;
654 }
655 // vres requires dt and gaps
656 if(vres > 0 && (dt == -1.0 || !dogaps))
657 {
658 ossx << "Warning - Option --vis <n> requires --gaps and --dt be given\n";
659 vres = 0;
660 }
661 if(vistab && vres == 0)
662 {
663 ossx << "Warning - Option --vtab requires that --vis <n> be given\n";
664 vistab = false;
665 }
666
667 // add new errors to the list
668 msg = oss.str();
669 if(!msg.empty())
670 errors += msg;
671 msg = ossx.str(); //stripTrailing(msg,'\n');
672 if(!msg.empty())
673 extras += msg;
674
675 return 0;
676
677 } // end Configuration::ExtraProcessing(string& errors) throw()
678
679 //-----------------------------------------------------------------------------
680 // Return 0 ok, >0 number of files successfully read, <0 fatal error
ProcessFiles()681 int ProcessFiles()
682 {
683 try
684 {
685 Configuration& C(Configuration::Instance());
686 int iret,ii,k;
687 size_t i,j,nfile,nfiles;
688 string tag;
689 CommonTime lastObsTime, prevObsTime, firstObsTime;
690 RinexSatID sat;
691 Rinex3ObsStream ostrm;
692 ostringstream oss;
693
694 // output in header.version, unless user requests otherwise
695 double outputVersion(C.Rversion);
696 if(C.doCurrRversion) outputVersion = Rinex3ObsBase::currentVersion;
697
698 // estimate time step
699 const size_t ndtmax=15;
700 double dt, bestdt[ndtmax];
701 int ndt[ndtmax];
702 // cache the out-of-time-order records
703 bool cacheon;
704 vector<CommonTime> cachetime;
705 vector<vector<Rinex3ObsData> > cache;
706
707 for(nfiles=0,nfile=0; nfile<C.InputObsFiles.size(); nfile++)
708 {
709 Rinex3ObsStream istrm;
710 Rinex3ObsHeader Rhead, Rheadout;
711 Rinex3ObsData Rdata;
712
713 // If command line specified P1/P2 are to be considered
714 // as Y-code, set the Rinex3ObsHeader flag to indicate such.
715 if (C.ycode)
716 {
717 Rhead.PisY = true;
718 Rheadout.PisY = true;
719 }
720
721 cacheon = false;
722 cache.clear();
723 cachetime.clear();
724
725 string filename(C.InputObsFiles[nfile]);
726
727 // iret is set to 0 ok, or could not: 1 open file, 2 read header, 3 read data
728 iret = 0;
729 for(i=0; i<ndtmax; i++)
730 ndt[i] = -1;
731
732 // open the file ------------------------------------------------
733 istrm.open(filename.c_str(),ios::in);
734 if(!istrm.is_open())
735 {
736 LOG(WARNING) << "Warning : could not open file " << filename;
737 iret = 1;
738 continue;
739 }
740 istrm.exceptions(ios::failbit);
741
742 // get file size - on windows its different b/c of CRs
743 //char ch;
744 //istrm.seekg(0,ios::end);
745 //streampos filesize(istrm.tellg());
746 //istrm.seekg(0,ios::beg);
747
748 // output file name
749 if(C.quiet)
750 {
751 std::string choppedFN(filename);
752 choppedFN.erase(0,1+filename.find_last_of("/\\"));
753 LOG(INFO) << "+++++++++++++ " << C.PrgmName
754 << " summary of Rinex obs file " << choppedFN
755 << " +++++++++++++";
756 }
757 else if(!C.brief)
758 {
759 LOG(INFO) << "+++++++++++++ " << C.PrgmName
760 << " summary of Rinex obs file " << filename
761 << " +++++++++++++";
762 }
763
764 // read the header ----------------------------------------------
765 try
766 {
767 istrm >> Rhead;
768 }
769 catch(Exception& e)
770 {
771 LOG(WARNING) << "Warning : Failed to read header: " << e.what()
772 << "\n Header dump follows.";
773 Rhead.dump(LOGstrm);
774 istrm.close();
775 iret = 2;
776 continue;
777 }
778 if(Rhead.lastObs.getTimeSystem() != Rhead.firstObs.getTimeSystem())
779 Rhead.lastObs.setTimeSystem(Rhead.firstObs.getTimeSystem());
780
781 // output file name and header
782 if(C.brief)
783 {
784 if(nfile > 0)
785 LOG(INFO) << "";
786 LOG(INFO) << "File name: " << filename
787 << " (RINEX ver. " << Rhead.version << ")";
788 LOG(INFO) << "Marker name: " << Rhead.markerName;
789 LOG(INFO) << "Antenna type: " << Rhead.antType;
790 LOG(INFO) << "Position (XYZ,m) : " << fixed << setprecision(4)
791 << Rhead.antennaPosition << ".";
792 LOG(INFO) << "Antenna offset (UEN,m) : " << fixed << setprecision(4)
793 << Rhead.antennaDeltaHEN << ".";
794 }
795 else if(!C.nohead)
796 {
797 LOG(DEBUG) << "RINEX header:";
798 if(outputVersion == 0.0)
799 outputVersion = Rhead.version;
800 else
801 LOG(INFO) << " (Header has version " << Rhead.version
802 << "; output as version " << outputVersion << ")";
803 Rhead.dump(LOGstrm,outputVersion);
804 }
805
806 if(!Rhead.isValid())
807 {
808 LOG(INFO) << "Abort: header is invalid.";
809 if(C.quiet)
810 {
811 std::string choppedFN(filename);
812 choppedFN.erase(0,1+filename.find_last_of("/\\"));
813 LOG(INFO) << "\n+++++++++++++ End of RinSum summary of "
814 << choppedFN << " +++++++++++++";
815 }
816 else if(!C.brief)
817 {
818 LOG(INFO) << "\n+++++++++++++ End of RinSum summary of "
819 << filename << " +++++++++++++";
820 }
821 continue;
822 }
823
824 // initialize counting -------------------------------------------
825 int nepochs(0), nauxheads(0), nmaxobs(0);
826 vector<TableData> table; // table of counts per sat,obs
827 map<char, vector<int> > totals; // totals per system,obs
828
829 prevObsTime = CommonTime::BEGINNING_OF_TIME;
830 firstObsTime = CommonTime::BEGINNING_OF_TIME;
831
832 // initialize for all systems in the header
833 map<std::string,vector<RinexObsID> >::const_iterator sit; // used below often
834 for(sit=Rhead.mapObsTypes.begin(); sit != Rhead.mapObsTypes.end(); ++sit)
835 {
836 // Initialize the vectors contained in the map
837 totals[(sit->first)[0]] = vector<int>((sit->second).size());
838
839 LOG(DEBUG) << "GNSS " << (sit->first) << " is present with "
840 << (sit->second).size() << " observations...";
841
842 // find the max size of obs list
843 if(int((sit->second).size()) > nmaxobs)
844 nmaxobs = (sit->second).size();
845 }
846
847 // initialize millisecond handler with obstypes and wavelengths
848 vector<string> msots;
849 if(C.doms)
850 {
851 vector<double> waves;
852 // get obs types from header
853 for(sit=Rhead.mapObsTypes.begin(); sit != Rhead.mapObsTypes.end(); ++sit)
854 {
855 // get the system
856 RinexSatID rsid;
857 rsid.fromString(sit->first);
858 SatID sid(rsid);
859 // TD support only GPS currently
860 if(rsid.systemChar() != 'G') continue;
861 // excluded satellites/systems
862 if(find(C.exSats.begin(), C.exSats.end(), rsid) != C.exSats.end())
863 continue;
864 // get the obstypes, prepend the system character
865 for(i=0; i<sit->second.size(); i++)
866 {
867 tag = sit->second[i].asString(); // 3-char obs type
868 if(tag[0] == 'C' || tag[0] == 'L')
869 {
870 // code and phase only
871 msots.push_back(string(1,rsid.systemChar())+tag);
872 // get wavelength ... NB TD Glonass frequency channel not supported
873 if(tag[0] == 'L')
874 {
875 ii = asInt(string(1,tag[1]));
876 waves.push_back(getWavelength(sid.system, ii));
877 }
878 else
879 waves.push_back(0.0);
880 }
881 }
882 }
883
884 C.msh.setObstypes(msots,waves);
885 LOG(DEBUG) << "Initialize millisecond handler with obs type, wavelength:";
886 for(i=0; i<msots.size(); i++) LOG(DEBUG) << " " << msots[i]
887 << fixed << setprecision(6) << " " << waves[i];
888 }
889
890 if(pLOGstrm == &cout && !C.brief)
891 LOG(INFO) << "\nReading the observation data...";
892
893 // loop over epochs ---------------------------------------------
894 while(1)
895 {
896 try
897 {
898 istrm >> Rdata;
899 }
900 catch(Exception& e)
901 {
902 LOG(WARNING) << " Warning : Failed to read obs data (Exception "
903 << e.getText(0) << "); dump follows.";
904 Rdata.dump(LOGstrm,Rhead);
905 istrm.close();
906 iret = 3;
907 break;
908 }
909 catch(std::exception& e)
910 {
911 Exception ge(string("Std excep: ") + e.what());
912 GPSTK_THROW(ge);
913 }
914 catch(...)
915 {
916 Exception ue("Unknown exception while reading RINEX data.");
917 GPSTK_THROW(ue);
918 }
919
920 // normal EOF
921 if(!istrm.good() || istrm.eof())
922 {
923 iret = 0;
924 break;
925 }
926
927 // stay within time limits
928 if(Rdata.time < C.beginTime)
929 {
930 LOG(DEBUG) << " RINEX data timetag " << printTime(C.beginTime,C.longfmt)
931 << " is before begin time.";
932 continue;
933 }
934 if(Rdata.time > C.endTime)
935 {
936 LOG(DEBUG) << " RINEX data timetag " << printTime(C.endTime,C.longfmt)
937 << " is after end time.";
938 break;
939 }
940
941 // fix time systems - only for data, not aux headers
942 if(Rdata.epochFlag == 0) {
943 if(nepochs == 0 &&
944 Rdata.time.getTimeSystem() != Rhead.lastObs.getTimeSystem())
945 {
946 Rhead.lastObs.setTimeSystem(Rdata.time.getTimeSystem());
947 Rhead.firstObs.setTimeSystem(Rdata.time.getTimeSystem());
948 }
949 lastObsTime = Rdata.time;
950 lastObsTime.setTimeSystem(Rhead.lastObs.getTimeSystem());
951 firstObsTime.setTimeSystem(Rhead.lastObs.getTimeSystem());
952 prevObsTime.setTimeSystem(Rhead.lastObs.getTimeSystem());
953 if(firstObsTime == CommonTime::BEGINNING_OF_TIME)
954 firstObsTime = lastObsTime;
955 }
956
957 LOG(DEBUG) << " Read RINEX data: flag " << Rdata.epochFlag
958 << ", timetag " << printTime(Rdata.time,C.longfmt);
959
960 // if aux header data, either output or skip
961 if(Rdata.epochFlag > 1)
962 {
963 //if(C.debug > -1)
964 // for(j=0; j<Rdata.auxHeader.commentList.size(); j++)
965 // LOG(DEBUG) << "Comment: " << Rdata.auxHeader.commentList[j];
966 if(C.verbose) {
967 LOG(INFO) << "\nDump auxiliary header information:";
968 Rdata.auxHeader.dump(LOGstrm);
969 }
970 nauxheads++;
971 continue;
972 }
973
974 // debug: dump the RINEX data object
975 if(C.debug > -1)
976 Rdata.dump(LOGstrm,Rhead);
977
978 // count this epoch
979 nepochs++;
980
981 // check for data out of time order
982 // use < 1.e-3 not < 0 b/c inline header info (epochFlag > 1) excluded
983 if(prevObsTime != CommonTime::BEGINNING_OF_TIME
984 && Rdata.time-prevObsTime < 1.e-3)
985 {
986 // save it
987 if(!cacheon)
988 {
989 // new block
990 cachetime.push_back(prevObsTime);
991 cacheon = true;
992 vector<Rinex3ObsData> v;
993 cache.push_back(v);
994 }
995 cache[cache.size()-1].push_back(Rdata);
996 continue;
997 }
998 cacheon = false;
999
1000 // look for gaps in the timetags
1001 int ncount;
1002 if(C.dt > 0.0)
1003 {
1004 ncount = int(0.5+(lastObsTime-firstObsTime)/C.dt);
1005 // update gap count
1006 if(C.gapcount.size() == 0)
1007 {
1008 // create the list
1009 C.gapcount.push_back(ncount); // start time
1010 C.gapcount.push_back(ncount-1); // end time
1011 }
1012 i = C.gapcount.size() - 1;
1013 if(ncount == C.gapcount[i] + 1) // no gap
1014 C.gapcount[i] = ncount;
1015 else
1016 {
1017 // found a gap
1018 C.gapcount.push_back(ncount); // start time
1019 C.gapcount.push_back(ncount); // end time
1020 }
1021
1022 // TD test after 50 epochs - wrong dt is disasterous
1023 }
1024
1025 // loop over satellites -------------------------------------
1026 Rinex3ObsData::DataMap::const_iterator it;
1027 for(it=Rdata.obs.begin(); it != Rdata.obs.end(); ++it)
1028 {
1029 const RinexSatID& sat(it->first);
1030
1031 // is sat included?
1032 if(C.onlySats.size() > 0 &&
1033 find(C.onlySats.begin(), C.onlySats.end(), sat) == C.onlySats.end()
1034 && find(C.onlySats.begin(), C.onlySats.end(),
1035 RinexSatID(-1,sat.system)) == C.onlySats.end())
1036 continue;
1037
1038 // is sat excluded?
1039 if(find(C.exSats.begin(), C.exSats.end(), sat) != C.exSats.end())
1040 continue;
1041 // check for all sats of this system
1042 else if(find(C.exSats.begin(), C.exSats.end(),
1043 RinexSatID(-1,sat.system)) != C.exSats.end())
1044 continue;
1045
1046 const vector<RinexDatum>& vecData(it->second);
1047
1048 // find this sat in the table; add it if necessary
1049 vector<TableData>::iterator ptab;
1050 ptab = find(table.begin(),table.end(),TableData(sat,nmaxobs));
1051 if(ptab == table.end())
1052 {
1053 // add it
1054 table.push_back(TableData(sat,nmaxobs));
1055 ptab = find(table.begin(),table.end(),TableData(sat,nmaxobs));
1056 ptab->begin = lastObsTime;
1057 if(C.dt > 0.0)
1058 {
1059 ptab->gapcount.push_back(ncount); // start time
1060 ptab->gapcount.push_back(ncount-1); // end time
1061 }
1062 }
1063
1064 // update list of gap times
1065 if(C.dt > 0.0)
1066 {
1067 i = ptab->gapcount.size() - 1; // index of curr end time
1068 if(ncount == ptab->gapcount[i] + 1) // no gap
1069 ptab->gapcount[i] = ncount;
1070 else
1071 {
1072 // found a gap
1073 ptab->gapcount.push_back(ncount); // start time
1074 ptab->gapcount.push_back(ncount); // end time
1075 }
1076 }
1077
1078 // set the end time for this satellite to the current epoch
1079 ptab->end = lastObsTime;
1080 if(C.debug > -1)
1081 {
1082 oss.str("");
1083 oss << "Sat " << setw(2) << sat;
1084 }
1085
1086 // first, find the current system...
1087 char sysCode = sat.systemChar();
1088 string sysStr(string(1,sysCode));
1089
1090 // update Obs data totals
1091 for(size_t index=0; index != vecData.size(); index++)
1092 {
1093 if(C.debug > -1)
1094 oss << " (" << index << ")";
1095
1096 // if this observation is not zero, update it's total count
1097 if(vecData[index].data != 0)
1098 {
1099 (ptab->nobs)[index]++; // per obs
1100 if(totals[sysCode].size() == 0)
1101 totals[sysCode] = vector<int>(vecData.size());
1102 totals[sysCode][index]++; // per system
1103 }
1104
1105 // if looking for milliseconds, update handler
1106 if(C.doms && vecData[index].data != 0)
1107 {
1108 tag = sysStr + Rhead.mapObsTypes[sysStr][index].asString();
1109 if(vectorindex(msots,tag) != -1)
1110 {
1111 C.msh.add(lastObsTime, sat, tag, vecData[index].data);
1112 }
1113 }
1114
1115 if(C.debug > -1)
1116 oss << fixed << setprecision(3)
1117 << " " << asString(Rhead.mapObsTypes[sysStr][index])
1118 << " " << setw(13) << vecData[index].data
1119 << " " << vecData[index].lli
1120 << " " << vecData[index].ssi;
1121
1122 } // end loop over observations
1123
1124 if(C.debug > -1)
1125 LOG(DEBUG) << oss.str();
1126
1127 } // end loop over satellites
1128
1129 if(prevObsTime != CommonTime::BEGINNING_OF_TIME)
1130 {
1131 dt = lastObsTime-prevObsTime;
1132 if(dt > 0.0)
1133 {
1134 for(i=0; i<ndtmax; i++)
1135 {
1136 if(ndt[i] <= 0)
1137 {
1138 bestdt[i]=dt;
1139 ndt[i]=1;
1140 break;
1141 }
1142 if(fabs(dt-bestdt[i]) < 0.0001)
1143 {
1144 ndt[i]++;
1145 break;
1146 }
1147 if(i == ndtmax-1)
1148 {
1149 k = 0;
1150 int nleast = ndt[k];
1151 for(j=1; j<ndtmax; j++)
1152 {
1153 if(ndt[j] <= nleast)
1154 {
1155 k = j;
1156 nleast = ndt[j];
1157 }
1158 }
1159 ndt[k] = 1;
1160 bestdt[k] = dt;
1161 }
1162 }
1163 }
1164 else if(dt == 0)
1165 {
1166 LOG(WARNING) << "Warning - repeated time tag at "
1167 << printTime(lastObsTime,C.longfmt);
1168 }
1169 else
1170 {
1171 LOG(WARNING) << "Warning - time tags out of order: "
1172 << printTime(prevObsTime,C.longfmt) << " > "
1173 << printTime(lastObsTime,C.longfmt);
1174 //<< " " << scientific << setprecision(4) << dt;
1175 }
1176 }
1177 prevObsTime = lastObsTime;
1178
1179 } // end while loop over epochs
1180
1181 istrm.close();
1182
1183 // check that we found some data
1184 if(nepochs <= 0)
1185 {
1186 LOG(INFO) << "File " << filename
1187 << " : no data found. Are time limits wrong?";
1188 continue;
1189 }
1190
1191 // Compute interval -------------------------------------------------
1192 for(i=1,j=0; i < ndtmax; i++)
1193 {
1194 if(ndt[i] > ndt[j])
1195 j = i;
1196 dt = bestdt[j];
1197 }
1198
1199 // Summary info -----------------------------------------------------
1200 LOG(INFO) << "Computed interval " << fixed << setw(5) << setprecision(2)
1201 << dt << " seconds.";
1202 LOG(INFO) << "Computed first epoch: " << printTime(firstObsTime,C.longfmt);
1203 LOG(INFO) << "Computed last epoch: " << printTime(lastObsTime,C.longfmt);
1204
1205 // compute time span of dataset in days/hours/minutes/seconds
1206 oss.str("");
1207 oss << "Computed time span: ";
1208 double secs = lastObsTime - firstObsTime;
1209 int remainder = int(secs);
1210 CivilTime delta(firstObsTime);
1211 delta.day = remainder / 86400; remainder %= 86400;
1212 delta.hour = remainder / 3600; remainder %= 3600;
1213 delta.minute = remainder / 60; remainder %= 60;
1214 delta.second = remainder;
1215 if(delta.day > 0)
1216 oss << delta.day << "d ";
1217
1218 LOG(INFO) << oss.str() << delta.hour << "h " << delta.minute << "m "
1219 << delta.second << "s = " << secs << " seconds.";
1220
1221 //LOG(INFO) << "Computed file size: " << filesize << " bytes.";
1222
1223 // Reusing secs, as it is equivalent to the original expression
1224 // i = 1+int(0.5+(lastObsTime-firstObsTime)/dt);
1225 i = 1+int(0.5 + secs / dt);
1226
1227 LOG(INFO) << "There were " << nepochs << " epochs ("
1228 << fixed << setprecision(2) << double(nepochs*100)/i
1229 << "% of " << i << " possible epochs in this timespan) and "
1230 << nauxheads << " inline header blocks.";
1231
1232 // Sort table
1233 if(C.sorttime)
1234 sort(table.begin(),table.end(),TableBegLessThan());
1235 else
1236 sort(table.begin(),table.end(),TableSATLessThan());
1237
1238 // output table
1239 // header
1240 vector<TableData>::iterator tabIt;
1241 if(table.size() > 0)
1242 table.begin()->sat.setfill('0');
1243
1244 if(!C.brief && !C.notab)
1245 {
1246 // non-brief output ------------
1247 LOG(INFO) << "\n Summary of data available in this file: "
1248 << "(Spans are based on times and interval)";
1249 string fmt(C.gpstime ? C.gpsfmt : C.calfmt);
1250 j = 0;
1251 for(sit=Rhead.mapObsTypes.begin(); sit != Rhead.mapObsTypes.end(); ++sit)
1252 {
1253 RinexSatID sat(sit->first);
1254
1255 map<char, vector<int> >::const_iterator totalsIter;
1256 // compute grand total first
1257 totalsIter = totals.find((sit->first)[0]);
1258 const vector<int>& vec = totalsIter->second;
1259 for(i=0,k=0; k<vec.size(); k++) i += vec[k];
1260 if(i == 0)
1261 continue;
1262
1263 // print the table
1264 if(++j > 1)
1265 LOG(INFO) << "";
1266 LOG(INFO) << "System " << sit->first <<" = "<< sat.systemString() << ":";
1267 oss.str("");
1268 oss << " Sat\\OT:";
1269
1270 // print line of RINEX 3 codes
1271 for(k=0; k < (sit->second).size(); k++) {
1272 RinexObsID rot((sit->second)[k]);
1273 oss << setw(k==0?4:7) << rot.asString(outputVersion);
1274 }
1275 LOG(INFO) << oss.str() << " Span Begin time - End time";
1276
1277 // print the table
1278 for(tabIt = table.begin(); tabIt != table.end(); ++tabIt)
1279 {
1280 std::string sysChar;
1281 sysChar += (tabIt->sat).systemChar();
1282 if((sit->first) == sysChar)
1283 {
1284 oss.str("");
1285 oss << " " << tabIt->sat << " ";
1286 size_t obsSize = (Rhead.mapObsTypes.find(sysChar)->second).size();
1287 for(k = 0; k < obsSize; k++)
1288 oss << setw(7) << tabIt->nobs[k];
1289
1290 oss << setw(7) << 1+int(0.5+(tabIt->end-tabIt->begin)/dt);
1291
1292 LOG(INFO) << oss.str() << " " << printTime(tabIt->begin,fmt)
1293 << " - " << printTime(tabIt->end,fmt);
1294 }
1295 }
1296
1297 oss.str("");
1298 oss << "TOTAL";
1299 for(k=0; k<vec.size(); k++) oss << setw(7) << vec[k];
1300 LOG(INFO) << oss.str();
1301 }
1302 LOG(INFO) << "";
1303 }
1304 else
1305 {
1306 // brief output ---------------
1307 // output satellites
1308 oss.str(""); oss << "SATs(" << table.size() << "):";
1309 i = 0;
1310 for(tabIt = table.begin(); tabIt != table.end(); ++tabIt)
1311 {
1312 oss << " " << tabIt->sat;
1313 if((++i % 20) == 0)
1314 {
1315 LOG(INFO) << oss.str();
1316 oss.str(""); i=0;
1317 oss << "SATs ...:";
1318 }
1319 }
1320 LOG(INFO) << oss.str();
1321
1322 // output obs types
1323 sit = Rhead.mapObsTypes.begin();
1324 for( ; sit != Rhead.mapObsTypes.end(); ++sit)
1325 {
1326 string sysCode = (sit->first);
1327 const vector<RinexObsID>& vec(sit->second);
1328
1329 // is this system found in the list of sats?
1330 map<char, vector<int> >::const_iterator totalsIter;
1331 totalsIter = totals.find(sysCode[0]);
1332 const vector<int>& vectot = totalsIter->second;
1333 for(i=0,k=0; k<vectot.size(); k++) i += vectot[k];
1334 if(i == 0)
1335 continue; // no, skip it
1336
1337 oss.str("");
1338 oss << "System " << RinexSatID(sysCode).systemString3()
1339 << " Obs types(" << vec.size() << "): ";
1340
1341 for(i=0; i<vec.size(); i++) {
1342 RinexObsID rot(vec[i]);
1343 oss << " " << rot.asString(outputVersion);
1344 }
1345
1346 // if RINEX ver. 2, then add ver 2 obstypes in parentheses
1347 //map<string, map<string, RinexObsID> > Rinex3ObsHeader::mapSysR2toR3ObsID
1348 //Rhead.mapSysR2toR3ObsID[sys][ot2] = OT3;
1349 if(Rhead.version < 3)
1350 {
1351 oss << " [v2:";
1352 for(i=0; i<vec.size(); i++)
1353 {
1354 map<string,RinexObsID>::iterator it;
1355 for(it = Rhead.mapSysR2toR3ObsID[sysCode].begin();
1356 it != Rhead.mapSysR2toR3ObsID[sysCode].end(); ++it)
1357 {
1358 if(it->second == vec[i])
1359 {
1360 oss << " " << it->first;
1361 break;
1362 }
1363 }
1364 }
1365 oss << "]";
1366 }
1367
1368 LOG(INFO) << oss.str();
1369 }
1370 }
1371
1372 // gaps
1373 if(C.dogaps)
1374 {
1375 // summary of gaps using count
1376 oss.str("");
1377 oss << "Summary of gaps (vs count) in the data in this file, "
1378 << "assuming dt = " << C.dt << " sec.\n";
1379 if(C.dt != dt)
1380 oss << " Warning - computed dt does not match input dt\n";
1381 oss << " First epoch = " << printTime(firstObsTime,C.longfmt)
1382 << " and last epoch = " << printTime(lastObsTime,C.longfmt) << endl;
1383 oss << " Sat beg - end (count,size) ... "
1384 << "[count = # of dt's from first epoch]\n";
1385 // print for timetags = all sats
1386 k = C.gapcount.size()-1; // size() is at least 2
1387 oss << "GAP ALL " << setw(5) << C.gapcount[0]
1388 << " - " << setw(5) << C.gapcount[k];
1389
1390 // NB DO NOT make ii size_t
1391 for(ii=1; ii<=k-2; ii+=2)
1392 oss << " (" << C.gapcount[ii]+1 // begin of gap
1393 << "," << C.gapcount[ii+1]-C.gapcount[ii]-1 << ")"; // size
1394 oss << endl;
1395
1396 // loop over sats
1397 for(tabIt = table.begin(); tabIt != table.end(); ++tabIt)
1398 {
1399 k = tabIt->gapcount.size() - 1;
1400 oss << "GAP " << tabIt->sat << " " << setw(5) << tabIt->gapcount[0]
1401 << " - " << setw(5) << tabIt->gapcount[k];
1402 // NB DO NOT make ii size_t
1403 for(ii=1; ii<=k-2; ii+=2)
1404 oss << " (" << tabIt->gapcount[ii]+1 << "," // begin count of gap
1405 << tabIt->gapcount[ii+1]-tabIt->gapcount[ii]-1 << ")"; // size
1406 oss << endl;
1407 }
1408
1409 tag = oss.str(); stripTrailing(tag,"\n");
1410 LOG(INFO) << tag;
1411
1412 // summary of gaps using sow
1413 oss.str("");
1414 double t(static_cast<GPSWeekSecond>(firstObsTime).sow), d(C.dt);
1415 oss << "\nSummary of gaps (vs SOW) in the data in this file, assuming dt = "
1416 << C.dt << " sec.\n";
1417 if(C.dt != dt)
1418 oss << " Warning - computed dt does not match input dt\n";
1419 oss << " First epoch = " << printTime(firstObsTime,C.longfmt)
1420 << " and last epoch = " << printTime(lastObsTime,C.longfmt) << endl;
1421 oss << " Sat beg - end (sow,number of missing points)\n";
1422
1423 // print for timetags = all sats
1424 k = C.gapcount.size()-1; // size() is at least 2
1425 oss << "GAP ALL " << fixed << setprecision(1) << setw(8) << t+d*C.gapcount[0]
1426 << " - " << setw(8) << t+d*C.gapcount[k];
1427 // NB DO NOT make ii size_t
1428 for(ii=1; ii<=k-2; ii+=2)
1429 oss << " (" << t+d*(C.gapcount[ii]+1) // begin of gap
1430 << "," << C.gapcount[ii+1]-C.gapcount[ii]-1 << ")"; // size
1431 oss << endl;
1432
1433 // loop over sats
1434 for(tabIt = table.begin(); tabIt != table.end(); ++tabIt)
1435 {
1436 k = tabIt->gapcount.size() - 1;
1437 oss << "GAP " << tabIt->sat << " " << fixed << setprecision(1)
1438 << setw(8) << t+d*tabIt->gapcount[0]
1439 << " - " << setw(8) << t+d*tabIt->gapcount[k];
1440 // NB DO NOT make ii size_t
1441 for(ii=1; ii<=k-2; ii+=2)
1442 oss << " (" << t+d*(tabIt->gapcount[ii]+1) << "," // begin sow of gap
1443 << tabIt->gapcount[ii+1]-tabIt->gapcount[ii]-1 << ")"; // size
1444 oss << endl;
1445 }
1446
1447 tag = oss.str(); stripTrailing(tag,"\n");
1448 LOG(INFO) << tag;
1449
1450 // visibility
1451 if(C.vres > 0)
1452 {
1453 // print visibility graphically, resolution C.vres = counts/character
1454 double dn(static_cast<double>(C.vres));
1455 oss.str("");
1456 oss << "\nVisibility - resolution is " << dn << " epochs = " << dn*C.dt
1457 << " seconds.\n";
1458 oss << " First epoch = " << printTime(firstObsTime,C.longfmt)
1459 << " and last epoch = " << printTime(lastObsTime,C.longfmt) << endl;
1460 oss << "VIS ALL ";
1461 bool isOn(false);
1462 for(k=0,i=0; i<C.gapcount.size()-1; i+=2)
1463 {
1464 ii = int(double(C.gapcount[i]/dn));
1465 if(ii-k > 0)
1466 {
1467 oss << string(ii-k,' ');
1468 k = ii;
1469 isOn = false;
1470 }
1471 ii = int(double(C.gapcount[i+1]/dn));
1472 if(ii-k > 0)
1473 {
1474 if(isOn)
1475 {
1476 oss << "x";
1477 ii--;
1478 }
1479 oss << string(ii-k,'X');
1480 k = ii;
1481 isOn = true;
1482 }
1483 }
1484 LOG(INFO) << oss.str();
1485
1486 // timetable of visibility, resolution dn epochs
1487 // to get resolution = 1 epoch, remove isOn, kk and //RES=1
1488 multimap<int,string> vtab;
1489
1490 // loop over sats
1491 //ostringstream ossvt;
1492 for(tabIt = table.begin(); tabIt != table.end(); ++tabIt)
1493 {
1494 oss.str("");
1495 oss << "VIS " << tabIt->sat << " ";
1496
1497 isOn = false;
1498 bool first(true);
1499 int jj,kk(static_cast<int>(tabIt->gapcount[0]/dn)); // + 0.5);
1500 for(k=0,i=0; i<tabIt->gapcount.size()-1; i+=2)
1501 {
1502 // satellite 'off'
1503 j = int(double(tabIt->gapcount[i]/dn));
1504 if(!first)
1505 {
1506 vtab.insert(multimap<int, string>::value_type(
1507 kk, string("-")+asString(tabIt->sat)));
1508 kk = j;
1509 }
1510 first = false;
1511 jj = j-k;
1512 if(jj > 0)
1513 {
1514 isOn = false;
1515 oss << string(jj,' ');
1516 k = j;
1517 }
1518 // satellite 'on'
1519 j = int(double(tabIt->gapcount[i+1]/dn));
1520 vtab.insert(multimap<int, string>::value_type(
1521 kk, string("+")+asString(tabIt->sat)));
1522 kk = j;
1523 jj = j-k;
1524 if(jj > 0)
1525 {
1526 if(!isOn)
1527 {
1528 isOn = true;
1529 }
1530 else
1531 {
1532 oss << "x";
1533 jj--;
1534 }
1535 oss << string(jj,'X');
1536 k = j;
1537 }
1538 }
1539 vtab.insert(multimap<int, string>::value_type(
1540 kk, string("-")+asString(tabIt->sat)));
1541 LOG(INFO) << oss.str();
1542 }
1543
1544 if(C.vistab)
1545 {
1546 LOG(INFO) << "\n Visibility Timetable - resolution is "
1547 << dn << " epochs = " << dn*C.dt << " seconds.\n"
1548 << " First epoch = " << printTime(firstObsTime,C.longfmt)
1549 << " and last epoch = " << printTime(lastObsTime,C.longfmt) << "\n"
1550 << " YYYY/MM/DD HH:MM:SS = week d secs-of-wk Xtot count nX "
1551 << "seconds nsats visible satellites";
1552 j = k = 0;
1553 CommonTime ttag(firstObsTime);
1554 vector<string> sats;
1555 multimap<int,string>::const_iterator vt;
1556 vt = vtab.begin();
1557 while(vt != vtab.end())
1558 {
1559 while(vt != vtab.end() && vt->first == k)
1560 {
1561 string str(vt->second);
1562 if(str[0] == '+')
1563 {
1564 //LOG(INFO) << "Add " << str.substr(1);
1565 sats.push_back(str.substr(1));
1566 }
1567 else
1568 {
1569 vector<string>::iterator fsat;
1570 fsat = find(sats.begin(),sats.end(),str.substr(1));
1571 if(fsat != sats.end())
1572 {
1573 sats.erase(fsat);
1574 }
1575 }
1576 ++vt;
1577 }
1578
1579 ttag += (k-j)*C.dt*dn;
1580
1581 if(vt == vtab.end())
1582 break;
1583
1584 sort(sats.begin(),sats.end());
1585
1586 oss.str("");
1587 oss << "VTAB " << setw(4) << printTime(ttag,C.longfmt)
1588 << " " << setw(4) << k
1589 << " " << setw(5) << k*C.vres
1590 << " " << setw(3) << vt->first - k
1591 << fixed << setprecision(1)
1592 << " " << setw(8) << (vt->first-k)*C.dt*dn
1593 << " " << setw(5) << sats.size();
1594 for(i=0; i<sats.size(); i++) oss << " " << sats[i];
1595 LOG(INFO) << oss.str();
1596
1597 j = k;
1598 k = vt->first;
1599 }
1600 LOG(INFO) << "VTAB " << setw(4) << printTime(ttag,C.longfmt)
1601 << " " << setw(4) << k
1602 << " " << setw(5) << int(0.5+(ttag-firstObsTime)/C.dt)
1603 << " END";
1604 }
1605
1606 } // end if C.vres > 0 (user chose vis output)
1607 }
1608
1609 // output milliseconds
1610 if(C.doms)
1611 {
1612 C.msh.afterAddbeforeFix();
1613
1614 // true b/c no fixing, but false b/c editing commands follow
1615 LOG(INFO) << C.msh.getFindMessage(false);
1616
1617 vector<string> cmds = C.msh.getEditCommands();
1618 for(i=0; i<cmds.size(); i++)
1619 LOG(INFO) << cmds[i] << " # edit cmd for millisecond clock adjust";
1620 LOG(INFO) << "";
1621 }
1622
1623 // Warnings ------------------------------------------------------------
1624 // there were records out of time order
1625 if(cache.size() > 0)
1626 {
1627 for(i=0; i<cache.size(); i++)
1628 LOG(INFO) << " Warning: " << setw(4) << cache[i].size()
1629 << " data records following epoch "
1630 << printTime(cachetime[i],C.calfmt) << " are out of time order,"
1631 << "\n with epochs " << printTime(cache[i][0].time,C.calfmt)
1632 << " to " << printTime(cache[i][cache[i].size()-1].time,C.calfmt)
1633 << endl;
1634 }
1635
1636 if((Rhead.valid & Rinex3ObsHeader::validInterval)
1637 && fabs(dt-Rhead.interval) > 1.e-3)
1638 LOG(INFO) << " Warning - Computed interval is " << setprecision(2)
1639 << dt << " sec, while input header has " << setprecision(2)
1640 << Rhead.interval << " sec.";
1641
1642 if(C.beginTime == CommonTime::BEGINNING_OF_TIME
1643 && fabs(firstObsTime-Rhead.firstObs) > 1.e-8)
1644 LOG(INFO) << " Warning - Computed first time does not agree with header";
1645
1646 if(C.endTime == CommonTime::END_OF_TIME
1647 && (Rhead.valid & Rinex3ObsHeader::validLastTime)
1648 && fabs(lastObsTime-Rhead.lastObs) > 1.e-8)
1649 LOG(INFO) << " Warning - Computed last time does not agree with header";
1650
1651 // look for empty systems
1652 for(sit=Rhead.mapObsTypes.begin(); sit != Rhead.mapObsTypes.end(); ++sit)
1653 {
1654 map<char,vector<int> >::const_iterator totIt(totals.find(sit->first[0]));
1655 const vector<int>& vec(totIt->second);
1656 for(i=0,k=0; k<vec.size(); k++)
1657 i += vec[k];
1658 if(i == 0)
1659 {
1660 RinexSatID sat(sit->first);
1661 if( (find(C.exSats.begin(), C.exSats.end(),
1662 RinexSatID(-1,sat.system)) == C.exSats.end()) // sys not excluded
1663 &&
1664 (C.onlySats.size() > 0 &&
1665 find(C.onlySats.begin(), C.onlySats.end(), // only system
1666 RinexSatID(-1,sat.system)) != C.onlySats.end()) )
1667 LOG(INFO) << " Warning - System " << sit->first << " = "
1668 << sat.systemString() << " should be removed from the header.";
1669 }
1670 }
1671
1672 // look for obs types that are completely empty
1673 // sit declared above map<std::string,vector<RinexObsID> >::const_iterator sit;
1674 for(sit=Rhead.mapObsTypes.begin(); sit != Rhead.mapObsTypes.end(); ++sit)
1675 {
1676 // loop over obs types in header - systems first
1677 RinexSatID sat(sit->first);
1678 map<char, vector<int> >::const_iterator totalsIter;
1679 totalsIter = totals.find((sit->first)[0]);
1680
1681 // this vector is printed after "TOTAL" above
1682 const vector<int>& totvec = totalsIter->second;
1683
1684 // compute grand total first - skip if this system has no data at all
1685 for(i=0,k=0; k<totvec.size(); k++) i += totvec[k];
1686 if(i == 0)
1687 continue;
1688
1689 for(k=0; k<totvec.size(); k++)
1690 {
1691 if(totvec[k] == 0)
1692 {
1693 tag = string();
1694 if(Rhead.version < 3)
1695 {
1696 map<string,RinexObsID>::iterator it;
1697 for(it = Rhead.mapSysR2toR3ObsID[sit->first].begin();
1698 it != Rhead.mapSysR2toR3ObsID[sit->first].end(); ++it)
1699 {
1700 if(it->second == sit->second[k])
1701 {
1702 tag = string(", ") + it->first + string(" in ver.2");
1703 break;
1704 }
1705 }
1706 }
1707 LOG(INFO) << " Warning - Obs type "
1708 << sit->first << asString((sit->second)[k])
1709 << " (" << sat.systemString()
1710 << " " << asString((sit->second)[k]) << tag
1711 << ") should be removed from header";
1712 }
1713 }
1714 }
1715
1716 // failure due to critical error
1717 if(iret < 0)
1718 break;
1719
1720 if(iret == 0)
1721 nfiles++;
1722
1723 } // end loop over files
1724
1725 if(iret < 0)
1726 return iret;
1727
1728 return nfiles;
1729 }
1730 catch(Exception& e)
1731 {
1732 GPSTK_RETHROW(e);
1733 }
1734 } // end ProcessFiles()
1735
1736 //-----------------------------------------------------------------------------
1737 //-----------------------------------------------------------------------------
1738 //-----------------------------------------------------------------------------
1739