source: trunk/src/Scantable.h@ 2019

Last change on this file since 2019 was 2012, checked in by WataruKawasaki, 14 years ago

New Development: Yes

JIRA Issue: Yes (CAS-2373, CAS-2620)

Ready for Test: Yes

Interface Changes: Yes

What Interface Changed: For Scantable::polyBaseline(), parameters and return value have been changed.

Test Programs:

Put in Release Notes: Yes

Module(s): sdbaseline, sd.linefinder

Description: (1) CAS-2373-related:

(1.1) Modified Scantable::polyBaseline() to have the row-based loop inside.

Now it fits and subtracts baseline for all rows and also output info
about the fitting result to logger and text file, while in the
previous version this method just did baseline fitting/subtraction
for one row only and had to be put inside a row-based loop at the
python side ("poly_baseline()" in scantable.py) and result output had
also to be controlled at the python side. Using a test data from NRO
45m telescope (348,000 rows, 512 channels), the processing time of
scantable.poly_baseline() has reduced from 130 minutes to 5-6 minutes.

(1.2) For accelerating the another polynomial fitting function, namely

scantable.auto_poly_baseline(), added a method
Scantable::autoPolyBaseline(). This basically does the same thing
with Scantable::polyBaseline(), but uses linefinfer also to
automatically flag the line regions.

(1.3) To make linefinder usable in Scantable class, added a method

linefinder.setdata(). This makes it possible to apply linefinder
for a float-list data given as a parameter, without setting scantable,
while it was indispensable to set scantable to use linefinger previously.

(2) CAS-2620-related:

Added Scantable::cubicSplineBaseline() and autoCubicSplineBaseline() for
fit baseline using the cubic spline function. Parameters include npiece
(number of polynomial pieces), clipthresh (clipping threshold), and
clipniter (maximum iteration number).



  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 19.9 KB
