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

Last change on this file since 12717 was 12687, checked in by lyard, 13 years ago
removed plotting in fitsdump svn version
File size: 35.7 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 <CCfits/CCfits>
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#include "Time.h"
30#endif
31
32using namespace std;
33
34class MyColumn : public CCfits::Column
35{
36public:
37 const string &comment() const { return CCfits::Column::comment(); }
38 long width() const
39 {//the width() method seems to be buggy (or empty ?) as it always return 1... redo it ourselves.
40 string inter = format();
41 inter = inter.substr(0, inter.size()-1);
42 return atoi(inter.c_str());
43 }
44};
45
46#ifdef PLOTTING_PLEASE
47class TimeScaleDraw: public QwtScaleDraw
48{
49public:
50 virtual QwtText label(double v) const
51 {
52 Time t(v);
53 string time = t.GetAsStr("%H:%M:%S%F");
54 while (time[time.size()-1] == '0' && time.size() > 2)
55 {
56 time = time.substr(0, time.size()-1);
57 }
58 return QwtText(time.c_str());
59 }
60};
61#endif
62
63class FitsDumper
64{
65public:
66 FitsDumper();
67 ~FitsDumper();
68
69private:
70
71 CCfits::FITS *fFile; /// FITS pointer
72 CCfits::Table *fTable; /// Table pointer
73 map<string, MyColumn*> fColMap; /// map between the column names and their CCfits objects
74
75 // Convert CCfits::ValueType into a human readable string
76 string ValueTypeToStr(CCfits::ValueType type) const;
77
78 // Convert CCfits::ValueType into a number of associated bytes
79 int ValueTypeToSize(CCfits::ValueType type) const;
80
81 /// Calculate the buffer size required to read a row of the fits table, as well as the offsets to each column
82 vector<int> CalculateOffsets() const;
83
84 template<class T>
85 T PtrToValue(const unsigned char* &ptr) const;
86// template<class T>
87// double PtrToDouble(const unsigned char *ptr) const;
88// double PtrToDouble(const unsigned char *ptr, CCfits::ValueType type) const;
89
90 /// Write a single row of the selected data
91 int WriteRow(ostream &, const vector<MyColumn*> &, const vector<int> &, unsigned char *, const vector<pair<int, int> >&) const;
92
93 bool OpenFile(const string &); /// Open a file
94 bool OpenTable(const string &); /// Open a table
95
96 /// Lists all columns of an open file
97 void List();
98 void ListHeader();
99 void ListKeywords(ostream &);
100
101 bool separateColumnsFromRanges(const vector<string>& list,
102 vector<pair<int, int> >& ranges,
103 vector<string>& listNamesOnly);
104 /// Perform the dumping, based on the current dump list
105 bool Dump(const string &, const vector<string> &list, int);
106 ///Display the selected columns values VS time
107#ifdef PLOTTING_PLEASE
108 int doCurvesDisplay( const vector<string> &list, const string& tableName);
109#endif
110 // bool Plot(const vector<string> &list);
111
112public:
113 ///Configures the fitsLoader from the config file and/or command arguments.
114 int ExecConfig(Configuration& conf);
115};
116
117// --------------------------------------------------------------------------
118//
119//! Constructor
120//! @param out
121//! the ostream where to redirect the outputs
122//
123FitsDumper::FitsDumper() : fFile(0), fTable(0)
124{
125}
126
127// --------------------------------------------------------------------------
128//
129//! Destructor
130//
131FitsDumper::~FitsDumper()
132{
133 if (fFile)
134 delete fFile;
135}
136
137
138string FitsDumper::ValueTypeToStr(CCfits::ValueType type) const
139{
140 switch (type)
141 {
142 case CCfits::Tbyte: return "uint8_t";
143 case CCfits::Tushort: return "uint16_t";
144 case CCfits::Tshort: return "int16_t";
145 case CCfits::Tuint: return "uint32_t";
146 case CCfits::Tint: return "int32_t";
147 case CCfits::Tulong: return "uint32_t";
148 case CCfits::Tlong: return "int32_t";
149 case CCfits::Tlonglong: return "int64_t";
150 case CCfits::Tfloat: return "float";
151 case CCfits::Tdouble: return "double";
152
153 default:
154 return "unknwown";
155 }
156}
157
158int FitsDumper::ValueTypeToSize(CCfits::ValueType type) const
159{
160 switch (type)
161 {
162 case CCfits::Tbyte: return sizeof(uint8_t);
163 case CCfits::Tushort:
164 case CCfits::Tshort: return sizeof(uint16_t);
165 case CCfits::Tuint:
166 case CCfits::Tint:
167 case CCfits::Tulong:
168 case CCfits::Tlong: return sizeof(uint32_t);
169 case CCfits::Tlonglong: return sizeof(uint64_t);
170 case CCfits::Tfloat: return sizeof(float);
171 case CCfits::Tdouble: return sizeof(double);
172
173 default:
174 return 0;
175 }
176}
177
178template<class T>
179T FitsDumper::PtrToValue(const unsigned char* &ptr) const
180{
181 T t;
182 reverse_copy(ptr, ptr+sizeof(T), reinterpret_cast<unsigned char*>(&t));
183 ptr += sizeof(T);
184
185 return t;
186}
187/*
188template<class T>
189double FitsDumper::PtrToDouble(const unsigned char *ptr) const
190{
191 T t;
192 reverse_copy(ptr, ptr+sizeof(T), reinterpret_cast<unsigned char*>(&t));
193 return t;
194}
195
196double FitsDumper::PtrToDouble(const unsigned char *ptr, CCfits::ValueType type) const
197{
198 switch (type)
199 {
200 case CCfits::Tbyte: return PtrToDouble<uint8_t> (ptr);
201 case CCfits::Tushort: return PtrToDouble<uint16_t>(ptr);
202 case CCfits::Tuint:
203 case CCfits::Tulong: return PtrToDouble<uint32_t>(ptr);
204 case CCfits::Tshort: return PtrToDouble<int16_t> (ptr);
205 case CCfits::Tint:
206 case CCfits::Tlong: return PtrToDouble<int32_t> (ptr);
207 case CCfits::Tlonglong: return PtrToDouble<int64_t> (ptr);
208 case CCfits::Tfloat: return PtrToDouble<float> (ptr);
209 case CCfits::Tdouble: return PtrToDouble<double> (ptr);
210
211 default:
212 throw runtime_error("Data type not implemented yet.");
213 }
214}
215*/
216
217// --------------------------------------------------------------------------
218//
219//! Writes a single row of the selected FITS data to the output file.
220//!
221//! Data type \b not yet implemented:
222//! CCfits::Tnull, CCfits::Tbit, CCfits::Tlogical, CCfits::Tstring,
223//! CCfits::Tcomplex, CCfits::Tdblcomplex, CCfits::VTbit,
224//! CCfits::VTbyte, CCfits::VTlogical, CCfits::VTushort,
225//! CCfits::VTshort, CCfits::VTuint, CCfits::VTint, CCfits::VTulong,
226//! CCfits::VTlong, CCfits::VTlonglong, CCfits::VTfloat,
227//! CCfits::VTdouble, CCfits::VTcomplex, CCfits::VTdblcomplex
228//!
229//! @param offsets
230//! a vector containing the offsets to the columns (in bytes)
231//! @param targetFile
232//! the ofstream where to write to
233//! @param fitsBuffer
234//! the memory were the row has been loaded by cfitsio
235//
236int FitsDumper::WriteRow(ostream &out, const vector<MyColumn*> &list, const vector<int> &offsets, unsigned char* fitsBuffer, const vector<pair<int, int> >& ranges) const
237{
238 int cnt = 0;
239 vector<pair<int, int> >::const_iterator jt = ranges.begin();
240 for (vector<MyColumn*>::const_iterator it=list.begin(); it!=list.end(); it++, jt++)
241 {
242 if (jt == ranges.end())
243 {
244 cout << "ERROR: END OF RANGE POINTER REACHED" << endl;
245 return false;
246 }
247 const MyColumn *col = *it;
248
249 // CCfits starts counting at 1 not 0
250 const int offset = offsets[col->index()-1];
251
252 // Get the pointer to the array which we are supposed to print
253 const unsigned char *ptr = fitsBuffer + offset;
254
255 // Loop over all array entries
256 int sizeToSkip = 0;
257 switch (col->type())
258 {
259 case CCfits::Tbyte: sizeToSkip = sizeof(uint8_t); break;
260 case CCfits::Tushort: sizeToSkip = sizeof(uint16_t); break;
261 case CCfits::Tuint:
262 case CCfits::Tulong: sizeToSkip = sizeof(uint32_t); break;
263 case CCfits::Tshort: sizeToSkip = sizeof(int16_t); break;
264 case CCfits::Tint:
265 case CCfits::Tlong: sizeToSkip = sizeof(int32_t); break;
266 case CCfits::Tlonglong: sizeToSkip = sizeof(int64_t); break;
267 case CCfits::Tfloat: sizeToSkip = sizeof(float); break;
268 case CCfits::Tdouble: sizeToSkip = sizeof(double); break;
269 default:
270 cerr << "Data type not implemented yet." << endl;
271 return 0;
272 }
273 ptr += sizeToSkip*jt->first;
274 for (int width=jt->first; width<jt->second; width++)
275 {
276 switch (col->type())
277 {
278 case CCfits::Tbyte: out << PtrToValue<uint8_t> (ptr); break;
279 case CCfits::Tushort: out << PtrToValue<uint16_t>(ptr); break;
280 case CCfits::Tuint:
281 case CCfits::Tulong: out << PtrToValue<uint32_t>(ptr); break;
282 case CCfits::Tshort: out << PtrToValue<int16_t> (ptr); break;
283 case CCfits::Tint:
284 case CCfits::Tlong: out << PtrToValue<int32_t> (ptr); break;
285 case CCfits::Tlonglong: out << PtrToValue<int64_t> (ptr); break;
286 case CCfits::Tfloat: out << PtrToValue<float> (ptr); break;
287 case CCfits::Tdouble: out << PtrToValue<double> (ptr); break;
288
289 default:
290 cerr << "Data type not implemented yet." << endl;
291 return 0;
292 }
293
294 out << " ";
295 cnt++;
296 }
297 }
298
299 if (cnt>0)
300 out << endl;
301
302 if (out.fail())
303 {
304 cerr << "ERROR - writing output: " << strerror(errno) << endl;
305 return -1;
306 }
307
308 return cnt;
309}
310
311// --------------------------------------------------------------------------
312//
313//! Calculates the required buffer size for reading one row of the current
314//! table. Also calculates the offsets to all the columns
315//
316vector<int> FitsDumper::CalculateOffsets() const
317{
318 map<int,int> sizes;
319
320 for (map<string, MyColumn*>::const_iterator it=fColMap.begin();
321 it!=fColMap.end(); it++)
322 {
323 const int &width = it->second->width();
324 const int &idx = it->second->index();
325
326 const int size = ValueTypeToSize(it->second->type());
327 if (size==0)
328 {
329 cerr << "Data type " << (int)it->second->type() << " not implemented yet." << endl;
330 return vector<int>();
331 }
332//cout << "data size: " << size << " width: " << width << " index: " << idx << endl;
333 int realwidth = (width == 0)? 1 : width;
334 sizes[idx] = size*realwidth;
335 }
336
337 //calculate the offsets in the vector.
338 vector<int> result(1, 0);
339
340 int size = 0;
341 int idx = 0;
342
343 for (map<int,int>::const_iterator it=sizes.begin(); it!=sizes.end(); it++)
344 {
345 size += it->second;
346 result.push_back(size);
347
348 if (it->first == ++idx)
349 continue;
350
351 cerr << "Expected index " << idx << ", but found " << it->first << endl;
352 return vector<int>();
353 }
354
355 return result;
356}
357
358// --------------------------------------------------------------------------
359//
360//! Loads the fits file based on the current parameters
361//
362bool FitsDumper::OpenFile(const string &filename)
363{
364 if (fFile)
365 delete fFile;
366
367 ostringstream str;
368 try
369 {
370 fFile = new CCfits::FITS(filename);
371 }
372 catch (CCfits::FitsException e)
373 {
374 cerr << "Could not open FITS file " << filename << " reason: " << e.message() << endl;
375 return false;
376 }
377
378 return true;
379}
380
381bool FitsDumper::OpenTable(const string &tablename)
382{
383 if (!fFile)
384 {
385 cerr << "No file open." << endl;
386 return false;
387 }
388
389 const multimap< string, CCfits::ExtHDU * > extMap = fFile->extension();
390 if (extMap.find(tablename) == extMap.end())
391 {
392 cerr << "Table '" << tablename << "' not found." << endl;
393 return false;
394 }
395
396 fTable = dynamic_cast<CCfits::Table*>(extMap.find(tablename)->second);
397 if (!fTable)
398 {
399 cerr << "Object '" << tablename << "' returned not a CCfits::Table." << endl;
400 return false;
401 }
402
403 map<string, CCfits::Column*>& tCols = fTable->column();
404 for (map<string, CCfits::Column*>::const_iterator it = tCols.begin(); it != tCols.end(); it++)
405 {
406 fColMap.insert(make_pair(it->first, static_cast<MyColumn*>(it->second)));
407 }
408// fColMap = fTable->column();
409
410 fTable->makeThisCurrent();
411
412 fTable->getComments();
413 fTable->getHistory();
414 fTable->readAllKeys();
415
416 return true;
417}
418
419
420void FitsDumper::List()
421{
422 if (!fFile)
423 {
424 cerr << "No file open." << endl;
425 return;
426 }
427
428 cout << "\nFile: " << fFile->name() << "\n";
429
430 const multimap< string, CCfits::ExtHDU * > extMap = fFile->extension();
431 for (std::multimap<string, CCfits::ExtHDU*>::const_iterator it=extMap.begin(); it != extMap.end(); it++)
432 {
433
434 CCfits::Table *table = dynamic_cast<CCfits::Table*>(extMap.find(it->first)->second);
435
436 table->makeThisCurrent();
437 table->readData();
438
439 cout << " " << it->first << " [" << table->rows() << "]\n";
440
441 const map<string, CCfits::Column*> &cols = table->column();
442
443// for (map<string, CCfits::Column*>::const_iterator id = cols.begin(); ib != cols.end(); ib++)
444// {
445// TFORM
446// }
447
448 for (map<string, CCfits::Column*>::const_iterator ic=cols.begin();
449 ic != cols.end(); ic++)
450 {
451 const MyColumn *col = static_cast<MyColumn*>(ic->second);
452
453 cout << " " << col->name() << "[" ;
454 if (col->width() == 0)
455 cout << 1;
456 else
457 cout << col->width();
458 cout << "] * " << col->scale() << " (" << col->unit() << ":" << ValueTypeToStr(col->type())<< ") " << col->comment() << "\n";
459 /*
460 inline size_t Column::repeat () const
461 inline bool Column::varLength () const
462 inline double Column::zero () const
463 inline const String& Column::display () const
464 inline const String& Column::dimen () const
465 inline const String& Column::TBCOL ()
466 inline const String& Column::TTYPE ()
467 inline const String& Column::TFORM ()
468 inline const String& Column::TDISP ()
469 inline const String& Column::TZERO ()
470 inline const String& Column::TDIM ()
471 inline const String& Column::TNULL ()
472 inline const String& Column::TLMIN ()
473 inline const String& Column::TLMAX ()
474 inline const String& Column::TDMAX ()
475 inline const String& Column::TDMIN ()
476 */
477 }
478 cout << '\n';
479 }
480 cout << flush;
481}
482
483class MyKeyword : public CCfits::Keyword
484{
485public:
486 CCfits::ValueType keytype() const { return CCfits::Keyword::keytype(); }
487};
488
489void FitsDumper::ListKeywords(ostream &out)
490{
491 map<string,CCfits::Keyword*> keys = fTable->keyWord();
492 for (map<string,CCfits::Keyword*>::const_iterator it=keys.begin();
493 it!=keys.end(); it++)
494 {
495 string str;
496 double d;
497 int l;
498
499 const MyKeyword *kw = static_cast<MyKeyword*>(it->second);
500 kw->keytype();
501 out << "## " << setw(8) << kw->name() << " = " << setw(10);
502 if (kw->keytype()==16)
503 out << ("'"+kw->value(str)+"'");
504 if (kw->keytype()==31)
505 out << kw->value(l);
506 if (kw->keytype()==82)
507 out << kw->value(d);
508 out << " / " << kw->comment() << endl;
509 }
510}
511
512void FitsDumper::ListHeader()
513{
514 if (!fTable)
515 {
516 cerr << "No table open." << endl;
517 return;
518 }
519
520 cout << "\nTable: " << fTable->name() << " (rows=" << fTable->rows() << ")\n";
521 if (!fTable->comment().empty())
522 cout << "Comment: \t" << fTable->comment() << '\n';
523 if (!fTable->history().empty())
524 cout << "History: \t" << fTable->history() << '\n';
525
526 ListKeywords(cout);
527 cout << endl;
528}
529
530bool FitsDumper::separateColumnsFromRanges(const vector<string>& list,
531 vector<pair<int, int> >& ranges,
532 vector<string>& listNamesOnly)
533{
534 for (vector<string>::const_iterator it=list.begin(); it!=list.end(); it++)
535 {
536 string columnNameOnly = *it;
537 unsigned long bracketIndex0 = columnNameOnly.find_first_of('[');
538 unsigned long bracketIndex1 = columnNameOnly.find_first_of(']');
539 unsigned long colonIndex = columnNameOnly.find_first_of(':');
540// cout << bracketIndex0 << " " << bracketIndex1 << " " << colonIndex << endl;
541 int columnStart = -1;
542 int columnEnd = -1;
543 if (bracketIndex0 != string::npos)
544 {//there is a range given. Extract the range
545 if (colonIndex != string::npos)
546 {//we have a range here
547 columnStart = atoi(columnNameOnly.substr(bracketIndex0+1, colonIndex-(bracketIndex0+1)).c_str());
548 columnEnd = atoi(columnNameOnly.substr(colonIndex+1, bracketIndex1-(colonIndex+1)).c_str());
549 columnEnd++;
550 }
551 else
552 {//only a single index there
553 columnStart = atoi(columnNameOnly.substr(bracketIndex0+1, bracketIndex1 - (bracketIndex0+1)).c_str());
554 columnEnd = columnStart+1;
555// cout << "Cstart " << columnStart << " end: " << columnEnd << endl;
556 }
557 columnNameOnly = columnNameOnly.substr(0, bracketIndex0);
558 }
559
560 if (fColMap.find(columnNameOnly) == fColMap.end())
561 {
562 cerr << "ERROR - Column '" << columnNameOnly << "' not found in table." << endl;
563 return false;
564 }
565// cout << "The column name is: " << columnNameOnly << endl;
566 MyColumn *cCol = static_cast<MyColumn*>(fColMap.find(columnNameOnly)->second);
567 if (bracketIndex0 == string::npos)
568 {//no range given: use the full range
569 ranges.push_back(make_pair(0, cCol->width()));
570 columnStart = 0;
571 columnEnd = 1;
572 }
573 else
574 {//use the range extracted earlier
575 if (columnStart < 0)
576 {
577 cerr << "ERROR - Start range for column " << columnNameOnly << " is less than zero (" << columnStart << "). Aborting" << endl;
578 return false;
579 }
580 if (columnEnd>1 && columnEnd > cCol->width())
581 {
582 cerr << "ERROR - End range for column " << columnNameOnly << " is greater than the last element (" << cCol->width() << " vs " << columnEnd << "). Aborting" << endl;
583 return false;
584 }
585 ranges.push_back(make_pair(columnStart, columnEnd));
586 }
587// cout << "Will be exporting from " << columnStart << " to " << columnEnd-1 << " for column " << columnNameOnly << endl;
588 listNamesOnly.push_back(columnNameOnly);
589 }
590 return true;
591}
592// --------------------------------------------------------------------------
593//
594//! Perform the actual dump, based on the current parameters
595//
596bool FitsDumper::Dump(const string &filename, const vector<string> &list, int precision)
597{
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 // FIXME: Maybe do this when opening a table?
610 const vector<int> offsets = CalculateOffsets();
611 // for (int i=0;i<offsets.size(); i++)
612 // cout << offsets[i] << " ";
613// cout << endl;
614 if (offsets.size()==0)
615 return false;
616
617 ofstream out(filename=="-"?"/dev/stdout":filename);
618 if (!out)
619 {
620 cerr << "Cannot open file " << filename << ": " << strerror(errno) << endl;
621 return false;
622 }
623
624 out.precision(precision);
625
626 out << "## --------------------------------------------------------------------------\n";
627 if (filename!="-")
628 out << "## File: \t" << filename << '\n';
629 out << "## Table: \t" << fTable->name() << '\n';
630 if (!fTable->comment().empty())
631 out << "## Comment: \t" << fTable->comment() << '\n';
632 if (!fTable->history().empty())
633 out << "## History: \t" << fTable->history() << '\n';
634 out << "## NumRows: \t" << fTable->rows() << '\n';
635 out << "## --------------------------------------------------------------------------\n";
636 ListKeywords(out);
637 out << "## --------------------------------------------------------------------------\n";
638 out << "#\n";
639 //vector<pair<int, int> >::const_iterator rangesIt = ranges.begin();
640 //the above is soooo yesterday ;) let's try the new auto feature of c++0x
641 auto rangesIt = ranges.begin();
642 for (vector<string>::const_iterator it=listNamesOnly.begin(); it!=listNamesOnly.end(); it++, rangesIt++)
643 {
644 const MyColumn *col = static_cast<MyColumn*>(fTable->column()[*it]);
645 if (rangesIt->first != 0 || rangesIt->second != col->width())
646 {
647 out << "#";
648 for (int i=rangesIt->first; i<rangesIt->second; i++)
649 out << " " << col->name() << "[" << i << "]";
650 out << ": " << col->unit();
651 }
652 else
653 out << "# " << col->name() << "[" << col->width() << "]: " << col->unit();
654
655 if (!col->comment().empty())
656 out << " (" <<col->comment() << ")";
657 out << '\n';
658 }
659 out << "#" << endl;
660
661
662 // Loop over all columns in our list of requested columns
663 vector<MyColumn*> columns;
664 for (vector<string>::const_iterator it=listNamesOnly.begin(); it!=listNamesOnly.end(); it++)
665 {
666 MyColumn *cCol = static_cast<MyColumn*>(fColMap.find(*it)->second);
667 columns.push_back(cCol);
668 }
669 const int size = offsets[offsets.size()-1];
670 unsigned char* fitsBuffer = new unsigned char[size];
671
672 int status = 0;
673 int endIndex = (fColMap.find("Time") == fColMap.end()) ? fTable->rows()-1 : fTable->rows();
674 for (int i=1; i<=endIndex; i++)
675 {
676 fits_read_tblbytes(fFile->fitsPointer(), i, 1, size, fitsBuffer, &status);
677 if (status)
678 {
679 char text[30];//max length of cfitsio error strings (from doc)
680 fits_get_errstatus(status, text);
681
682 cerr << "Reading row " << i << " failed: " << text << " (fits_insert_rows,rc=" << status << ")" << endl;
683 break;
684 }
685 if (WriteRow(out, columns, offsets, fitsBuffer, ranges)<0)
686 {
687 status=1;
688 break;
689 }
690 }
691 delete[] fitsBuffer;
692
693 return status==0;
694}
695
696// --------------------------------------------------------------------------
697//
698//! Perform the actual dump, based on the current parameters
699//
700/*
701bool FitsDumper::Plot(const vector<string> &list)
702{
703 for (vector<string>::const_iterator it=list.begin(); it!=list.end(); it++)
704 if (fColMap.find(*it) == fColMap.end())
705 {
706 cerr << "WARNING - Column '" << *it << "' not found in table." << endl;
707 return false;
708 }
709
710 // FIXME: Maybe do this when opening a table?
711 const vector<int> offsets = CalculateOffsets();
712 if (offsets.size()==0)
713 return false;
714
715 // Loop over all columns in our list of requested columns
716 const CCfits::Column *col[3] =
717 {
718 list.size()>0 ? fColMap.find(list[0])->second : 0,
719 list.size()>1 ? fColMap.find(list[1])->second : 0,
720 list.size()>2 ? fColMap.find(list[2])->second : 0
721 };
722
723 const int size = offsets[offsets.size()-1];
724 unsigned char* fitsBuffer = new unsigned char[size];
725
726 const int idx = 0;
727
728 // CCfits starts counting at 1 not 0
729 const size_t pos[3] =
730 {
731 col[0] ? offsets[col[0]->index()-1] + idx*col[0]->width()*ValueTypeToSize(col[0]->type()) : 0,
732 col[1] ? offsets[col[1]->index()-1] + idx*col[1]->width()*ValueTypeToSize(col[1]->type()) : 0,
733 col[2] ? offsets[col[2]->index()-1] + idx*col[2]->width()*ValueTypeToSize(col[2]->type()) : 0
734 };
735
736 int status = 0;
737 for (int i=1; i<=fTable->rows(); i++)
738 {
739 fits_read_tblbytes(fFile->fitsPointer(), i, 1, size, fitsBuffer, &status);
740 if (status)
741 {
742 cerr << "An error occurred while reading fits row #" << i << " error code: " << status << endl;
743 break;
744 }
745
746 try
747 {
748 const double x[3] =
749 {
750 col[0] ? PtrToDouble(fitsBuffer+pos[0], col[0]->type()) : 0,
751 col[1] ? PtrToDouble(fitsBuffer+pos[1], col[1]->type()) : 0,
752 col[2] ? PtrToDouble(fitsBuffer+pos[2], col[2]->type()) : 0
753 };
754 }
755 catch (const runtime_error &e)
756 {
757 cerr << e.what() << endl;
758 status=1;
759 break;
760 }
761 }
762 delete[] fitsBuffer;
763
764 return status==0;
765}
766*/
767
768// --------------------------------------------------------------------------
769//
770//! Retrieves the configuration parameters
771//! @param conf
772//! the configuration object
773//
774int FitsDumper::ExecConfig(Configuration& conf)
775{
776 if (conf.Has("fitsfile"))
777 {
778 if (!OpenFile(conf.Get<string>("fitsfile")))
779 return -1;
780 }
781
782 if (conf.Get<bool>("list"))
783 List();
784
785 if (conf.Has("tablename"))
786 {
787 if (!OpenTable(conf.Get<string>("tablename")))
788 return -1;
789 }
790
791#ifdef PLOTTING_PLEASE
792 if (conf.Get<bool>("graph"))
793 {
794 if (!conf.Has("col"))
795 {
796 cout << "Please specify the columns that should be dumped as arguments. Aborting" << endl;
797 return 0;
798 }
799 doCurvesDisplay(conf.Get<vector<string>>("col"),
800 conf.Get<string>("tablename"));
801 return 1;
802 }
803#endif
804
805 if (conf.Get<bool>("header"))
806 ListHeader();
807
808 if (conf.Get<bool>("header") || conf.Get<bool>("list"))
809 return 1;
810
811 if (conf.Has("outfile"))
812 {
813 if (!conf.Has("col"))
814 {
815 cout << "Please specify the columns that should be dumped as arguments. Aborting" << endl;
816 return 0;
817 }
818 if (!Dump(conf.Get<string>("outfile"),
819 conf.Get<vector<string>>("col"),
820 conf.Get<int>("precision")))
821 return -1;
822 }
823
824
825 return 0;
826}
827
828void PrintUsage()
829{
830 cout <<
831 "fitsdump is a tool to dump data from a FITS table as ascii.\n"
832 "\n"
833 "Usage: fitsdump [OPTIONS] fitsfile col col ... \n"
834 " or: fitsdump [OPTIONS]\n";
835 cout << endl;
836}
837
838void PrintHelp()
839{
840 //
841}
842#ifdef PLOTTING_PLEASE
843int FitsDumper::doCurvesDisplay( const vector<string> &list, const string& tableName)
844{
845 //first of all, let's separate the columns from their ranges and check that the requested columns are indeed part of the file
846 vector<pair<int, int> > ranges;
847 vector<string> listNamesOnly;
848 if (!separateColumnsFromRanges(list, ranges, listNamesOnly))
849 {
850 cerr << "Something went wrong while extracting the columns names from parameters. Aborting" << endl;
851 return false;
852 }
853 vector<string> curvesNames;
854 stringstream str;
855 for (auto it=ranges.begin(), jt=listNamesOnly.begin(); it != ranges.end(); it++, jt++)
856 {
857 for (int i=it->first; i<it->second;i++)
858 {
859 str.str("");
860 str << *jt << "[" << i << "]";
861 curvesNames.push_back(str.str());
862 }
863 }
864 char* handle = new char[17];
865 sprintf(handle,"FitsDump Display");
866// Qt::HANDLE h = *handle;//NULL
867 int argc = 1;
868 char ** argv = &handle;
869 QApplication a(argc, argv);
870
871
872
873 QwtPlot* plot = new QwtPlot();
874 QwtPlotGrid* grid = new QwtPlotGrid;
875 grid->enableX(false);
876 grid->enableY(true);
877 grid->enableXMin(false);
878 grid->enableYMin(false);
879 grid->setMajPen(QPen(Qt::black, 0, Qt::DotLine));
880 grid->attach(plot);
881 plot->setAutoReplot(true);
882 string title = tableName;
883 plot->setAxisScaleDraw( QwtPlot::xBottom, new TimeScaleDraw());
884
885 QWidget window;
886 QHBoxLayout* layout = new QHBoxLayout(&window);
887 layout->setContentsMargins(0,0,0,0);
888 layout->addWidget(plot);
889
890 QwtPlotZoomer zoom(plot->canvas());
891 zoom.setRubberBandPen(QPen(Qt::gray, 2, Qt::DotLine));
892 zoom.setTrackerPen(QPen(Qt::gray));
893 int totalSize = 0;
894 for (unsigned int i=0;i<list.size();i++)
895 totalSize += ranges[i].second - ranges[i].first;
896 // cout << "Total size: " << totalSize << endl;
897 vector<QwtPlotCurve*> curves(totalSize);
898 int ii=0;
899 for (auto it = curves.begin(), jt=curvesNames.begin(); it != curves.end(); it++, jt++)
900 {
901 *it = new QwtPlotCurve(jt->c_str());
902 switch (ii%6)
903 {
904 case 0:
905 (*it)->setPen(QColor(255,0,0));
906 break;
907 case 1:
908 (*it)->setPen(QColor(0,255,0));
909 break;
910 case 2:
911 (*it)->setPen(QColor(0,0,255));
912 break;
913 case 3:
914 (*it)->setPen(QColor(255,255,0));
915 break;
916 case 4:
917 (*it)->setPen(QColor(0,255,255));
918 break;
919 case 5:
920 (*it)->setPen(QColor(255,0,255));
921 break;
922 default:
923 (*it)->setPen(QColor(0,0,0));
924 };
925 ii++;
926 (*it)->setStyle(QwtPlotCurve::Lines);
927 (*it)->attach(plot);
928 }
929 plot->insertLegend(new QwtLegend(), QwtPlot::RightLegend);
930
931 const vector<int> offsets = CalculateOffsets();
932 if (offsets.size()==0)
933 return false;
934
935
936 // Loop over all columns in our list of requested columns
937 vector<MyColumn*> columns;
938 for (vector<string>::const_iterator it=listNamesOnly.begin(); it!=listNamesOnly.end(); it++)
939 {
940 MyColumn *cCol = static_cast<MyColumn*>(fColMap.find(*it)->second);
941 columns.push_back(cCol);
942 }
943 //add the time column to the given columns
944 if (fColMap.find("Time") == fColMap.end() && fColMap.find("UnixTimeUTC") == fColMap.end())
945 {
946 cerr << "Error: time column could not be found in given table. Aborting" << endl;
947 return false;
948 }
949 MyColumn* timeCol;
950 bool unixTime = false;
951 int endIndex = 0;
952 if (fColMap.find("Time") != fColMap.end())
953 {
954 timeCol = static_cast<MyColumn*>(fColMap.find("Time")->second);
955 ranges.push_back(make_pair(0,1));
956 endIndex = fTable->rows();
957 }
958 if (fColMap.find("UnixTimeUTC") != fColMap.end())
959 {
960 timeCol = static_cast<MyColumn*>(fColMap.find("UnixTimeUTC")->second);
961 ranges.push_back(make_pair(0,2));
962 endIndex = fTable->rows()-1;
963 unixTime = true;
964 }
965 columns.push_back(timeCol);
966 /////
967 const int size = offsets[offsets.size()-1];
968 unsigned char* fitsBuffer = new unsigned char[size];
969
970// stringstream str;
971 str.str("");
972 int status = 0;
973
974 vector<double*> xValues(totalSize);
975 double* yValues;
976 cout.precision(10);
977 str.precision(20);
978 for (auto it=xValues.begin(); it!=xValues.end(); it++)
979 *it = new double[fTable->rows()];
980
981 yValues = new double[fTable->rows()];
982
983 cout.precision(3);
984 for (int i=1; i<=endIndex; i++)
985 {
986 cout << "/r" << "Constructing graph " << ((float)(i)/(float)(endIndex))*100.0 << "%";
987 fits_read_tblbytes(fFile->fitsPointer(), i, 1, size, fitsBuffer, &status);
988 if (status)
989 {
990 cerr << "An error occurred while reading fits row #" << i << " error code: " << status << endl;
991 break;
992 }
993 if (WriteRow(str, columns, offsets, fitsBuffer, ranges)<0)
994 {
995 status=1;
996 cerr << "An Error occured while reading the fits row " << i << endl;
997 return -1;
998 }
999// yValues[i-1] = i;
1000 for (auto it=xValues.begin(); it!= xValues.end(); it++)
1001 {
1002 str >> (*it)[i-1];
1003// cout << (*it)[i-1] << " ";
1004 }
1005 if (unixTime)
1006 {
1007 long u1, u2;
1008 str >> u1 >> u2;
1009 // cout << u1 << " " << u2;
1010 boost::posix_time::ptime unixTimeT( boost::gregorian::date(1970, boost::gregorian::Jan, 1),
1011 boost::posix_time::seconds(u1) + boost::posix_time::microsec(u2));
1012
1013 Time mjdTime(unixTimeT);
1014 yValues[i-1] = mjdTime.Mjd();
1015 // cout << " " << mjdTime.Mjd() << endl;
1016 }
1017 else
1018 {
1019 str >> yValues[i-1];
1020 if (yValues[i-1] < 40587)
1021 yValues[i-1] += 40587;
1022 }
1023 if (i==1)
1024 {
1025 Time t(yValues[0]);
1026 title += " - " + t.GetAsStr("%Y-%m-%d");
1027 plot->setTitle(title.c_str());
1028 }
1029// cout << yValues[i-1] << " ";
1030// cout << endl;
1031 }
1032 //set the actual data.
1033 auto jt = xValues.begin();
1034 for (auto it=curves.begin(); it != curves.end(); it++, jt++)
1035 (*it)->setRawData(yValues, *jt, endIndex);
1036
1037 QStack<QRectF> stack;
1038 double minX, minY, maxX, maxY;
1039 minX = minY = 1e10;
1040 maxX = maxY = -1e10;
1041 QRectF rect;
1042 QPointF point;
1043 for (auto it=curves.begin(); it!= curves.end(); it++)
1044 {
1045 rect = (*it)->boundingRect();
1046 point = rect.bottomRight();
1047 if (point.x() < minX) minX = point.x();
1048 if (point.y() < minY) minY = point.y();
1049 if (point.x() > maxX) maxX = point.x();
1050 if (point.y() > maxY) maxY = point.y();
1051 point = rect.topLeft();
1052 if (point.x() < minX) minX = point.x();
1053 if (point.y() < minY) minY = point.y();
1054 if (point.x() > maxX) maxX = point.x();
1055 if (point.y() > maxY) maxY = point.y();
1056 }
1057 QPointF bottomRight(maxX, minY);
1058 QPointF topLeft(minX, maxY);
1059 QPointF center((bottomRight+topLeft)/2.f);
1060 stack.push(QRectF(topLeft + (topLeft-center)*(.5f),bottomRight + (bottomRight-center)*(.5f)));
1061 zoom.setZoomStack(stack);
1062
1063 delete[] fitsBuffer;
1064 window.resize(600, 400);
1065 window.show();
1066
1067 a.exec();
1068
1069
1070 for (auto it = curves.begin(); it != curves.end(); it++)
1071 {
1072 (*it)->detach();
1073 delete *it;
1074 }
1075 grid->detach();
1076 for (auto it = xValues.begin(); it != xValues.end(); it++)
1077 delete[] *it;
1078 delete[] yValues;
1079 delete[] handle;
1080 return 0;
1081}
1082#endif
1083
1084void SetupConfiguration(Configuration& conf)
1085{
1086 po::options_description configs("Fitsdump options");
1087 configs.add_options()
1088 ("fitsfile,f", var<string>()
1089#if BOOST_VERSION >= 104200
1090 ->required()
1091#endif
1092 , "Name of FITS file")
1093 ("tablename,t", var<string>("DATA")
1094#if BOOST_VERSION >= 104200
1095 ->required()
1096#endif
1097 , "Name of input table")
1098 ("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")
1099 ("outfile,o", var<string>("/dev/stdout"), "Name of output file (-:/dev/stdout)")
1100 ("precision,p", var<int>(20), "Precision of ofstream")
1101 ("list,l", po_switch(), "List all tables and columns in file")
1102 ("header,h", po_switch(), "Dump header of given table")
1103#ifdef PLOTTING_PLEASE
1104 ("graph,g", po_switch(), "Plot the columns instead of dumping them")
1105#endif
1106 ;
1107
1108 po::positional_options_description p;
1109 p.add("fitsfile", 1); // The first positional options
1110 p.add("col", -1); // All others
1111
1112 conf.AddOptions(configs);
1113 conf.SetArgumentPositions(p);
1114}
1115
1116int main(int argc, const char** argv)
1117{
1118 Configuration conf(argv[0]);
1119 conf.SetPrintUsage(PrintUsage);
1120 SetupConfiguration(conf);
1121
1122 if (!conf.DoParse(argc, argv, PrintHelp))
1123 return -1;
1124
1125 FitsDumper loader;
1126 return loader.ExecConfig(conf);
1127}
Note: See TracBrowser for help on using the repository browser.