source: trunk/src/Scantable.cpp@ 1081

Last change on this file since 1081 was 1068, checked in by mar637, 18 years ago

added Scantable::getDirectionString; added Scantable::setSourceType
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 30.0 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
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
340int Scantable::nscan() const {
341 Vector<uInt> scannos(scanCol_.getColumn());
342 uInt nout = GenSort<uInt>::sort( scannos, Sort::Ascending,
343 Sort::QuickSort|Sort::NoDuplicates );
344 return int(nout);
345}
346
347std::string Scantable::formatSec(Double x) const
348{
349 Double xcop = x;
350 MVTime mvt(xcop/24./3600.); // make days
351
352 if (x < 59.95)
353 return String(" ") + mvt.string(MVTime::TIME_CLEAN_NO_HM, 7)+"s";
354 else if (x < 3599.95)
355 return String(" ") + mvt.string(MVTime::TIME_CLEAN_NO_H,7)+" ";
356 else {
357 ostringstream oss;
358 oss << setw(2) << std::right << setprecision(1) << mvt.hour();
359 oss << ":" << mvt.string(MVTime::TIME_CLEAN_NO_H,7) << " ";
360 return String(oss);
361 }
362};
363
364std::string Scantable::formatDirection(const MDirection& md) const
365{
366 Vector<Double> t = md.getAngle(Unit(String("rad"))).getValue();
367 Int prec = 7;
368
369 MVAngle mvLon(t[0]);
370 String sLon = mvLon.string(MVAngle::TIME,prec);
371 uInt tp = md.getRef().getType();
372 if (tp == MDirection::GALACTIC ||
373 tp == MDirection::SUPERGAL ) {
374 sLon = mvLon(0.0).string(MVAngle::ANGLE_CLEAN,prec);
375 }
376 MVAngle mvLat(t[1]);
377 String sLat = mvLat.string(MVAngle::ANGLE+MVAngle::DIG2,prec);
378 return sLon + String(" ") + sLat;
379}
380
381
382std::string Scantable::getFluxUnit() const
383{
384 return table_.keywordSet().asString("FluxUnit");
385}
386
387void Scantable::setFluxUnit(const std::string& unit)
388{
389 String tmp(unit);
390 Unit tU(tmp);
391 if (tU==Unit("K") || tU==Unit("Jy")) {
392 table_.rwKeywordSet().define(String("FluxUnit"), tmp);
393 } else {
394 throw AipsError("Illegal unit - must be compatible with Jy or K");
395 }
396}
397
398void Scantable::setInstrument(const std::string& name)
399{
400 bool throwIt = true;
401 // create an Instrument to see if this is valid
402 STAttr::convertInstrument(name, throwIt);
403 String nameU(name);
404 nameU.upcase();
405 table_.rwKeywordSet().define(String("AntennaName"), nameU);
406}
407
408MPosition Scantable::getAntennaPosition () const
409{
410 Vector<Double> antpos;
411 table_.keywordSet().get("AntennaPosition", antpos);
412 MVPosition mvpos(antpos(0),antpos(1),antpos(2));
413 return MPosition(mvpos);
414}
415
416void Scantable::makePersistent(const std::string& filename)
417{
418 String inname(filename);
419 Path path(inname);
420 inname = path.expandedName();
421 table_.deepCopy(inname, Table::New);
422}
423
424int Scantable::nbeam( int scanno ) const
425{
426 if ( scanno < 0 ) {
427 Int n;
428 table_.keywordSet().get("nBeam",n);
429 return int(n);
430 } else {
431 // take the first POLNO,IFNO,CYCLENO as nbeam shouldn't vary with these
432 Table t = table_(table_.col("SCANNO") == scanno);
433 ROTableRow row(t);
434 const TableRecord& rec = row.get(0);
435 Table subt = t( t.col("IFNO") == Int(rec.asuInt("IFNO"))
436 && t.col("POLNO") == Int(rec.asuInt("POLNO"))
437 && t.col("CYCLENO") == Int(rec.asuInt("CYCLENO")) );
438 ROTableVector<uInt> v(subt, "BEAMNO");
439 return int(v.nelements());
440 }
441 return 0;
442}
443
444int Scantable::nif( int scanno ) const
445{
446 if ( scanno < 0 ) {
447 Int n;
448 table_.keywordSet().get("nIF",n);
449 return int(n);
450 } else {
451 // take the first POLNO,BEAMNO,CYCLENO as nbeam shouldn't vary with these
452 Table t = table_(table_.col("SCANNO") == scanno);
453 ROTableRow row(t);
454 const TableRecord& rec = row.get(0);
455 Table subt = t( t.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
456 && t.col("POLNO") == Int(rec.asuInt("POLNO"))
457 && t.col("CYCLENO") == Int(rec.asuInt("CYCLENO")) );
458 if ( subt.nrow() == 0 ) return 0;
459 ROTableVector<uInt> v(subt, "IFNO");
460 return int(v.nelements());
461 }
462 return 0;
463}
464
465int Scantable::npol( int scanno ) const
466{
467 if ( scanno < 0 ) {
468 Int n;
469 table_.keywordSet().get("nPol",n);
470 return n;
471 } else {
472 // take the first POLNO,IFNO,CYCLENO as nbeam shouldn't vary with these
473 Table t = table_(table_.col("SCANNO") == scanno);
474 ROTableRow row(t);
475 const TableRecord& rec = row.get(0);
476 Table subt = t( t.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
477 && t.col("IFNO") == Int(rec.asuInt("IFNO"))
478 && t.col("CYCLENO") == Int(rec.asuInt("CYCLENO")) );
479 if ( subt.nrow() == 0 ) return 0;
480 ROTableVector<uInt> v(subt, "POLNO");
481 return int(v.nelements());
482 }
483 return 0;
484}
485
486int Scantable::ncycle( int scanno ) const
487{
488 if ( scanno < 0 ) {
489 Block<String> cols(2);
490 cols[0] = "SCANNO";
491 cols[1] = "CYCLENO";
492 TableIterator it(table_, cols);
493 int n = 0;
494 while ( !it.pastEnd() ) {
495 ++n;
496 ++it;
497 }
498 return n;
499 } else {
500 Table t = table_(table_.col("SCANNO") == scanno);
501 ROTableRow row(t);
502 const TableRecord& rec = row.get(0);
503 Table subt = t( t.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
504 && t.col("POLNO") == Int(rec.asuInt("POLNO"))
505 && t.col("IFNO") == Int(rec.asuInt("IFNO")) );
506 if ( subt.nrow() == 0 ) return 0;
507 return int(subt.nrow());
508 }
509 return 0;
510}
511
512
513int Scantable::nrow( int scanno ) const
514{
515 return int(table_.nrow());
516}
517
518int Scantable::nchan( int ifno ) const
519{
520 if ( ifno < 0 ) {
521 Int n;
522 table_.keywordSet().get("nChan",n);
523 return int(n);
524 } else {
525 // take the first SCANNO,POLNO,BEAMNO,CYCLENO as nbeam shouldn't vary with these
526 Table t = table_(table_.col("IFNO") == ifno);
527 if ( t.nrow() == 0 ) return 0;
528 ROArrayColumn<Float> v(t, "SPECTRA");
529 return v.shape(0)(0);
530 }
531 return 0;
532}
533
534int Scantable::getChannels(int whichrow) const
535{
536 return specCol_.shape(whichrow)(0);
537}
538
539int Scantable::getBeam(int whichrow) const
540{
541 return beamCol_(whichrow);
542}
543
544int Scantable::getIF(int whichrow) const
545{
546 return ifCol_(whichrow);
547}
548
549int Scantable::getPol(int whichrow) const
550{
551 return polCol_(whichrow);
552}
553
554std::string Scantable::formatTime(const MEpoch& me, bool showdate) const
555{
556 MVTime mvt(me.getValue());
557 if (showdate)
558 mvt.setFormat(MVTime::YMD);
559 else
560 mvt.setFormat(MVTime::TIME);
561 ostringstream oss;
562 oss << mvt;
563 return String(oss);
564}
565
566void Scantable::calculateAZEL()
567{
568 MPosition mp = getAntennaPosition();
569 MEpoch::ROScalarColumn timeCol(table_, "TIME");
570 ostringstream oss;
571 oss << "Computed azimuth/elevation using " << endl
572 << mp << endl;
573 for (Int i=0; i<nrow(); ++i) {
574 MEpoch me = timeCol(i);
575 MDirection md = getDirection(i);
576 oss << " Time: " << formatTime(me,False) << " Direction: " << formatDirection(md)
577 << endl << " => ";
578 MeasFrame frame(mp, me);
579 Vector<Double> azel =
580 MDirection::Convert(md, MDirection::Ref(MDirection::AZEL,
581 frame)
582 )().getAngle("rad").getValue();
583 azCol_.put(i,Float(azel[0]));
584 elCol_.put(i,Float(azel[1]));
585 oss << "azel: " << azel[0]/C::pi*180.0 << " "
586 << azel[1]/C::pi*180.0 << " (deg)" << endl;
587 }
588 pushLog(String(oss));
589}
590
591void Scantable::flag(const std::vector<bool>& msk)
592{
593 if ( selector_.empty() && msk.size() == 0 )
594 throw(AipsError("Trying to flag whole scantable."));
595 if ( msk.size() == 0 ) {
596 uChar userflag = 1 << 7;
597 for ( uInt i=0; i<table_.nrow(); ++i) {
598 Vector<uChar> flgs = flagsCol_(i);
599 flgs = userflag;
600 flagsCol_.put(i, flgs);
601 }
602 return;
603 }
604 if ( int(msk.size()) != nchan() ) {
605 throw(AipsError("Mask has incorrect number of channels."));
606 }
607 for ( uInt i=0; i<table_.nrow(); ++i) {
608 Vector<uChar> flgs = flagsCol_(i);
609 if ( flgs.nelements() != msk.size() ) {
610 throw(AipsError("Mask has incorrect number of channels."
611 " Probably varying with IF. Please flag per IF"));
612 }
613 std::vector<bool>::const_iterator it;
614 uInt j = 0;
615 uChar userflag = 1 << 7;
616 for (it = msk.begin(); it != msk.end(); ++it) {
617 if ( *it ) {
618 flgs(j) = userflag;
619 }
620 ++j;
621 }
622 flagsCol_.put(i, flgs);
623 }
624}
625
626std::vector<bool> Scantable::getMask(int whichrow) const
627{
628 Vector<uChar> flags;
629 flagsCol_.get(uInt(whichrow), flags);
630 Vector<Bool> bflag(flags.shape());
631 convertArray(bflag, flags);
632 bflag = !bflag;
633 std::vector<bool> mask;
634 bflag.tovector(mask);
635 return mask;
636}
637
638std::vector<float> Scantable::getSpectrum( int whichrow,
639 const std::string& poltype ) const
640{
641 String ptype = poltype;
642 if (poltype == "" ) ptype = getPolType();
643 if ( whichrow < 0 || whichrow >= nrow() )
644 throw(AipsError("Illegal row number."));
645 std::vector<float> out;
646 Vector<Float> arr;
647 uInt requestedpol = polCol_(whichrow);
648 String basetype = getPolType();
649 if ( ptype == basetype ) {
650 specCol_.get(whichrow, arr);
651 } else {
652 STPol* stpol = 0;
653 stpol =STPol::getPolClass(Scantable::factories_, basetype);
654 try {
655 uInt row = uInt(whichrow);
656 stpol->setSpectra(getPolMatrix(row));
657 Float fang,fhand,parang;
658 fang = focusTable_.getTotalFeedAngle(mfocusidCol_(row));
659 fhand = focusTable_.getFeedHand(mfocusidCol_(row));
660 parang = paraCol_(row);
661 /// @todo re-enable this
662 // disable total feed angle to support paralactifying Caswell style
663 stpol->setPhaseCorrections(parang, -parang, fhand);
664 arr = stpol->getSpectrum(requestedpol, ptype);
665 delete stpol;
666 } catch (AipsError& e) {
667 delete stpol;
668 throw(e);
669 }
670 }
671 if ( arr.nelements() == 0 )
672 pushLog("Not enough polarisations present to do the conversion.");
673 arr.tovector(out);
674 return out;
675}
676
677void asap::Scantable::setSpectrum( const std::vector<float>& spec,
678 int whichrow )
679{
680 Vector<Float> spectrum(spec);
681 Vector<Float> arr;
682 specCol_.get(whichrow, arr);
683 if ( spectrum.nelements() != arr.nelements() )
684 throw AipsError("The spectrum has incorrect number of channels.");
685 specCol_.put(whichrow, spectrum);
686}
687
688
689String Scantable::generateName()
690{
691 return (File::newUniqueName("./","temp")).baseName();
692}
693
694const casa::Table& Scantable::table( ) const
695{
696 return table_;
697}
698
699casa::Table& Scantable::table( )
700{
701 return table_;
702}
703
704std::string Scantable::getPolType() const
705{
706 return table_.keywordSet().asString("POLTYPE");
707}
708
709void Scantable::unsetSelection()
710{
711 table_ = originalTable_;
712 attach();
713 selector_.reset();
714}
715
716void Scantable::setSelection( const STSelector& selection )
717{
718 Table tab = const_cast<STSelector&>(selection).apply(originalTable_);
719 if ( tab.nrow() == 0 ) {
720 throw(AipsError("Selection contains no data. Not applying it."));
721 }
722 table_ = tab;
723 attach();
724 selector_ = selection;
725}
726
727std::string Scantable::summary( bool verbose )
728{
729 // Format header info
730 ostringstream oss;
731 oss << endl;
732 oss << asap::SEPERATOR << endl;
733 oss << " Scan Table Summary" << endl;
734 oss << asap::SEPERATOR << endl;
735 oss.flags(std::ios_base::left);
736 oss << setw(15) << "Beams:" << setw(4) << nbeam() << endl
737 << setw(15) << "IFs:" << setw(4) << nif() << endl
738 << setw(15) << "Polarisations:" << setw(4) << npol()
739 << "(" << getPolType() << ")" << endl
740 << setw(15) << "Channels:" << setw(4) << nchan() << endl;
741 oss << endl;
742 String tmp;
743 oss << setw(15) << "Observer:"
744 << table_.keywordSet().asString("Observer") << endl;
745 oss << setw(15) << "Obs Date:" << getTime(-1,true) << endl;
746 table_.keywordSet().get("Project", tmp);
747 oss << setw(15) << "Project:" << tmp << endl;
748 table_.keywordSet().get("Obstype", tmp);
749 oss << setw(15) << "Obs. Type:" << tmp << endl;
750 table_.keywordSet().get("AntennaName", tmp);
751 oss << setw(15) << "Antenna Name:" << tmp << endl;
752 table_.keywordSet().get("FluxUnit", tmp);
753 oss << setw(15) << "Flux Unit:" << tmp << endl;
754 Vector<Double> vec(moleculeTable_.getRestFrequencies());
755 oss << setw(15) << "Rest Freqs:";
756 if (vec.nelements() > 0) {
757 oss << setprecision(10) << vec << " [Hz]" << endl;
758 } else {
759 oss << "none" << endl;
760 }
761
762 oss << setw(15) << "Abcissa:" << getAbcissaLabel(0) << endl;
763 oss << selector_.print() << endl;
764 oss << endl;
765 // main table
766 String dirtype = "Position ("
767 + getDirectionRefString()
768 + ")";
769 oss << setw(5) << "Scan" << setw(15) << "Source"
770 << setw(10) << "Time" << setw(18) << "Integration" << endl;
771 oss << setw(5) << "" << setw(5) << "Beam" << setw(3) << "" << dirtype << endl;
772 oss << setw(10) << "" << setw(3) << "IF" << setw(6) << ""
773 << setw(8) << "Frame" << setw(16)
774 << "RefVal" << setw(10) << "RefPix" << setw(12) << "Increment" <<endl;
775 oss << asap::SEPERATOR << endl;
776 TableIterator iter(table_, "SCANNO");
777 while (!iter.pastEnd()) {
778 Table subt = iter.table();
779 ROTableRow row(subt);
780 MEpoch::ROScalarColumn timeCol(subt,"TIME");
781 const TableRecord& rec = row.get(0);
782 oss << setw(4) << std::right << rec.asuInt("SCANNO")
783 << std::left << setw(1) << ""
784 << setw(15) << rec.asString("SRCNAME")
785 << setw(10) << formatTime(timeCol(0), false);
786 // count the cycles in the scan
787 TableIterator cyciter(subt, "CYCLENO");
788 int nint = 0;
789 while (!cyciter.pastEnd()) {
790 ++nint;
791 ++cyciter;
792 }
793 oss << setw(3) << std::right << nint << setw(3) << " x " << std::left
794 << setw(6) << formatSec(rec.asFloat("INTERVAL")) << endl;
795
796 TableIterator biter(subt, "BEAMNO");
797 while (!biter.pastEnd()) {
798 Table bsubt = biter.table();
799 ROTableRow brow(bsubt);
800 const TableRecord& brec = brow.get(0);
801 uInt row0 = bsubt.rowNumbers(table_)[0];
802 oss << setw(5) << "" << setw(4) << std::right << brec.asuInt("BEAMNO")<< std::left;
803 oss << setw(4) << "" << formatDirection(getDirection(row0)) << endl;
804 TableIterator iiter(bsubt, "IFNO");
805 while (!iiter.pastEnd()) {
806 Table isubt = iiter.table();
807 ROTableRow irow(isubt);
808 const TableRecord& irec = irow.get(0);
809 oss << setw(10) << "";
810 oss << setw(3) << std::right << irec.asuInt("IFNO") << std::left
811 << setw(2) << "" << frequencies().print(irec.asuInt("FREQ_ID"));
812
813 ++iiter;
814 }
815 ++biter;
816 }
817 ++iter;
818 }
819 /// @todo implement verbose mode
820 return String(oss);
821}
822
823std::string Scantable::getTime(int whichrow, bool showdate) const
824{
825 MEpoch::ROScalarColumn timeCol(table_, "TIME");
826 MEpoch me;
827 if (whichrow > -1) {
828 me = timeCol(uInt(whichrow));
829 } else {
830 Double tm;
831 table_.keywordSet().get("UTC",tm);
832 me = MEpoch(MVEpoch(tm));
833 }
834 return formatTime(me, showdate);
835}
836
837std::string Scantable::getDirectionString(int whichrow) const
838{
839 return formatDirection(getDirection(uInt(whichrow)));
840}
841
842std::vector< double > asap::Scantable::getAbcissa( int whichrow ) const
843{
844 if ( whichrow > int(table_.nrow()) ) throw(AipsError("Illegal ro number"));
845 std::vector<double> stlout;
846 int nchan = specCol_(whichrow).nelements();
847 String us = freqTable_.getUnitString();
848 if ( us == "" || us == "pixel" || us == "channel" ) {
849 for (int i=0; i<nchan; ++i) {
850 stlout.push_back(double(i));
851 }
852 return stlout;
853 }
854
855 const MPosition& mp = getAntennaPosition();
856 const MDirection& md = getDirection(whichrow);
857 const MEpoch& me = timeCol_(whichrow);
858 Double rf = moleculeTable_.getRestFrequency(mmolidCol_(whichrow));
859 SpectralCoordinate spc =
860 freqTable_.getSpectralCoordinate(md, mp, me, rf, mfreqidCol_(whichrow));
861 Vector<Double> pixel(nchan);
862 Vector<Double> world;
863 indgen(pixel);
864 if ( Unit(us) == Unit("Hz") ) {
865 for ( int i=0; i < nchan; ++i) {
866 Double world;
867 spc.toWorld(world, pixel[i]);
868 stlout.push_back(double(world));
869 }
870 } else if ( Unit(us) == Unit("km/s") ) {
871 Vector<Double> world;
872 spc.pixelToVelocity(world, pixel);
873 world.tovector(stlout);
874 }
875 return stlout;
876}
877void asap::Scantable::setDirectionRefString( const std::string & refstr )
878{
879 MDirection::Types mdt;
880 if (refstr != "" && !MDirection::getType(mdt, refstr)) {
881 throw(AipsError("Illegal Direction frame."));
882 }
883 if ( refstr == "" ) {
884 String defaultstr = MDirection::showType(dirCol_.getMeasRef().getType());
885 table_.rwKeywordSet().define("DIRECTIONREF", defaultstr);
886 } else {
887 table_.rwKeywordSet().define("DIRECTIONREF", String(refstr));
888 }
889}
890
891std::string asap::Scantable::getDirectionRefString( ) const
892{
893 return table_.keywordSet().asString("DIRECTIONREF");
894}
895
896MDirection Scantable::getDirection(int whichrow ) const
897{
898 String usertype = table_.keywordSet().asString("DIRECTIONREF");
899 String type = MDirection::showType(dirCol_.getMeasRef().getType());
900 if ( usertype != type ) {
901 MDirection::Types mdt;
902 if (!MDirection::getType(mdt, usertype)) {
903 throw(AipsError("Illegal Direction frame."));
904 }
905 return dirCol_.convert(uInt(whichrow), mdt);
906 } else {
907 return dirCol_(uInt(whichrow));
908 }
909}
910
911std::string Scantable::getAbcissaLabel( int whichrow ) const
912{
913 if ( whichrow > int(table_.nrow()) ) throw(AipsError("Illegal ro number"));
914 const MPosition& mp = getAntennaPosition();
915 const MDirection& md = getDirection(whichrow);
916 const MEpoch& me = timeCol_(whichrow);
917 const Double& rf = mmolidCol_(whichrow);
918 SpectralCoordinate spc =
919 freqTable_.getSpectralCoordinate(md, mp, me, rf, mfreqidCol_(whichrow));
920
921 String s = "Channel";
922 Unit u = Unit(freqTable_.getUnitString());
923 if (u == Unit("km/s")) {
924 s = CoordinateUtil::axisLabel(spc,0,True,True,True);
925 } else if (u == Unit("Hz")) {
926 Vector<String> wau(1);wau = u.getName();
927 spc.setWorldAxisUnits(wau);
928 s = CoordinateUtil::axisLabel(spc,0,True,True,False);
929 }
930 return s;
931
932}
933
934void asap::Scantable::setRestFrequencies( double rf, const std::string& unit )
935{
936 ///@todo lookup in line table to fill in name and formattedname
937 Unit u(unit);
938 Quantum<Double> urf(rf, u);
939 uInt id = moleculeTable_.addEntry(urf.getValue("Hz"), "", "");
940 TableVector<uInt> tabvec(table_, "MOLECULE_ID");
941 tabvec = id;
942}
943
944void asap::Scantable::setRestFrequencies( const std::string& name )
945{
946 throw(AipsError("setRestFrequencies( const std::string& name ) NYI"));
947 ///@todo implement
948}
949
950std::vector< unsigned int > asap::Scantable::rownumbers( ) const
951{
952 std::vector<unsigned int> stlout;
953 Vector<uInt> vec = table_.rowNumbers();
954 vec.tovector(stlout);
955 return stlout;
956}
957
958
959Matrix<Float> asap::Scantable::getPolMatrix( uInt whichrow ) const
960{
961 ROTableRow row(table_);
962 const TableRecord& rec = row.get(whichrow);
963 Table t =
964 originalTable_( originalTable_.col("SCANNO") == Int(rec.asuInt("SCANNO"))
965 && originalTable_.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
966 && originalTable_.col("IFNO") == Int(rec.asuInt("IFNO"))
967 && originalTable_.col("CYCLENO") == Int(rec.asuInt("CYCLENO")) );
968 ROArrayColumn<Float> speccol(t, "SPECTRA");
969 return speccol.getColumn();
970}
971
972std::vector< std::string > asap::Scantable::columnNames( ) const
973{
974 Vector<String> vec = table_.tableDesc().columnNames();
975 return mathutil::tovectorstring(vec);
976}
977
978casa::MEpoch::Types asap::Scantable::getTimeReference( ) const
979{
980 return MEpoch::castType(timeCol_.getMeasRef().getType());
981}
982
983void asap::Scantable::addFit( const STFitEntry & fit, int row )
984{
985 cout << mfitidCol_(uInt(row)) << endl;
986 uInt id = fitTable_.addEntry(fit, mfitidCol_(uInt(row)));
987 mfitidCol_.put(uInt(row), id);
988}
989
990
991}
992 //namespace asap
Note: See TracBrowser for help on using the repository browser.