Line 
1//
2// C++ Interface: Scantable
3//
4// Description:
5//
6//
7// Author: Malte Marquarding <asap@atnf.csiro.au>, (C) 2005
8//
9// Copyright: See COPYING file that comes with this distribution
10//
11//
12#ifndef ASAPSCANTABLE_H
13#define ASAPSCANTABLE_H
14
15// STL
16#include <string>
17#include <vector>
18#include <iostream>
19#include <fstream>
20// AIPS++
21#include <casa/aips.h>
22#include <casa/Containers/Record.h>
23#include <casa/Arrays/MaskedArray.h>
24#include <casa/BasicSL/String.h>
25#include <casa/Utilities/CountedPtr.h>
26
27#include <tables/Tables/Table.h>
28#include <tables/Tables/ArrayColumn.h>
29#include <tables/Tables/ScalarColumn.h>
30
31#include <measures/TableMeasures/ScalarMeasColumn.h>
32
33#include <coordinates/Coordinates/SpectralCoordinate.h>
34
35#include <casa/Arrays/Vector.h>
36#include <casa/Quanta/Quantum.h>
37
38#include <casa/Exceptions/Error.h>
39
40#include "Logger.h"
41#include "STHeader.h"
42#include "STFrequencies.h"
43#include "STWeather.h"
44#include "STFocus.h"
45#include "STTcal.h"
46#include "STMolecules.h"
47#include "STSelector.h"
48#include "STHistory.h"
49#include "STPol.h"
50#include "STFit.h"
51#include "STFitEntry.h"
52#include "STFitter.h"
53
54namespace asap {
55
56/**
57 * This class contains and wraps a casa::Table, which is used to store
58 * all the information. This can be either a MemoryTable or a
59 * disk based Table.
60 * It provides access functions to the underlying table
61 * It contains n subtables:
62 * @li weather
63 * @li frequencies
64 * @li molecules
65 * @li tcal
66 * @li focus
67 * @li fits
68 *
69 * @brief The main ASAP data container
70 * @author Malte Marquarding
71 * @date
72 * @version
73*/
74class Scantable : private Logger
75{
76
77friend class STMath;
78
79public:
80 /**
81 * Default constructor
82 */
83 explicit Scantable(casa::Table::TableType ttype = casa::Table::Memory);
84
85 /**
86 * Create a Scantable object form an existing table on disk
87 * @param[in] name the name of the existing Scantable
88 */
89 explicit Scantable(const std::string& name, casa::Table::TableType ttype = casa::Table::Memory);
90
91 /// @fixme this is only sensible for MemoryTables....
92 Scantable(const Scantable& other, bool clear=true);
93
94 /**
95 * Destructor
96 */
97 virtual ~Scantable();
98
99 /**
100 * get a const reference to the underlying casa::Table
101 * @return const \ref casa::Table reference
102 */
103 const casa::Table& table() const;
104
105 /**
106 * get a reference to the underlying casa::Table with the Selection
107 * object applied if set
108 * @return casa::Table reference
109 */
110 casa::Table& table();
111
112
113 /**
114 * Get a handle to the selection object
115 * @return constant STSelector reference
116 */
117 const STSelector& getSelection() const { return selector_; }
118
119 /**
120 * Set the data to be a subset as defined by the STSelector
121 * @param selection a STSelector object
122 */
123 void setSelection(const STSelector& selection);
124
125 /**
126 * unset the selection of the data
127 */
128 void unsetSelection();
129 /**
130 * set the header
131 * @param[in] sth an STHeader object
132 */
133 void setHeader( const STHeader& sth );
134
135 /**
136 * get the header information
137 * @return an STHeader object
138 */
139 STHeader getHeader( ) const;
140
141 /**
142 * Checks if the "other" Scantable is conformant with this,
143 * i.e. if header values are the same.
144 * @param[in] other another Scantable
145 * @return true or false
146 */
147 bool conformant( const Scantable& other);
148
149 /**
150 *
151 * @param stype The type of the source, 0 = on, 1 = off
152 */
153 void setSourceType(int stype);
154
155
156 /**
157 * The number of scans in the table
158 * @return number of scans in the table
159 */
160 int nscan() const;
161
162 casa::MEpoch::Types getTimeReference() const;
163
164
165 casa::MEpoch getEpoch(int whichrow) const;
166
167 /**
168 * Get global antenna position
169 * @return casa::MPosition
170 */
171 casa::MPosition getAntennaPosition() const;
172
173 /**
174 * the @ref casa::MDirection for a specific row
175 * @param[in] whichrow the row number
176 * return casa::MDirection
177 */
178 casa::MDirection getDirection( int whichrow ) const;
179
180 /**
181 * get the direction type as a string, e.g. "J2000"
182 * @param[in] whichrow the row number
183 * return the direction string
184 */
185 std::string getDirectionString( int whichrow ) const;
186
187 /**
188 * set the direction type as a string, e.g. "J2000"
189 * @param[in] refstr the direction type
190 */
191 void setDirectionRefString(const std::string& refstr="");
192
193 /**
194 * get the direction reference string
195 * @return a string describing the direction reference
196 */
197 std::string getDirectionRefString() const;
198
199 /**
200 * Return the Flux unit of the data, e.g. "Jy" or "K"
201 * @return string
202 */
203 std::string getFluxUnit() const;
204
205 /**
206 * Set the Flux unit of the data
207 * @param unit a string representing the unit, e.g "Jy" or "K"
208 */
209 void setFluxUnit( const std::string& unit );
210
211 /**
212 * Set the Stokes type of the data
213 * @param feedtype a string representing the type, e.g "circular" or "linear"
214 */
215 void setFeedType( const std::string& feedtype );
216
217 /**
218 *
219 * @param instrument a string representing an insturment. see xxx
220 */
221 void setInstrument( const std::string& instrument );
222
223 /**
224 * (Re)calculate the azimuth and elevationnfor all rows
225 */
226 void calculateAZEL();
227
228 /**
229 * "hard" flag the data, this flags everything selected in setSelection()
230 * param[in] msk a boolean mask of length nchan describing the points to
231 * to be flagged
232 */
233 //void flag( const std::vector<bool>& msk = std::vector<bool>());
234 //void flag( const std::vector<bool>& msk = std::vector<bool>(), bool unflag=false);
235
236 void flag( int whichrow = -1, const std::vector<bool>& msk = std::vector<bool>(), bool unflag=false);
237
238 /**
239 * Flag the data in a row-based manner. (CAS-1433 Wataru Kawasaki)
240 * param[in] rows list of row numbers to be flagged
241 */
242 void flagRow( const std::vector<casa::uInt>& rows = std::vector<casa::uInt>(), bool unflag=false);
243
244 /**
245 * Get flagRow info at the specified row. If true, the whole data
246 * at the row should be flagged.
247 */
248 bool getFlagRow(int whichrow) const
249 { return (flagrowCol_(whichrow) > 0); }
250
251 /**
252 * Flag the data outside a specified range (in a channel-based manner).
253 * (CAS-1807 Wataru Kawasaki)
254 */
255 void clip(const casa::Float uthres, const casa::Float dthres, bool clipoutside, bool unflag);
256
257 /**
258 * Return a list of booleans with the size of nchan for a specified row, to get info
259 * about which channel is clipped.
260 */
261 std::vector<bool> getClipMask(int whichrow, const casa::Float uthres, const casa::Float dthres, bool clipoutside, bool unflag);
262 void srchChannelsToClip(casa::uInt whichrow, const casa::Float uthres, const casa::Float dthres, bool clipoutside, bool unflag,
263 casa::Vector<casa::uChar> flgs);
264
265 /**
266 * Return a list of row numbers with respect to the original table.
267 * @return a list of unsigned ints
268 */
269 std::vector<unsigned int> rownumbers() const;
270
271
272 /**
273 * Get the number of beams in the data or a specific scan
274 * @param scanno the scan number to get the number of beams for.
275 * If scanno<0 the number is retrieved from the header.
276 * @return an integer number
277 */
278 int nbeam(int scanno=-1) const;
279 /**
280 * Get the number of IFs in the data or a specific scan
281 * @param scanno the scan number to get the number of IFs for.
282 * If scanno<0 the number is retrieved from the header.
283 * @return an integer number
284 */
285 int nif(int scanno=-1) const;
286 /**
287 * Get the number of polarizations in the data or a specific scan
288 * @param scanno the scan number to get the number of polarizations for.
289 * If scanno<0 the number is retrieved from the header.
290 * @return an integer number
291 */
292 int npol(int scanno=-1) const;
293
294 std::string getPolType() const;
295
296
297 /**
298 * Get the number of integartion cycles
299 * @param scanno the scan number to get the number of rows for.
300 * If scanno<0 the number is retrieved from the header.
301 * @return the number of rows (for the specified scanno)
302 */
303 int nrow(int scanno=-1) const;
304
305 int getBeam(int whichrow) const;
306 std::vector<uint> getBeamNos() const { return getNumbers(beamCol_); }
307
308 int getIF(int whichrow) const;
309 std::vector<uint> getIFNos() const { return getNumbers(ifCol_); }
310
311 int getPol(int whichrow) const;
312 std::vector<uint> getPolNos() const { return getNumbers(polCol_); }
313
314 std::vector<uint> getScanNos() const { return getNumbers(scanCol_); }
315 int getScan(int whichrow) const { return scanCol_(whichrow); }
316
317 //TT addition
318 std::vector<uint> getMolNos() {return getNumbers(mmolidCol_); }
319
320 /**
321 * Get the number of channels in the data or a specific IF. This currently
322 * varies only with IF number
323 * @param ifno the IF number to get the number of channels for.
324 * If ifno<0 the number is retrieved from the header.
325 * @return an integer number
326 */
327 int nchan(int ifno=-1) const;
328 int getChannels(int whichrow) const;
329
330 int ncycle(int scanno=-1) const;
331 int getCycle(int whichrow) const { return cycleCol_(whichrow); }
332
333 double getInterval(int whichrow) const
334 { return integrCol_(whichrow); }
335
336 float getTsys(int whichrow) const
337 { return casa::Vector<casa::Float>(tsysCol_(whichrow))(0); }
338 float getElevation(int whichrow) const
339 { return elCol_(whichrow); }
340 float getAzimuth(int whichrow) const
341 { return azCol_(whichrow); }
342 float getParAngle(int whichrow) const
343 { return focus().getParAngle(mfocusidCol_(whichrow)); }
344 int getTcalId(int whichrow) const
345 { return mtcalidCol_(whichrow); }
346
347 std::string getSourceName(int whichrow) const
348 { return srcnCol_(whichrow); }
349
350 std::vector<bool> getMask(int whichrow) const;
351 std::vector<float> getSpectrum(int whichrow,
352 const std::string& poltype = "" ) const;
353
354 void setSpectrum(const std::vector<float>& spec, int whichrow);
355
356 std::string getPolarizationLabel(int index, const std::string& ptype) const
357 { return STPol::getPolLabel(index, ptype ); }
358
359 /**
360 * Write the Scantable to disk
361 * @param filename the output file name
362 */
363 void makePersistent(const std::string& filename);
364
365 std::vector<std::string> getHistory() const
366 { return historyTable_.getHistory(); };
367
368 void addHistory(const std::string& hist) { historyTable_.addEntry(hist); }
369
370 void appendToHistoryTable(const STHistory& otherhist)
371 { historyTable_.append(otherhist); }
372
373 std::string summary(bool verbose=false);
374 //std::string getTime(int whichrow=-1, bool showdate=true) const;
375 std::string getTime(int whichrow=-1, bool showdate=true, casa::uInt prec=0) const;
376 double getIntTime(int whichrow) const { return integrCol_(whichrow); }
377
378 // returns unit, conversion frame, doppler, base-frame
379
380 /**
381 * Get the frequency set up
382 * This is forwarded to the STFrequencies subtable
383 * @return unit, frame, doppler
384 */
385 std::vector<std::string> getCoordInfo() const
386 { return freqTable_.getInfo(); };
387
388 void setCoordInfo(std::vector<string> theinfo)
389 { return freqTable_.setInfo(theinfo); };
390
391
392 std::vector<double> getAbcissa(int whichrow) const;
393
394 std::vector<float> getWeather(int whichrow) const;
395
396 std::string getAbcissaLabel(int whichrow) const;
397 std::vector<double> getRestFrequencies() const
398 { return moleculeTable_.getRestFrequencies(); }
399 std::vector<double> getRestFrequency(int id) const
400 { return moleculeTable_.getRestFrequency(id); }
401
402 /**
403 void setRestFrequencies(double rf, const std::string& name = "",
404 const std::string& = "Hz");
405 **/
406 // Modified by Takeshi Nakazato 05/09/2008
407 /***
408 void setRestFrequencies(vector<double> rf, const vector<std::string>& name = "",
409 const std::string& = "Hz");
410 ***/
411 void setRestFrequencies(vector<double> rf,
412 const vector<std::string>& name = vector<std::string>(1,""),
413 const std::string& = "Hz");
414
415 //void setRestFrequencies(const std::string& name);
416 void setRestFrequencies(const vector<std::string>& name);
417
418 void shift(int npix);
419
420 casa::SpectralCoordinate getSpectralCoordinate(int whichrow) const;
421
422 void convertDirection(const std::string& newframe);
423
424 STFrequencies& frequencies() { return freqTable_; }
425 const STFrequencies& frequencies() const { return freqTable_; }
426 STWeather& weather() { return weatherTable_; }
427 const STWeather& weather() const { return weatherTable_; }
428 STFocus& focus() { return focusTable_; }
429 const STFocus& focus() const { return focusTable_; }
430 STTcal& tcal() { return tcalTable_; }
431 const STTcal& tcal() const { return tcalTable_; }
432 STMolecules& molecules() { return moleculeTable_; }
433 const STMolecules& molecules() const { return moleculeTable_; }
434 STHistory& history() { return historyTable_; }
435 const STHistory& history() const { return historyTable_; }
436 STFit& fit() { return fitTable_; }
437 const STFit& fit() const { return fitTable_; }
438
439 std::vector<std::string> columnNames() const;
440
441 void addFit(const STFitEntry& fit, int row);
442 STFitEntry getFit(int row) const
443 { STFitEntry fe; fitTable_.getEntry(fe, mfitidCol_(row)); return fe; }
444
445 //Added by TT
446 /**
447 * Get the antenna name
448 * @return antenna name string
449 */
450 casa::String getAntennaName() const;
451
452 /**
453 * For GBT MS data only. check a scan list
454 * against the information found in GBT_GO table for
455 * scan number orders to get correct pairs.
456 * @param[in] scan list
457 * @return status
458 */
459 int checkScanInfo(const std::vector<int>& scanlist) const;
460
461 /**
462 * Get the direction as a vector, for a specific row
463 * @param[in] whichrow the row numbyyer
464 * @return the direction in a vector
465 */
466 std::vector<double> getDirectionVector(int whichrow) const;
467
468 /**
469 * Set a flag indicating whether the data was parallactified
470 * @param[in] flag true or false
471 */
472 void parallactify(bool flag)
473 { focus().setParallactify(flag); }
474
475 /**
476 * Reshape spectrum
477 * @param[in] nmin, nmax minimum and maximum channel
478 * @param[in] irow row number
479 *
480 * 30/07/2008 Takeshi Nakazato
481 **/
482 void reshapeSpectrum( int nmin, int nmax ) throw( casa::AipsError );
483 void reshapeSpectrum( int nmin, int nmax, int irow ) ;
484
485 /**
486 * Change channel number under fixed bandwidth
487 * @param[in] nchan, dnu new channel number and spectral resolution
488 * @param[in] irow row number
489 *
490 * 27/08/2008 Takeshi Nakazato
491 **/
492 void regridChannel( int nchan, double dnu ) ;
493 void regridChannel( int nchan, double dnu, int irow ) ;
494
495 bool getFlagtraFast(int whichrow);
496
497 void polyBaseline(const std::vector<bool>& mask,
498 int order,
499 bool outLogger=false,
500 const std::string& blfile="");
501 void autoPolyBaseline(const std::vector<bool>& mask,
502 int order,
503 const std::vector<int>& edge,
504 float threshold=3.0,
505 int chanAvgLimit=1,
506 bool outLogger=false,
507 const std::string& blfile="");
508 void cubicSplineBaseline(const std::vector<bool>& mask,
509 int nPiece,
510 float thresClip,
511 int nIterClip,
512 bool outLogger=false,
513 const std::string& blfile="");
514 void autoCubicSplineBaseline(const std::vector<bool>& mask,
515 int nPiece,
516 float thresClip,
517 int nIterClip,
518 const std::vector<int>& edge,
519 float threshold=3.0,
520 int chanAvgLimit=1,
521 bool outLogger=false,
522 const std::string& blfile="");
523 float getRms(const std::vector<bool>& mask, int whichrow);
524 std::string formatBaselineParams(const std::vector<float>& params,
525 const std::vector<bool>& fixed,
526 float rms,
527 const std::string& masklist,
528 int whichrow,
529 bool verbose=false) const;
530 std::string formatPiecewiseBaselineParams(const std::vector<int>& ranges,
531 const std::vector<float>& params,
532 const std::vector<bool>& fixed,
533 float rms,
534 const std::string& masklist,
535 int whichrow,
536 bool verbose=false) const;
537
538
539private:
540
541 casa::Matrix<casa::Float> getPolMatrix( casa::uInt whichrow ) const;
542
543 /**
544 * Turns a time value into a formatted string
545 * @param x
546 * @return
547 */
548 std::string formatSec(casa::Double x) const;
549
550 std::string formatTime(const casa::MEpoch& me, bool showdate)const;
551 std::string formatTime(const casa::MEpoch& me, bool showdate, casa::uInt prec)const;
552
553 /**
554 * Turns a casa::MDirection into a nicely formatted string
555 * @param md an casa::MDirection
556 * @return
557 */
558 std::string formatDirection(const casa::MDirection& md) const;
559
560 /**
561 * Create a unique file name for the paged (temporary) table
562 * @return just the name
563 */
564 static casa::String generateName();
565
566 /**
567 * attach to cached columns
568 */
569 void attach();
570
571 /**
572 * Set up the main casa::Table
573 */
574 void setupMainTable();
575
576 void attachSubtables();
577 void copySubtables(const Scantable& other);
578
579 /**
580 * Convert an "old" asap1 style row index into a new index
581 * @param[in] therow
582 * @return and index into @table_
583 */
584 int rowToScanIndex(int therow);
585
586 std::vector<uint> getNumbers(const casa::ScalarColumn<casa::uInt>& col) const;
587
588 static const casa::uInt version_ = 3;
589
590 STSelector selector_;
591
592 casa::Table::TableType type_;
593
594 // the actual data
595 casa::Table table_;
596 casa::Table originalTable_;
597
598 STTcal tcalTable_;
599 STFrequencies freqTable_;
600 STWeather weatherTable_;
601 STFocus focusTable_;
602 STMolecules moleculeTable_;
603 STHistory historyTable_;
604 STFit fitTable_;
605
606 // Cached Columns to avoid reconstructing them for each row get/put
607 casa::ScalarColumn<casa::Double> integrCol_;
608 casa::MDirection::ScalarColumn dirCol_;
609 casa::MEpoch::ScalarColumn timeCol_;
610 casa::ScalarColumn<casa::Float> azCol_;
611 casa::ScalarColumn<casa::Float> elCol_;
612 casa::ScalarColumn<casa::String> srcnCol_, fldnCol_;
613 casa::ScalarColumn<casa::uInt> scanCol_, beamCol_, ifCol_, polCol_, cycleCol_, flagrowCol_;
614 casa::ScalarColumn<casa::Int> rbeamCol_, srctCol_;
615 casa::ArrayColumn<casa::Float> specCol_, tsysCol_;
616 casa::ArrayColumn<casa::uChar> flagsCol_;
617
618 // id in frequencies table
619 casa::ScalarColumn<casa::uInt> mfreqidCol_;
620 // id in tcal table
621 casa::ScalarColumn<casa::uInt> mtcalidCol_;
622
623 casa::ArrayColumn<casa::String> histitemCol_;
624 casa::ScalarColumn<casa::Int> mfitidCol_;
625 casa::ScalarColumn<casa::uInt> mweatheridCol_;
626
627 casa::ScalarColumn<casa::uInt> mfocusidCol_;
628
629 casa::ScalarColumn<casa::uInt> mmolidCol_;
630
631 static std::map<std::string, STPol::STPolFactory *> factories_;
632 void initFactories();
633
634 /**
635 * Add an auxiliary column to the main table and attach it to a
636 * cached column. Use for adding new columns that the original asap2
637 * tables do not have.
638 * @param[in] col reference to the cached column to be attached
639 * @param[in] colName column name in asap table
640 * @param[in] defValue default value to fill in the column
641 *
642 * 25/10/2009 Wataru Kawasaki
643 */
644 template<class T, class T2> void attachAuxColumnDef(casa::ScalarColumn<T>&,
645 const casa::String&,
646 const T2&);
647 template<class T, class T2> void attachAuxColumnDef(casa::ArrayColumn<T>&,
648 const casa::String&,
649 const casa::Array<T2>&);
650
651 void fitBaseline(const std::vector<bool>& mask, int whichrow, Fitter& fitter);
652 std::vector<float> doCubicSplineFitting(const std::vector<float>& data,
653 const std::vector<bool>& mask,
654 std::vector<int>& sectionRanges,
655 std::vector<float>& params,
656 int nPiece=2,
657 float thresClip=3.0,
658 int nIterClip=1);
659 bool hasSameNchanOverIFs();
660 std::string getMaskRangeList(const std::vector<bool>& mask,
661 int whichrow,
662 const casa::String& coordInfo,
663 bool hasSameNchan,
664 int firstIF,
665 bool silent=false);
666 std::vector<int> getMaskEdgeIndices(const std::vector<bool>& mask, bool getStartIndices=true);
667 std::string formatBaselineParamsHeader(int whichrow, const std::string& masklist, bool verbose) const;
668 std::string formatBaselineParamsFooter(float rms, bool verbose) const;
669 std::vector<bool> getCompositeChanMask(int whichrow, const std::vector<bool>& inMask);
670 //std::vector<bool> getCompositeChanMask(int whichrow, const std::vector<bool>& inMask, const std::vector<int>& edge, const int minEdgeSize, STLineFinder& lineFinder);
671
672 void applyChanFlag( casa::uInt whichrow, const std::vector<bool>& msk, casa::uChar flagval);
673
674
675};
676
677} // namespace
678
679#endif
Note: See TracBrowser for help on using the repository browser.