source: trunk/src/Scantable.cpp@ 1176

Last change on this file since 1176 was 1170, checked in by mar637, 18 years ago

added line name to setRestFrequency

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 30.4 KB
Line 
1//
2// C++ Implementation: 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#include <map>
13
14#include <casa/aips.h>
15#include <casa/iostream.h>
16#include <casa/iomanip.h>
17#include <casa/OS/Path.h>
18#include <casa/OS/File.h>
19#include <casa/Arrays/Array.h>
20#include <casa/Arrays/ArrayMath.h>
21#include <casa/Arrays/MaskArrMath.h>
22#include <casa/Arrays/ArrayLogical.h>
23#include <casa/Arrays/ArrayAccessor.h>
24#include <casa/Arrays/VectorSTLIterator.h>
25#include <casa/Arrays/Vector.h>
26#include <casa/BasicMath/Math.h>
27#include <casa/BasicSL/Constants.h>
28#include <casa/Quanta/MVAngle.h>
29#include <casa/Containers/RecordField.h>
30#include <casa/Utilities/GenSort.h>
31
32#include <tables/Tables/TableParse.h>
33#include <tables/Tables/TableDesc.h>
34#include <tables/Tables/TableCopy.h>
35#include <tables/Tables/SetupNewTab.h>
36#include <tables/Tables/ScaColDesc.h>
37#include <tables/Tables/ArrColDesc.h>
38#include <tables/Tables/TableRow.h>
39#include <tables/Tables/TableVector.h>
40#include <tables/Tables/TableIter.h>
41
42#include <tables/Tables/ExprNode.h>
43#include <tables/Tables/TableRecord.h>
44#include <measures/Measures/MFrequency.h>
45#include <measures/Measures/MEpoch.h>
46#include <measures/Measures/MeasTable.h>
47#include <measures/Measures/MeasRef.h>
48#include <measures/TableMeasures/TableMeasRefDesc.h>
49#include <measures/TableMeasures/TableMeasValueDesc.h>
50#include <measures/TableMeasures/TableMeasDesc.h>
51#include <measures/TableMeasures/ScalarMeasColumn.h>
52#include <coordinates/Coordinates/CoordinateUtil.h>
53#include <casa/Quanta/MVTime.h>
54#include <casa/Quanta/MVAngle.h>
55
56#include "Scantable.h"
57#include "STPolLinear.h"
58#include "STPolStokes.h"
59#include "STAttr.h"
60#include "MathUtils.h"
61
62using namespace casa;
63
64namespace asap {
65
66std::map<std::string, STPol::STPolFactory *> Scantable::factories_;
67
68void Scantable::initFactories() {
69 if ( factories_.empty() ) {
70 Scantable::factories_["linear"] = &STPolLinear::myFactory;
71 Scantable::factories_["stokes"] = &STPolStokes::myFactory;
72 }
73}
74
75Scantable::Scantable(Table::TableType ttype) :
76 type_(ttype)
77{
78 initFactories();
79 setupMainTable();
80 freqTable_ = STFrequencies(*this);
81 table_.rwKeywordSet().defineTable("FREQUENCIES", freqTable_.table());
82 weatherTable_ = STWeather(*this);
83 table_.rwKeywordSet().defineTable("WEATHER", weatherTable_.table());
84 focusTable_ = STFocus(*this);
85 table_.rwKeywordSet().defineTable("FOCUS", focusTable_.table());
86 tcalTable_ = STTcal(*this);
87 table_.rwKeywordSet().defineTable("TCAL", tcalTable_.table());
88 moleculeTable_ = STMolecules(*this);
89 table_.rwKeywordSet().defineTable("MOLECULES", moleculeTable_.table());
90 historyTable_ = STHistory(*this);
91 table_.rwKeywordSet().defineTable("HISTORY", historyTable_.table());
92 fitTable_ = STFit(*this);
93 table_.rwKeywordSet().defineTable("FIT", fitTable_.table());
94 originalTable_ = table_;
95 attach();
96}
97
98Scantable::Scantable(const std::string& name, Table::TableType ttype) :
99 type_(ttype)
100{
101 initFactories();
102 Table tab(name, Table::Update);
103 uInt version = tab.keywordSet().asuInt("VERSION");
104 if (version != version_) {
105 throw(AipsError("Unsupported version of ASAP file."));
106 }
107 if ( type_ == Table::Memory ) {
108 table_ = tab.copyToMemoryTable(generateName());
109 } else {
110 table_ = tab;
111 }
112
113 attachSubtables();
114 originalTable_ = table_;
115 attach();
116}
117
118Scantable::Scantable( const Scantable& other, bool clear )
119{
120 // with or without data
121 String newname = String(generateName());
122 type_ = other.table_.tableType();
123 if ( other.table_.tableType() == Table::Memory ) {
124 if ( clear ) {
125 table_ = TableCopy::makeEmptyMemoryTable(newname,
126 other.table_, True);
127 } else
128 table_ = other.table_.copyToMemoryTable(newname);
129 } else {
130 other.table_.deepCopy(newname, Table::New, False,
131 other.table_.endianFormat(),
132 Bool(clear));
133 table_ = Table(newname, Table::Update);
134 table_.markForDelete();
135 }
136 /// @todo reindex SCANNO, recompute nbeam, nif, npol
137 if ( clear ) copySubtables(other);
138 attachSubtables();
139 originalTable_ = table_;
140 attach();
141}
142
143void Scantable::copySubtables(const Scantable& other) {
144 Table t = table_.rwKeywordSet().asTable("FREQUENCIES");
145 TableCopy::copyRows(t, other.freqTable_.table());
146 t = table_.rwKeywordSet().asTable("FOCUS");
147 TableCopy::copyRows(t, other.focusTable_.table());
148 t = table_.rwKeywordSet().asTable("WEATHER");
149 TableCopy::copyRows(t, other.weatherTable_.table());
150 t = table_.rwKeywordSet().asTable("TCAL");
151 TableCopy::copyRows(t, other.tcalTable_.table());
152 t = table_.rwKeywordSet().asTable("MOLECULES");
153 TableCopy::copyRows(t, other.moleculeTable_.table());
154 t = table_.rwKeywordSet().asTable("HISTORY");
155 TableCopy::copyRows(t, other.historyTable_.table());
156 t = table_.rwKeywordSet().asTable("FIT");
157 TableCopy::copyRows(t, other.fitTable_.table());
158}
159
160void Scantable::attachSubtables()
161{
162 freqTable_ = STFrequencies(table_);
163 focusTable_ = STFocus(table_);
164 weatherTable_ = STWeather(table_);
165 tcalTable_ = STTcal(table_);
166 moleculeTable_ = STMolecules(table_);
167 historyTable_ = STHistory(table_);
168 fitTable_ = STFit(table_);
169}
170
171Scantable::~Scantable()
172{
173 //cout << "~Scantable() " << this << endl;
174}
175
176void Scantable::setupMainTable()
177{
178 TableDesc td("", "1", TableDesc::Scratch);
179 td.comment() = "An ASAP Scantable";
180 td.rwKeywordSet().define("VERSION", uInt(version_));
181
182 // n Cycles
183 td.addColumn(ScalarColumnDesc<uInt>("SCANNO"));
184 // new index every nBeam x nIF x nPol
185 td.addColumn(ScalarColumnDesc<uInt>("CYCLENO"));
186
187 td.addColumn(ScalarColumnDesc<uInt>("BEAMNO"));
188 td.addColumn(ScalarColumnDesc<uInt>("IFNO"));
189 // linear, circular, stokes
190 td.rwKeywordSet().define("POLTYPE", String("linear"));
191 td.addColumn(ScalarColumnDesc<uInt>("POLNO"));
192
193 td.addColumn(ScalarColumnDesc<uInt>("FREQ_ID"));
194 td.addColumn(ScalarColumnDesc<uInt>("MOLECULE_ID"));
195 td.addColumn(ScalarColumnDesc<Int>("REFBEAMNO"));
196
197 td.addColumn(ScalarColumnDesc<Double>("TIME"));
198 TableMeasRefDesc measRef(MEpoch::UTC); // UTC as default
199 TableMeasValueDesc measVal(td, "TIME");
200 TableMeasDesc<MEpoch> mepochCol(measVal, measRef);
201 mepochCol.write(td);
202
203 td.addColumn(ScalarColumnDesc<Double>("INTERVAL"));
204
205 td.addColumn(ScalarColumnDesc<String>("SRCNAME"));
206 // Type of source (on=0, off=1, other=-1)
207 td.addColumn(ScalarColumnDesc<Int>("SRCTYPE", Int(-1)));
208 td.addColumn(ScalarColumnDesc<String>("FIELDNAME"));
209
210 //The actual Data Vectors
211 td.addColumn(ArrayColumnDesc<Float>("SPECTRA"));
212 td.addColumn(ArrayColumnDesc<uChar>("FLAGTRA"));
213 td.addColumn(ArrayColumnDesc<Float>("TSYS"));
214
215 td.addColumn(ArrayColumnDesc<Double>("DIRECTION",
216 IPosition(1,2),
217 ColumnDesc::Direct));
218 TableMeasRefDesc mdirRef(MDirection::J2000); // default
219 TableMeasValueDesc tmvdMDir(td, "DIRECTION");
220 // the TableMeasDesc gives the column a type
221 TableMeasDesc<MDirection> mdirCol(tmvdMDir, mdirRef);
222 // a uder set table type e.g. GALCTIC, B1950 ...
223 td.rwKeywordSet().define("DIRECTIONREF", String("J2000"));
224 // writing create the measure column
225 mdirCol.write(td);
226 td.addColumn(ScalarColumnDesc<Float>("AZIMUTH"));
227 td.addColumn(ScalarColumnDesc<Float>("ELEVATION"));
228 td.addColumn(ScalarColumnDesc<Float>("PARANGLE"));
229 td.addColumn(ScalarColumnDesc<Float>("OPACITY"));
230
231 td.addColumn(ScalarColumnDesc<uInt>("TCAL_ID"));
232 ScalarColumnDesc<Int> fitColumn("FIT_ID");
233 fitColumn.setDefault(Int(-1));
234 td.addColumn(fitColumn);
235
236 td.addColumn(ScalarColumnDesc<uInt>("FOCUS_ID"));
237 td.addColumn(ScalarColumnDesc<uInt>("WEATHER_ID"));
238
239 // columns which just get dragged along, as they aren't used in asap
240 td.addColumn(ScalarColumnDesc<Double>("SRCVELOCITY"));
241 td.addColumn(ArrayColumnDesc<Double>("SRCPROPERMOTION"));
242 td.addColumn(ArrayColumnDesc<Double>("SRCDIRECTION"));
243 td.addColumn(ArrayColumnDesc<Double>("SCANRATE"));
244
245 td.rwKeywordSet().define("OBSMODE", String(""));
246
247 // Now create Table SetUp from the description.
248 SetupNewTable aNewTab(generateName(), td, Table::Scratch);
249 table_ = Table(aNewTab, type_, 0);
250 originalTable_ = table_;
251}
252
253
254void Scantable::attach()
255{
256 timeCol_.attach(table_, "TIME");
257 srcnCol_.attach(table_, "SRCNAME");
258 srctCol_.attach(table_, "SRCTYPE");
259 specCol_.attach(table_, "SPECTRA");
260 flagsCol_.attach(table_, "FLAGTRA");
261 tsysCol_.attach(table_, "TSYS");
262 cycleCol_.attach(table_,"CYCLENO");
263 scanCol_.attach(table_, "SCANNO");
264 beamCol_.attach(table_, "BEAMNO");
265 ifCol_.attach(table_, "IFNO");
266 polCol_.attach(table_, "POLNO");
267 integrCol_.attach(table_, "INTERVAL");
268 azCol_.attach(table_, "AZIMUTH");
269 elCol_.attach(table_, "ELEVATION");
270 dirCol_.attach(table_, "DIRECTION");
271 paraCol_.attach(table_, "PARANGLE");
272 fldnCol_.attach(table_, "FIELDNAME");
273 rbeamCol_.attach(table_, "REFBEAMNO");
274
275 mfitidCol_.attach(table_,"FIT_ID");
276 mfreqidCol_.attach(table_, "FREQ_ID");
277 mtcalidCol_.attach(table_, "TCAL_ID");
278 mfocusidCol_.attach(table_, "FOCUS_ID");
279 mmolidCol_.attach(table_, "MOLECULE_ID");
280}
281
282void Scantable::setHeader(const STHeader& sdh)
283{
284 table_.rwKeywordSet().define("nIF", sdh.nif);
285 table_.rwKeywordSet().define("nBeam", sdh.nbeam);
286 table_.rwKeywordSet().define("nPol", sdh.npol);
287 table_.rwKeywordSet().define("nChan", sdh.nchan);
288 table_.rwKeywordSet().define("Observer", sdh.observer);
289 table_.rwKeywordSet().define("Project", sdh.project);
290 table_.rwKeywordSet().define("Obstype", sdh.obstype);
291 table_.rwKeywordSet().define("AntennaName", sdh.antennaname);
292 table_.rwKeywordSet().define("AntennaPosition", sdh.antennaposition);
293 table_.rwKeywordSet().define("Equinox", sdh.equinox);
294 table_.rwKeywordSet().define("FreqRefFrame", sdh.freqref);
295 table_.rwKeywordSet().define("FreqRefVal", sdh.reffreq);
296 table_.rwKeywordSet().define("Bandwidth", sdh.bandwidth);
297 table_.rwKeywordSet().define("UTC", sdh.utc);
298 table_.rwKeywordSet().define("FluxUnit", sdh.fluxunit);
299 table_.rwKeywordSet().define("Epoch", sdh.epoch);
300 table_.rwKeywordSet().define("POLTYPE", sdh.poltype);
301}
302
303STHeader Scantable::getHeader() const
304{
305 STHeader sdh;
306 table_.keywordSet().get("nBeam",sdh.nbeam);
307 table_.keywordSet().get("nIF",sdh.nif);
308 table_.keywordSet().get("nPol",sdh.npol);
309 table_.keywordSet().get("nChan",sdh.nchan);
310 table_.keywordSet().get("Observer", sdh.observer);
311 table_.keywordSet().get("Project", sdh.project);
312 table_.keywordSet().get("Obstype", sdh.obstype);
313 table_.keywordSet().get("AntennaName", sdh.antennaname);
314 table_.keywordSet().get("AntennaPosition", sdh.antennaposition);
315 table_.keywordSet().get("Equinox", sdh.equinox);
316 table_.keywordSet().get("FreqRefFrame", sdh.freqref);
317 table_.keywordSet().get("FreqRefVal", sdh.reffreq);
318 table_.keywordSet().get("Bandwidth", sdh.bandwidth);
319 table_.keywordSet().get("UTC", sdh.utc);
320 table_.keywordSet().get("FluxUnit", sdh.fluxunit);
321 table_.keywordSet().get("Epoch", sdh.epoch);
322 table_.keywordSet().get("POLTYPE", sdh.poltype);
323 return sdh;
324}
325
326void asap::Scantable::setSourceType( int stype )
327{
328 if ( stype < 0 || stype > 1 )
329 throw(AipsError("Illegal sourcetype."));
330 TableVector<Int> tabvec(table_, "SRCTYPE");
331 tabvec = Int(stype);
332}
333
334bool Scantable::conformant( const Scantable& other )
335{
336 return this->getHeader().conformant(other.getHeader());
337}
338
339
340
341std::string Scantable::formatSec(Double x) const
342{
343 Double xcop = x;
344 MVTime mvt(xcop/24./3600.); // make days
345
346 if (x < 59.95)
347 return String(" ") + mvt.string(MVTime::TIME_CLEAN_NO_HM, 7)+"s";
348 else if (x < 3599.95)
349 return String(" ") + mvt.string(MVTime::TIME_CLEAN_NO_H,7)+" ";
350 else {
351 ostringstream oss;
352 oss << setw(2) << std::right << setprecision(1) << mvt.hour();
353 oss << ":" << mvt.string(MVTime::TIME_CLEAN_NO_H,7) << " ";
354 return String(oss);
355 }
356};
357
358std::string Scantable::formatDirection(const MDirection& md) const
359{
360 Vector<Double> t = md.getAngle(Unit(String("rad"))).getValue();
361 Int prec = 7;
362
363 MVAngle mvLon(t[0]);
364 String sLon = mvLon.string(MVAngle::TIME,prec);
365 uInt tp = md.getRef().getType();
366 if (tp == MDirection::GALACTIC ||
367 tp == MDirection::SUPERGAL ) {
368 sLon = mvLon(0.0).string(MVAngle::ANGLE_CLEAN,prec);
369 }
370 MVAngle mvLat(t[1]);
371 String sLat = mvLat.string(MVAngle::ANGLE+MVAngle::DIG2,prec);
372 return sLon + String(" ") + sLat;
373}
374
375
376std::string Scantable::getFluxUnit() const
377{
378 return table_.keywordSet().asString("FluxUnit");
379}
380
381void Scantable::setFluxUnit(const std::string& unit)
382{
383 String tmp(unit);
384 Unit tU(tmp);
385 if (tU==Unit("K") || tU==Unit("Jy")) {
386 table_.rwKeywordSet().define(String("FluxUnit"), tmp);
387 } else {
388 throw AipsError("Illegal unit - must be compatible with Jy or K");
389 }
390}
391
392void Scantable::setInstrument(const std::string& name)
393{
394 bool throwIt = true;
395 // create an Instrument to see if this is valid
396 STAttr::convertInstrument(name, throwIt);
397 String nameU(name);
398 nameU.upcase();
399 table_.rwKeywordSet().define(String("AntennaName"), nameU);
400}
401
402MPosition Scantable::getAntennaPosition () const
403{
404 Vector<Double> antpos;
405 table_.keywordSet().get("AntennaPosition", antpos);
406 MVPosition mvpos(antpos(0),antpos(1),antpos(2));
407 return MPosition(mvpos);
408}
409
410void Scantable::makePersistent(const std::string& filename)
411{
412 String inname(filename);
413 Path path(inname);
414 /// @todo reindex SCANNO, recompute nbeam, nif, npol
415 inname = path.expandedName();
416 table_.deepCopy(inname, Table::New);
417}
418
419int Scantable::nbeam( int scanno ) const
420{
421 if ( scanno < 0 ) {
422 Int n;
423 table_.keywordSet().get("nBeam",n);
424 return int(n);
425 } else {
426 // take the first POLNO,IFNO,CYCLENO as nbeam shouldn't vary with these
427 Table t = table_(table_.col("SCANNO") == scanno);
428 ROTableRow row(t);
429 const TableRecord& rec = row.get(0);
430 Table subt = t( t.col("IFNO") == Int(rec.asuInt("IFNO"))
431 && t.col("POLNO") == Int(rec.asuInt("POLNO"))
432 && t.col("CYCLENO") == Int(rec.asuInt("CYCLENO")) );
433 ROTableVector<uInt> v(subt, "BEAMNO");
434 return int(v.nelements());
435 }
436 return 0;
437}
438
439int Scantable::nif( int scanno ) const
440{
441 if ( scanno < 0 ) {
442 Int n;
443 table_.keywordSet().get("nIF",n);
444 return int(n);
445 } else {
446 // take the first POLNO,BEAMNO,CYCLENO as nbeam shouldn't vary with these
447 Table t = table_(table_.col("SCANNO") == scanno);
448 ROTableRow row(t);
449 const TableRecord& rec = row.get(0);
450 Table subt = t( t.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
451 && t.col("POLNO") == Int(rec.asuInt("POLNO"))
452 && t.col("CYCLENO") == Int(rec.asuInt("CYCLENO")) );
453 if ( subt.nrow() == 0 ) return 0;
454 ROTableVector<uInt> v(subt, "IFNO");
455 return int(v.nelements());
456 }
457 return 0;
458}
459
460int Scantable::npol( int scanno ) const
461{
462 if ( scanno < 0 ) {
463 Int n;
464 table_.keywordSet().get("nPol",n);
465 return n;
466 } else {
467 // take the first POLNO,IFNO,CYCLENO as nbeam shouldn't vary with these
468 Table t = table_(table_.col("SCANNO") == scanno);
469 ROTableRow row(t);
470 const TableRecord& rec = row.get(0);
471 Table subt = t( t.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
472 && t.col("IFNO") == Int(rec.asuInt("IFNO"))
473 && t.col("CYCLENO") == Int(rec.asuInt("CYCLENO")) );
474 if ( subt.nrow() == 0 ) return 0;
475 ROTableVector<uInt> v(subt, "POLNO");
476 return int(v.nelements());
477 }
478 return 0;
479}
480
481int Scantable::ncycle( int scanno ) const
482{
483 if ( scanno < 0 ) {
484 Block<String> cols(2);
485 cols[0] = "SCANNO";
486 cols[1] = "CYCLENO";
487 TableIterator it(table_, cols);
488 int n = 0;
489 while ( !it.pastEnd() ) {
490 ++n;
491 ++it;
492 }
493 return n;
494 } else {
495 Table t = table_(table_.col("SCANNO") == scanno);
496 ROTableRow row(t);
497 const TableRecord& rec = row.get(0);
498 Table subt = t( t.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
499 && t.col("POLNO") == Int(rec.asuInt("POLNO"))
500 && t.col("IFNO") == Int(rec.asuInt("IFNO")) );
501 if ( subt.nrow() == 0 ) return 0;
502 return int(subt.nrow());
503 }
504 return 0;
505}
506
507
508int Scantable::nrow( int scanno ) const
509{
510 return int(table_.nrow());
511}
512
513int Scantable::nchan( int ifno ) const
514{
515 if ( ifno < 0 ) {
516 Int n;
517 table_.keywordSet().get("nChan",n);
518 return int(n);
519 } else {
520 // take the first SCANNO,POLNO,BEAMNO,CYCLENO as nbeam shouldn't vary with these
521 Table t = table_(table_.col("IFNO") == ifno);
522 if ( t.nrow() == 0 ) return 0;
523 ROArrayColumn<Float> v(t, "SPECTRA");
524 return v.shape(0)(0);
525 }
526 return 0;
527}
528
529int Scantable::nscan() const {
530 Vector<uInt> scannos(scanCol_.getColumn());
531 uInt nout = genSort( scannos, Sort::Ascending,
532 Sort::QuickSort|Sort::NoDuplicates );
533 return int(nout);
534}
535
536int Scantable::getChannels(int whichrow) const
537{
538 return specCol_.shape(whichrow)(0);
539}
540
541int Scantable::getBeam(int whichrow) const
542{
543 return beamCol_(whichrow);
544}
545
546std::vector<uint> Scantable::getNumbers(ScalarColumn<uInt>& col)
547{
548 Vector<uInt> nos(col.getColumn());
549 uInt n = genSort( nos, Sort::Ascending, Sort::QuickSort|Sort::NoDuplicates );
550 nos.resize(n, True);
551 std::vector<uint> stlout;
552 nos.tovector(stlout);
553 return stlout;
554}
555
556int Scantable::getIF(int whichrow) const
557{
558 return ifCol_(whichrow);
559}
560
561int Scantable::getPol(int whichrow) const
562{
563 return polCol_(whichrow);
564}
565
566std::string Scantable::formatTime(const MEpoch& me, bool showdate) const
567{
568 MVTime mvt(me.getValue());
569 if (showdate)
570 mvt.setFormat(MVTime::YMD);
571 else
572 mvt.setFormat(MVTime::TIME);
573 ostringstream oss;
574 oss << mvt;
575 return String(oss);
576}
577
578void Scantable::calculateAZEL()
579{
580 MPosition mp = getAntennaPosition();
581 MEpoch::ROScalarColumn timeCol(table_, "TIME");
582 ostringstream oss;
583 oss << "Computed azimuth/elevation using " << endl
584 << mp << endl;
585 for (Int i=0; i<nrow(); ++i) {
586 MEpoch me = timeCol(i);
587 MDirection md = getDirection(i);
588 oss << " Time: " << formatTime(me,False) << " Direction: " << formatDirection(md)
589 << endl << " => ";
590 MeasFrame frame(mp, me);
591 Vector<Double> azel =
592 MDirection::Convert(md, MDirection::Ref(MDirection::AZEL,
593 frame)
594 )().getAngle("rad").getValue();
595 azCol_.put(i,Float(azel[0]));
596 elCol_.put(i,Float(azel[1]));
597 oss << "azel: " << azel[0]/C::pi*180.0 << " "
598 << azel[1]/C::pi*180.0 << " (deg)" << endl;
599 }
600 pushLog(String(oss));
601}
602
603void Scantable::flag(const std::vector<bool>& msk)
604{
605 if ( selector_.empty() && msk.size() == 0 )
606 throw(AipsError("Trying to flag whole scantable."));
607 if ( msk.size() == 0 ) {
608 uChar userflag = 1 << 7;
609 for ( uInt i=0; i<table_.nrow(); ++i) {
610 Vector<uChar> flgs = flagsCol_(i);
611 flgs = userflag;
612 flagsCol_.put(i, flgs);
613 }
614 return;
615 }
616 if ( int(msk.size()) != nchan() ) {
617 throw(AipsError("Mask has incorrect number of channels."));
618 }
619 for ( uInt i=0; i<table_.nrow(); ++i) {
620 Vector<uChar> flgs = flagsCol_(i);
621 if ( flgs.nelements() != msk.size() ) {
622 throw(AipsError("Mask has incorrect number of channels."
623 " Probably varying with IF. Please flag per IF"));
624 }
625 std::vector<bool>::const_iterator it;
626 uInt j = 0;
627 uChar userflag = 1 << 7;
628 for (it = msk.begin(); it != msk.end(); ++it) {
629 if ( *it ) {
630 flgs(j) = userflag;
631 }
632 ++j;
633 }
634 flagsCol_.put(i, flgs);
635 }
636}
637
638std::vector<bool> Scantable::getMask(int whichrow) const
639{
640 Vector<uChar> flags;
641 flagsCol_.get(uInt(whichrow), flags);
642 Vector<Bool> bflag(flags.shape());
643 convertArray(bflag, flags);
644 bflag = !bflag;
645 std::vector<bool> mask;
646 bflag.tovector(mask);
647 return mask;
648}
649
650std::vector<float> Scantable::getSpectrum( int whichrow,
651 const std::string& poltype ) const
652{
653 String ptype = poltype;
654 if (poltype == "" ) ptype = getPolType();
655 if ( whichrow < 0 || whichrow >= nrow() )
656 throw(AipsError("Illegal row number."));
657 std::vector<float> out;
658 Vector<Float> arr;
659 uInt requestedpol = polCol_(whichrow);
660 String basetype = getPolType();
661 if ( ptype == basetype ) {
662 specCol_.get(whichrow, arr);
663 } else {
664 STPol* stpol = 0;
665 stpol =STPol::getPolClass(Scantable::factories_, basetype);
666 try {
667 uInt row = uInt(whichrow);
668 stpol->setSpectra(getPolMatrix(row));
669 Float fang,fhand,parang;
670 fang = focusTable_.getTotalFeedAngle(mfocusidCol_(row));
671 fhand = focusTable_.getFeedHand(mfocusidCol_(row));
672 parang = paraCol_(row);
673 /// @todo re-enable this
674 // disable total feed angle to support paralactifying Caswell style
675 stpol->setPhaseCorrections(parang, -parang, fhand);
676 arr = stpol->getSpectrum(requestedpol, ptype);
677 delete stpol;
678 } catch (AipsError& e) {
679 delete stpol;
680 throw(e);
681 }
682 }
683 if ( arr.nelements() == 0 )
684 pushLog("Not enough polarisations present to do the conversion.");
685 arr.tovector(out);
686 return out;
687}
688
689void asap::Scantable::setSpectrum( const std::vector<float>& spec,
690 int whichrow )
691{
692 Vector<Float> spectrum(spec);
693 Vector<Float> arr;
694 specCol_.get(whichrow, arr);
695 if ( spectrum.nelements() != arr.nelements() )
696 throw AipsError("The spectrum has incorrect number of channels.");
697 specCol_.put(whichrow, spectrum);
698}
699
700
701String Scantable::generateName()
702{
703 return (File::newUniqueName("./","temp")).baseName();
704}
705
706const casa::Table& Scantable::table( ) const
707{
708 return table_;
709}
710
711casa::Table& Scantable::table( )
712{
713 return table_;
714}
715
716std::string Scantable::getPolType() const
717{
718 return table_.keywordSet().asString("POLTYPE");
719}
720
721void Scantable::unsetSelection()
722{
723 table_ = originalTable_;
724 attach();
725 selector_.reset();
726}
727
728void Scantable::setSelection( const STSelector& selection )
729{
730 Table tab = const_cast<STSelector&>(selection).apply(originalTable_);
731 if ( tab.nrow() == 0 ) {
732 throw(AipsError("Selection contains no data. Not applying it."));
733 }
734 table_ = tab;
735 attach();
736 selector_ = selection;
737}
738
739std::string Scantable::summary( bool verbose )
740{
741 // Format header info
742 ostringstream oss;
743 oss << endl;
744 oss << asap::SEPERATOR << endl;
745 oss << " Scan Table Summary" << endl;
746 oss << asap::SEPERATOR << endl;
747 oss.flags(std::ios_base::left);
748 oss << setw(15) << "Beams:" << setw(4) << nbeam() << endl
749 << setw(15) << "IFs:" << setw(4) << nif() << endl
750 << setw(15) << "Polarisations:" << setw(4) << npol()
751 << "(" << getPolType() << ")" << endl
752 << setw(15) << "Channels:" << setw(4) << nchan() << endl;
753 oss << endl;
754 String tmp;
755 oss << setw(15) << "Observer:"
756 << table_.keywordSet().asString("Observer") << endl;
757 oss << setw(15) << "Obs Date:" << getTime(-1,true) << endl;
758 table_.keywordSet().get("Project", tmp);
759 oss << setw(15) << "Project:" << tmp << endl;
760 table_.keywordSet().get("Obstype", tmp);
761 oss << setw(15) << "Obs. Type:" << tmp << endl;
762 table_.keywordSet().get("AntennaName", tmp);
763 oss << setw(15) << "Antenna Name:" << tmp << endl;
764 table_.keywordSet().get("FluxUnit", tmp);
765 oss << setw(15) << "Flux Unit:" << tmp << endl;
766 Vector<Double> vec(moleculeTable_.getRestFrequencies());
767 oss << setw(15) << "Rest Freqs:";
768 if (vec.nelements() > 0) {
769 oss << setprecision(10) << vec << " [Hz]" << endl;
770 } else {
771 oss << "none" << endl;
772 }
773
774 oss << setw(15) << "Abcissa:" << getAbcissaLabel(0) << endl;
775 oss << selector_.print() << endl;
776 oss << endl;
777 // main table
778 String dirtype = "Position ("
779 + getDirectionRefString()
780 + ")";
781 oss << setw(5) << "Scan" << setw(15) << "Source"
782 << setw(10) << "Time" << setw(18) << "Integration" << endl;
783 oss << setw(5) << "" << setw(5) << "Beam" << setw(3) << "" << dirtype << endl;
784 oss << setw(10) << "" << setw(3) << "IF" << setw(6) << ""
785 << setw(8) << "Frame" << setw(16)
786 << "RefVal" << setw(10) << "RefPix" << setw(12) << "Increment" <<endl;
787 oss << asap::SEPERATOR << endl;
788 TableIterator iter(table_, "SCANNO");
789 while (!iter.pastEnd()) {
790 Table subt = iter.table();
791 ROTableRow row(subt);
792 MEpoch::ROScalarColumn timeCol(subt,"TIME");
793 const TableRecord& rec = row.get(0);
794 oss << setw(4) << std::right << rec.asuInt("SCANNO")
795 << std::left << setw(1) << ""
796 << setw(15) << rec.asString("SRCNAME")
797 << setw(10) << formatTime(timeCol(0), false);
798 // count the cycles in the scan
799 TableIterator cyciter(subt, "CYCLENO");
800 int nint = 0;
801 while (!cyciter.pastEnd()) {
802 ++nint;
803 ++cyciter;
804 }
805 oss << setw(3) << std::right << nint << setw(3) << " x " << std::left
806 << setw(6) << formatSec(rec.asFloat("INTERVAL")) << endl;
807
808 TableIterator biter(subt, "BEAMNO");
809 while (!biter.pastEnd()) {
810 Table bsubt = biter.table();
811 ROTableRow brow(bsubt);
812 const TableRecord& brec = brow.get(0);
813 uInt row0 = bsubt.rowNumbers(table_)[0];
814 oss << setw(5) << "" << setw(4) << std::right << brec.asuInt("BEAMNO")<< std::left;
815 oss << setw(4) << "" << formatDirection(getDirection(row0)) << endl;
816 TableIterator iiter(bsubt, "IFNO");
817 while (!iiter.pastEnd()) {
818 Table isubt = iiter.table();
819 ROTableRow irow(isubt);
820 const TableRecord& irec = irow.get(0);
821 oss << setw(10) << "";
822 oss << setw(3) << std::right << irec.asuInt("IFNO") << std::left
823 << setw(2) << "" << frequencies().print(irec.asuInt("FREQ_ID"));
824
825 ++iiter;
826 }
827 ++biter;
828 }
829 ++iter;
830 }
831 /// @todo implement verbose mode
832 return String(oss);
833}
834
835std::string Scantable::getTime(int whichrow, bool showdate) const
836{
837 MEpoch::ROScalarColumn timeCol(table_, "TIME");
838 MEpoch me;
839 if (whichrow > -1) {
840 me = timeCol(uInt(whichrow));
841 } else {
842 Double tm;
843 table_.keywordSet().get("UTC",tm);
844 me = MEpoch(MVEpoch(tm));
845 }
846 return formatTime(me, showdate);
847}
848
849std::string Scantable::getDirectionString(int whichrow) const
850{
851 return formatDirection(getDirection(uInt(whichrow)));
852}
853
854std::vector< double > asap::Scantable::getAbcissa( int whichrow ) const
855{
856 if ( whichrow > int(table_.nrow()) ) throw(AipsError("Illegal ro number"));
857 std::vector<double> stlout;
858 int nchan = specCol_(whichrow).nelements();
859 String us = freqTable_.getUnitString();
860 if ( us == "" || us == "pixel" || us == "channel" ) {
861 for (int i=0; i<nchan; ++i) {
862 stlout.push_back(double(i));
863 }
864 return stlout;
865 }
866
867 const MPosition& mp = getAntennaPosition();
868 const MDirection& md = getDirection(whichrow);
869 const MEpoch& me = timeCol_(whichrow);
870 Double rf = moleculeTable_.getRestFrequency(mmolidCol_(whichrow));
871 SpectralCoordinate spc =
872 freqTable_.getSpectralCoordinate(md, mp, me, rf, mfreqidCol_(whichrow));
873 Vector<Double> pixel(nchan);
874 Vector<Double> world;
875 indgen(pixel);
876 if ( Unit(us) == Unit("Hz") ) {
877 for ( int i=0; i < nchan; ++i) {
878 Double world;
879 spc.toWorld(world, pixel[i]);
880 stlout.push_back(double(world));
881 }
882 } else if ( Unit(us) == Unit("km/s") ) {
883 Vector<Double> world;
884 spc.pixelToVelocity(world, pixel);
885 world.tovector(stlout);
886 }
887 return stlout;
888}
889void asap::Scantable::setDirectionRefString( const std::string & refstr )
890{
891 MDirection::Types mdt;
892 if (refstr != "" && !MDirection::getType(mdt, refstr)) {
893 throw(AipsError("Illegal Direction frame."));
894 }
895 if ( refstr == "" ) {
896 String defaultstr = MDirection::showType(dirCol_.getMeasRef().getType());
897 table_.rwKeywordSet().define("DIRECTIONREF", defaultstr);
898 } else {
899 table_.rwKeywordSet().define("DIRECTIONREF", String(refstr));
900 }
901}
902
903std::string asap::Scantable::getDirectionRefString( ) const
904{
905 return table_.keywordSet().asString("DIRECTIONREF");
906}
907
908MDirection Scantable::getDirection(int whichrow ) const
909{
910 String usertype = table_.keywordSet().asString("DIRECTIONREF");
911 String type = MDirection::showType(dirCol_.getMeasRef().getType());
912 if ( usertype != type ) {
913 MDirection::Types mdt;
914 if (!MDirection::getType(mdt, usertype)) {
915 throw(AipsError("Illegal Direction frame."));
916 }
917 return dirCol_.convert(uInt(whichrow), mdt);
918 } else {
919 return dirCol_(uInt(whichrow));
920 }
921}
922
923std::string Scantable::getAbcissaLabel( int whichrow ) const
924{
925 if ( whichrow > int(table_.nrow()) ) throw(AipsError("Illegal ro number"));
926 const MPosition& mp = getAntennaPosition();
927 const MDirection& md = getDirection(whichrow);
928 const MEpoch& me = timeCol_(whichrow);
929 const Double& rf = mmolidCol_(whichrow);
930 SpectralCoordinate spc =
931 freqTable_.getSpectralCoordinate(md, mp, me, rf, mfreqidCol_(whichrow));
932
933 String s = "Channel";
934 Unit u = Unit(freqTable_.getUnitString());
935 if (u == Unit("km/s")) {
936 s = CoordinateUtil::axisLabel(spc, 0, True,True, True);
937 } else if (u == Unit("Hz")) {
938 Vector<String> wau(1);wau = u.getName();
939 spc.setWorldAxisUnits(wau);
940 s = CoordinateUtil::axisLabel(spc, 0, True, True, False);
941 }
942 return s;
943
944}
945
946void asap::Scantable::setRestFrequencies( double rf, const std::string& name,
947 const std::string& unit )
948{
949 ///@todo lookup in line table to fill in name and formattedname
950 Unit u(unit);
951 Quantum<Double> urf(rf, u);
952 uInt id = moleculeTable_.addEntry(urf.getValue("Hz"), name, "");
953 TableVector<uInt> tabvec(table_, "MOLECULE_ID");
954 tabvec = id;
955}
956
957void asap::Scantable::setRestFrequencies( const std::string& name )
958{
959 throw(AipsError("setRestFrequencies( const std::string& name ) NYI"));
960 ///@todo implement
961}
962
963std::vector< unsigned int > asap::Scantable::rownumbers( ) const
964{
965 std::vector<unsigned int> stlout;
966 Vector<uInt> vec = table_.rowNumbers();
967 vec.tovector(stlout);
968 return stlout;
969}
970
971
972Matrix<Float> asap::Scantable::getPolMatrix( uInt whichrow ) const
973{
974 ROTableRow row(table_);
975 const TableRecord& rec = row.get(whichrow);
976 Table t =
977 originalTable_( originalTable_.col("SCANNO") == Int(rec.asuInt("SCANNO"))
978 && originalTable_.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
979 && originalTable_.col("IFNO") == Int(rec.asuInt("IFNO"))
980 && originalTable_.col("CYCLENO") == Int(rec.asuInt("CYCLENO")) );
981 ROArrayColumn<Float> speccol(t, "SPECTRA");
982 return speccol.getColumn();
983}
984
985std::vector< std::string > asap::Scantable::columnNames( ) const
986{
987 Vector<String> vec = table_.tableDesc().columnNames();
988 return mathutil::tovectorstring(vec);
989}
990
991casa::MEpoch::Types asap::Scantable::getTimeReference( ) const
992{
993 return MEpoch::castType(timeCol_.getMeasRef().getType());
994}
995
996void Scantable::addFit( const STFitEntry & fit, int row )
997{
998 cout << mfitidCol_(uInt(row)) << endl;
999 uInt id = fitTable_.addEntry(fit, mfitidCol_(uInt(row)));
1000 mfitidCol_.put(uInt(row), id);
1001}
1002
1003
1004}
1005 //namespace asap
Note: See TracBrowser for help on using the repository browser.