source: trunk/src/Scantable.cpp@ 1533

Last change on this file since 1533 was 1507, checked in by Max Voronkov, 16 years ago

fixed typo in the exception text

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