source: trunk/src/Scantable.cpp@ 1799

Last change on this file since 1799 was 1743, checked in by Malte Marquarding, 14 years ago

Cosmetics

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