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

Last change on this file since 12900 was 12887, checked in by lyard, 14 years ago
fixed problem with sub-column selection
File size: 19.1 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 <float.h>
11
12#include <map>
13#include <fstream>
14
15#include <boost/regex.hpp>
16
17#include "Time.h"
18#include "externals/fits.h"
19
20using namespace std;
21
22struct MyColumn
23{
24 string name;
25
26 fits::Table::Column col;
27
28 uint32_t first;
29 uint32_t last;
30
31 void *ptr;
32};
33
34struct minMaxStruct
35{
36 double min;
37 double max;
38 long double average;
39 long double squared;
40 long numValues;
41 minMaxStruct() : min(FLT_MAX), max(-FLT_MAX), average(0), squared(0), numValues(0) { }
42
43 void add(double val)
44 {
45 average += val;
46 squared += val*val;
47
48 if (val<min)
49 min = val;
50
51 if (val>max)
52 max = val;
53
54 numValues++;
55 }
56};
57
58
59class FitsDumper : public fits
60{
61private:
62 string fFilename;
63
64 // Convert CCfits::ValueType into a human readable string
65 string ValueTypeToStr(char type) const;
66
67 /// Lists all columns of an open file
68 void List();
69 void ListHeader(const string& filename);
70 void ListKeywords(ostream &);
71
72 vector<MyColumn> InitColumns(vector<string> list);
73
74 ///Display the selected columns values VS time
75 void Dump(ofstream &, const vector<MyColumn> &, const string &);
76 void DumpMinMax(ofstream &, const vector<MyColumn> &, bool);
77 void DumpStats(ofstream &, const vector<MyColumn> &);
78
79public:
80 FitsDumper(const string &fname);
81
82 ///Configures the fitsLoader from the config file and/or command arguments.
83 int Exec(Configuration& conf);
84};
85
86// --------------------------------------------------------------------------
87//
88//! Constructor
89//! @param out
90//! the ostream where to redirect the outputs
91//
92FitsDumper::FitsDumper(const string &fname) : fits(fname), fFilename(fname)
93{
94}
95
96string FitsDumper::ValueTypeToStr(char type) const
97{
98 switch (type)
99 {
100 case 'L': return "bool(8)";
101 case 'A': return "char(8)";
102 case 'B': return "byte(8)";
103 case 'I': return "short(16)";
104 case 'J': return "int(32)";
105 case 'K': return "int(64)";
106 case 'E': return "float(32)";
107 case 'D': return "double(64)";
108 default:
109 return "unknown";
110 }
111}
112
113void FitsDumper::List()
114{
115 const fits::Table::Keys &fKeyMap = GetKeys();
116 const fits::Table::Columns &fColMap = GetColumns();
117
118 cout << "\nFile: " << fFilename << "\n";
119
120 cout << " " << fKeyMap.find("EXTNAME")->second.value << " [";
121 cout << fKeyMap.find("NAXIS2")->second.value << "]\n";
122
123 for (auto it = fColMap.begin(); it != fColMap.end(); it++)
124 {
125 cout << " " << it->first << "[" << it->second.num << "] (" << it->second.unit << ":" << ValueTypeToStr(it->second.type) << ") ";
126 for (auto jt = fKeyMap.begin(); jt != fKeyMap.end(); jt++)
127 if (jt->second.value == it->first)
128 cout << "/ " << jt->second.comment << endl;
129 }
130
131 cout << endl;
132}
133
134void FitsDumper::ListKeywords(ostream &fout)
135{
136 const fits::Table::Keys &fKeyMap = GetKeys();
137
138 for (auto it=fKeyMap.begin(); it != fKeyMap.end(); it++)
139 {
140 fout << "## " << ::left << setw(8) << it->first << "= ";
141
142 if (it->second.type=='T')
143 fout << ::left << setw(20) << ("'"+it->second.value+"'");
144 else
145 fout << ::right << setw(20) << it->second.value;
146
147 if (it->second.comment.size()>0)
148 fout << " / " << it->second.comment;
149 fout << '\n';
150 }
151
152 fout << flush;
153}
154
155void FitsDumper::ListHeader(const string& filename)
156{
157 ofstream fout(filename=="-"?"/dev/stdout":filename);
158 if (!fout)
159 {
160 cerr << "Cannot open file " << filename << ": " << strerror(errno) << endl;
161 return;
162 }
163
164 const fits::Table::Keys &fKeyMap = GetKeys();
165
166 fout << "\nTable: " << fKeyMap.find("EXTNAME")->second.value << " (rows=" << fKeyMap.find("NAXIS2")->second.value << ")\n";
167 if (fKeyMap.find("COMMENT") != fKeyMap.end())
168 fout << "Comment: \t" << fKeyMap.find("COMMENT")->second.value << "\n";
169
170 ListKeywords(fout);
171 fout << endl;
172}
173
174vector<MyColumn> FitsDumper::InitColumns(vector<string> names)
175{
176 static const boost::regex expr("([[:word:].]+)(\\[([[:digit:]]+)?(:)?([[:digit:]]+)?\\])?");
177
178 const fits::Table::Columns &fColMap = GetColumns();
179
180 if (names.size()==0)
181 for (auto it=fColMap.begin(); it!=fColMap.end(); it++)
182 if (it->second.num>0)
183 names.push_back(it->first);
184
185 vector<MyColumn> vec;
186
187 for (auto it=names.begin(); it!=names.end(); it++)
188 {
189 boost::smatch what;
190 if (!boost::regex_match(*it, what, expr, boost::match_extra))
191 {
192 cerr << "Couldn't parse expression '" << *it << "' " << endl;
193 return vector<MyColumn>();
194 }
195
196 const string name = what[1];
197
198 const auto iter = fColMap.find(name);
199 if (iter==fColMap.end())
200 {
201 cerr << "ERROR - Column '" << name << "' not found in table." << endl;
202 return vector<MyColumn>();
203 }
204
205 const fits::Table::Column &col = iter->second;
206
207 const string val0 = what[3];
208 const string delim = what[4];
209 const string val1 = what[5];
210
211 const uint32_t first = val0.empty() ? 0 : atoi(val0.c_str());
212 const uint32_t last = (val0.empty() && delim.empty()) ? col.num-1 : (val1.empty() ? first : atoi(val1.c_str()));
213
214 if (first>=col.num)
215 {
216 cerr << "ERROR - First index " << first << " for column " << name << " exceeds number of elements " << col.num << endl;
217 return vector<MyColumn>();
218 }
219
220 if (last>=col.num)
221 {
222 cerr << "ERROR - Last index " << last << " for column " << name << " exceeds number of elements " << col.num << endl;
223 return vector<MyColumn>();
224 }
225
226 if (first>last)
227 {
228 cerr << "ERROR - Last index " << last << " for column " << name << " exceeds first index " << first << endl;
229 return vector<MyColumn>();
230 }
231
232 MyColumn mycol;
233
234 mycol.name = name;
235 mycol.col = col;
236 mycol.first = first;
237 mycol.last = last;
238
239 vec.push_back(mycol);
240 }
241
242 for (auto it=vec.begin(); it!=vec.end(); it++)
243 it->ptr = SetPtrAddress(it->name);
244
245 return vec;
246}
247
248// --------------------------------------------------------------------------
249//
250//! Perform the actual dump, based on the current parameters
251//
252void FitsDumper::Dump(ofstream &fout, const vector<MyColumn> &cols, const string &filename)
253{
254 const fits::Table::Keys &fKeyMap = GetKeys();
255
256 fout << "## --------------------------------------------------------------------------\n";
257 fout << "## Fits file:\t" << fFilename << '\n';
258 if (filename!="-")
259 fout << "## File: \t" << filename << '\n';
260 fout << "## Table: \t" << fKeyMap.find("EXTNAME")->second.value << '\n';
261 fout << "## NumRows: \t" << GetInt("NAXIS2") << '\n';
262 fout << "## Comment: \t" << ((fKeyMap.find("COMMENT") != fKeyMap.end()) ? fKeyMap.find("COMMENT")->second.value : "") << '\n';
263 fout << "## --------------------------------------------------------------------------\n";
264 ListKeywords(fout);
265 fout << "## --------------------------------------------------------------------------\n";
266 fout << "#\n";
267
268 for (auto it=cols.begin(); it!=cols.end(); it++)
269 {
270 fout << "# " << it->name;
271
272 if (it->first==it->last)
273 {
274 if (it->first!=0)
275 fout << "[" << it->first << "]";
276 }
277 else
278 fout << "[" << it->first << ":" << it->last << "]";
279
280 fout << ": " << it->col.unit << '\n';
281 }
282 fout << "#" << endl;
283
284 // -----------------------------------------------------------------
285
286 while (GetNextRow())
287 {
288 const size_t row = GetRow();
289 if (row==GetNumRows())
290 break;
291
292 for (auto it=cols.begin(); it!=cols.end(); it++)
293 {
294 string msg;
295 for (uint32_t i=it->first; i<=it->last; i++)
296 {
297 switch (it->col.type)
298 {
299 case 'A':
300 msg += reinterpret_cast<const char*>(it->ptr)[i];
301 break;
302 case 'B':
303 fout << (unsigned int)reinterpret_cast<const unsigned char*>(it->ptr)[i] << " ";
304 break;
305 case 'L':
306 fout << reinterpret_cast<const bool*>(it->ptr)[i] << " ";
307 break;
308 case 'I':
309 fout << reinterpret_cast<const int16_t*>(it->ptr)[i] << " ";
310 break;
311 case 'J':
312 fout << reinterpret_cast<const int32_t*>(it->ptr)[i] << " ";
313 break;
314 case 'K':
315 fout << reinterpret_cast<const int64_t*>(it->ptr)[i] << " ";
316 break;
317 case 'E':
318 fout << reinterpret_cast<const float*>(it->ptr)[i] << " ";
319 break;
320 case 'D':
321 fout << reinterpret_cast<const double*>(it->ptr)[i] << " ";
322 break;
323 default:
324 ;
325 }
326 }
327
328 if (it->col.type=='A')
329 fout << "'" << msg << "' ";
330 }
331 fout << endl;
332 }
333}
334
335void FitsDumper::DumpMinMax(ofstream &fout, const vector<MyColumn> &cols, bool fNoZeroPlease)
336{
337 vector<minMaxStruct> statData(cols.size());
338
339 // Loop over all columns in our list of requested columns
340 while (GetNextRow())
341 {
342 const size_t row = GetRow();
343 if (row==GetNumRows())
344 break;
345
346 auto statsIt = statData.begin();
347
348 for (auto it=cols.begin(); it!=cols.end(); it++, statsIt++)
349 {
350 if ((it->name=="UnixTimeUTC" || it->name=="PCTime") && it->first==0 && it->last==1)
351 {
352 const uint32_t *val = reinterpret_cast<const uint32_t*>(it->ptr);
353 if (fNoZeroPlease && val[0]==0 && val[1]==0)
354 continue;
355
356 statsIt->add(Time(val[0], val[1]).Mjd());
357 continue;
358 }
359
360 for (uint32_t i=it->first; i<=it->last; i++)
361 {
362 double cValue = 0;
363 switch (it->col.type)
364 {
365 case 'L':
366 cValue = reinterpret_cast<const bool*>(it->ptr)[i];
367 break;
368 case 'B':
369 cValue = reinterpret_cast<const int8_t*>(it->ptr)[i];
370 break;
371 case 'I':
372 cValue = reinterpret_cast<const int16_t*>(it->ptr)[i];
373 break;
374 case 'J':
375 cValue = reinterpret_cast<const int32_t*>(it->ptr)[i];
376 break;
377 case 'K':
378 cValue = reinterpret_cast<const int64_t*>(it->ptr)[i];
379 break;
380 case 'E':
381 cValue = reinterpret_cast<const float*>(it->ptr)[i];
382 break;
383 case 'D':
384 cValue = reinterpret_cast<const double*>(it->ptr)[i];
385 break;
386 default:
387 ;
388 }
389
390 if (fNoZeroPlease && cValue == 0)
391 continue;
392
393 statsIt->add(cValue);
394 }
395 }
396 }
397
398 // okay. So now I've got ALL the data, loaded.
399 // let's do the summing and averaging in a safe way (i.e. avoid overflow
400 // of variables as much as possible)
401 auto statsIt = statData.begin();
402 for (auto it=cols.begin(); it!=cols.end(); it++, statsIt++)
403 {
404 fout << "\n[" << it->name << ':' << it->first;
405 if (it->first!=it->last)
406 fout << ':' << it->last;
407 fout << "]\n";
408
409 if (statsIt->numValues == 0)
410 {
411 fout << "Min: -\nMax: -\nAvg: -\nRms: -" << endl;
412 continue;
413 }
414
415 const long &num = statsIt->numValues;
416
417 long double &avg = statsIt->average;
418 long double &rms = statsIt->squared;
419
420 avg /= num;
421 rms = sqrt(rms/num - avg*avg);
422
423 fout << "Min: " << statsIt->min << '\n';
424 fout << "Max: " << statsIt->max << '\n';
425 fout << "Avg: " << avg << '\n';
426 fout << "Rms: " << rms << endl;
427 }
428}
429
430template<typename T>
431void displayStats(vector<char> &array, ofstream& out)
432{
433 const size_t numElems = array.size()/sizeof(T);
434 if (numElems == 0)
435 {
436 out << "Min: -\nMax: -\nMed: -\nAvg: -\nRms: -" << endl;
437 return;
438 }
439
440 T *val = reinterpret_cast<T*>(array.data());
441
442 sort(val, val+numElems);
443
444 out << "Min: " << double(val[0]) << '\n';
445 out << "Max: " << double(val[numElems-1]) << '\n';
446
447 if (numElems>2)
448 {
449 if (numElems%2 == 0)
450 out << "Med: " << (double(val[numElems/2]) + double(val[numElems/2+1]))/2 << '\n';
451 else
452 out << "Med: " << double(val[numElems/2+1]) << '\n';
453 }
454
455 long double avg = 0;
456 long double rms = 0;
457 for (uint32_t i=0;i<numElems;i++)
458 {
459 avg += double(val[i]);
460 rms += double(val[i])*double(val[i]);
461 }
462
463 avg /= numElems;
464 rms = sqrt(rms/numElems - avg*avg);
465
466 out << "Avg: " << avg << '\n';
467 out << "Rms: " << rms << endl;
468
469}
470
471void FitsDumper::DumpStats(ofstream &fout, const vector<MyColumn> &cols)
472{
473 // Loop over all columns in our list of requested columns
474 vector<vector<char>> statData;
475
476 for (auto it=cols.begin(); it!=cols.end(); it++)
477 statData.push_back(vector<char>(it->col.size*GetNumRows()*(it->last-it->first+1)));
478
479 while (GetNextRow())
480 {
481 const size_t row = GetRow();
482 if (row==GetNumRows())
483 break;
484
485 auto statsIt = statData.begin();
486 for (auto it=cols.begin(); it!=cols.end(); it++, statsIt++)
487 {
488 const char *src = reinterpret_cast<const char*>(it->ptr);
489 const size_t sz = (it->last-it->first+1)*it->col.size;
490 memcpy(statsIt->data()+row*sz, src+it->first*it->col.size, sz);
491 }
492 }
493
494 auto statsIt = statData.begin();
495 for (auto it=cols.begin(); it!=cols.end(); it++, statsIt++)
496 {
497 fout << "\n[" << it->name << ':' << it->first;
498 if (it->last!=it->first)
499 fout << ':' << it->last;
500 fout << "]\n";
501
502 switch (it->col.type)
503 {
504 case 'L':
505 displayStats<bool>(*statsIt, fout);
506 break;
507 case 'B':
508 displayStats<char>(*statsIt, fout);
509 break;
510 case 'I':
511 displayStats<int16_t>(*statsIt, fout);
512 break;
513 case 'J':
514 displayStats<int32_t>(*statsIt, fout);
515 break;
516 case 'K':
517 displayStats<int64_t>(*statsIt, fout);
518 break;
519 case 'E':
520 displayStats<float>(*statsIt, fout);
521 break;
522 case 'D':
523 displayStats<double>(*statsIt, fout);
524 break;
525 default:
526 ;
527 }
528 }
529}
530
531// --------------------------------------------------------------------------
532//
533//! Retrieves the configuration parameters
534//! @param conf
535//! the configuration object
536//
537int FitsDumper::Exec(Configuration& conf)
538{
539 if (conf.Get<bool>("list"))
540 List();
541
542 if (conf.Get<bool>("header"))
543 ListHeader(conf.Get<string>("outfile"));
544
545
546 if (conf.Get<bool>("header") || conf.Get<bool>("list"))
547 return 1;
548
549 // ------------------------------------------------------------
550
551 if (conf.Get<bool>("minmax") && conf.Get<bool>("stat"))
552 {
553 cerr << "Invalid combination of options: cannot do stats and minmax. Aborting" << endl;
554 return -1;
555 }
556 if (conf.Get<bool>("stat") && conf.Get<bool>("nozero"))
557 {
558 cerr << "Invalid combination of options: nozero only works with minmax. Aborting" << endl;
559 return -1;
560 }
561
562 // ------------------------------------------------------------
563
564 const string filename = conf.Get<string>("outfile");
565
566 ofstream fout(filename=="-"?"/dev/stdout":filename);
567 if (!fout)
568 {
569 cerr << "Cannot open file " << filename << ": " << strerror(errno) << endl;
570 return false;
571 }
572 fout.precision(conf.Get<int>("precision"));
573
574 const vector<MyColumn> cols = InitColumns(conf.Vec<string>("col"));
575 if (cols.size()==0)
576 return false;
577
578 if (conf.Get<bool>("minmax"))
579 {
580 DumpMinMax(fout, cols, conf.Get<bool>("nozero"));
581 return 0;
582 }
583
584 if (conf.Get<bool>("stat"))
585 {
586 DumpStats(fout, cols);
587 return 0;
588 }
589
590 Dump(fout, cols, filename);
591
592 return 0;
593}
594
595void PrintUsage()
596{
597 cout <<
598 "fitsdump is a tool to dump data from a FITS table as ascii.\n"
599 "\n"
600 "Usage: fitsdump [OPTIONS] fitsfile col col ... \n"
601 " or: fitsdump [OPTIONS]\n";
602 cout << endl;
603}
604
605void PrintHelp()
606{
607 //
608}
609
610void SetupConfiguration(Configuration& conf)
611{
612 po::options_description configs("Fitsdump options");
613 configs.add_options()
614 ("fitsfile,f", var<string>()
615#if BOOST_VERSION >= 104200
616 ->required()
617#endif
618 , "Name of FITS file")
619 ("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")
620 ("outfile,o", var<string>("/dev/stdout"), "Name of output file (-:/dev/stdout)")
621 ("precision,p", var<int>(20), "Precision of ofstream")
622 ("list,l", po_switch(), "List all tables and columns in file")
623 ("header,h", po_switch(), "Dump header of given table")
624 ("stat,s", po_switch(), "Perform statistics instead of dump")
625 ("minmax,m", po_switch(), "Calculates min and max of data")
626 ("nozero,z", po_switch(), "skip 0 values for stats")
627 ("force", po_switch(), "Force reading the fits file even if END key is missing")
628 ;
629
630 po::positional_options_description p;
631 p.add("fitsfile", 1); // The first positional options
632 p.add("col", -1); // All others
633
634 conf.AddOptions(configs);
635 conf.SetArgumentPositions(p);
636}
637
638int main(int argc, const char** argv)
639{
640 Configuration conf(argv[0]);
641 conf.SetPrintUsage(PrintUsage);
642 SetupConfiguration(conf);
643
644 if (!conf.DoParse(argc, argv, PrintHelp))
645 return -1;
646
647 if (!conf.Has("fitsfile"))
648 {
649 cerr << "Filename required." << endl;
650 return -1;
651 }
652
653 FitsDumper loader(conf.Get<string>("fitsfile"));
654 if (!loader)
655 {
656 cerr << "ERROR - Opening " << conf.Get<string>("fitsfile");
657 cerr << " failed: " << strerror(errno) << endl;
658 return -1;
659 }
660
661 return loader.Exec(conf);
662}
Note: See TracBrowser for help on using the repository browser.