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 /**
40 * @file ReadRawData.cpp
41 * Read RINEX observation file data, all epochs and all files; part of program DDBase.
42 */
43
44 //------------------------------------------------------------------------------------
45 // TD ReadRawData.cpp error msg for failure to open output RAW file
46 // TD ReadRawData.cpp is LastEpoch ever used?
47
48 //------------------------------------------------------------------------------------
49 // includes
50 // system
51 #include <fstream>
52 #include "TimeString.hpp"
53 #include "Epoch.hpp"
54 #include "GPSWeekSecond.hpp"
55
56 // GPSTk
57
58 // DDBase
59 #include "DDBase.hpp"
60
61 //------------------------------------------------------------------------------------
62 using namespace std;
63 using namespace gpstk;
64
65 //------------------------------------------------------------------------------------
66 // local data
67 static CommonTime EarliestTime;// earliest timetag among newly-input observation epochs
68 static int ngood; // number of good data points, this epoch
69 static double sow; // GPS seconds of week of current epoch
70 ofstream ofprs; // output file for PRS solution
71 ofstream *pofs=NULL; // pointer to output file stream (&ofprs)
72
73 //------------------------------------------------------------------------------------
74 // prototypes -- others
75 int OutputClockData(void); // DataOutput.cpp
76 int ReadNextObs(ObsFile& of); // ReadObsFiles.cpp
77 int ProcessRawData(ObsFile& obsfile, CommonTime& timetag, ofstream *pofs)
78 ; // ProcessRawData.cpp
79 // prototypes -- this module only
80 int FindEarliestTime(void);
81 void ComputeSolutionEpoch(void);
82
83 //------------------------------------------------------------------------------------
ReadAndProcessRawData(void)84 int ReadAndProcessRawData(void)
85 {
86 try {
87 int iret,ntotal;
88 size_t nfile;
89
90 if(CI.Verbose) oflog << "BEGIN ReadAndProcessRawData()"
91 << " at total time " << fixed << setprecision(3)
92 << double(clock()-totaltime)/double(CLOCKS_PER_SEC) << " seconds."
93 << endl;
94 if(CI.Screen) cout << "Reading raw data and computing PR solution ..." << endl;
95
96 if(!CI.OutputPRSFile.empty()) {
97 ofprs.open(CI.OutputPRSFile.c_str(),ios::out);
98 if(!ofprs.is_open()) {
99 // TD error msg
100 CI.OutputPRSFile = string("");
101 }
102 else {
103 oflog << "Opened file " << CI.OutputPRSFile << " for PRS output." << endl;
104 ofprs << "# " << Title << endl;
105 ofprs << "PRS site ns week sec wk dX(m) dY(m)"
106 << " dZ(m) clk(m) rms(m) slope PRNs..."
107 //<< " (ret) Valid/NotValid"
108 << endl;
109 }
110 pofs = &ofprs;
111 }
112
113 // loop over all epochs in all files
114 do {
115
116 // find earliest time among open, active files, and synchronize reading
117 iret = FindEarliestTime();
118 if(iret == 1) {
119 if(CI.Debug) oflog << "End of data reached in ReadAndProcessRawData."
120 << endl;
121 iret = 0;
122 break;
123 }
124 if(iret == 2) {
125 if(CI.Verbose) oflog << "After end time (quit) : "
126 << printTime(EarliestTime,"%Y/%02m/%02d %2H:%02M:%6.3f=%F/%10.3g") << endl;
127 iret = 0;
128 break;
129 }
130 if(iret == 3) {
131 if(CI.Debug) oflog << "Before begin time : "
132 << printTime(EarliestTime,"%Y/%02m/%02d %2H:%02M:%6.3f=%F/%10.3g") << endl;
133 iret = 0;
134 continue;
135 }
136
137 if(CI.Debug) oflog << "Found " << ngood << " stations with data at epoch "
138 << printTime(EarliestTime,"%Y/%m/%d %H:%M:%6.3f=%F/%10.3g") << endl;
139
140 // round receiver epoch to even multiple of data interval, else even second
141 ComputeSolutionEpoch();
142
143 // preprocess at this epoch
144 for(nfile=0; nfile<ObsFileList.size(); nfile++) {
145
146 // skip files that are 'dead' or out of synch
147 if(!ObsFileList[nfile].valid) continue;
148 if(fabs(ObsFileList[nfile].Robs.time - EarliestTime) >= 0.5) continue;
149
150 // process at the nominal receive time
151 iret = ProcessRawData(ObsFileList[nfile],ObsFileList[nfile].Robs.time,pofs);
152 if(iret) break;
153
154 } // end loop over observation files
155
156 } while(iret == 0); // end loop over all epochs
157
158 if(!CI.OutputPRSFile.empty()) ofprs.close();
159
160 if(iret) return iret;
161
162 if(CI.Screen) cout << "Last epoch is "
163 << printTime(SolutionEpoch,"%Y/%02m/%02d %2H:%02M:%6.3f = %F/%10.3g") << endl;
164 if(CI.Verbose) oflog << "Last epoch is "
165 << printTime(SolutionEpoch,"%Y/%02m/%02d %2H:%02M:%6.3f = %F/%10.3g") << endl;
166
167 // was there any data?
168 for(ntotal=0,nfile=0; nfile<ObsFileList.size(); nfile++) {
169 //if(!ObsFileList[nfile].valid) continue;
170 if(ObsFileList[nfile].nread <= 0)
171 ObsFileList[nfile].valid = false;
172 else
173 ntotal += ObsFileList[nfile].nread;
174 }
175 if(CI.Verbose)
176 oflog << "Total: " << ObsFileList.size() << " files, "
177 << ntotal << " epochs were read." << endl;
178 if(CI.Screen)
179 cout << "Total: " << ObsFileList.size() << " files, "
180 << ntotal << " epochs were read." << endl;
181
182 if(ntotal == 0) {
183 oflog << "No data found. Abort." << endl;
184 if(CI.Screen)
185 cout << "No data found. Abort." << endl;
186 return -3;
187 }
188
189 // average PR solution
190 {
191 bool ok=true;
192 map<string,Station>::const_iterator it;
193 for(it=Stations.begin(); it != Stations.end(); it++) {
194 Position& pos=Stations[it->first].pos;
195
196 if(CI.Verbose)
197 oflog << "For station " << it->first << " read "
198 << it->second.PRSXstats.N() << " good data epochs." << endl;
199
200 if(it->second.PRSXstats.N() <= 0) {
201 oflog << "Warning - No good data found for station " << it->first << endl;
202 ok = false;
203 continue;
204 }
205
206 Position PRsol;
207 PRsol.setECEF(it->second.PRSXstats.Average(),
208 it->second.PRSYstats.Average(),
209 it->second.PRSZstats.Average());
210 if(CI.Verbose) {
211 oflog << "Average PR solution for site " << it->first
212 << fixed << setprecision(5)
213 << " " << setw(15) << it->second.PRSXstats.Average()
214 << " " << setw(15) << it->second.PRSYstats.Average()
215 << " " << setw(15) << it->second.PRSZstats.Average()
216 << endl;
217 oflog << "Std-dev PR solution for site " << it->first
218 << fixed << setprecision(5)
219 << " " << setw(15) << it->second.PRSXstats.StdDev()
220 << " " << setw(15) << it->second.PRSYstats.StdDev()
221 << " " << setw(15) << it->second.PRSZstats.StdDev()
222 << endl;
223 }
224 if(CI.Screen) {
225 cout << "Average PR solution for site " << it->first
226 << fixed << setprecision(5)
227 << " " << setw(15) << it->second.PRSXstats.Average()
228 << " " << setw(15) << it->second.PRSYstats.Average()
229 << " " << setw(15) << it->second.PRSZstats.Average()
230 << endl;
231 cout << "Std-dev PR solution for site " << it->first
232 << fixed << setprecision(5)
233 << " " << setw(15) << it->second.PRSXstats.StdDev()
234 << " " << setw(15) << it->second.PRSYstats.StdDev()
235 << " " << setw(15) << it->second.PRSZstats.StdDev()
236 << endl;
237 }
238
239 // use PR solution if apriori position not given
240 if(it->second.usePRS) {
241 pos = PRsol;
242 oflog << "Adopting average pseudorange solution for "
243 << it->first << " position"
244 //<< pos.printf(" : %.4x %.4y %.4z meters Position")
245 << endl;
246 if(CI.Screen)
247 cout << "Adopting average pseudorange solution for "
248 << it->first << " position"
249 //<< pos.printf(" : %.4x %.4y %.4z meters Position")
250 << endl;
251 }
252 //else if(pos.getRadius() < 1.) {
253 //}
254 else {
255 // sanity check...
256 // keep this low! large position errors have enduring effects in editing!
257 if(range(pos,PRsol) > 50.0) {
258 oflog << "Warning - Pseudorange solution is far from input "
259 << "position for station " << it->first << " : delta = "
260 << fixed << setprecision(3) << range(pos,PRsol)
261 << " meters. Abort." << endl;
262 cerr << "Warning - Pseudorange solution is far from input "
263 << "position for station " << it->first << " : delta = "
264 << fixed << setprecision(3) << range(pos,PRsol)
265 << " meters. Abort." << endl;
266 iret = -1;
267 OutputClockData(); // usually done in ClockModel() later...
268 }
269 }
270 }
271
272 if(!ok) {
273 oflog << "One or more stations have no data. Abort." << endl;
274 cerr << "One or more stations have no data. Abort." << endl;
275 iret = -3;
276 }
277 }
278
279 return iret;
280 }
281 catch(Exception& e) { GPSTK_RETHROW(e); }
282 catch(std::exception& e) { Exception E("std except: "+string(e.what())); GPSTK_THROW(E); }
283 catch(...) { Exception e("Unknown exception"); GPSTK_THROW(e); }
284 } // end ReadAndProcessRawData()
285
286 //------------------------------------------------------------------------------------
287 // read the data for the next (earliest in future) observation epoch
FindEarliestTime(void)288 int FindEarliestTime(void)
289 {
290 try {
291 int iret;
292 size_t nfile;
293
294 EarliestTime = CommonTime::END_OF_TIME;
295
296 // loop over all (open) obs files
297 for(nfile=0; nfile<ObsFileList.size(); nfile++) {
298
299 // is this a valid, active file?
300 if(!ObsFileList[nfile].valid) continue;
301
302 iret = ReadNextObs(ObsFileList[nfile]);
303 if(iret < 0) { // error or EOF -- set file 'dead'
304 ObsFileList[nfile].valid = false;
305 continue;
306 }
307 // success - file is active
308 else if(ObsFileList[nfile].Robs.time < EarliestTime)
309 EarliestTime = ObsFileList[nfile].Robs.time;
310
311 } // end loop over all obs files
312
313 // if no more data is available, EarliestTime will never get set
314 if(EarliestTime == CommonTime::END_OF_TIME) return 1;
315
316 // if past end time, quit
317 if(EarliestTime > CI.EndTime) return 2;
318
319 // synchronize reading at EarliestTime
320 for(ngood=0,nfile=0; nfile<ObsFileList.size(); nfile++) {
321 if(!ObsFileList[nfile].valid) continue;
322 // if this data time == EarliestTime, process and set flag to read again
323 if(fabs(ObsFileList[nfile].Robs.time - EarliestTime) < 1.) {
324 ngood++;
325 ObsFileList[nfile].getNext = true;
326 }
327 else ObsFileList[nfile].getNext = false;
328 }
329
330 // apply time limits
331 if(EarliestTime < CI.BegTime) return 3;
332
333 return 0;
334 }
335 catch(Exception& e) { GPSTK_RETHROW(e); }
336 catch(std::exception& e) { Exception E("std except: "+string(e.what())); GPSTK_THROW(E); }
337 catch(...) { Exception e("Unknown exception"); GPSTK_THROW(e); }
338 }
339
340 //------------------------------------------------------------------------------------
ComputeSolutionEpoch(void)341 void ComputeSolutionEpoch(void)
342 {
343 try {
344 double dt;
345
346 // round receiver epoch to even multiple of data interval, else even second
347 SolutionEpoch = EarliestTime;
348 sow = static_cast<GPSWeekSecond>(SolutionEpoch).sow;
349 sow = CI.DataInterval * double(int(sow/CI.DataInterval + 0.5));
350 SolutionEpoch += (sow - static_cast<GPSWeekSecond>(SolutionEpoch).sow);
351 if(CI.Debug) oflog << "Solution epoch is "
352 << printTime(SolutionEpoch,"%Y/%02m/%02d %2H:%02M:%6.3f = %F/%10.3g") << endl;
353
354 // save first and last epochs
355 if(fabs(FirstEpoch-CommonTime::BEGINNING_OF_TIME) < 0.1) {
356 FirstEpoch = SolutionEpoch;
357 if(CI.Screen)
358 cout << "First epoch is "
359 << printTime(FirstEpoch,"%Y/%02m/%02d %2H:%02M:%6.3f = %F/%10.3g") << endl;
360 if(CI.Verbose)
361 oflog << "First epoch is "
362 << printTime(FirstEpoch,"%Y/%02m/%02d %2H:%02M:%6.3f = %F/%10.3g") << endl;
363
364 // compute rotation matrix that corrects for earth orientation
365 //PNS = ident<double>(3);
366 } // end if first epoch
367
368 LastEpoch = SolutionEpoch; // TD use LastEpoch?
369
370 // compute the current count
371 dt = SolutionEpoch-FirstEpoch;
372 Count = int(dt/CI.DataInterval + 0.5);
373
374 }
375 catch(Exception& e) { GPSTK_RETHROW(e); }
376 catch(std::exception& e) { Exception E("std except: "+string(e.what())); GPSTK_THROW(E); }
377 catch(...) { Exception e("Unknown exception"); GPSTK_THROW(e); }
378 }
379
380 //------------------------------------------------------------------------------------
381 //------------------------------------------------------------------------------------
382