Changeset 805 for trunk/src/Scantable.cpp
- Timestamp:
- 02/16/06 12:02:18 (18 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/src/Scantable.cpp
r804 r805 1 //#--------------------------------------------------------------------------- 2 //# SDMemTable.cc: A MemoryTable container for single dish integrations 3 //#--------------------------------------------------------------------------- 4 //# Copyright (C) 2004 5 //# ATNF 6 //# 7 //# This program is free software; you can redistribute it and/or modify it 8 //# under the terms of the GNU General Public License as published by the Free 9 //# Software Foundation; either version 2 of the License, or (at your option) 10 //# any later version. 11 //# 12 //# This program is distributed in the hope that it will be useful, but 13 //# WITHOUT ANY WARRANTY; without even the implied warranty of 14 //# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General 15 //# Public License for more details. 16 //# 17 //# You should have received a copy of the GNU General Public License along 18 //# with this program; if not, write to the Free Software Foundation, Inc., 19 //# 675 Massachusetts Ave, Cambridge, MA 02139, USA. 20 //# 21 //# Correspondence concerning this software should be addressed as follows: 22 //# Internet email: Malte.Marquarding@csiro.au 23 //# Postal address: Malte Marquarding, 24 //# Australia Telescope National Facility, 25 //# P.O. Box 76, 26 //# Epping, NSW, 2121, 27 //# AUSTRALIA 28 //# 29 //# $Id: 30 //#--------------------------------------------------------------------------- 31 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 // 32 12 #include <map> 33 13 … … 35 15 #include <casa/iostream.h> 36 16 #include <casa/iomanip.h> 17 #include <casa/OS/Path.h> 18 #include <casa/OS/File.h> 37 19 #include <casa/Arrays/Array.h> 38 20 #include <casa/Arrays/ArrayMath.h> … … 45 27 #include <casa/BasicSL/Constants.h> 46 28 #include <casa/Quanta/MVAngle.h> 29 #include <casa/Containers/RecordField.h> 47 30 48 31 #include <tables/Tables/TableParse.h> … … 52 35 #include <tables/Tables/ScaColDesc.h> 53 36 #include <tables/Tables/ArrColDesc.h> 37 #include <tables/Tables/TableRow.h> 38 #include <tables/Tables/TableVector.h> 39 #include <tables/Tables/TableIter.h> 54 40 55 41 #include <tables/Tables/ExprNode.h> 56 42 #include <tables/Tables/TableRecord.h> 57 43 #include <measures/Measures/MFrequency.h> 44 #include <measures/Measures/MEpoch.h> 58 45 #include <measures/Measures/MeasTable.h> 46 #include <measures/Measures/MeasRef.h> 47 #include <measures/TableMeasures/TableMeasRefDesc.h> 48 #include <measures/TableMeasures/TableMeasValueDesc.h> 49 #include <measures/TableMeasures/TableMeasDesc.h> 50 #include <measures/TableMeasures/ScalarMeasColumn.h> 59 51 #include <coordinates/Coordinates/CoordinateUtil.h> 60 52 #include <casa/Quanta/MVTime.h> 61 53 #include <casa/Quanta/MVAngle.h> 62 54 63 #include "SDDefs.h" 64 #include "SDContainer.h" 65 #include "MathUtils.h" 66 #include "SDPol.h" 55 #include "Scantable.h" 67 56 #include "SDAttr.h" 68 57 69 #include "SDMemTable.h"70 71 58 using namespace casa; 72 using namespace asap; 73 74 SDMemTable::SDMemTable() : 75 IFSel_(0), 76 beamSel_(0), 77 polSel_(0) 78 { 79 setup(); 59 60 namespace asap { 61 62 Scantable::Scantable(Table::TableType ttype) : 63 type_(ttype), 64 freqTable_(ttype), 65 focusTable_(ttype), 66 weatherTable_(ttype), 67 tcalTable_(ttype), 68 moleculeTable_(ttype) 69 { 70 setupMainTable(); 71 table_.rwKeywordSet().defineTable("FREQUENCIES", freqTable_.table()); 72 table_.rwKeywordSet().defineTable("WEATHER", weatherTable_.table()); 73 table_.rwKeywordSet().defineTable("FOCUS", focusTable_.table()); 74 table_.rwKeywordSet().defineTable("TCAL", tcalTable_.table()); 75 table_.rwKeywordSet().defineTable("MOLECULES", moleculeTable_.table()); 76 setupHistoryTable(); 77 setupFitTable(); 78 79 originalTable_ = table_; 80 80 attach(); 81 81 } 82 82 83 SDMemTable::SDMemTable(const std::string& name) : 84 IFSel_(0), 85 beamSel_(0), 86 polSel_(0) 83 Scantable::Scantable(const std::string& name, Table::TableType ttype) : 84 type_(ttype), 85 freqTable_(ttype) 87 86 { 88 87 Table tab(name); … … 92 91 throw(AipsError("Unsupported version of ASAP file.")); 93 92 } 94 table_ = tab.copyToMemoryTable("dummy"); 93 if ( type_ == Table::Memory ) 94 table_ = tab.copyToMemoryTable("dummy"); 95 else 96 table_ = tab; 97 98 originalTable_ = table_; 95 99 attach(); 96 100 } 97 101 98 SDMemTable::SDMemTable(const SDMemTable& other, Bool clear) 99 { 100 table_ = other.table_.copyToMemoryTable(String("dummy")); 101 // clear all rows()102 103 Scantable::Scantable( const Scantable& other, bool clear ) 104 { 105 // with or without data 102 106 if (clear) { 103 table_.removeRow(this->table_.rowNumbers()); 104 IFSel_= 0; 105 beamSel_= 0; 106 polSel_= 0; 107 } else { 108 IFSel_= other.IFSel_; 109 beamSel_= other.beamSel_; 110 polSel_= other.polSel_; 111 } 112 107 table_ = TableCopy::makeEmptyMemoryTable(String("dummy"), other.table_, 108 True); 109 } else { 110 table_ = other.table_.copyToMemoryTable(String("dummy")); 111 } 112 113 originalTable_ = table_; 113 114 attach(); 114 115 } 115 116 116 SDMemTable::SDMemTable(const Table& tab, const std::string& exprs) : 117 IFSel_(0), 118 beamSel_(0), 119 polSel_(0) 120 { 121 Table t = tableCommand(exprs,tab); 122 if (t.nrow() == 0) 123 throw(AipsError("Query unsuccessful.")); 124 table_ = t.copyToMemoryTable("dummy"); 125 attach(); 126 renumber(); 127 } 128 129 SDMemTable::~SDMemTable() 130 { 131 //cerr << "goodbye from SDMemTable @ " << this << endl; 132 } 133 134 SDMemTable SDMemTable::getScan(Int scanID) const 135 { 136 String cond("SELECT * from $1 WHERE SCANID == "); 137 cond += String::toString(scanID); 138 return SDMemTable(table_, cond); 139 } 140 141 SDMemTable &SDMemTable::operator=(const SDMemTable& other) 142 { 143 if (this != &other) { 144 IFSel_= other.IFSel_; 145 beamSel_= other.beamSel_; 146 polSel_= other.polSel_; 147 table_ = other.table_.copyToMemoryTable(String("dummy")); 148 attach(); 149 } 150 return *this; 151 } 152 153 SDMemTable SDMemTable::getSource(const std::string& source) const 154 { 155 String cond("SELECT * from $1 WHERE SRCNAME == "); 156 cond += source; 157 return SDMemTable(table_, cond); 158 } 159 160 void SDMemTable::setup() 117 Scantable::~Scantable() 118 { 119 cout << "~Scantable() " << this << endl; 120 } 121 122 void Scantable::setupMainTable() 161 123 { 162 124 TableDesc td("", "1", TableDesc::Scratch); 163 td.comment() = "A SDMemTable";125 td.comment() = "An ASAP Scantable"; 164 126 td.rwKeywordSet().define("VERSION", Int(version_)); 165 127 128 // n Cycles 129 td.addColumn(ScalarColumnDesc<uInt>("SCANNO")); 130 // new index every nBeam x nIF x nPol 131 td.addColumn(ScalarColumnDesc<uInt>("CYCLENO")); 132 133 td.addColumn(ScalarColumnDesc<uInt>("BEAMNO")); 134 td.addColumn(ScalarColumnDesc<uInt>("IFNO")); 135 td.rwKeywordSet().define("POLTYPE", String("linear")); 136 td.addColumn(ScalarColumnDesc<uInt>("POLNO")); 137 138 td.addColumn(ScalarColumnDesc<uInt>("FREQ_ID")); 139 td.addColumn(ScalarColumnDesc<uInt>("MOLECULE_ID")); 140 // linear, circular, stokes [I Q U V], stokes1 [I Plinear Pangle V] 141 td.addColumn(ScalarColumnDesc<Int>("REFBEAMNO")); 142 166 143 td.addColumn(ScalarColumnDesc<Double>("TIME")); 144 TableMeasRefDesc measRef(MEpoch::UTC); // UTC as default 145 TableMeasValueDesc measVal(td, "TIME"); 146 TableMeasDesc<MEpoch> mepochCol(measVal, measRef); 147 mepochCol.write(td); 148 149 td.addColumn(ScalarColumnDesc<Double>("INTERVAL")); 150 167 151 td.addColumn(ScalarColumnDesc<String>("SRCNAME")); 152 // Type of source (on=0, off=1, other=-1) 153 td.addColumn(ScalarColumnDesc<Int>("SRCTYPE", Int(-1))); 154 td.addColumn(ScalarColumnDesc<String>("FIELDNAME")); 155 156 //The actual Data Vectors 168 157 td.addColumn(ArrayColumnDesc<Float>("SPECTRA")); 169 158 td.addColumn(ArrayColumnDesc<uChar>("FLAGTRA")); 170 159 td.addColumn(ArrayColumnDesc<Float>("TSYS")); 171 td.addColumn(ArrayColumnDesc<Float>("STOKES")); 172 td.addColumn(ScalarColumnDesc<Int>("SCANID")); 173 td.addColumn(ScalarColumnDesc<Double>("INTERVAL")); 174 td.addColumn(ArrayColumnDesc<uInt>("FREQID")); 175 td.addColumn(ArrayColumnDesc<uInt>("RESTFREQID")); 176 td.addColumn(ArrayColumnDesc<Double>("DIRECTION")); 177 td.addColumn(ScalarColumnDesc<String>("FIELDNAME")); 178 td.addColumn(ScalarColumnDesc<String>("TCALTIME")); 179 td.addColumn(ArrayColumnDesc<Float>("TCAL")); 180 td.addColumn(ScalarColumnDesc<Float>("AZIMUTH")); 181 td.addColumn(ScalarColumnDesc<Float>("ELEVATION")); 160 161 td.addColumn(ArrayColumnDesc<Double>("DIRECTION", 162 IPosition(1,2), 163 ColumnDesc::Direct)); 164 TableMeasRefDesc mdirRef(MDirection::J2000); // default 165 TableMeasValueDesc tmvdMDir(td, "DIRECTION"); 166 // the TableMeasDesc gives the column a type 167 TableMeasDesc<MDirection> mdirCol(tmvdMDir, mdirRef); 168 // writing create the measure column 169 mdirCol.write(td); 170 td.addColumn(ScalarColumnDesc<Double>("AZIMUTH")); 171 td.addColumn(ScalarColumnDesc<Double>("ELEVATION")); 182 172 td.addColumn(ScalarColumnDesc<Float>("PARANGLE")); 183 td.addColumn(ScalarColumnDesc<Int>("REFBEAM")); 184 td.addColumn(ArrayColumnDesc<Int>("FITID")); 173 174 td.addColumn(ScalarColumnDesc<uInt>("TCAL_ID")); 175 td.addColumn(ScalarColumnDesc<uInt>("FIT_ID")); 176 177 td.addColumn(ScalarColumnDesc<uInt>("FOCUS_ID")); 178 td.addColumn(ScalarColumnDesc<uInt>("WEATHER_ID")); 179 180 td.rwKeywordSet().define("OBSMODE", String("")); 185 181 186 182 // Now create Table SetUp from the description. 187 183 SetupNewTable aNewTab("dummy", td, Table::New); 188 189 // Bind the Stokes Virtual machine to the STOKES column Because we190 // don't know how many polarizations will be in the data at this191 // point, we must bind the Virtual Engine regardless. The STOKES192 // column won't be accessed if not appropriate (nPol=4)193 SDStokesEngine::registerClass();194 SDStokesEngine stokesEngine(String("STOKES"), String("SPECTRA"));195 aNewTab.bindColumn("STOKES", stokesEngine);196 197 // Create Table198 184 table_ = Table(aNewTab, Table::Memory, 0); 199 // add subtable 200 TableDesc tdf("", "1", TableDesc::Scratch); 201 tdf.addColumn(ArrayColumnDesc<String>("FUNCTIONS")); 202 tdf.addColumn(ArrayColumnDesc<Int>("COMPONENTS")); 203 tdf.addColumn(ArrayColumnDesc<Double>("PARAMETERS")); 204 tdf.addColumn(ArrayColumnDesc<Bool>("PARMASK")); 205 tdf.addColumn(ArrayColumnDesc<String>("FRAMEINFO")); 206 SetupNewTable fittab("fits", tdf, Table::New); 207 Table fitTable(fittab, Table::Memory); 208 table_.rwKeywordSet().defineTable("FITS", fitTable); 209 185 186 originalTable_ = table_; 187 188 } 189 190 void asap::Scantable::setupHistoryTable( ) 191 { 210 192 TableDesc tdh("", "1", TableDesc::Scratch); 211 193 tdh.addColumn(ScalarColumnDesc<String>("ITEM")); 212 SetupNewTable histtab("hist ", tdh, Table::New);194 SetupNewTable histtab("history", tdh, Table::New); 213 195 Table histTable(histtab, Table::Memory); 214 196 table_.rwKeywordSet().defineTable("HISTORY", histTable); 215 197 } 216 198 217 void SDMemTable::attach() 218 { 219 timeCol_.attach(table_, "TIME"); 220 srcnCol_.attach(table_, "SRCNAME"); 221 specCol_.attach(table_, "SPECTRA"); 222 flagsCol_.attach(table_, "FLAGTRA"); 223 tsCol_.attach(table_, "TSYS"); 224 stokesCol_.attach(table_, "STOKES"); 225 scanCol_.attach(table_, "SCANID"); 226 integrCol_.attach(table_, "INTERVAL"); 227 freqidCol_.attach(table_, "FREQID"); 228 restfreqidCol_.attach(table_, "RESTFREQID"); 229 dirCol_.attach(table_, "DIRECTION"); 230 fldnCol_.attach(table_, "FIELDNAME"); 231 tcaltCol_.attach(table_, "TCALTIME"); 232 tcalCol_.attach(table_, "TCAL"); 233 azCol_.attach(table_, "AZIMUTH"); 234 elCol_.attach(table_, "ELEVATION"); 235 paraCol_.attach(table_, "PARANGLE"); 236 rbeamCol_.attach(table_, "REFBEAM"); 237 fitCol_.attach(table_,"FITID"); 238 } 239 240 241 std::string SDMemTable::getSourceName(Int whichRow) const 242 { 243 String name; 244 srcnCol_.get(whichRow, name); 245 return name; 246 } 247 248 float SDMemTable::getElevation(Int whichRow) const 249 { 250 float elevation; 251 elCol_.get(whichRow, elevation); 252 return elevation; 253 } 254 255 float SDMemTable::getAzimuth(Int whichRow) const 256 { 257 float azimuth; 258 azCol_.get(whichRow, azimuth); 259 return azimuth; 260 } 261 262 float SDMemTable::getParAngle(Int whichRow) const 263 { 264 float parangle; 265 paraCol_.get(whichRow, parangle); 266 return parangle; 267 } 268 269 std::string SDMemTable::getTime(Int whichRow, Bool showDate) const 270 { 271 Double tm; 272 if (whichRow > -1) { 273 timeCol_.get(whichRow, tm); 274 } else { 275 table_.keywordSet().get("UTC",tm); 276 } 277 MVTime mvt(tm); 278 if (showDate) 279 mvt.setFormat(MVTime::YMD); 280 else 281 mvt.setFormat(MVTime::TIME); 282 ostringstream oss; 283 oss << mvt; 284 return String(oss); 285 } 286 287 double SDMemTable::getInterval(Int whichRow) const 288 { 289 Double intval; 290 integrCol_.get(whichRow, intval); 291 return intval; 292 } 293 294 bool SDMemTable::setIF(Int whichIF) 295 { 296 if ( whichIF >= 0 && whichIF < nIF()) { 297 IFSel_ = whichIF; 298 return true; 299 } 300 return false; 301 } 302 303 bool SDMemTable::setBeam(Int whichBeam) 304 { 305 if ( whichBeam >= 0 && whichBeam < nBeam()) { 306 beamSel_ = whichBeam; 307 return true; 308 } 309 return false; 310 } 311 312 bool SDMemTable::setPol(Int whichPol) 313 { 314 if ( whichPol >= 0 && whichPol < nPol()) { 315 polSel_ = whichPol; 316 return true; 317 } 318 return false; 319 } 320 321 void SDMemTable::resetCursor() 322 { 323 polSel_ = 0; 324 IFSel_ = 0; 325 beamSel_ = 0; 326 } 327 328 std::vector<bool> SDMemTable::getMask(Int whichRow) const 329 { 330 331 std::vector<bool> mask; 332 333 Array<uChar> arr; 334 flagsCol_.get(whichRow, arr); 335 336 ArrayAccessor<uChar, Axis<asap::BeamAxis> > aa0(arr); 337 aa0.reset(aa0.begin(uInt(beamSel_)));//go to beam 338 ArrayAccessor<uChar, Axis<asap::IFAxis> > aa1(aa0); 339 aa1.reset(aa1.begin(uInt(IFSel_)));// go to IF 340 ArrayAccessor<uChar, Axis<asap::PolAxis> > aa2(aa1); 341 aa2.reset(aa2.begin(uInt(polSel_)));// go to pol 342 343 for (ArrayAccessor<uChar, Axis<asap::ChanAxis> > i(aa2); i != i.end(); ++i) { 344 bool out =!static_cast<bool>(*i); 345 mask.push_back(out); 346 } 347 return mask; 348 } 349 350 351 352 std::vector<float> SDMemTable::getSpectrum(Int whichRow) const 353 { 354 Array<Float> arr; 355 specCol_.get(whichRow, arr); 356 return getFloatSpectrum(arr); 357 } 358 359 360 int SDMemTable::nStokes() const 361 { 362 return stokesCol_.shape(0).nelements(); // All rows same shape 363 } 364 365 366 std::vector<float> SDMemTable::getStokesSpectrum(Int whichRow, 367 Bool doPol) const 368 // 369 // Gets one STokes parameter depending on cursor polSel location 370 // doPol=False : I,Q,U,V 371 // doPol=True : I,P,PA,V ; P = sqrt(Q**2+U**2), PA = 0.5*atan2(Q,U) 372 // 373 { 374 AlwaysAssert(asap::nAxes==4,AipsError); 375 if (nPol()!=1 && nPol()!=2 && nPol()!=4) { 376 throw (AipsError("You must have 1,2 or 4 polarizations to get the Stokes parameters")); 377 } 378 379 // For full conversion we are only supporting linears at the moment 380 381 if (nPol() > 2) { 382 String antName; 383 table_.keywordSet().get("AntennaName", antName); 384 Instrument inst = SDAttr::convertInstrument (antName, True); 385 SDAttr sdAtt; 386 if (sdAtt.feedPolType(inst) != LINEAR) { 387 throw(AipsError("Only linear polarizations are supported")); 388 } 389 } 390 391 Array<Float> arr; 392 stokesCol_.get(whichRow, arr); 393 394 if (doPol && (polSel_==1 || polSel_==2)) { // Q,U --> P, P.A. 395 396 // Set current cursor location 397 398 const IPosition& shape = arr.shape(); 399 IPosition start, end; 400 getCursorSlice(start, end, shape); 401 402 // Get Q and U slices 403 404 Array<Float> Q = SDPolUtil::getStokesSlice(arr,start,end,"Q"); 405 Array<Float> U = SDPolUtil::getStokesSlice(arr,start,end,"U"); 406 407 // Compute output 408 409 Array<Float> out; 410 if (polSel_==1) { // P 411 out = SDPolUtil::polarizedIntensity(Q,U); 412 } else if (polSel_==2) { // P.A. 413 out = SDPolUtil::positionAngle(Q,U); 414 } 415 416 // Copy to output 417 418 IPosition vecShape(1,shape(asap::ChanAxis)); 419 Vector<Float> outV = out.reform(vecShape); 420 std::vector<float> stlout; 421 outV.tovector(stlout); 422 return stlout; 423 424 } else { 425 // Selects at the cursor location 426 return getFloatSpectrum(arr); 427 } 428 } 429 430 std::string SDMemTable::getPolarizationLabel(Bool linear, Bool stokes, 431 Bool linPol, Int polIdx) const 432 { 433 uInt idx = polSel_; 434 if (polIdx >=0) idx = polIdx; 435 return SDPolUtil::polarizationLabel(idx, linear, stokes, linPol); 436 } 437 438 439 440 std::vector<float> SDMemTable::stokesToPolSpectrum(Int whichRow, 441 Bool toLinear, 442 Int polIdx) const 443 // 444 // polIdx 445 // 0:3 -> RR,LL,Real(RL),Imag(RL) 446 // XX,YY,Real(XY),Image(XY) 447 // 448 // Gets only 449 // RR = I + V 450 // LL = I - V 451 // at the moment 452 // 453 { 454 AlwaysAssert(asap::nAxes==4,AipsError); 455 if (nStokes()!=4) { 456 throw (AipsError("You must have 4 Stokes to convert to linear or circular")); 457 } 458 // 459 Array<Float> arr, out; 460 stokesCol_.get(whichRow, arr); 461 462 // Set current cursor location 463 464 const IPosition& shape = arr.shape(); 465 IPosition start, end; 466 getCursorSlice(start, end, shape); 467 468 // Get the slice 469 470 if (toLinear) { 471 throw(AipsError("Conversion to linears not yet supported")); 472 } else { 473 uInt selection = polSel_; 474 if (polIdx > -1) selection = polIdx; 475 Bool doRR = (selection==0); 476 if (selection>1) { 477 throw(AipsError("Only conversion to RR & LL is currently supported")); 478 } 479 480 // Get I and V slices 481 Array<Float> I = SDPolUtil::getStokesSlice(arr,start,end,"I"); 482 Array<Float> V = SDPolUtil::getStokesSlice(arr,start,end,"V"); 483 484 // Compute output 485 out = SDPolUtil::circularPolarizationFromStokes(I, V, doRR); 486 } 487 488 // Copy to output 489 IPosition vecShape(1,shape(asap::ChanAxis)); 490 Vector<Float> outV = out.reform(vecShape); 491 std::vector<float> stlout; 492 outV.tovector(stlout); 493 // 494 return stlout; 495 } 496 497 498 499 500 Array<Float> SDMemTable::getStokesSpectrum(Int whichRow, Int iBeam, Int iIF) const 501 { 502 503 // Get data 504 505 Array<Float> arr; 506 stokesCol_.get(whichRow, arr); 507 508 // Set current cursor location and overwrite polarization axis 509 510 const IPosition& shape = arr.shape(); 511 IPosition start(shape.nelements(),0); 512 IPosition end(shape-1); 513 if (iBeam!=-1) { 514 start(asap::BeamAxis) = iBeam; 515 end(asap::BeamAxis) = iBeam; 516 } 517 if (iIF!=-1) { 518 start(asap::IFAxis) = iIF; 519 end(asap::IFAxis) = iIF; 520 } 521 522 // Get slice 523 524 return arr(start,end); 525 } 526 527 528 std::vector<string> SDMemTable::getCoordInfo() const 529 { 530 String un; 531 Table t = table_.keywordSet().asTable("FREQUENCIES"); 532 String sunit; 533 t.keywordSet().get("UNIT",sunit); 534 String dpl; 535 t.keywordSet().get("DOPPLER",dpl); 536 if (dpl == "") dpl = "RADIO"; 537 String rfrm; 538 t.keywordSet().get("REFFRAME",rfrm); 539 std::vector<string> inf; 540 inf.push_back(sunit); 541 inf.push_back(rfrm); 542 inf.push_back(dpl); 543 t.keywordSet().get("BASEREFFRAME",rfrm); 544 inf.push_back(rfrm); 545 return inf; 546 } 547 548 void SDMemTable::setCoordInfo(std::vector<string> theinfo) 549 { 550 std::vector<string>::iterator it; 551 String un,rfrm, brfrm,dpl; 552 un = theinfo[0]; // Abcissa unit 553 rfrm = theinfo[1]; // Active (or conversion) frame 554 dpl = theinfo[2]; // Doppler 555 brfrm = theinfo[3]; // Base frame 556 Table t = table_.rwKeywordSet().asTable("FREQUENCIES"); 557 558 Vector<Double> rstf; 559 t.keywordSet().get("RESTFREQS",rstf); 560 561 Bool canDo = True; 562 Unit u1("km/s");Unit u2("Hz"); 563 if (Unit(un) == u1) { 564 Vector<Double> rstf; 565 t.keywordSet().get("RESTFREQS",rstf); 566 if (rstf.nelements() == 0) { 567 throw(AipsError("Can't set unit to km/s if no restfrequencies are specified")); 568 } 569 } else if (Unit(un) != u2 && un != "") { 570 throw(AipsError("Unit not conformant with Spectral Coordinates")); 571 } 572 t.rwKeywordSet().define("UNIT", un); 573 574 MFrequency::Types mdr; 575 if (!MFrequency::getType(mdr, rfrm)) { 576 577 Int a,b;const uInt* c; 578 const String* valid = MFrequency::allMyTypes(a, b, c); 579 String pfix = "Please specify a legal frame type. Types are\n"; 580 throw(AipsError(pfix+(*valid))); 581 } else { 582 t.rwKeywordSet().define("REFFRAME",rfrm); 583 } 584 585 MDoppler::Types dtype; 586 dpl.upcase(); 587 if (!MDoppler::getType(dtype, dpl)) { 588 throw(AipsError("Doppler type unknown")); 589 } else { 590 t.rwKeywordSet().define("DOPPLER",dpl); 591 } 592 593 if (!MFrequency::getType(mdr, brfrm)) { 594 Int a,b;const uInt* c; 595 const String* valid = MFrequency::allMyTypes(a, b, c); 596 String pfix = "Please specify a legal frame type. Types are\n"; 597 throw(AipsError(pfix+(*valid))); 598 } else { 599 t.rwKeywordSet().define("BASEREFFRAME",brfrm); 600 } 601 } 602 603 604 std::vector<double> SDMemTable::getAbcissa(Int whichRow) const 605 { 606 std::vector<double> abc(nChan()); 607 608 // Get header units keyword 609 Table t = table_.keywordSet().asTable("FREQUENCIES"); 610 String sunit; 611 t.keywordSet().get("UNIT",sunit); 612 if (sunit == "") sunit = "pixel"; 613 Unit u(sunit); 614 615 // Easy if just wanting pixels 616 if (sunit==String("pixel")) { 617 // assume channels/pixels 618 std::vector<double>::iterator it; 619 uInt i=0; 620 for (it = abc.begin(); it != abc.end(); ++it) { 621 (*it) = Double(i++); 622 } 623 return abc; 624 } 625 626 // Continue with km/s or Hz. Get FreqIDs 627 Vector<uInt> freqIDs; 628 freqidCol_.get(whichRow, freqIDs); 629 uInt freqID = freqIDs(IFSel_); 630 restfreqidCol_.get(whichRow, freqIDs); 631 uInt restFreqID = freqIDs(IFSel_); 632 633 // Get SpectralCoordinate, set reference frame conversion, 634 // velocity conversion, and rest freq state 635 636 SpectralCoordinate spc = getSpectralCoordinate(freqID, restFreqID, whichRow); 637 Vector<Double> pixel(nChan()); 638 indgen(pixel); 639 640 if (u == Unit("km/s")) { 641 Vector<Double> world; 642 spc.pixelToVelocity(world,pixel); 643 std::vector<double>::iterator it; 644 uInt i = 0; 645 for (it = abc.begin(); it != abc.end(); ++it) { 646 (*it) = world[i]; 647 i++; 648 } 649 } else if (u == Unit("Hz")) { 650 651 // Set world axis units 652 Vector<String> wau(1); wau = u.getName(); 653 spc.setWorldAxisUnits(wau); 654 655 std::vector<double>::iterator it; 656 Double tmp; 657 uInt i = 0; 658 for (it = abc.begin(); it != abc.end(); ++it) { 659 spc.toWorld(tmp,pixel[i]); 660 (*it) = tmp; 661 i++; 662 } 663 } 664 return abc; 665 } 666 667 std::string SDMemTable::getAbcissaString(Int whichRow) const 668 { 669 Table t = table_.keywordSet().asTable("FREQUENCIES"); 670 671 String sunit; 672 t.keywordSet().get("UNIT",sunit); 673 if (sunit == "") sunit = "pixel"; 674 Unit u(sunit); 675 676 Vector<uInt> freqIDs; 677 freqidCol_.get(whichRow, freqIDs); 678 uInt freqID = freqIDs(IFSel_); 679 restfreqidCol_.get(whichRow, freqIDs); 680 uInt restFreqID = freqIDs(IFSel_); 681 682 // Get SpectralCoordinate, with frame, velocity, rest freq state set 683 SpectralCoordinate spc = getSpectralCoordinate(freqID, restFreqID, whichRow); 684 685 String s = "Channel"; 686 if (u == Unit("km/s")) { 687 s = CoordinateUtil::axisLabel(spc,0,True,True,True); 688 } else if (u == Unit("Hz")) { 689 Vector<String> wau(1);wau = u.getName(); 690 spc.setWorldAxisUnits(wau); 691 692 s = CoordinateUtil::axisLabel(spc,0,True,True,False); 693 } 694 return s; 695 } 696 697 void SDMemTable::setSpectrum(std::vector<float> spectrum, int whichRow) 698 { 699 Array<Float> arr; 700 specCol_.get(whichRow, arr); 701 if (spectrum.size() != arr.shape()(asap::ChanAxis)) { 702 throw(AipsError("Attempting to set spectrum with incorrect length.")); 703 } 704 705 // Setup accessors 706 ArrayAccessor<Float, Axis<asap::BeamAxis> > aa0(arr); 707 aa0.reset(aa0.begin(uInt(beamSel_))); // Beam selection 708 ArrayAccessor<Float, Axis<asap::IFAxis> > aa1(aa0); 709 aa1.reset(aa1.begin(uInt(IFSel_))); // IF selection 710 ArrayAccessor<Float, Axis<asap::PolAxis> > aa2(aa1); 711 aa2.reset(aa2.begin(uInt(polSel_))); // Pol selection 712 713 // Iterate 714 std::vector<float>::iterator it = spectrum.begin(); 715 for (ArrayAccessor<Float, Axis<asap::ChanAxis> > i(aa2); i != i.end(); ++i) { 716 (*i) = Float(*it); 717 it++; 718 } 719 specCol_.put(whichRow, arr); 720 } 721 722 void SDMemTable::getSpectrum(Vector<Float>& spectrum, Int whichRow) const 723 { 724 Array<Float> arr; 725 specCol_.get(whichRow, arr); 726 727 // Iterate and extract 728 spectrum.resize(arr.shape()(3)); 729 ArrayAccessor<Float, Axis<asap::BeamAxis> > aa0(arr); 730 aa0.reset(aa0.begin(uInt(beamSel_)));//go to beam 731 ArrayAccessor<Float, Axis<asap::IFAxis> > aa1(aa0); 732 aa1.reset(aa1.begin(uInt(IFSel_)));// go to IF 733 ArrayAccessor<Float, Axis<asap::PolAxis> > aa2(aa1); 734 aa2.reset(aa2.begin(uInt(polSel_)));// go to pol 735 736 ArrayAccessor<Float, Axis<asap::BeamAxis> > va(spectrum); 737 for (ArrayAccessor<Float, Axis<asap::ChanAxis> > i(aa2); i != i.end(); ++i) { 738 (*va) = (*i); 739 va++; 740 } 741 } 742 743 744 /* 745 void SDMemTable::getMask(Vector<Bool>& mask, Int whichRow) const { 746 Array<uChar> arr; 747 flagsCol_.get(whichRow, arr); 748 mask.resize(arr.shape()(3)); 749 750 ArrayAccessor<uChar, Axis<asap::BeamAxis> > aa0(arr); 751 aa0.reset(aa0.begin(uInt(beamSel_)));//go to beam 752 ArrayAccessor<uChar, Axis<asap::IFAxis> > aa1(aa0); 753 aa1.reset(aa1.begin(uInt(IFSel_)));// go to IF 754 ArrayAccessor<uChar, Axis<asap::PolAxis> > aa2(aa1); 755 aa2.reset(aa2.begin(uInt(polSel_)));// go to pol 756 757 Bool useUserMask = ( chanMask_.size() == arr.shape()(3) ); 758 759 ArrayAccessor<Bool, Axis<asap::BeamAxis> > va(mask); 760 std::vector<bool> tmp; 761 tmp = chanMask_; // WHY the fxxx do I have to make a copy here. The 762 // iterator should work on chanMask_?? 763 std::vector<bool>::iterator miter; 764 miter = tmp.begin(); 765 766 for (ArrayAccessor<uChar, Axis<asap::ChanAxis> > i(aa2); i != i.end(); ++i) { 767 bool out =!static_cast<bool>(*i); 768 if (useUserMask) { 769 out = out && (*miter); 770 miter++; 771 } 772 (*va) = out; 773 va++; 774 } 775 } 776 */ 777 778 MaskedArray<Float> SDMemTable::rowAsMaskedArray(uInt whichRow, 779 Bool toStokes) const 780 { 781 // Get flags 782 Array<uChar> farr; 783 flagsCol_.get(whichRow, farr); 784 785 // Get data and convert mask to Bool 786 Array<Float> arr; 787 Array<Bool> mask; 788 if (toStokes) { 789 stokesCol_.get(whichRow, arr); 790 791 Array<Bool> tMask(farr.shape()); 792 convertArray(tMask, farr); 793 mask = SDPolUtil::stokesData (tMask, True); 794 } else { 795 specCol_.get(whichRow, arr); 796 mask.resize(farr.shape()); 797 convertArray(mask, farr); 798 } 799 800 return MaskedArray<Float>(arr,!mask); 801 } 802 803 Float SDMemTable::getTsys(Int whichRow) const 804 { 805 Array<Float> arr; 806 tsCol_.get(whichRow, arr); 807 Float out; 808 809 IPosition ip(arr.shape()); 810 ip(asap::BeamAxis) = beamSel_; 811 ip(asap::IFAxis) = IFSel_; 812 ip(asap::PolAxis) = polSel_; 813 ip(asap::ChanAxis)=0; // First channel only 814 815 out = arr(ip); 816 return out; 817 } 818 819 MDirection SDMemTable::getDirection(Int whichRow, Bool refBeam) const 820 { 821 MDirection::Types mdr = getDirectionReference(); 822 Array<Double> posit; 823 dirCol_.get(whichRow,posit); 824 Vector<Double> wpos(2); 825 Int rb; 826 rbeamCol_.get(whichRow,rb); 827 wpos[0] = posit(IPosition(2,beamSel_,0)); 828 wpos[1] = posit(IPosition(2,beamSel_,1)); 829 if (refBeam && rb > -1) { // use refBeam instead if it exists 830 wpos[0] = posit(IPosition(2,rb,0)); 831 wpos[1] = posit(IPosition(2,rb,1)); 832 } 833 834 Quantum<Double> lon(wpos[0],Unit(String("rad"))); 835 Quantum<Double> lat(wpos[1],Unit(String("rad"))); 836 return MDirection(lon, lat, mdr); 837 } 838 839 MEpoch SDMemTable::getEpoch(Int whichRow) const 840 { 841 MEpoch::Types met = getTimeReference(); 842 843 Double obstime; 844 timeCol_.get(whichRow,obstime); 845 MVEpoch tm2(Quantum<Double>(obstime, Unit(String("d")))); 846 return MEpoch(tm2, met); 847 } 848 849 MPosition SDMemTable::getAntennaPosition () const 850 { 851 Vector<Double> antpos; 852 table_.keywordSet().get("AntennaPosition", antpos); 853 MVPosition mvpos(antpos(0),antpos(1),antpos(2)); 854 return MPosition(mvpos); 855 } 856 857 858 SpectralCoordinate SDMemTable::getSpectralCoordinate(uInt freqID) const 859 { 860 861 Table t = table_.keywordSet().asTable("FREQUENCIES"); 862 if (freqID> t.nrow() ) { 863 throw(AipsError("SDMemTable::getSpectralCoordinate - freqID out of range")); 864 } 865 866 Double rp,rv,inc; 867 String rf; 868 ROScalarColumn<Double> rpc(t, "REFPIX"); 869 ROScalarColumn<Double> rvc(t, "REFVAL"); 870 ROScalarColumn<Double> incc(t, "INCREMENT"); 871 t.keywordSet().get("BASEREFFRAME",rf); 872 873 // Create SpectralCoordinate (units Hz) 874 MFrequency::Types mft; 875 if (!MFrequency::getType(mft, rf)) { 876 ostringstream oss; 877 pushLog("WARNING: Frequency type unknown assuming TOPO"); 878 mft = MFrequency::TOPO; 879 } 880 rpc.get(freqID, rp); 881 rvc.get(freqID, rv); 882 incc.get(freqID, inc); 883 884 SpectralCoordinate spec(mft,rv,inc,rp); 885 return spec; 886 } 887 888 889 SpectralCoordinate SDMemTable::getSpectralCoordinate(uInt freqID, 890 uInt restFreqID, 891 uInt whichRow) const 892 { 893 894 // Create basic SC 895 SpectralCoordinate spec = getSpectralCoordinate (freqID); 896 897 Table t = table_.keywordSet().asTable("FREQUENCIES"); 898 899 // Set rest frequency 900 Vector<Double> restFreqIDs; 901 t.keywordSet().get("RESTFREQS",restFreqIDs); 902 if (restFreqID < restFreqIDs.nelements()) { // Should always be true 903 spec.setRestFrequency(restFreqIDs(restFreqID),True); 904 } 905 906 // Set up frame conversion layer 907 String frm; 908 t.keywordSet().get("REFFRAME",frm); 909 if (frm == "") frm = "TOPO"; 910 MFrequency::Types mtype; 911 if (!MFrequency::getType(mtype, frm)) { 912 // Should never happen 913 pushLog("WARNING: Frequency type unknown assuming TOPO"); 914 mtype = MFrequency::TOPO; 915 } 916 917 // Set reference frame conversion (requires row) 918 MDirection dir = getDirection(whichRow); 919 MEpoch epoch = getEpoch(whichRow); 920 MPosition pos = getAntennaPosition(); 921 922 if (!spec.setReferenceConversion(mtype,epoch,pos,dir)) { 923 throw(AipsError("Couldn't convert frequency frame.")); 924 } 925 926 // Now velocity conversion if appropriate 927 String unitStr; 928 t.keywordSet().get("UNIT",unitStr); 929 930 String dpl; 931 t.keywordSet().get("DOPPLER",dpl); 932 MDoppler::Types dtype; 933 MDoppler::getType(dtype, dpl); 934 935 // Only set velocity unit if non-blank and non-Hz 936 if (!unitStr.empty()) { 937 Unit unitU(unitStr); 938 if (unitU==Unit("Hz")) { 939 } else { 940 spec.setVelocity(unitStr, dtype); 941 } 942 } 943 944 return spec; 945 } 946 947 948 Bool SDMemTable::setCoordinate(const SpectralCoordinate& speccord, 949 uInt freqID) { 950 Table t = table_.rwKeywordSet().asTable("FREQUENCIES"); 951 if (freqID > t.nrow() ) { 952 throw(AipsError("SDMemTable::setCoordinate - coord no out of range")); 953 } 954 ScalarColumn<Double> rpc(t, "REFPIX"); 955 ScalarColumn<Double> rvc(t, "REFVAL"); 956 ScalarColumn<Double> incc(t, "INCREMENT"); 957 958 rpc.put(freqID, speccord.referencePixel()[0]); 959 rvc.put(freqID, speccord.referenceValue()[0]); 960 incc.put(freqID, speccord.increment()[0]); 961 962 return True; 963 } 964 965 Int SDMemTable::nCoordinates() const 966 { 967 return table_.keywordSet().asTable("FREQUENCIES").nrow(); 968 } 969 970 971 std::vector<double> SDMemTable::getRestFreqs() const 972 { 973 Table t = table_.keywordSet().asTable("FREQUENCIES"); 974 Vector<Double> tvec; 975 t.keywordSet().get("RESTFREQS",tvec); 976 std::vector<double> stlout; 977 tvec.tovector(stlout); 978 return stlout; 979 } 980 981 bool SDMemTable::putSDFreqTable(const SDFrequencyTable& sdft) 199 void Scantable::setupFitTable() 982 200 { 983 201 TableDesc td("", "1", TableDesc::Scratch); 984 td.addColumn(ScalarColumnDesc<Double>("REFPIX")); 985 td.addColumn(ScalarColumnDesc<Double>("REFVAL")); 986 td.addColumn(ScalarColumnDesc<Double>("INCREMENT")); 987 SetupNewTable aNewTab("freqs", td, Table::New); 988 Table aTable (aNewTab, Table::Memory, sdft.length()); 989 ScalarColumn<Double> sc0(aTable, "REFPIX"); 990 ScalarColumn<Double> sc1(aTable, "REFVAL"); 991 ScalarColumn<Double> sc2(aTable, "INCREMENT"); 992 for (uInt i=0; i < sdft.length(); ++i) { 993 sc0.put(i,sdft.referencePixel(i)); 994 sc1.put(i,sdft.referenceValue(i)); 995 sc2.put(i,sdft.increment(i)); 996 } 997 String rf = sdft.refFrame(); 998 if (rf.contains("TOPO")) rf = "TOPO"; 999 String brf = sdft.baseRefFrame(); 1000 if (brf.contains("TOPO")) brf = "TOPO"; 1001 1002 aTable.rwKeywordSet().define("BASEREFFRAME", brf); 1003 aTable.rwKeywordSet().define("REFFRAME", rf); 1004 aTable.rwKeywordSet().define("EQUINOX", sdft.equinox()); 1005 aTable.rwKeywordSet().define("UNIT", sdft.unit()); 1006 aTable.rwKeywordSet().define("DOPPLER", String("RADIO")); 1007 Vector<Double> rfvec; 1008 String rfunit; 1009 sdft.restFrequencies(rfvec,rfunit); 1010 Quantum<Vector<Double> > q(rfvec, rfunit); 1011 rfvec.resize(); 1012 rfvec = q.getValue("Hz"); 1013 aTable.rwKeywordSet().define("RESTFREQS", rfvec); 1014 table_.rwKeywordSet().defineTable("FREQUENCIES", aTable); 1015 return true; 1016 } 1017 1018 bool SDMemTable::putSDFitTable(const SDFitTable& sdft) 1019 { 1020 TableDesc td("", "1", TableDesc::Scratch); 202 td.addColumn(ScalarColumnDesc<uInt>("FIT_ID")); 1021 203 td.addColumn(ArrayColumnDesc<String>("FUNCTIONS")); 1022 204 td.addColumn(ArrayColumnDesc<Int>("COMPONENTS")); … … 1026 208 SetupNewTable aNewTab("fits", td, Table::New); 1027 209 Table aTable(aNewTab, Table::Memory); 1028 ArrayColumn<String> sc0(aTable, "FUNCTIONS");1029 ArrayColumn<Int> sc1(aTable, "COMPONENTS");1030 ArrayColumn<Double> sc2(aTable, "PARAMETERS");1031 ArrayColumn<Bool> sc3(aTable, "PARMASK");1032 ArrayColumn<String> sc4(aTable, "FRAMEINFO");1033 for (uInt i; i<sdft.length(); ++i) {1034 const Vector<Double>& parms = sdft.getParameters(i);1035 const Vector<Bool>& parmask = sdft.getParameterMask(i);1036 const Vector<String>& funcs = sdft.getFunctions(i);1037 const Vector<Int>& comps = sdft.getComponents(i);1038 const Vector<String>& finfo = sdft.getFrameInfo(i);1039 sc0.put(i,funcs);1040 sc1.put(i,comps);1041 sc3.put(i,parmask);1042 sc2.put(i,parms);1043 sc4.put(i,finfo);1044 }1045 210 table_.rwKeywordSet().defineTable("FITS", aTable); 1046 return true; 1047 } 1048 1049 SDFitTable SDMemTable::getSDFitTable() const 1050 { 1051 const Table& t = table_.keywordSet().asTable("FITS"); 1052 Vector<Double> parms; 1053 Vector<Bool> parmask; 1054 Vector<String> funcs; 1055 Vector<String> finfo; 1056 Vector<Int> comps; 1057 ROArrayColumn<Double> parmsCol(t, "PARAMETERS"); 1058 ROArrayColumn<Bool> parmaskCol(t, "PARMASK"); 1059 ROArrayColumn<Int> compsCol(t, "COMPONENTS"); 1060 ROArrayColumn<String> funcsCol(t, "FUNCTIONS"); 1061 ROArrayColumn<String> finfoCol(t, "FRAMEINFO"); 1062 uInt n = t.nrow(); 1063 SDFitTable sdft; 1064 for (uInt i=0; i<n; ++i) { 1065 parmaskCol.get(i, parmask); 1066 parmsCol.get(i, parms); 1067 funcsCol.get(i, funcs); 1068 compsCol.get(i, comps); 1069 finfoCol.get(i, finfo); 1070 sdft.addFit(parms, parmask, funcs, comps, finfo); 1071 } 1072 return sdft; 1073 } 1074 1075 SDFitTable SDMemTable::getSDFitTable(uInt whichRow) const { 1076 const Table& t = table_.keywordSet().asTable("FITS"); 1077 if (t.nrow() == 0 || whichRow >= t.nrow()) return SDFitTable(); 1078 Array<Int> fitid; 1079 fitCol_.get(whichRow, fitid); 1080 if (fitid.nelements() == 0) return SDFitTable(); 1081 1082 IPosition shp = fitid.shape(); 1083 IPosition start(4, beamSel_, IFSel_, polSel_,0); 1084 IPosition end(4, beamSel_, IFSel_, polSel_, shp[3]-1); 1085 1086 // reform the output array slice to be of dim=1 1087 Vector<Int> tmp = (fitid(start, end)).reform(IPosition(1,shp[3])); 1088 1089 Vector<Double> parms; 1090 Vector<Bool> parmask; 1091 Vector<String> funcs; 1092 Vector<String> finfo; 1093 Vector<Int> comps; 1094 ROArrayColumn<Double> parmsCol(t, "PARAMETERS"); 1095 ROArrayColumn<Bool> parmaskCol(t, "PARMASK"); 1096 ROArrayColumn<Int> compsCol(t, "COMPONENTS"); 1097 ROArrayColumn<String> funcsCol(t, "FUNCTIONS"); 1098 ROArrayColumn<String> finfoCol(t, "FRAMEINFO"); 1099 SDFitTable sdft; 1100 Int k=-1; 1101 for (uInt i=0; i< tmp.nelements(); ++i) { 1102 k = tmp[i]; 1103 if ( k > -1 && k < t.nrow() ) { 1104 parms.resize(); 1105 parmsCol.get(k, parms); 1106 parmask.resize(); 1107 parmaskCol.get(k, parmask); 1108 funcs.resize(); 1109 funcsCol.get(k, funcs); 1110 comps.resize(); 1111 compsCol.get(k, comps); 1112 finfo.resize(); 1113 finfoCol.get(k, finfo); 1114 sdft.addFit(parms, parmask, funcs, comps, finfo); 1115 } 1116 } 1117 return sdft; 1118 } 1119 1120 void SDMemTable::addFit(uInt whichRow, 1121 const Vector<Double>& p, const Vector<Bool>& m, 1122 const Vector<String>& f, const Vector<Int>& c) 1123 { 1124 if (whichRow >= nRow()) { 1125 throw(AipsError("Specified row out of range")); 1126 } 1127 Table t = table_.keywordSet().asTable("FITS"); 1128 uInt nrow = t.nrow(); 1129 t.addRow(); 1130 ArrayColumn<Double> parmsCol(t, "PARAMETERS"); 1131 ArrayColumn<Bool> parmaskCol(t, "PARMASK"); 1132 ArrayColumn<Int> compsCol(t, "COMPONENTS"); 1133 ArrayColumn<String> funcsCol(t, "FUNCTIONS"); 1134 ArrayColumn<String> finfoCol(t, "FRAMEINFO"); 1135 parmsCol.put(nrow, p); 1136 parmaskCol.put(nrow, m); 1137 compsCol.put(nrow, c); 1138 funcsCol.put(nrow, f); 1139 Vector<String> fi = mathutil::toVectorString(getCoordInfo()); 1140 finfoCol.put(nrow, fi); 1141 1142 Array<Int> fitarr; 1143 fitCol_.get(whichRow, fitarr); 1144 1145 Array<Int> newarr; // The new Array containing the fitid 1146 Int pos =-1; // The fitid position in the array 1147 if ( fitarr.nelements() == 0 ) { // no fits at all in this row 1148 Array<Int> arr(IPosition(4,nBeam(),nIF(),nPol(),1)); 1149 arr = -1; 1150 newarr.reference(arr); 1151 pos = 0; 1152 } else { 1153 IPosition shp = fitarr.shape(); 1154 IPosition start(4, beamSel_, IFSel_, polSel_,0); 1155 IPosition end(4, beamSel_, IFSel_, polSel_, shp[3]-1); 1156 // reform the output array slice to be of dim=1 1157 Array<Int> tmp = (fitarr(start, end)).reform(IPosition(1,shp[3])); 1158 const Vector<Int>& fits = tmp; 1159 VectorSTLIterator<Int> it(fits); 1160 Int i = 0; 1161 while (it != fits.end()) { 1162 if (*it == -1) { 1163 pos = i; 1164 break; 1165 } 1166 ++i; 1167 ++it; 1168 }; 1169 } 1170 if (pos == -1) { 1171 mathutil::extendLastArrayAxis(newarr,fitarr, -1); 1172 pos = fitarr.shape()[3]; // the new element position 1173 } else { 1174 if (fitarr.nelements() > 0) 1175 newarr = fitarr; 1176 } 1177 newarr(IPosition(4, beamSel_, IFSel_, polSel_, pos)) = Int(nrow); 1178 fitCol_.put(whichRow, newarr); 1179 1180 } 1181 1182 SDFrequencyTable SDMemTable::getSDFreqTable() const 1183 { 1184 const Table& t = table_.keywordSet().asTable("FREQUENCIES"); 1185 SDFrequencyTable sdft; 1186 1187 // Add refpix/refval/incr. What are the units ? Hz I suppose 1188 // but it's nowhere described... 1189 Vector<Double> refPix, refVal, incr; 1190 ScalarColumn<Double> refPixCol(t, "REFPIX"); 1191 ScalarColumn<Double> refValCol(t, "REFVAL"); 1192 ScalarColumn<Double> incrCol(t, "INCREMENT"); 1193 refPix = refPixCol.getColumn(); 1194 refVal = refValCol.getColumn(); 1195 incr = incrCol.getColumn(); 1196 1197 uInt n = refPix.nelements(); 1198 for (uInt i=0; i<n; i++) { 1199 sdft.addFrequency(refPix[i], refVal[i], incr[i]); 1200 } 1201 1202 // Frequency reference frame. I don't know if this 1203 // is the correct frame. It might be 'REFFRAME' 1204 // rather than 'BASEREFFRAME' ? 1205 String frame; 1206 t.keywordSet().get("REFFRAME",frame); 1207 sdft.setRefFrame(frame); 1208 t.keywordSet().get("BASEREFFRAME",frame); 1209 sdft.setBaseRefFrame(frame); 1210 1211 // Equinox 1212 Float equinox; 1213 t.keywordSet().get("EQUINOX", equinox); 1214 sdft.setEquinox(equinox); 1215 1216 String unit; 1217 t.keywordSet().get("UNIT", unit); 1218 sdft.setUnit(unit); 1219 1220 // Rest Frequency 1221 Vector<Double> restFreqs; 1222 t.keywordSet().get("RESTFREQS", restFreqs); 1223 for (uInt i=0; i<restFreqs.nelements(); i++) { 1224 sdft.addRestFrequency(restFreqs[i]); 1225 } 1226 sdft.setRestFrequencyUnit(String("Hz")); 1227 1228 return sdft; 1229 } 1230 1231 bool SDMemTable::putSDContainer(const SDContainer& sdc) 1232 { 1233 uInt rno = table_.nrow(); 1234 table_.addRow(); 1235 1236 timeCol_.put(rno, sdc.timestamp); 1237 srcnCol_.put(rno, sdc.sourcename); 1238 fldnCol_.put(rno, sdc.fieldname); 1239 specCol_.put(rno, sdc.getSpectrum()); 1240 flagsCol_.put(rno, sdc.getFlags()); 1241 tsCol_.put(rno, sdc.getTsys()); 1242 scanCol_.put(rno, sdc.scanid); 1243 integrCol_.put(rno, sdc.interval); 1244 freqidCol_.put(rno, sdc.getFreqMap()); 1245 restfreqidCol_.put(rno, sdc.getRestFreqMap()); 1246 dirCol_.put(rno, sdc.getDirection()); 1247 rbeamCol_.put(rno, sdc.refbeam); 1248 tcalCol_.put(rno, sdc.tcal); 1249 tcaltCol_.put(rno, sdc.tcaltime); 1250 azCol_.put(rno, sdc.azimuth); 1251 elCol_.put(rno, sdc.elevation); 1252 paraCol_.put(rno, sdc.parangle); 1253 fitCol_.put(rno, sdc.getFitMap()); 1254 return true; 1255 } 1256 1257 SDContainer SDMemTable::getSDContainer(uInt whichRow) const 1258 { 1259 SDContainer sdc(nBeam(),nIF(),nPol(),nChan()); 1260 timeCol_.get(whichRow, sdc.timestamp); 1261 srcnCol_.get(whichRow, sdc.sourcename); 1262 integrCol_.get(whichRow, sdc.interval); 1263 scanCol_.get(whichRow, sdc.scanid); 1264 fldnCol_.get(whichRow, sdc.fieldname); 1265 rbeamCol_.get(whichRow, sdc.refbeam); 1266 azCol_.get(whichRow, sdc.azimuth); 1267 elCol_.get(whichRow, sdc.elevation); 1268 paraCol_.get(whichRow, sdc.parangle); 1269 Vector<Float> tc; 1270 tcalCol_.get(whichRow, tc); 1271 sdc.tcal[0] = tc[0];sdc.tcal[1] = tc[1]; 1272 tcaltCol_.get(whichRow, sdc.tcaltime); 1273 1274 Array<Float> spectrum; 1275 Array<Float> tsys; 1276 Array<uChar> flagtrum; 1277 Vector<uInt> fmap; 1278 Array<Double> direction; 1279 Array<Int> fits; 1280 1281 specCol_.get(whichRow, spectrum); 1282 sdc.putSpectrum(spectrum); 1283 flagsCol_.get(whichRow, flagtrum); 1284 sdc.putFlags(flagtrum); 1285 tsCol_.get(whichRow, tsys); 1286 sdc.putTsys(tsys); 1287 freqidCol_.get(whichRow, fmap); 1288 sdc.putFreqMap(fmap); 1289 restfreqidCol_.get(whichRow, fmap); 1290 sdc.putRestFreqMap(fmap); 1291 dirCol_.get(whichRow, direction); 1292 sdc.putDirection(direction); 1293 fitCol_.get(whichRow, fits); 1294 sdc.putFitMap(fits); 1295 return sdc; 1296 } 1297 1298 bool SDMemTable::putSDHeader(const SDHeader& sdh) 211 } 212 213 void Scantable::attach() 214 { 215 historyTable_ = table_.keywordSet().asTable("HISTORY"); 216 fitTable_ = table_.keywordSet().asTable("FITS"); 217 218 timeCol_.attach(table_, "TIME"); 219 srcnCol_.attach(table_, "SRCNAME"); 220 specCol_.attach(table_, "SPECTRA"); 221 flagsCol_.attach(table_, "FLAGTRA"); 222 tsCol_.attach(table_, "TSYS"); 223 cycleCol_.attach(table_,"CYCLENO"); 224 scanCol_.attach(table_, "SCANNO"); 225 beamCol_.attach(table_, "BEAMNO"); 226 integrCol_.attach(table_, "INTERVAL"); 227 azCol_.attach(table_, "AZIMUTH"); 228 elCol_.attach(table_, "ELEVATION"); 229 dirCol_.attach(table_, "DIRECTION"); 230 paraCol_.attach(table_, "PARANGLE"); 231 fldnCol_.attach(table_, "FIELDNAME"); 232 rbeamCol_.attach(table_, "REFBEAMNO"); 233 234 mfitidCol_.attach(table_,"FIT_ID"); 235 fitidCol_.attach(fitTable_,"FIT_ID"); 236 237 mfreqidCol_.attach(table_, "FREQ_ID"); 238 239 mtcalidCol_.attach(table_, "TCAL_ID"); 240 241 mfocusidCol_.attach(table_, "FOCUS_ID"); 242 243 mmolidCol_.attach(table_, "MOLECULE_ID"); 244 } 245 246 void Scantable::putSDHeader(const SDHeader& sdh) 1299 247 { 1300 248 table_.rwKeywordSet().define("nIF", sdh.nif); … … 1314 262 table_.rwKeywordSet().define("FluxUnit", sdh.fluxunit); 1315 263 table_.rwKeywordSet().define("Epoch", sdh.epoch); 1316 return true; 1317 } 1318 1319 SDHeader SDMemTable::getSDHeader() const 264 } 265 266 SDHeader Scantable::getSDHeader() const 1320 267 { 1321 268 SDHeader sdh; … … 1338 285 return sdh; 1339 286 } 1340 void SDMemTable::makePersistent(const std::string& filename) 1341 { 1342 table_.deepCopy(filename,Table::New); 1343 1344 } 1345 1346 Int SDMemTable::nScan() const { 1347 Int n = 0; 1348 Int previous = -1;Int current=0; 287 288 int Scantable::rowToScanIndex( int therow ) 289 { 290 int therealrow = -1; 291 292 return therealrow; 293 } 294 295 int Scantable::nScan() const { 296 int n = 0; 297 int previous = -1; int current = 0; 1349 298 for (uInt i=0; i< scanCol_.nrow();i++) { 1350 299 scanCol_.getScalar(i,current); … … 1357 306 } 1358 307 1359 String SDMemTable::formatSec(Double x) const308 std::string Scantable::formatSec(Double x) const 1360 309 { 1361 310 Double xcop = x; … … 1374 323 }; 1375 324 1376 String SDMemTable::formatDirection(const MDirection& md) const325 std::string Scantable::formatDirection(const MDirection& md) const 1377 326 { 1378 327 Vector<Double> t = md.getAngle(Unit(String("rad"))).getValue(); … … 1387 336 1388 337 1389 std::string S DMemTable::getFluxUnit() const338 std::string Scantable::getFluxUnit() const 1390 339 { 1391 340 String tmp; … … 1394 343 } 1395 344 1396 void S DMemTable::setFluxUnit(const std::string& unit)345 void Scantable::setFluxUnit(const std::string& unit) 1397 346 { 1398 347 String tmp(unit); … … 1405 354 } 1406 355 1407 1408 void SDMemTable::setInstrument(const std::string& name) 1409 { 1410 Bool throwIt = True; 356 void Scantable::setInstrument(const std::string& name) 357 { 358 bool throwIt = true; 1411 359 Instrument ins = SDAttr::convertInstrument(name, throwIt); 1412 360 String nameU(name); … … 1415 363 } 1416 364 1417 std::string SDMemTable::summary(bool verbose) const { 1418 1419 // Format header info 1420 ostringstream oss; 1421 oss << endl; 1422 oss << "--------------------------------------------------------------------------------" << endl; 1423 oss << " Scan Table Summary" << endl; 1424 oss << "--------------------------------------------------------------------------------" << endl; 1425 oss.flags(std::ios_base::left); 1426 oss << setw(15) << "Beams:" << setw(4) << nBeam() << endl 1427 << setw(15) << "IFs:" << setw(4) << nIF() << endl 1428 << setw(15) << "Polarisations:" << setw(4) << nPol() << endl 1429 << setw(15) << "Channels:" << setw(4) << nChan() << endl; 1430 oss << endl; 1431 String tmp; 1432 table_.keywordSet().get("Observer", tmp); 1433 oss << setw(15) << "Observer:" << tmp << endl; 1434 oss << setw(15) << "Obs Date:" << getTime(-1,True) << endl; 1435 table_.keywordSet().get("Project", tmp); 1436 oss << setw(15) << "Project:" << tmp << endl; 1437 table_.keywordSet().get("Obstype", tmp); 1438 oss << setw(15) << "Obs. Type:" << tmp << endl; 1439 table_.keywordSet().get("AntennaName", tmp); 1440 oss << setw(15) << "Antenna Name:" << tmp << endl; 1441 table_.keywordSet().get("FluxUnit", tmp); 1442 oss << setw(15) << "Flux Unit:" << tmp << endl; 1443 Table t = table_.keywordSet().asTable("FREQUENCIES"); 1444 Vector<Double> vec; 1445 t.keywordSet().get("RESTFREQS",vec); 1446 oss << setw(15) << "Rest Freqs:"; 1447 if (vec.nelements() > 0) { 1448 oss << setprecision(10) << vec << " [Hz]" << endl; 1449 } else { 1450 oss << "None set" << endl; 1451 } 1452 oss << setw(15) << "Abcissa:" << getAbcissaString() << endl; 1453 oss << setw(15) << "Cursor:" << "Beam[" << getBeam() << "] " 1454 << "IF[" << getIF() << "] " << "Pol[" << getPol() << "]" << endl; 1455 oss << endl; 1456 1457 String dirtype ="Position ("+ MDirection::showType(getDirectionReference()) + ")"; 1458 oss << setw(5) << "Scan" 1459 << setw(15) << "Source" 1460 << setw(24) << dirtype 1461 << setw(10) << "Time" 1462 << setw(18) << "Integration" 1463 << setw(7) << "FreqIDs" << endl; 1464 oss << "--------------------------------------------------------------------------------" << endl; 1465 1466 // Generate list of scan start and end integrations 1467 Vector<Int> scanIDs = scanCol_.getColumn(); 1468 Vector<uInt> startInt, endInt; 1469 mathutil::scanBoundaries(startInt, endInt, scanIDs); 1470 1471 const uInt nScans = startInt.nelements(); 1472 String name; 1473 Vector<uInt> freqIDs, listFQ; 1474 Vector<uInt> restFreqIDs, listRestFQ; 1475 uInt nInt; 1476 1477 for (uInt i=0; i<nScans; i++) { 1478 // Get things from first integration of scan 1479 String time = getTime(startInt(i),False); 1480 String tInt = formatSec(Double(getInterval(startInt(i)))); 1481 String posit = formatDirection(getDirection(startInt(i),True)); 1482 srcnCol_.getScalar(startInt(i),name); 1483 1484 // Find all the FreqIDs in this scan 1485 listFQ.resize(0); 1486 listRestFQ.resize(0); 1487 for (uInt j=startInt(i); j<endInt(i)+1; j++) { 1488 freqidCol_.get(j, freqIDs); 1489 for (uInt k=0; k<freqIDs.nelements(); k++) { 1490 mathutil::addEntry(listFQ, freqIDs(k)); 1491 } 1492 // 1493 restfreqidCol_.get(j, restFreqIDs); 1494 for (uInt k=0; k<restFreqIDs.nelements(); k++) { 1495 mathutil::addEntry(listRestFQ, restFreqIDs(k)); 1496 } 1497 } 1498 1499 nInt = endInt(i) - startInt(i) + 1; 1500 oss << setw(3) << std::right << i << std::left << setw(2) << " " 1501 << setw(15) << name 1502 << setw(24) << posit 1503 << setw(10) << time 1504 << setw(3) << std::right << nInt << setw(3) << " x " << std::left 1505 << setw(6) << tInt 1506 << " " << listFQ << " " << listRestFQ << endl; 1507 } 1508 oss << endl; 1509 oss << "Table contains " << table_.nrow() << " integration(s) in " 1510 << nScans << " scan(s)." << endl; 1511 1512 // Frequency Table 1513 if (verbose) { 1514 std::vector<string> info = getCoordInfo(); 1515 SDFrequencyTable sdft = getSDFreqTable(); 1516 oss << endl << endl; 1517 oss << "FreqID Frame RefFreq(Hz) RefPix Increment(Hz)" << endl; 1518 oss << "--------------------------------------------------------------------------------" << endl; 1519 for (uInt i=0; i<sdft.length(); i++) { 1520 oss << setw(8) << i << setw(8) 1521 << info[3] << setw(16) << setprecision(8) 1522 << sdft.referenceValue(i) << setw(10) 1523 << sdft.referencePixel(i) << setw(12) 1524 << sdft.increment(i) << endl; 1525 } 1526 oss << "--------------------------------------------------------------------------------" << endl; 1527 } 1528 return String(oss); 1529 } 1530 /* 1531 std::string SDMemTable::scanSummary(const std::vector<int>& whichScans) { 1532 ostringstream oss; 1533 Vector<Int> scanIDs = scanCol_.getColumn(); 1534 Vector<uInt> startInt, endInt; 1535 mathutil::scanBoundaries(startInt, endInt, scanIDs); 1536 const uInt nScans = startInt.nelements(); 1537 std::vector<int>::const_iterator it(whichScans); 1538 return String(oss); 1539 } 1540 */ 1541 Int SDMemTable::nBeam() const 1542 { 1543 Int n; 1544 table_.keywordSet().get("nBeam",n); 1545 return n; 1546 } 1547 1548 Int SDMemTable::nIF() const { 1549 Int n; 1550 table_.keywordSet().get("nIF",n); 1551 return n; 1552 } 1553 1554 Int SDMemTable::nPol() const { 1555 Int n; 1556 table_.keywordSet().get("nPol",n); 1557 return n; 1558 } 1559 1560 Int SDMemTable::nChan() const { 1561 Int n; 1562 table_.keywordSet().get("nChan",n); 1563 return n; 1564 } 1565 1566 Table SDMemTable::getHistoryTable() const 365 MPosition Scantable::getAntennaPosition () const 366 { 367 Vector<Double> antpos; 368 table_.keywordSet().get("AntennaPosition", antpos); 369 MVPosition mvpos(antpos(0),antpos(1),antpos(2)); 370 return MPosition(mvpos); 371 } 372 373 void Scantable::makePersistent(const std::string& filename) 374 { 375 String inname(filename); 376 Path path(inname); 377 inname = path.expandedName(); 378 table_.deepCopy(inname, Table::New); 379 } 380 381 int asap::Scantable::nbeam( int scanno ) const 382 { 383 if ( scanno < 0 ) { 384 Int n; 385 table_.keywordSet().get("nBeam",n); 386 return int(n); 387 } else { 388 // take the first POLNO,IFNO,CYCLENO as nbeam shouldn't vary with these 389 Table tab = table_(table_.col("SCANNO") == scanno 390 && table_.col("POLNO") == 0 391 && table_.col("IFNO") == 0 392 && table_.col("CYCLENO") == 0 ); 393 ROTableVector<uInt> v(tab, "BEAMNO"); 394 return int(v.nelements()); 395 } 396 return 0; 397 } 398 399 int asap::Scantable::nif( int scanno ) const 400 { 401 if ( scanno < 0 ) { 402 Int n; 403 table_.keywordSet().get("nIF",n); 404 return int(n); 405 } else { 406 // take the first POLNO,BEAMNO,CYCLENO as nbeam shouldn't vary with these 407 Table tab = table_(table_.col("SCANNO") == scanno 408 && table_.col("POLNO") == 0 409 && table_.col("BEAMNO") == 0 410 && table_.col("CYCLENO") == 0 ); 411 ROTableVector<uInt> v(tab, "IFNO"); 412 return int(v.nelements()); 413 } 414 return 0; 415 } 416 417 int asap::Scantable::npol( int scanno ) const 418 { 419 if ( scanno < 0 ) { 420 Int n; 421 table_.keywordSet().get("nPol",n); 422 return n; 423 } else { 424 // take the first POLNO,IFNO,CYCLENO as nbeam shouldn't vary with these 425 Table tab = table_(table_.col("SCANNO") == scanno 426 && table_.col("IFNO") == 0 427 && table_.col("BEAMNO") == 0 428 && table_.col("CYCLENO") == 0 ); 429 ROTableVector<uInt> v(tab, "POLNO"); 430 return int(v.nelements()); 431 } 432 return 0; 433 } 434 435 int asap::Scantable::nrow( int scanno ) const 436 { 437 if ( scanno < 0 ) { 438 return int(table_.nrow()); 439 } else { 440 // take the first POLNO,IFNO,CYCLENO as nbeam shouldn't vary with these 441 Table tab = table_(table_.col("SCANNO") == scanno 442 && table_.col("BEAMNO") == 0 443 && table_.col("IFNO") == 0 444 && table_.col("POLNO") == 0 ); 445 return int(tab.nrow()); 446 } 447 return 0; 448 } 449 450 451 int asap::Scantable::nchan( int scanno, int ifno ) const 452 { 453 if ( scanno < 0 && ifno < 0 ) { 454 Int n; 455 table_.keywordSet().get("nChan",n); 456 return int(n); 457 } else { 458 // take the first POLNO,IFNO,CYCLENO as nbeam shouldn't vary with these 459 Table tab = table_(table_.col("SCANNO") == scanno 460 && table_.col("IFNO") == ifno 461 && table_.col("BEAMNO") == 0 462 && table_.col("POLNO") == 0 463 && table_.col("CYCLENO") == 0 ); 464 ROArrayColumn<Float> v(tab, "SPECTRA"); 465 cout << v.shape(0) << endl; 466 return 0; 467 } 468 return 0; 469 } 470 471 Table Scantable::getHistoryTable() const 1567 472 { 1568 473 return table_.keywordSet().asTable("HISTORY"); 1569 474 } 1570 475 1571 void S DMemTable::appendToHistoryTable(const Table& otherHist)476 void Scantable::appendToHistoryTable(const Table& otherHist) 1572 477 { 1573 478 Table t = table_.rwKeywordSet().asTable("HISTORY"); 1574 const String sep = "--------------------------------------------------------------------------------"; 1575 addHistory( sep);479 480 addHistory(asap::SEPERATOR); 1576 481 TableCopy::copyRows(t, otherHist, t.nrow(), 0, otherHist.nrow()); 1577 addHistory( sep);1578 } 1579 1580 void S DMemTable::addHistory(const std::string& hist)482 addHistory(asap::SEPERATOR); 483 } 484 485 void Scantable::addHistory(const std::string& hist) 1581 486 { 1582 487 Table t = table_.rwKeywordSet().asTable("HISTORY"); … … 1587 492 } 1588 493 1589 std::vector<std::string> S DMemTable::getHistory() const494 std::vector<std::string> Scantable::getHistory() const 1590 495 { 1591 496 Vector<String> history; … … 1602 507 } 1603 508 1604 1605 /* 1606 void SDMemTable::maskChannels(const std::vector<Int>& whichChans ) { 1607 1608 std::vector<int>::iterator it; 1609 ArrayAccessor<uChar, Axis<asap::PolAxis> > j(flags_); 1610 for (it = whichChans.begin(); it != whichChans.end(); it++) { 1611 j.reset(j.begin(uInt(*it))); 1612 for (ArrayAccessor<uChar, Axis<asap::BeamAxis> > i(j); i != i.end(); ++i) { 1613 for (ArrayAccessor<uChar, Axis<asap::IFAxis> > ii(i); ii != ii.end(); ++ii) { 1614 for (ArrayAccessor<uChar, Axis<asap::ChanAxis> > iii(ii); 1615 iii != iii.end(); ++iii) { 1616 (*iii) = 1617 } 1618 } 1619 } 1620 } 1621 1622 } 1623 */ 1624 void SDMemTable::flag(int whichRow) 1625 { 1626 Array<uChar> arr; 1627 flagsCol_.get(whichRow, arr); 1628 1629 ArrayAccessor<uChar, Axis<asap::BeamAxis> > aa0(arr); 1630 aa0.reset(aa0.begin(uInt(beamSel_)));//go to beam 1631 ArrayAccessor<uChar, Axis<asap::IFAxis> > aa1(aa0); 1632 aa1.reset(aa1.begin(uInt(IFSel_)));// go to IF 1633 ArrayAccessor<uChar, Axis<asap::PolAxis> > aa2(aa1); 1634 aa2.reset(aa2.begin(uInt(polSel_)));// go to pol 1635 1636 for (ArrayAccessor<uChar, Axis<asap::ChanAxis> > i(aa2); i != i.end(); ++i) { 1637 (*i) = uChar(True); 1638 } 1639 1640 flagsCol_.put(whichRow, arr); 1641 } 1642 1643 MDirection::Types SDMemTable::getDirectionReference() const 1644 { 1645 Float eq; 1646 table_.keywordSet().get("Equinox",eq); 1647 std::map<float,string> mp; 1648 mp[2000.0] = "J2000"; 1649 mp[1950.0] = "B1950"; 1650 MDirection::Types mdr; 1651 if (!MDirection::getType(mdr, mp[eq])) { 1652 mdr = MDirection::J2000; 1653 pushLog("WARNING: Unknown equinox using J2000"); 1654 } 1655 1656 return mdr; 1657 } 1658 1659 MEpoch::Types SDMemTable::getTimeReference() const 1660 { 1661 MEpoch::Types met; 1662 String ep; 1663 table_.keywordSet().get("Epoch",ep); 1664 if (!MEpoch::getType(met, ep)) { 1665 pushLog("WARNING: Epoch type unknown - using UTC"); 1666 met = MEpoch::UTC; 1667 } 1668 1669 return met; 1670 } 1671 1672 1673 Bool SDMemTable::setRestFreqs(const Vector<Double>& restFreqsIn, 1674 const String& sUnit, 1675 const vector<string>& lines, 1676 const String& source, 1677 Int whichIF) 1678 { 1679 const Int nIFs = nIF(); 1680 if (whichIF>=nIFs) { 1681 throw(AipsError("Illegal IF index")); 1682 } 1683 1684 // Find vector of restfrequencies 1685 // Double takes precedence over String 1686 Unit unit; 1687 Vector<Double> restFreqs; 1688 if (restFreqsIn.nelements()>0) { 1689 restFreqs.resize(restFreqsIn.nelements()); 1690 restFreqs = restFreqsIn; 1691 unit = Unit(sUnit); 1692 } else if (lines.size()>0) { 1693 const uInt nLines = lines.size(); 1694 unit = Unit("Hz"); 1695 restFreqs.resize(nLines); 1696 MFrequency lineFreq; 1697 for (uInt i=0; i<nLines; i++) { 1698 String tS(lines[i]); 1699 tS.upcase(); 1700 if (MeasTable::Line(lineFreq, tS)) { 1701 restFreqs[i] = lineFreq.getValue().getValue(); // Hz 1702 } else { 1703 String s = String(lines[i]) + 1704 String(" is an unrecognized spectral line"); 1705 throw(AipsError(s)); 1706 } 1707 } 1708 } else { 1709 throw(AipsError("You have not specified any rest frequencies or lines")); 1710 } 1711 1712 // If multiple restfreqs, must be length nIF. In this 1713 // case we will just replace the rest frequencies 1714 // We can't disinguish scalar and vector of length 1 1715 const uInt nRestFreqs = restFreqs.nelements(); 1716 Int idx = -1; 1717 SDFrequencyTable sdft = getSDFreqTable(); 1718 1719 ostringstream oss; 1720 if (nRestFreqs>1) { 1721 // Replace restFreqs, one per IF 1722 if (nRestFreqs != nIFs) { 1723 throw (AipsError("Number of rest frequencies must be equal to the number of IFs")); 1724 } 1725 ostringstream oss; 1726 oss << "Replaced rest frequencies, one per IF, with given list : " << restFreqs; 1727 sdft.deleteRestFrequencies(); 1728 for (uInt i=0; i<nRestFreqs; i++) { 1729 Quantum<Double> rf(restFreqs[i], unit); 1730 sdft.addRestFrequency(rf.getValue("Hz")); 1731 } 1732 } else { 1733 1734 // Add new rest freq 1735 Quantum<Double> rf(restFreqs[0], unit); 1736 idx = sdft.addRestFrequency(rf.getValue("Hz")); 1737 if (whichIF>=0) { 1738 oss << "Selected given rest frequency (" << restFreqs[0] << ") for IF " << whichIF << endl; 1739 } else { 1740 oss << "Selected given rest frequency (" << restFreqs[0] << ") for all IFs" << endl; 1741 } 1742 } 1743 pushLog(String(oss)); 1744 // Replace 1745 Bool empty = source.empty(); 1746 Bool ok = False; 1747 if (putSDFreqTable(sdft)) { 1748 const uInt nRow = table_.nrow(); 1749 String srcName; 1750 Vector<uInt> restFreqIDs; 1751 for (uInt i=0; i<nRow; i++) { 1752 srcnCol_.get(i, srcName); 1753 restfreqidCol_.get(i,restFreqIDs); 1754 if (idx==-1) { 1755 // Replace vector of restFreqs; one per IF. 1756 // No selection possible 1757 for (uInt i=0; i<nIFs; i++) restFreqIDs[i] = i; 1758 } else { 1759 // Set RestFreqID for selected data 1760 if (empty || source==srcName) { 1761 if (whichIF<0) { 1762 restFreqIDs = idx; 1763 } else { 1764 restFreqIDs[whichIF] = idx; 1765 } 1766 } 1767 } 1768 restfreqidCol_.put(i,restFreqIDs); 1769 } 1770 ok = True; 1771 } else { 1772 ok = False; 1773 } 1774 1775 return ok; 1776 } 1777 1778 std::string SDMemTable::spectralLines() const 1779 { 1780 Vector<String> lines = MeasTable::Lines(); 1781 MFrequency lineFreq; 1782 Double freq; 1783 ostringstream oss; 1784 1785 oss.flags(std::ios_base::left); 1786 oss << "Line Frequency (Hz)" << endl; 1787 oss << "-----------------------" << endl; 1788 for (uInt i=0; i<lines.nelements(); i++) { 1789 MeasTable::Line(lineFreq, lines[i]); 1790 freq = lineFreq.getValue().getValue(); // Hz 1791 oss << setw(11) << lines[i] << setprecision(10) << freq << endl; 1792 } 1793 return String(oss); 1794 } 1795 1796 void SDMemTable::renumber() 1797 { 1798 uInt nRow = scanCol_.nrow(); 1799 Int newscanid = 0; 1800 Int cIdx;// the current scanid 1801 // get the first scanid 1802 scanCol_.getScalar(0,cIdx); 1803 Int pIdx = cIdx;// the scanid of the previous row 1804 for (uInt i=0; i<nRow;++i) { 1805 scanCol_.getScalar(i,cIdx); 1806 if (pIdx == cIdx) { 1807 // renumber 1808 scanCol_.put(i,newscanid); 1809 } else { 1810 ++newscanid; 1811 pIdx = cIdx; // store scanid 1812 --i; // don't increment next loop 1813 } 1814 } 1815 } 1816 1817 1818 void SDMemTable::getCursorSlice(IPosition& start, IPosition& end, 1819 const IPosition& shape) const 1820 { 1821 const uInt nDim = shape.nelements(); 1822 start.resize(nDim); 1823 end.resize(nDim); 1824 1825 start(asap::BeamAxis) = beamSel_; 1826 end(asap::BeamAxis) = beamSel_; 1827 start(asap::IFAxis) = IFSel_; 1828 end(asap::IFAxis) = IFSel_; 1829 1830 start(asap::PolAxis) = polSel_; 1831 end(asap::PolAxis) = polSel_; 1832 1833 start(asap::ChanAxis) = 0; 1834 end(asap::ChanAxis) = shape(asap::ChanAxis) - 1; 1835 } 1836 1837 1838 std::vector<float> SDMemTable::getFloatSpectrum(const Array<Float>& arr) const 1839 // Get spectrum at cursor location 1840 { 1841 1842 // Setup accessors 1843 ArrayAccessor<Float, Axis<asap::BeamAxis> > aa0(arr); 1844 aa0.reset(aa0.begin(uInt(beamSel_))); // Beam selection 1845 1846 ArrayAccessor<Float, Axis<asap::IFAxis> > aa1(aa0); 1847 aa1.reset(aa1.begin(uInt(IFSel_))); // IF selection 1848 1849 ArrayAccessor<Float, Axis<asap::PolAxis> > aa2(aa1); 1850 aa2.reset(aa2.begin(uInt(polSel_))); // Pol selection 1851 1852 std::vector<float> spectrum; 1853 for (ArrayAccessor<Float, Axis<asap::ChanAxis> > i(aa2); i != i.end(); ++i) { 1854 spectrum.push_back(*i); 1855 } 1856 return spectrum; 1857 } 1858 1859 void SDMemTable::calculateAZEL() 509 std::string Scantable::formatTime(const MEpoch& me, bool showdate) const 510 { 511 MVTime mvt(me.getValue()); 512 if (showdate) 513 mvt.setFormat(MVTime::YMD); 514 else 515 mvt.setFormat(MVTime::TIME); 516 ostringstream oss; 517 oss << mvt; 518 return String(oss); 519 } 520 521 void Scantable::calculateAZEL() 1860 522 { 1861 523 MPosition mp = getAntennaPosition(); 524 MEpoch::ROScalarColumn timeCol(table_, "TIME"); 1862 525 ostringstream oss; 1863 526 oss << "Computed azimuth/elevation using " << endl 1864 527 << mp << endl; 1865 for (uInt i=0; i<nRow();++i) { 1866 MEpoch me = getEpoch(i); 1867 MDirection md = getDirection(i,False); 1868 oss << " Time: " << getTime(i,False) << " Direction: " << formatDirection(md) 528 for (uInt i=0; i<nrow(); ++i) { 529 MEpoch me = timeCol(i); 530 MDirection md = dirCol_(i); 531 dirCol_.get(i,md); 532 oss << " Time: " << formatTime(me,False) << " Direction: " << formatDirection(md) 1869 533 << endl << " => "; 1870 534 MeasFrame frame(mp, me); … … 1880 544 pushLog(String(oss)); 1881 545 } 546 547 double Scantable::getInterval(int whichrow) const 548 { 549 if (whichrow < 0) return 0.0; 550 Double intval; 551 integrCol_.get(Int(whichrow), intval); 552 return intval; 553 } 554 555 std::vector<bool> Scantable::getMask(int whichrow) const 556 { 557 Vector<uChar> flags; 558 flagsCol_.get(uInt(whichrow), flags); 559 Vector<Bool> bflag(flags.shape()); 560 convertArray(bflag, flags); 561 bflag = !bflag; 562 std::vector<bool> mask; 563 bflag.tovector(mask); 564 return mask; 565 } 566 567 std::vector<float> Scantable::getSpectrum(int whichrow) const 568 { 569 Vector<Float> arr; 570 specCol_.get(whichrow, arr); 571 std::vector<float> out; 572 arr.tovector(out); 573 return out; 574 } 575 576 String Scantable::generateName() 577 { 578 return (File::newUniqueName("./","temp")).baseName(); 579 } 580 581 const casa::Table& Scantable::table( ) const 582 { 583 return table_; 584 } 585 586 casa::Table& Scantable::table( ) 587 { 588 return table_; 589 } 590 591 std::string Scantable::getPolarizationLabel(bool linear, bool stokes, 592 bool linpol, int polidx) const 593 { 594 uInt idx = 0; 595 if (polidx >=0) idx = polidx; 596 return ""; 597 //return SDPolUtil::polarizationLabel(idx, linear, stokes, linpol); 598 } 599 600 void Scantable::unsetSelection() 601 { 602 table_ = originalTable_; 603 selector_.reset(); 604 } 605 606 void Scantable::setSelection( const STSelector& selection ) 607 { 608 Table tab = const_cast<STSelector&>(selection).apply(originalTable_); 609 if ( tab.nrow() == 0 ) { 610 throw(AipsError("Selection contains no data. Not applying it.")); 611 } 612 table_ = tab; 613 selector_ = selection; 614 } 615 616 std::string asap::Scantable::summary( bool verbose ) 617 { 618 // Format header info 619 ostringstream oss; 620 oss << endl; 621 oss << asap::SEPERATOR << endl; 622 oss << " Scan Table Summary" << endl; 623 oss << asap::SEPERATOR << endl; 624 oss.flags(std::ios_base::left); 625 oss << setw(15) << "Beams:" << setw(4) << nbeam() << endl 626 << setw(15) << "IFs:" << setw(4) << nif() << endl 627 << setw(15) << "Polarisations:" << setw(4) << npol() << endl 628 << setw(15) << "Channels:" << setw(4) << nchan() << endl; 629 oss << endl; 630 String tmp; 631 //table_.keywordSet().get("Observer", tmp); 632 oss << setw(15) << "Observer:" << table_.keywordSet().asString("Observer") << endl; 633 oss << setw(15) << "Obs Date:" << getTime(-1,true) << endl; 634 table_.keywordSet().get("Project", tmp); 635 oss << setw(15) << "Project:" << tmp << endl; 636 table_.keywordSet().get("Obstype", tmp); 637 oss << setw(15) << "Obs. Type:" << tmp << endl; 638 table_.keywordSet().get("AntennaName", tmp); 639 oss << setw(15) << "Antenna Name:" << tmp << endl; 640 table_.keywordSet().get("FluxUnit", tmp); 641 oss << setw(15) << "Flux Unit:" << tmp << endl; 642 Vector<Float> vec; 643 oss << setw(15) << "Rest Freqs:"; 644 if (vec.nelements() > 0) { 645 oss << setprecision(10) << vec << " [Hz]" << endl; 646 } else { 647 oss << "none" << endl; 648 } 649 oss << setw(15) << "Abcissa:" << "channel" << endl; 650 oss << selector_.print() << endl; 651 oss << endl; 652 // main table 653 String dirtype = "Position (" 654 + MDirection::showType(dirCol_.getMeasRef().getType()) 655 + ")"; 656 oss << setw(5) << "Scan" 657 << setw(15) << "Source" 658 // << setw(24) << dirtype 659 << setw(10) << "Time" 660 << setw(18) << "Integration" << endl 661 << setw(5) << "" << setw(10) << "Beam" << dirtype << endl 662 << setw(15) << "" << setw(5) << "IF" 663 << setw(8) << "Frame" << setw(16) 664 << "RefVal" << setw(10) << "RefPix" << setw(12) << "Increment" <<endl; 665 oss << asap::SEPERATOR << endl; 666 TableIterator iter(table_, "SCANNO"); 667 while (!iter.pastEnd()) { 668 Table subt = iter.table(); 669 ROTableRow row(subt); 670 MEpoch::ROScalarColumn timeCol(subt,"TIME"); 671 const TableRecord& rec = row.get(0); 672 oss << setw(4) << std::right << rec.asuInt("SCANNO") 673 << std::left << setw(1) << "" 674 << setw(15) << rec.asString("SRCNAME") 675 << setw(10) << formatTime(timeCol(0), false); 676 // count the cycles in the scan 677 TableIterator cyciter(subt, "CYCLENO"); 678 int nint = 0; 679 while (!cyciter.pastEnd()) { 680 ++nint; 681 ++cyciter; 682 } 683 oss << setw(3) << std::right << nint << setw(3) << " x " << std::left 684 << setw(6) << formatSec(rec.asFloat("INTERVAL")) << endl; 685 686 TableIterator biter(subt, "BEAMNO"); 687 while (!biter.pastEnd()) { 688 Table bsubt = biter.table(); 689 ROTableRow brow(bsubt); 690 MDirection::ROScalarColumn bdirCol(bsubt,"DIRECTION"); 691 const TableRecord& brec = brow.get(0); 692 oss << setw(6) << "" << setw(10) << brec.asuInt("BEAMNO"); 693 oss << setw(24) << formatDirection(bdirCol(0)) << endl; 694 TableIterator iiter(bsubt, "IFNO"); 695 while (!iiter.pastEnd()) { 696 Table isubt = iiter.table(); 697 ROTableRow irow(isubt); 698 const TableRecord& irec = irow.get(0); 699 oss << std::right <<setw(8) << "" << std::left << irec.asuInt("IFNO"); 700 oss << frequencies().print(irec.asuInt("FREQ_ID")); 701 702 ++iiter; 703 } 704 ++biter; 705 } 706 ++iter; 707 } 708 /// @todo implement verbose mode 709 return String(oss); 710 } 711 712 std::string Scantable::getTime(int whichrow, bool showdate) const 713 { 714 MEpoch::ROScalarColumn timeCol(table_, "TIME"); 715 MEpoch me; 716 if (whichrow > -1) { 717 me = timeCol(uInt(whichrow)); 718 } else { 719 Double tm; 720 table_.keywordSet().get("UTC",tm); 721 me = MEpoch(MVEpoch(tm)); 722 } 723 return formatTime(me, showdate); 724 } 725 726 }//namespace asap
Note: See TracChangeset
for help on using the changeset viewer.