source: trunk/src/Scantable.cpp@ 1369

Last change on this file since 1369 was 1360, checked in by mar637, 18 years ago

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