source: trunk/src/Scantable.cpp@ 1008

Last change on this file since 1008 was 1000, checked in by mar637, 19 years ago

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