source: trunk/FACT++/src/fitsdump.cc@ 12804

Last change on this file since 12804 was 12766, checked in by lyard, 14 years ago
fixed end of line issue with fitsdump
File size: 43.5 KB
Line 
1//****************************************************************
2/** @class FitsDumper
3
4 @brief Dumps contents of fits tables to stdout or a file
5
6 */
7 //****************************************************************
8#include "Configuration.h"
9
10#include <map>
11#include <fstream>
12
13#include "externals/fits.h"
14
15//#define PLOTTING_PLEASE
16
17#ifdef PLOTTING_PLEASE
18#include <QPen>
19#include <QtGui>
20#include <QApplication>
21
22#include <qwt_plot.h>
23#include <qwt_plot_grid.h>
24#include <qwt_plot_curve.h>
25#include <qwt_plot_zoomer.h>
26#include <qwt_legend.h>
27#include <qwt_scale_draw.h>
28
29#endif
30#include "Time.h"
31
32using namespace std;
33
34
35#ifdef PLOTTING_PLEASE
36class TimeScaleDraw: public QwtScaleDraw
37{
38public:
39 virtual QwtText label(double v) const
40 {
41 Time t(v);
42 string time = t.GetAsStr("%H:%M:%S%F");
43 while (time[time.size()-1] == '0' && time.size() > 2)
44 {
45 time = time.substr(0, time.size()-1);
46 }
47 return QwtText(time.c_str());
48 }
49};
50#endif
51
52class FitsDumper
53{
54public:
55 FitsDumper();
56 ~FitsDumper();
57
58private:
59 fits* fFile;
60 bool fDotsPlease;
61 bool fNoZeroPlease;
62 string fFilename;
63
64 fits::Table::Columns fColMap;
65 fits::Table::Keys fKeyMap;
66
67 // Convert CCfits::ValueType into a human readable string
68 string ValueTypeToStr(char type) const;
69
70 // Convert CCfits::ValueType into a number of associated bytes
71 int ValueTypeToSize(char type) const;
72
73 /// Calculate the buffer size required to read a row of the fits table, as well as the offsets to each column
74// vector<int> CalculateOffsets() const;
75
76 template<class T>
77 T PtrToValue(const unsigned char* &ptr) const;
78// template<class T>
79// double PtrToDouble(const unsigned char *ptr) const;
80// double PtrToDouble(const unsigned char *ptr, CCfits::ValueType type) const;
81
82 /// Write a single row of the selected data
83 // int WriteRow(ostream &, const vector<MyColumn*> &, const vector<int> &, unsigned char *, const vector<pair<int, int> >&) const;
84
85 bool OpenFile(const string &); /// Open a file
86 bool OpenTable(const string &); /// Open a table
87
88 /// Lists all columns of an open file
89 void List();
90 void ListHeader(const string& filename);
91 void ListKeywords(ostream &);
92
93 bool separateColumnsFromRanges(const vector<string>& list,
94 vector<pair<int, int> >& ranges,
95 vector<string>& listNamesOnly);
96 /// Perform the dumping, based on the current dump list
97 bool Dump(const string &, const vector<string> &list, int);
98 ///Display the selected columns values VS time
99#ifdef PLOTTING_PLEASE
100 int doCurvesDisplay( const vector<string> &list, const string& tableName);
101#endif
102 int doMinMaxPlease(const string& filename, const vector<string>& list, int precision);
103 int doStatsPlease(const string &filename, const vector<string>& list, int precision);
104// void doTBoundary(conf.Get<string>("outfile"), conf.Get<int>("precision"), true);
105 // bool Plot(const vector<string> &list);
106
107public:
108 ///Configures the fitsLoader from the config file and/or command arguments.
109 int ExecConfig(Configuration& conf);
110};
111
112// --------------------------------------------------------------------------
113//
114//! Constructor
115//! @param out
116//! the ostream where to redirect the outputs
117//
118FitsDumper::FitsDumper() : fFile(0), fDotsPlease(false), fNoZeroPlease(false)
119{
120}
121
122// --------------------------------------------------------------------------
123//
124//! Destructor
125//
126FitsDumper::~FitsDumper()
127{
128 if (fFile)
129 delete fFile;
130}
131
132
133string FitsDumper::ValueTypeToStr(char type) const
134{
135 switch (type)
136 {
137 case 'L': return "bool(8)";
138 case 'B': return "byte(8)";
139 case 'I': return "short(16)";
140 case 'J': return "int(32)";
141 case 'K': return "int(64)";
142 case 'E': return "float(32)";
143 case 'D': return "double(64)";
144 default:
145 return "unknwown";
146 }
147}
148
149int FitsDumper::ValueTypeToSize(char type) const
150{
151 switch (type)
152 {
153 case 'L': return sizeof(uint8_t);
154 case 'B': return sizeof(int8_t);
155 case 'I': return sizeof(int16_t);
156 case 'J': return sizeof(int32_t);
157 case 'K': return sizeof(int64_t);
158 case 'E': return sizeof(float);
159 case 'D': return sizeof(double);
160 default:
161 return 0;
162 }
163}
164
165template<class T>
166T FitsDumper::PtrToValue(const unsigned char* &ptr) const
167{
168 T t;
169 reverse_copy(ptr, ptr+sizeof(T), reinterpret_cast<unsigned char*>(&t));
170 ptr += sizeof(T);
171
172 return t;
173}
174// --------------------------------------------------------------------------
175//
176//! Loads the fits file based on the current parameters
177//
178bool FitsDumper::OpenFile(const string &filename)
179{
180 if (fFile)
181 {
182 fFile->close();
183 delete fFile;
184 }
185
186 try {
187 fFile = new fits(filename);
188 }
189 catch (std::runtime_error e)
190 {
191 cout << "Something went wrong while trying to open " << filename;
192 cout << ": " << e.what() << " Aborting dump." << endl;
193 return false;
194 }
195 fFilename = filename;
196
197 const fits::Table::Columns& tCols = fFile->GetColumns();
198
199 for (auto it=tCols.begin(); it != tCols.end(); it++)
200 fColMap.insert(*it);
201
202 const fits::Table::Keys& tkeys = fFile->GetKeys();
203
204 for (auto it=tkeys.begin(); it != tkeys.end(); it++)
205 fKeyMap.insert(*it);
206
207 return true;
208}
209
210bool FitsDumper::OpenTable(const string &)
211{
212 if (!fFile)
213 {
214 cerr << "No file open." << endl;
215 return false;
216 }
217
218
219 return true;
220}
221
222
223void FitsDumper::List()
224{
225 if (!fFile)
226 {
227 cerr << "No file open." << endl;
228 return;
229 }
230
231 cout << "\nFile: " << fFilename << "\n";
232
233 cout << " " << fKeyMap.find("EXTNAME")->second.value << " [";
234 cout << fKeyMap.find("NAXIS2")->second.value << "]\n";
235
236 for (auto it = fColMap.begin(); it != fColMap.end(); it++)
237 {
238 cout << " " << it->first << "[" << it->second.num << "] (" << it->second.unit << ":" << ValueTypeToStr(it->second.type) << ") ";
239 for (auto jt = fKeyMap.begin(); jt != fKeyMap.end(); jt++)
240 if (jt->second.value == it->first)
241 cout << jt->second.comment << endl;
242 }
243
244 cout << endl;
245 cout << flush;
246}
247
248void FitsDumper::ListKeywords(ostream &out)
249{
250 for (auto it=fKeyMap.begin(); it != fKeyMap.end(); it++) {
251 out << "## " << setw(8) << it->first << " = " << setw(10);
252 out << "'" << it->second.value << "'" << " / " << it->second.comment << endl;
253 }
254}
255
256void FitsDumper::ListHeader(const string& filename)
257{
258 if (!fFile)
259 {
260 cerr << "No table open." << endl;
261 return;
262 }
263 ofstream out(filename=="-"?"/dev/stdout":filename);
264 if (!out)
265 {
266 cerr << "Cannot open file " << filename << ": " << strerror(errno) << endl;
267 return;
268 }
269
270 out << "\nTable: " << fKeyMap.find("EXTNAME")->second.value << " (rows=" << fKeyMap.find("NAXIS2")->second.value << ")\n";
271 if (fKeyMap.find("COMMENT") != fKeyMap.end())
272 out << "Comment: \t" << fKeyMap.find("COMMENT")->second.value << "\n";
273
274 ListKeywords(out);
275 out << endl;
276
277}
278
279bool FitsDumper::separateColumnsFromRanges(const vector<string>& list,
280 vector<pair<int, int> >& ranges,
281 vector<string>& listNamesOnly)
282{
283 for (vector<string>::const_iterator it=list.begin(); it!=list.end(); it++)
284 {
285 string columnNameOnly = *it;
286 unsigned long bracketIndex0 = columnNameOnly.find_first_of('[');
287 unsigned long bracketIndex1 = columnNameOnly.find_first_of(']');
288 unsigned long colonIndex = columnNameOnly.find_first_of(':');
289 int columnStart = -1;
290 int columnEnd = -1;
291 if (bracketIndex0 != string::npos)
292 {//there is a range given. Extract the range
293 if (colonIndex != string::npos)
294 {//we have a range here
295 columnStart = atoi(columnNameOnly.substr(bracketIndex0+1, colonIndex-(bracketIndex0+1)).c_str());
296 columnEnd = atoi(columnNameOnly.substr(colonIndex+1, bracketIndex1-(colonIndex+1)).c_str());
297 columnEnd++;
298 }
299 else
300 {//only a single index there
301 columnStart = atoi(columnNameOnly.substr(bracketIndex0+1, bracketIndex1 - (bracketIndex0+1)).c_str());
302 columnEnd = columnStart+1;
303 }
304 columnNameOnly = columnNameOnly.substr(0, bracketIndex0);
305 }
306
307 if (fColMap.find(columnNameOnly) == fColMap.end())
308 {
309 cerr << "ERROR - Column '" << columnNameOnly << "' not found in table." << endl;
310 return false;
311 }
312 fits::Table::Column& cCol = fColMap.find(columnNameOnly)->second;
313 if (bracketIndex0 == string::npos)
314 {//no range given: use the full range
315 ranges.push_back(make_pair(0, cCol.num));
316 columnStart = 0;
317 columnEnd = 1;
318 }
319 else
320 {//use the range extracted earlier
321 if (columnStart < 0)
322 {
323 cerr << "ERROR - Start range for column " << columnNameOnly << " is less than zero (" << columnStart << "). Aborting" << endl;
324 return false;
325 }
326 if (columnEnd>1 && columnEnd > (int)(cCol.num))
327 {
328 cerr << "ERROR - End range for column " << columnNameOnly << " is greater than the last element (" << cCol.num << " vs " << columnEnd << "). Aborting" << endl;
329 return false;
330 }
331 ranges.push_back(make_pair(columnStart, columnEnd));
332 }
333// cout << "Will be exporting from " << columnStart << " to " << columnEnd-1 << " for column " << columnNameOnly << endl;
334 listNamesOnly.push_back(columnNameOnly);
335 }
336 return true;
337}
338// --------------------------------------------------------------------------
339//
340//! Perform the actual dump, based on the current parameters
341//
342bool FitsDumper::Dump(const string &filename, const vector<string> &list, int precision)
343{
344
345 //first of all, let's separate the columns from their ranges and check that the requested columns are indeed part of the file
346 vector<pair<int, int> > ranges;
347 vector<string> listNamesOnly;
348
349 if (!separateColumnsFromRanges(list, ranges, listNamesOnly))
350 {
351 cerr << "Something went wrong while extracting the columns names from parameters. Aborting" << endl;
352 return false;
353 }
354
355 ofstream out(filename=="-"?"/dev/stdout":filename);
356 if (!out)
357 {
358 cerr << "Cannot open file " << filename << ": " << strerror(errno) << endl;
359 return false;
360 }
361
362 out.precision(precision);
363
364 out << "## --------------------------------------------------------------------------\n";
365 if (filename!="-")
366 out << "## File: \t" << filename << '\n';
367 out << "## Table: \t" << fKeyMap.find("EXTNAME")->second.value << '\n';
368 out << "## Comment: \t" << ((fKeyMap.find("COMMENT") != fKeyMap.end()) ? fKeyMap.find("COMMENT")->second.value : "") << '\n';
369
370 out << "## NumRows: \t" << fFile->GetInt("NAXIS2") << '\n';
371 out << "## --------------------------------------------------------------------------\n";
372 ListKeywords(out);
373 out << "## --------------------------------------------------------------------------\n";
374 out << "#\n";
375
376 auto rangesIt = ranges.begin();
377 for (vector<string>::const_iterator it=listNamesOnly.begin(); it!=listNamesOnly.end(); it++, rangesIt++)
378 {
379 const fits::Table::Column& col = fColMap[*it];
380// const MyColumn *col = static_cast<MyColumn*>(fTable->column()[*it]);
381 if (rangesIt->first != 0 || rangesIt->second != (int)(col.num))
382 {
383 out << "#";
384 for (int i=rangesIt->first; i<rangesIt->second; i++)
385 out << " " << *it << "[" << i << "]";
386 out << ": " << col.unit;
387 }
388 else
389 out << "# " << *it << "[" << col.num << "]: " << col.unit;
390
391// FIXME: retrive the column comment
392// if (!col->comment().empty())
393// out << " (" <<col->comment() << ")";
394 out << '\n';
395 }
396 out << "#" << endl;
397
398
399 // Loop over all columns in our list of requested columns
400 vector<pair<char, char*> > columnsData;
401
402 for (vector<string>::const_iterator it=listNamesOnly.begin(); it!=listNamesOnly.end(); it++)
403 {
404 fits::Table::Column& cCol = fColMap.find(*it)->second;
405 columnsData.push_back(make_pair(cCol.type, new char[cCol.num*cCol.size]));
406 fFile->SetPtrAddress(*it, columnsData[columnsData.size()-1].second);
407 }
408 int numRows = fFile->GetInt("NAXIS2");
409 int row = 0;
410 while (fFile->GetNextRow() && row < numRows)
411 {
412 row++;
413 rangesIt = ranges.begin();
414 for (auto it=columnsData.begin(); it != columnsData.end(); it++, rangesIt++)
415 {
416 for (int i=rangesIt->first; i<rangesIt->second; i++)
417 {
418 switch (it->first) {
419 case 'L':
420 out << reinterpret_cast<bool*>(it->second)[i] << " ";
421 break;
422 case 'B':
423 out << reinterpret_cast<char*>(it->second)[i] << " ";
424 break;
425 case 'I':
426 out << reinterpret_cast<int16_t*>(it->second)[i] << " ";
427 break;
428 case 'J':
429 out << reinterpret_cast<int32_t*>(it->second)[i] << " ";
430 break;
431 case 'K':
432 out << reinterpret_cast<int64_t*>(it->second)[i] << " ";
433 break;
434 case 'E':
435 out << reinterpret_cast<float*>(it->second)[i] << " ";
436 break;
437 case 'D':
438 out << reinterpret_cast<double*>(it->second)[i] << " ";
439 break;
440 default:
441 ;
442 }
443 }
444 }
445 out << endl;
446 }
447 for (auto it = columnsData.begin(); it != columnsData.end(); it++)
448 delete[] it->second;
449 return true;
450}
451
452// --------------------------------------------------------------------------
453//
454//! Retrieves the configuration parameters
455//! @param conf
456//! the configuration object
457//
458int FitsDumper::ExecConfig(Configuration& conf)
459{
460 if (conf.Has("fitsfile"))
461 {
462 if (!OpenFile(conf.Get<string>("fitsfile")))
463 return -1;
464 }
465#ifdef PLOTTING_PLEASE
466 if (conf.Get<bool>("stat") && conf.Get<bool>("graph"))
467 {
468 cout << "Invalid conbination of options: cannot graph stats. Aborting" << endl;
469 return -1;
470 }
471#endif
472 if (conf.Get<bool>("minmax") && conf.Get<bool>("stat"))
473 {
474 cout << "Invalid combination of options: cannot do stats and minmax. Aborting" << endl;
475 return -1;
476 }
477#ifdef PLOTTING_PLEASE
478 if (conf.Get<bool>("minmax") && conf.Get<bool>("graph"))
479 {
480 cout << "Invalid combination of options: cannot graph minmax. Aborting" << endl;
481 return -1;
482 }
483#endif
484 if (conf.Get<bool>("stat") && conf.Get<bool>("nozero"))
485 {
486 cout << "Invalid combination of options: nozero only works with minmax. Aborting" << endl;
487 return -1;
488 }
489
490 if (conf.Get<bool>("nozero"))
491 {
492 fNoZeroPlease = true;
493 }
494
495 if (conf.Has("tablename"))
496 {
497 if (!OpenTable(conf.Get<string>("tablename")))
498 return -1;
499 }
500
501 if (conf.Get<bool>("list"))
502 List();
503 /*
504 if (conf.Get<bool>("tstart"))
505 {
506 doTBoundary(conf.Get<string>("outfile"), conf.Get<int>("precision"), true);
507 }
508 if (conf.Get<bool>("tstop"))
509 {
510 doTBoundary(conf.Get<string>("outfile"), conf.Get<int>("precision"), false);
511 }
512 */
513 if (conf.Get<bool>("minmax"))
514 {
515 if (!conf.Has("col"))
516 {
517 cout << "Please specify the columns that should be dumped as arguments. Aborting" << endl;
518 return 0;
519 }
520 doMinMaxPlease(conf.Get<string>("outfile"), conf.Get<vector<string>>("col"), conf.Get<int>("precision"));
521 return 0;
522 }
523
524 if (conf.Get<bool>("stat"))
525 {
526 if (!conf.Has("col"))
527 {
528 cout << "Please specify the columns that should be dumped as arguments. Aborting" << endl;
529 return 0;
530 }
531 doStatsPlease(conf.Get<string>("outfile"), conf.Get<vector<string>>("col"), conf.Get<int>("precision"));
532 return 0;
533 }
534
535#ifdef PLOTTING_PLEASE
536 if (conf.Get<bool>("graph"))
537 {
538 if (!conf.Has("col"))
539 {
540 cout << "Please specify the columns that should be dumped as arguments. Aborting" << endl;
541 return 0;
542 }
543 if (conf.Get<bool>("dots"))
544 fDotsPlease = true;
545 doCurvesDisplay(conf.Get<vector<string>>("col"),
546 conf.Get<string>("tablename"));
547 return 1;
548 }
549#endif
550
551 if (conf.Get<bool>("header"))
552 ListHeader(conf.Get<string>("outfile"));
553
554 if (conf.Get<bool>("header") || conf.Get<bool>("list"))
555 return 1;
556
557 if (conf.Has("outfile"))
558 {
559 if (!conf.Has("col"))
560 {
561 cout << "Please specify the columns that should be dumped as arguments. Aborting" << endl;
562 return 0;
563 }
564 if (!Dump(conf.Get<string>("outfile"),
565 conf.Get<vector<string>>("col"),
566 conf.Get<int>("precision")))
567 return -1;
568 }
569
570
571 return 0;
572}
573
574void PrintUsage()
575{
576 cout <<
577 "fitsdump is a tool to dump data from a FITS table as ascii.\n"
578 "\n"
579 "Usage: fitsdump [OPTIONS] fitsfile col col ... \n"
580 " or: fitsdump [OPTIONS]\n";
581 cout << endl;
582}
583
584void PrintHelp()
585{
586 //
587}
588
589struct minMaxStruct {
590 double min;
591 double max;
592 double average;
593 long numValues;
594 minMaxStruct() : min(1e10), max(-1e10), average(0), numValues(0) {}
595};
596
597int FitsDumper::doMinMaxPlease(const string& filename, const vector<string>& list, int precision)
598{
599 //first of all, let's separate the columns from their ranges and check that the requested columns are indeed part of the file
600 vector<pair<int, int> > ranges;
601 vector<string> listNamesOnly;
602
603 if (!separateColumnsFromRanges(list, ranges, listNamesOnly))
604 {
605 cerr << "Something went wrong while extracting the columns names from parameters. Aborting" << endl;
606 return false;
607 }
608
609 ofstream out(filename=="-"?"/dev/stdout":filename);
610 if (!out)
611 {
612 cerr << "Cannot open file " << filename << ": " << strerror(errno) << endl;
613 return false;
614 }
615
616 out.precision(precision);
617
618 // Loop over all columns in our list of requested columns
619 vector<pair<char, char*> > columnsData;
620 vector<minMaxStruct> statData;
621 int numRows = fFile->GetInt("NAXIS2");
622 auto rangesIt = ranges.begin();
623 for (vector<string>::const_iterator it=listNamesOnly.begin(); it!=listNamesOnly.end(); it++, rangesIt++)
624 {
625 fits::Table::Column& cCol = fColMap.find(*it)->second;
626 columnsData.push_back(make_pair(cCol.type, new char[cCol.num*cCol.size]));
627// minMaxStuct initData;
628 statData.push_back(minMaxStruct());
629 fFile->SetPtrAddress(*it, columnsData[columnsData.size()-1].second);
630 }
631
632 int row = 0;
633 long UTCvalue0=0;
634 long UTCvalue1=0;
635 while (fFile->GetNextRow() && row < numRows)
636 {
637 rangesIt = ranges.begin();
638 auto statsIt = statData.begin();
639 for (auto it=columnsData.begin(); it != columnsData.end(); it++, rangesIt++, statsIt++)
640 {
641 double cValue = 0;
642 for (int i=rangesIt->first; i<rangesIt->second; i++)
643 {
644 switch (it->first) {
645 case 'L':
646 cValue = reinterpret_cast<bool*>(it->second)[i];
647 break;
648 case 'B':
649 cValue = reinterpret_cast<bool*>(it->second)[i];
650 break;
651 case 'I':
652 cValue = reinterpret_cast<int16_t*>(it->second)[i];
653 break;
654 case 'J':
655 cValue = reinterpret_cast<int32_t*>(it->second)[i];
656 break;
657 case 'K':
658 cValue = reinterpret_cast<int64_t*>(it->second)[i];
659 break;
660 case 'E':
661 cValue = reinterpret_cast<float*>(it->second)[i];
662 break;
663 case 'D':
664 cValue = reinterpret_cast<double*>(it->second)[i];
665 break;
666 default:
667 ;
668 }
669 if (list.size() == 1 && (list[0] == "UnixTimeUTC" || list[0] == "PCTime"))
670 {
671 if (i==0)
672 {
673 UTCvalue0 = cValue;
674 }
675 else
676 {
677 UTCvalue1 = cValue;
678 boost::posix_time::ptime unixTimeT( boost::gregorian::date(1970, boost::gregorian::Jan, 1),
679 boost::posix_time::seconds(UTCvalue0) + boost::posix_time::microsec(UTCvalue1));
680
681 Time mjdTime(unixTimeT);
682 cValue = mjdTime.Mjd();
683 if (!fNoZeroPlease || cValue != 0)
684 {
685 statsIt->average += cValue;
686 if (cValue < statsIt->min)
687 statsIt->min = cValue;
688 if (cValue > statsIt->max)
689 statsIt->max = cValue;
690 statsIt->numValues++;
691 }
692 }
693 }
694 else {
695 if (!fNoZeroPlease || cValue != 0)
696 {
697 statsIt->average += cValue;
698 if (cValue < statsIt->min)
699 statsIt->min = cValue;
700 if (cValue > statsIt->max)
701 statsIt->max = cValue;
702 statsIt->numValues++;
703 }
704 }
705 }
706 }
707 row++;
708 }
709 for (auto it = columnsData.begin(); it != columnsData.end(); it++)
710 delete[] it->second;
711
712 //okay. So now I've got ALL the data, loaded.
713 //let's do the summing and averaging in a safe way (i.e. avoid overflow of variables as much as possible)
714 rangesIt = ranges.begin();
715 auto statsIt = statData.begin();
716
717 auto nameIt = listNamesOnly.begin();
718 for (auto it=columnsData.begin(); it != columnsData.end(); it++, rangesIt++, statsIt++, nameIt++)
719 {
720// int span = rangesIt->second - rangesIt->first;
721 cout << *nameIt << ": " << endl;
722 if (statsIt->numValues != 0)
723 {
724 statsIt->average /= statsIt->numValues;
725 out << "min: " << statsIt->min << endl;
726 out << "max: " << statsIt->max << endl;
727 out << "mea: " << statsIt->average << endl;
728 }
729 else
730 {
731 out << "min: 0" << endl << "max: 0" << endl << "mea: " << endl;
732 }
733
734 }
735 return true;
736}
737
738/*
739void FitsDumper::doTBoundary(const string& filename, int precision, bool tStop)
740{
741
742 //first of all, let's separate the columns from their ranges and check that the requested columns are indeed part of the file
743 vector<pair<int, int> > ranges;
744 vector<string> listNamesOnly;
745
746 if (!separateColumnsFromRanges(list, ranges, listNamesOnly))
747 {
748 cerr << "Something went wrong while extracting the columns names from parameters. Aborting" << endl;
749 return false;
750 }
751
752 ofstream out(filename=="-"?"/dev/stdout":filename);
753 if (!out)
754 {
755 cerr << "Cannot open file " << filename << ": " << strerror(errno) << endl;
756 return false;
757 }
758
759 out.precision(precision);
760
761 // Loop over all columns in our list of requested columns
762 vector<pair<char, char*> > columnsData;
763 vector<minMaxStruct> statData;
764 int numRows = fFile->GetInt("NAXIS2");
765 auto rangesIt = ranges.begin();
766 for (vector<string>::const_iterator it=listNamesOnly.begin(); it!=listNamesOnly.end(); it++, rangesIt++)
767 {
768 fits::Table::Column& cCol = fColMap.find(*it)->second;
769 columnsData.push_back(make_pair(cCol.type, new char[cCol.num*cCol.size]));
770 // minMaxStuct initData;
771 statData.push_back(minMaxStruct());
772 fFile->SetPtrAddress(*it, columnsData[columnsData.size()-1].second);
773 }
774
775 int row = 0;
776 while (fFile->GetNextRow() && row < numRows)
777 {
778 rangesIt = ranges.begin();
779 auto statsIt = statData.begin();
780 for (auto it=columnsData.begin(); it != columnsData.end(); it++, rangesIt++, statsIt++)
781 {
782 double cValue = 0;
783 for (int i=rangesIt->first; i<rangesIt->second; i++)
784 {
785 switch (it->first) {
786 case 'L':
787 cValue = reinterpret_cast<bool*>(it->second)[i];
788 break;
789 case 'B':
790 cValue = reinterpret_cast<bool*>(it->second)[i];
791 break;
792 case 'I':
793 cValue = reinterpret_cast<int16_t*>(it->second)[i];
794 break;
795 case 'J':
796 cValue = reinterpret_cast<int32_t*>(it->second)[i];
797 break;
798 case 'K':
799 cValue = reinterpret_cast<int64_t*>(it->second)[i];
800 break;
801 case 'E':
802 cValue = reinterpret_cast<float*>(it->second)[i];
803 break;
804 case 'D':
805 cValue = reinterpret_cast<double*>(it->second)[i];
806 break;
807 default:
808 ;
809 }
810 if (!fNoZeroPlease || cValue != 0)
811 {
812 statsIt->average += cValue;
813 if (cValue < statsIt->min)
814 statsIt->min = cValue;
815 if (cValue > statsIt->max)
816 statsIt->max = cValue;
817 statsIt->numValues++;
818 }
819 }
820 }
821 row++;
822 }
823 for (auto it = columnsData.begin(); it != columnsData.end(); it++)
824 delete[] it->second;
825
826 //okay. So now I've got ALL the data, loaded.
827 //let's do the summing and averaging in a safe way (i.e. avoid overflow of variables as much as possible)
828 rangesIt = ranges.begin();
829 auto statsIt = statData.begin();
830
831 auto nameIt = listNamesOnly.begin();
832 for (auto it=columnsData.begin(); it != columnsData.end(); it++, rangesIt++, statsIt++, nameIt++)
833 {
834 int span = rangesIt->second - rangesIt->first;
835 cout << *nameIt << ": " << endl;
836 if (statsIt->numValues != 0)
837 {
838 statsIt->average /= statsIt->numValues;
839 out << "min: " << statsIt->min << endl;
840 out << "max: " << statsIt->max << endl;
841 out << "mea: " << statsIt->average << endl;
842 }
843 else
844 {
845 out << "min: 0" << endl << "max: 0" << endl << "mea: " << endl;
846 }
847
848 }
849 return true;
850
851}
852*/
853template<typename T>
854void displayStats(char* array, int numElems, ofstream& out)
855{
856 if (numElems == 0)
857 {
858 out << "Min: 0" << endl;
859 out << "Max: 0" << endl;
860 out << "Med: 0" << endl;
861 out << "Mea: 0" << endl;
862 return;
863 }
864 sort(reinterpret_cast<T*>(array), &reinterpret_cast<T*>(array)[numElems]);
865 out << "Min: " << reinterpret_cast<T*>(array)[0] << endl;
866 out << "Max: " << reinterpret_cast<T*>(array)[numElems-1] << endl;
867 if (numElems%2 == 0)
868 out << "Med: " << (reinterpret_cast<T*>(array)[numElems/2] + reinterpret_cast<int16_t*>(array)[numElems/2+1])/2.f << endl;
869 else
870 out << "Med: " << reinterpret_cast<T*>(array)[numElems/2+1] << endl;
871
872 double average = 0;
873 for (int i=0;i<numElems;i++)
874 average += reinterpret_cast<T*>(array)[i];
875 out << "Mea: " << average/(double)numElems << endl;
876
877}
878int FitsDumper::doStatsPlease(const string &filename, const vector<string>& list, int precision)
879{
880
881 //first of all, let's separate the columns from their ranges and check that the requested columns are indeed part of the file
882 vector<pair<int, int> > ranges;
883 vector<string> listNamesOnly;
884
885 if (!separateColumnsFromRanges(list, ranges, listNamesOnly))
886 {
887 cerr << "Something went wrong while extracting the columns names from parameters. Aborting" << endl;
888 return false;
889 }
890
891 ofstream out(filename=="-"?"/dev/stdout":filename);
892 if (!out)
893 {
894 cerr << "Cannot open file " << filename << ": " << strerror(errno) << endl;
895 return false;
896 }
897
898 out.precision(precision);
899
900 // Loop over all columns in our list of requested columns
901 vector<pair<char, char*> > columnsData;
902 vector<char*> statData;
903 int numRows = fFile->GetInt("NAXIS2");
904 auto rangesIt = ranges.begin();
905 for (vector<string>::const_iterator it=listNamesOnly.begin(); it!=listNamesOnly.end(); it++, rangesIt++)
906 {
907 fits::Table::Column& cCol = fColMap.find(*it)->second;
908 columnsData.push_back(make_pair(cCol.type, new char[cCol.num*cCol.size]));
909 statData.push_back(new char[(rangesIt->second - rangesIt->first)*cCol.size*numRows]);
910 fFile->SetPtrAddress(*it, columnsData[columnsData.size()-1].second);
911 }
912
913 int row = 0;
914 while (fFile->GetNextRow() && row < numRows)
915 {
916 rangesIt = ranges.begin();
917 auto statsIt = statData.begin();
918 for (auto it=columnsData.begin(); it != columnsData.end(); it++, rangesIt++, statsIt++)
919 {
920 int span = rangesIt->second - rangesIt->first;
921 for (int i=rangesIt->first; i<rangesIt->second; i++)
922 {
923 switch (it->first) {
924 case 'L':
925 reinterpret_cast<bool*>(*statsIt)[i - rangesIt->first + row*span] = reinterpret_cast<bool*>(it->second)[i];
926 break;
927 case 'B':
928 reinterpret_cast<char*>(*statsIt)[i - rangesIt->first + row*span] = reinterpret_cast<char*>(it->second)[i];
929 break;
930 case 'I':
931 reinterpret_cast<int16_t*>(*statsIt)[i - rangesIt->first + row*span] = reinterpret_cast<int16_t*>(it->second)[i];
932 break;
933 case 'J':
934 reinterpret_cast<int32_t*>(*statsIt)[i - rangesIt->first + row*span] = reinterpret_cast<int32_t*>(it->second)[i];
935 break;
936 case 'K':
937 reinterpret_cast<int64_t*>(*statsIt)[i - rangesIt->first + row*span] = reinterpret_cast<int64_t*>(it->second)[i];
938 break;
939 case 'E':
940 reinterpret_cast<float*>(*statsIt)[i - rangesIt->first + row*span] = reinterpret_cast<float*>(it->second)[i];
941 break;
942 case 'D':
943 reinterpret_cast<double*>(*statsIt)[i - rangesIt->first + row*span] = reinterpret_cast<double*>(it->second)[i];
944 break;
945 default:
946 ;
947 }
948 }
949 }
950 row++;
951 }
952 for (auto it = columnsData.begin(); it != columnsData.end(); it++)
953 delete[] it->second;
954
955 //okay. So now I've got ALL the data, loaded.
956 //let's do the summing and averaging in a safe way (i.e. avoid overflow of variables as much as possible)
957 rangesIt = ranges.begin();
958 auto statsIt = statData.begin();
959
960 auto nameIt = list.begin();
961 for (auto it=columnsData.begin(); it != columnsData.end(); it++, rangesIt++, statsIt++, nameIt++)
962 {
963 int span = rangesIt->second - rangesIt->first;
964 out << *nameIt << ": " << endl;
965 switch (it->first) {
966 case 'L':
967 displayStats<bool>(*statsIt, numRows*span, out);
968 break;
969 case 'B':
970 displayStats<char>(*statsIt, numRows*span, out);
971 break;
972 case 'I':
973 displayStats<int16_t>(*statsIt, numRows*span, out);
974 break;
975 case 'J':
976 displayStats<int32_t>(*statsIt, numRows*span, out);
977 break;
978 case 'K':
979 displayStats<int64_t>(*statsIt, numRows*span, out);
980 break;
981 case 'E':
982 displayStats<float>(*statsIt, numRows*span, out);
983 break;
984 case 'D':
985 displayStats<double>(*statsIt, numRows*span, out);
986 break;
987 default:
988 ;
989 }
990 }
991 return true;
992}
993#ifdef PLOTTING_PLEASE
994int FitsDumper::doCurvesDisplay( const vector<string> &list, const string& tableName)
995{
996 //first of all, let's separate the columns from their ranges and check that the requested columns are indeed part of the file
997 vector<pair<int, int> > ranges;
998 vector<string> listNamesOnly;
999 if (!separateColumnsFromRanges(list, ranges, listNamesOnly))
1000 {
1001 cerr << "Something went wrong while extracting the columns names from parameters. Aborting" << endl;
1002 return false;
1003 }
1004 vector<string> curvesNames;
1005 stringstream str;
1006 for (auto it=ranges.begin(), jt=listNamesOnly.begin(); it != ranges.end(); it++, jt++)
1007 {
1008 for (int i=it->first; i<it->second;i++)
1009 {
1010 str.str("");
1011 str << *jt << "[" << i << "]";
1012 curvesNames.push_back(str.str());
1013 }
1014 }
1015 char* handle = new char[17];
1016 sprintf(handle,"FitsDump Display");
1017// Qt::HANDLE h = *handle;//NULL
1018 int argc = 1;
1019 char ** argv = &handle;
1020 QApplication a(argc, argv);
1021
1022
1023
1024 QwtPlot* plot = new QwtPlot();
1025 QwtPlotGrid* grid = new QwtPlotGrid;
1026 grid->enableX(false);
1027 grid->enableY(true);
1028 grid->enableXMin(false);
1029 grid->enableYMin(false);
1030 grid->setMajPen(QPen(Qt::black, 0, Qt::DotLine));
1031 grid->attach(plot);
1032 plot->setAutoReplot(true);
1033 string title = tableName;
1034 plot->setAxisScaleDraw( QwtPlot::xBottom, new TimeScaleDraw());
1035
1036 QWidget window;
1037 QHBoxLayout* layout = new QHBoxLayout(&window);
1038 layout->setContentsMargins(0,0,0,0);
1039 layout->addWidget(plot);
1040
1041 QwtPlotZoomer zoom(plot->canvas());
1042 zoom.setRubberBandPen(QPen(Qt::gray, 2, Qt::DotLine));
1043 zoom.setTrackerPen(QPen(Qt::gray));
1044 int totalSize = 0;
1045 for (unsigned int i=0;i<list.size();i++)
1046 totalSize += ranges[i].second - ranges[i].first;
1047
1048 vector<QwtPlotCurve*> curves(totalSize);
1049 int ii=0;
1050 for (auto it = curves.begin(), jt=curvesNames.begin(); it != curves.end(); it++, jt++)
1051 {
1052 *it = new QwtPlotCurve(jt->c_str());
1053 switch (ii%6)
1054 {
1055 case 0:
1056 (*it)->setPen(QColor(255,0,0));
1057 break;
1058 case 1:
1059 (*it)->setPen(QColor(0,255,0));
1060 break;
1061 case 2:
1062 (*it)->setPen(QColor(0,0,255));
1063 break;
1064 case 3:
1065 (*it)->setPen(QColor(255,255,0));
1066 break;
1067 case 4:
1068 (*it)->setPen(QColor(0,255,255));
1069 break;
1070 case 5:
1071 (*it)->setPen(QColor(255,0,255));
1072 break;
1073 default:
1074 (*it)->setPen(QColor(0,0,0));
1075 };
1076 ii++;
1077 if (fDotsPlease)
1078 (*it)->setStyle(QwtPlotCurve::Dots);
1079 else
1080 (*it)->setStyle(QwtPlotCurve::Lines);
1081 (*it)->attach(plot);
1082 }
1083 plot->insertLegend(new QwtLegend(), QwtPlot::RightLegend);
1084
1085
1086 vector<pair<char, char*> > columnsData;
1087
1088 for (vector<string>::const_iterator it=listNamesOnly.begin(); it!=listNamesOnly.end(); it++)
1089 {
1090 fits::Table::Column& cCol = fColMap.find(*it)->second;
1091 columnsData.push_back(make_pair(cCol.type, new char[cCol.num*cCol.size]));
1092 fFile->SetPtrAddress(*it, columnsData[columnsData.size()-1].second);
1093 }
1094
1095 //add the time column to the given columns
1096 if (fColMap.find("Time") == fColMap.end() &&
1097 fColMap.find("UnixTimeUTC") == fColMap.end())
1098 {
1099 cerr << "Error: time column could not be found in given table. Aborting" << endl;
1100 return false;
1101 }
1102 const fits::Table::Column& timeCol = (fColMap.find("Time") != fColMap.end()) ? fColMap.find("Time")->second : fColMap.find("UnixTimeUTC")->second;
1103 bool unixTime = (fColMap.find("Time") == fColMap.end());
1104 if (unixTime)
1105 ranges.push_back(make_pair(0,2));
1106 else
1107 ranges.push_back(make_pair(0,1));
1108 columnsData.push_back(make_pair(timeCol.type, new char[timeCol.num*timeCol.size]));
1109 fFile->SetPtrAddress(unixTime ? "UnixTimeUTC" : "Time", columnsData[columnsData.size()-1].second);
1110
1111// stringstream str;
1112 str.str("");
1113
1114
1115 vector<double*> xValues(totalSize);
1116 double* yValues;
1117 cout.precision(10);
1118 str.precision(20);
1119 for (auto it=xValues.begin(); it!=xValues.end(); it++)
1120 *it = new double[fFile->GetInt("NAXIS2")];
1121
1122 yValues = new double[fFile->GetInt("NAXIS2")];
1123
1124 cout.precision(3);
1125 int endIndex = 0;
1126 int numRows = fFile->GetInt("NAXIS2");
1127 for (int i=1;i<numRows;i++)
1128 {
1129 fFile->GetNextRow();
1130 cout << "\r" << "Constructing graph " << ((float)(endIndex)/(float)(fFile->GetInt("NAXIS2")))*100.0 << "%";
1131 endIndex++;
1132 auto rangesIt = ranges.begin();
1133 for (auto it=columnsData.begin(); it != columnsData.end(); it++, rangesIt++)
1134 {
1135 for (int j=rangesIt->first; j<rangesIt->second; j++)
1136 {
1137 switch (it->first) {
1138 case 'L':
1139 str << reinterpret_cast<bool*>(it->second)[j] << " ";
1140 break;
1141 case 'B':
1142 str << reinterpret_cast<char*>(it->second)[j] << " ";
1143 break;
1144 case 'I':
1145 str << reinterpret_cast<int16_t*>(it->second)[j] << " ";
1146 break;
1147 case 'J':
1148 str << reinterpret_cast<int32_t*>(it->second)[j] << " ";
1149 break;
1150 case 'K':
1151 str << reinterpret_cast<int64_t*>(it->second)[j] << " ";
1152 break;
1153 case 'E':
1154 str << reinterpret_cast<float*>(it->second)[j] << " ";
1155 break;
1156 case 'D':
1157 str << reinterpret_cast<double*>(it->second)[j] << " ";
1158 break;
1159 default:
1160 ;
1161 }
1162 }
1163 }
1164
1165 for (auto it=xValues.begin(); it!= xValues.end(); it++)
1166 {
1167 str >> (*it)[i-1];
1168 }
1169 if (unixTime)
1170 {
1171 long u1, u2;
1172 str >> u1 >> u2;
1173
1174 boost::posix_time::ptime unixTimeT( boost::gregorian::date(1970, boost::gregorian::Jan, 1),
1175 boost::posix_time::seconds(u1) + boost::posix_time::microsec(u2));
1176
1177 Time mjdTime(unixTimeT);
1178 yValues[i-1] = mjdTime.Mjd();
1179 if (yValues[i-1] < 40587)
1180 yValues[i-1] += 40587;
1181 }
1182 else
1183 {
1184 str >> yValues[i-1];
1185 if (yValues[i-1] < 40587)
1186 yValues[i-1] += 40587;
1187
1188 Time t(yValues[i-1]);
1189 string time = t.GetAsStr("%H:%M:%S%F");
1190 while (time[time.size()-1] == '0' && time.size() > 2)
1191 {
1192 time = time.substr(0, time.size()-1);
1193 }
1194 }
1195 if (i==1)
1196 {
1197 Time t(yValues[0]);
1198 title += " - " + t.GetAsStr("%Y-%m-%d");
1199 plot->setTitle(title.c_str());
1200 }
1201 }
1202 //set the actual data.
1203 auto jt = xValues.begin();
1204 for (auto it=curves.begin(); it != curves.end(); it++, jt++)
1205 (*it)->setRawData(yValues, *jt, endIndex-1);
1206
1207 QStack<QRectF> stack;
1208 double minX, minY, maxX, maxY;
1209 minX = minY = 1e10;
1210 maxX = maxY = -1e10;
1211 QRectF rect;
1212 QPointF point;
1213 for (auto it=curves.begin(); it!= curves.end(); it++)
1214 {
1215 rect = (*it)->boundingRect();
1216 point = rect.bottomRight();
1217 if (point.x() < minX) minX = point.x();
1218 if (point.y() < minY) minY = point.y();
1219 if (point.x() > maxX) maxX = point.x();
1220 if (point.y() > maxY) maxY = point.y();
1221 point = rect.topLeft();
1222 if (point.x() < minX) minX = point.x();
1223 if (point.y() < minY) minY = point.y();
1224 if (point.x() > maxX) maxX = point.x();
1225 if (point.y() > maxY) maxY = point.y();
1226 }
1227 QPointF bottomRight(maxX, minY);
1228 QPointF topLeft(minX, maxY);
1229 QPointF center((bottomRight+topLeft)/2.f);
1230 stack.push(QRectF(topLeft + (topLeft-center)*(.5f),bottomRight + (bottomRight-center)*(.5f)));
1231 zoom.setZoomStack(stack);
1232
1233// delete[] fitsBuffer;
1234 for (auto it = columnsData.begin(); it != columnsData.end(); it++)
1235 delete[] it->second;
1236 window.resize(600, 400);
1237 window.show();
1238
1239 a.exec();
1240
1241
1242 for (auto it = curves.begin(); it != curves.end(); it++)
1243 {
1244 (*it)->detach();
1245 delete *it;
1246 }
1247 grid->detach();
1248 for (auto it = xValues.begin(); it != xValues.end(); it++)
1249 delete[] *it;
1250
1251 delete[] yValues;
1252
1253 delete[] handle;
1254
1255 return 0;
1256}
1257#endif
1258
1259void SetupConfiguration(Configuration& conf)
1260{
1261 po::options_description configs("Fitsdump options");
1262 configs.add_options()
1263 ("fitsfile,f", var<string>()
1264#if BOOST_VERSION >= 104200
1265 ->required()
1266#endif
1267 , "Name of FITS file")
1268 ("tablename,t", var<string>("DATA")
1269#if BOOST_VERSION >= 104200
1270 ->required()
1271#endif
1272 , "Name of input table")
1273 ("col,c", vars<string>(), "List of columns to dump\narg is a list of columns, separated by a space.\nAdditionnally, a list of sub-columns can be added\ne.g. Data[3] will dump sub-column 3 of column Data\nData[3:4] will dump sub-columns 3 and 4\nOmitting this argument dump the entire column\nnota: all indices start at zero")
1274 ("outfile,o", var<string>("/dev/stdout"), "Name of output file (-:/dev/stdout)")
1275 ("precision,p", var<int>(20), "Precision of ofstream")
1276 ("list,l", po_switch(), "List all tables and columns in file")
1277 ("header,h", po_switch(), "Dump header of given table")
1278 ("stat,s", po_switch(), "Perform statistics instead of dump")
1279 ("minmax,m", po_switch(), "Calculates min and max of data")
1280 ("nozero,z", po_switch(), "skip 0 values for stats")
1281 ("tstart,a", po_switch(), "Give the mjdStart from reading the file data")
1282 ("tstop,b", po_switch(), "Give the mjdStop from reading the file data")
1283#ifdef PLOTTING_PLEASE
1284 ("graph,g", po_switch(), "Plot the columns instead of dumping them")
1285 ("dots,d", po_switch(), "Plot using dots instead of lines")
1286#endif
1287 ;
1288
1289 po::positional_options_description p;
1290 p.add("fitsfile", 1); // The first positional options
1291 p.add("col", -1); // All others
1292
1293 conf.AddOptions(configs);
1294 conf.SetArgumentPositions(p);
1295}
1296
1297int main(int argc, const char** argv)
1298{
1299 Configuration conf(argv[0]);
1300 conf.SetPrintUsage(PrintUsage);
1301 SetupConfiguration(conf);
1302
1303 if (!conf.DoParse(argc, argv, PrintHelp))
1304 return -1;
1305
1306 FitsDumper loader;
1307 return loader.ExecConfig(conf);
1308}
Note: See TracBrowser for help on using the repository browser.