source: trunk/src/Scantable.cpp@ 1141

Last change on this file since 1141 was 1111, checked in by mar637, 18 years ago

add getNumbers utility function, which is now called by all function which retrieve numbers from the Scantable

  • 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<uInt>::sort( 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 GenSort<uInt>::sort( nos, Sort::Ascending,
550 Sort::QuickSort|Sort::NoDuplicates );
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& unit )
947{
948 ///@todo lookup in line table to fill in name and formattedname
949 Unit u(unit);
950 Quantum<Double> urf(rf, u);
951 uInt id = moleculeTable_.addEntry(urf.getValue("Hz"), "", "");
952 TableVector<uInt> tabvec(table_, "MOLECULE_ID");
953 tabvec = id;
954}
955
956void asap::Scantable::setRestFrequencies( const std::string& name )
957{
958 throw(AipsError("setRestFrequencies( const std::string& name ) NYI"));
959 ///@todo implement
960}
961
962std::vector< unsigned int > asap::Scantable::rownumbers( ) const
963{
964 std::vector<unsigned int> stlout;
965 Vector<uInt> vec = table_.rowNumbers();
966 vec.tovector(stlout);
967 return stlout;
968}
969
970
971Matrix<Float> asap::Scantable::getPolMatrix( uInt whichrow ) const
972{
973 ROTableRow row(table_);
974 const TableRecord& rec = row.get(whichrow);
975 Table t =
976 originalTable_( originalTable_.col("SCANNO") == Int(rec.asuInt("SCANNO"))
977 && originalTable_.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
978 && originalTable_.col("IFNO") == Int(rec.asuInt("IFNO"))
979 && originalTable_.col("CYCLENO") == Int(rec.asuInt("CYCLENO")) );
980 ROArrayColumn<Float> speccol(t, "SPECTRA");
981 return speccol.getColumn();
982}
983
984std::vector< std::string > asap::Scantable::columnNames( ) const
985{
986 Vector<String> vec = table_.tableDesc().columnNames();
987 return mathutil::tovectorstring(vec);
988}
989
990casa::MEpoch::Types asap::Scantable::getTimeReference( ) const
991{
992 return MEpoch::castType(timeCol_.getMeasRef().getType());
993}
994
995void Scantable::addFit( const STFitEntry & fit, int row )
996{
997 cout << mfitidCol_(uInt(row)) << endl;
998 uInt id = fitTable_.addEntry(fit, mfitidCol_(uInt(row)));
999 mfitidCol_.put(uInt(row), id);
1000}
1001
1002
1003}
1004 //namespace asap
Note: See TracBrowser for help on using the repository browser.