source: trunk/src/SDMemTable.cc@ 277

Last change on this file since 277 was 275, checked in by kil064, 20 years ago

finalize implement DOPPLER setting and using. Before it was
just set to RADIO and never changed. Now you can set it
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 34.2 KB
Line 
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
32#include <map>
33
34#include <casa/aips.h>
35#include <casa/iostream.h>
36#include <casa/iomanip.h>
37#include <casa/Arrays/Array.h>
38#include <casa/Arrays/ArrayMath.h>
39#include <casa/Arrays/MaskArrMath.h>
40#include <casa/Arrays/ArrayLogical.h>
41#include <casa/Arrays/ArrayAccessor.h>
42#include <casa/Arrays/Vector.h>
43
44#include <tables/Tables/TableParse.h>
45#include <tables/Tables/TableDesc.h>
46#include <tables/Tables/SetupNewTab.h>
47#include <tables/Tables/ScaColDesc.h>
48#include <tables/Tables/ArrColDesc.h>
49
50#include <tables/Tables/ExprNode.h>
51#include <tables/Tables/ScalarColumn.h>
52#include <tables/Tables/ArrayColumn.h>
53#include <tables/Tables/TableRecord.h>
54#include <measures/Measures/MFrequency.h>
55#include <measures/Measures/MeasTable.h>
56#include <coordinates/Coordinates/CoordinateUtil.h>
57#include <casa/Quanta/MVTime.h>
58
59#include "SDDefs.h"
60#include "SDMemTable.h"
61#include "SDContainer.h"
62
63
64using namespace casa;
65using namespace asap;
66
67SDMemTable::SDMemTable() :
68 IFSel_(0),
69 beamSel_(0),
70 polSel_(0)
71{
72 setup();
73}
74
75SDMemTable::SDMemTable(const std::string& name) :
76 IFSel_(0),
77 beamSel_(0),
78 polSel_(0)
79{
80 Table tab(name);
81 table_ = tab.copyToMemoryTable("dummy");
82 //cerr << "hello from C SDMemTable @ " << this << endl;
83}
84
85SDMemTable::SDMemTable(const SDMemTable& other, Bool clear)
86{
87 IFSel_= other.IFSel_;
88 beamSel_= other.beamSel_;
89 polSel_= other.polSel_;
90 chanMask_ = other.chanMask_;
91 table_ = other.table_.copyToMemoryTable(String("dummy"));
92 // clear all rows()
93 if (clear) {
94 table_.removeRow(this->table_.rowNumbers());
95 } else {
96 IFSel_ = other.IFSel_;
97 beamSel_ = other.beamSel_;
98 polSel_ = other.polSel_;
99 }
100 //cerr << "hello from CC SDMemTable @ " << this << endl;
101}
102
103SDMemTable::SDMemTable(const Table& tab, const std::string& exprs) :
104 IFSel_(0),
105 beamSel_(0),
106 polSel_(0)
107{
108 Table t = tableCommand(exprs,tab);
109 if (t.nrow() == 0)
110 throw(AipsError("Query unsuccessful."));
111 table_ = t.copyToMemoryTable("dummy");
112}
113
114SDMemTable::~SDMemTable()
115{
116 //cerr << "goodbye from SDMemTable @ " << this << endl;
117}
118
119SDMemTable SDMemTable::getScan(Int scanID) const
120{
121 String cond("SELECT * from $1 WHERE SCANID == ");
122 cond += String::toString(scanID);
123 return SDMemTable(table_, cond);
124}
125
126SDMemTable &SDMemTable::operator=(const SDMemTable& other)
127{
128 if (this != &other) {
129 IFSel_= other.IFSel_;
130 beamSel_= other.beamSel_;
131 polSel_= other.polSel_;
132 chanMask_.resize(0);
133 chanMask_ = other.chanMask_;
134 table_ = other.table_.copyToMemoryTable(String("dummy"));
135 }
136 //cerr << "hello from ASS SDMemTable @ " << this << endl;
137 return *this;
138}
139
140SDMemTable SDMemTable::getSource(const std::string& source) const
141{
142 String cond("SELECT * from $1 WHERE SRCNAME == ");
143 cond += source;
144 return SDMemTable(table_, cond);
145}
146
147void SDMemTable::setup()
148{
149 TableDesc td("", "1", TableDesc::Scratch);
150 td.comment() = "A SDMemTable";
151 td.addColumn(ScalarColumnDesc<Double>("TIME"));
152 td.addColumn(ScalarColumnDesc<String>("SRCNAME"));
153 td.addColumn(ArrayColumnDesc<Float>("SPECTRA"));
154 td.addColumn(ArrayColumnDesc<uChar>("FLAGTRA"));
155 td.addColumn(ArrayColumnDesc<Float>("TSYS"));
156 td.addColumn(ScalarColumnDesc<Int>("SCANID"));
157 td.addColumn(ScalarColumnDesc<Double>("INTERVAL"));
158 td.addColumn(ArrayColumnDesc<uInt>("FREQID"));
159 td.addColumn(ArrayColumnDesc<Double>("DIRECTION"));
160 td.addColumn(ScalarColumnDesc<String>("FIELDNAME"));
161 td.addColumn(ScalarColumnDesc<String>("TCALTIME"));
162 td.addColumn(ArrayColumnDesc<Float>("TCAL"));
163 td.addColumn(ScalarColumnDesc<Float>("AZIMUTH"));
164 td.addColumn(ScalarColumnDesc<Float>("ELEVATION"));
165 td.addColumn(ScalarColumnDesc<Float>("PARANGLE"));
166 td.addColumn(ScalarColumnDesc<Int>("REFBEAM"));
167 td.addColumn(ArrayColumnDesc<String>("HISTORY"));
168
169 // Now create a new table from the description.
170
171 SetupNewTable aNewTab("dummy", td, Table::New);
172 table_ = Table(aNewTab, Table::Memory, 0);
173}
174
175std::string SDMemTable::getSourceName(Int whichRow) const
176{
177 ROScalarColumn<String> src(table_, "SRCNAME");
178 String name;
179 src.get(whichRow, name);
180 return name;
181}
182
183std::string SDMemTable::getTime(Int whichRow) const
184{
185 ROScalarColumn<Double> src(table_, "TIME");
186 Double tm;
187 src.get(whichRow, tm);
188 MVTime mvt(tm);
189 mvt.setFormat(MVTime::TIME);
190 ostringstream oss;
191 oss << mvt;
192 String str(oss);
193 return str;
194}
195double SDMemTable::getInterval(Int whichRow) const
196{
197 ROScalarColumn<Double> src(table_, "INTERVAL");
198 Double intval;
199 src.get(whichRow, intval);
200 return intval;
201}
202
203bool SDMemTable::setIF(Int whichIF)
204{
205 if ( whichIF >= 0 && whichIF < nIF()) {
206 IFSel_ = whichIF;
207 return true;
208 }
209 return false;
210}
211
212bool SDMemTable::setBeam(Int whichBeam)
213{
214 if ( whichBeam >= 0 && whichBeam < nBeam()) {
215 beamSel_ = whichBeam;
216 return true;
217 }
218 return false;
219}
220
221bool SDMemTable::setPol(Int whichPol)
222{
223 if ( whichPol >= 0 && whichPol < nPol()) {
224 polSel_ = whichPol;
225 return true;
226 }
227 return false;
228}
229
230bool SDMemTable::setMask(std::vector<int> whichChans)
231{
232 ROArrayColumn<uChar> spec(table_, "FLAGTRA");
233 std::vector<int>::iterator it;
234 uInt n = spec.shape(0)(3);
235 if (whichChans.empty()) {
236 chanMask_ = std::vector<bool>(n,true);
237 return true;
238 }
239 chanMask_.resize(n,true);
240 for (it = whichChans.begin(); it != whichChans.end(); ++it) {
241 if (*it < n) {
242 chanMask_[*it] = false;
243 }
244 }
245 return true;
246}
247
248std::vector<bool> SDMemTable::getMask(Int whichRow) const {
249 std::vector<bool> mask;
250 ROArrayColumn<uChar> spec(table_, "FLAGTRA");
251 Array<uChar> arr;
252 spec.get(whichRow, arr);
253 ArrayAccessor<uChar, Axis<asap::BeamAxis> > aa0(arr);
254 aa0.reset(aa0.begin(uInt(beamSel_)));//go to beam
255 ArrayAccessor<uChar, Axis<asap::IFAxis> > aa1(aa0);
256 aa1.reset(aa1.begin(uInt(IFSel_)));// go to IF
257 ArrayAccessor<uChar, Axis<asap::PolAxis> > aa2(aa1);
258 aa2.reset(aa2.begin(uInt(polSel_)));// go to pol
259
260 Bool useUserMask = ( chanMask_.size() == arr.shape()(3) );
261
262 std::vector<bool> tmp;
263 tmp = chanMask_; // WHY the fxxx do I have to make a copy here
264 std::vector<bool>::iterator miter;
265 miter = tmp.begin();
266
267 for (ArrayAccessor<uChar, Axis<asap::ChanAxis> > i(aa2); i != i.end(); ++i) {
268 bool out =!static_cast<bool>(*i);
269 if (useUserMask) {
270 out = out && (*miter);
271 miter++;
272 }
273 mask.push_back(out);
274 }
275 return mask;
276}
277std::vector<float> SDMemTable::getSpectrum(Int whichRow) const
278{
279 std::vector<float> spectrum;
280 ROArrayColumn<Float> spec(table_, "SPECTRA");
281 Array<Float> arr;
282 spec.get(whichRow, arr);
283 ArrayAccessor<Float, Axis<asap::BeamAxis> > aa0(arr);
284 aa0.reset(aa0.begin(uInt(beamSel_)));//go to beam
285 ArrayAccessor<Float, Axis<asap::IFAxis> > aa1(aa0);
286 aa1.reset(aa1.begin(uInt(IFSel_)));// go to IF
287 ArrayAccessor<Float, Axis<asap::PolAxis> > aa2(aa1);
288 aa2.reset(aa2.begin(uInt(polSel_)));// go to pol
289 for (ArrayAccessor<Float, Axis<asap::ChanAxis> > i(aa2); i != i.end(); ++i) {
290 spectrum.push_back(*i);
291 }
292 return spectrum;
293}
294std::vector<string> SDMemTable::getCoordInfo() const
295{
296 String un;
297 Table t = table_.keywordSet().asTable("FREQUENCIES");
298 String sunit;
299 t.keywordSet().get("UNIT",sunit);
300 String dpl;
301 t.keywordSet().get("DOPPLER",dpl);
302 if (dpl == "") dpl = "RADIO";
303 String rfrm;
304 t.keywordSet().get("REFFRAME",rfrm);
305 std::vector<string> inf;
306 inf.push_back(sunit);
307 inf.push_back(rfrm);
308 inf.push_back(dpl);
309 t.keywordSet().get("BASEREFFRAME",rfrm);
310 inf.push_back(rfrm);
311 return inf;
312}
313
314void SDMemTable::setCoordInfo(std::vector<string> theinfo)
315{
316 std::vector<string>::iterator it;
317 String un,rfrm,dpl;
318 un = theinfo[0];
319 rfrm = theinfo[1];
320 dpl = theinfo[2];
321
322 Table t = table_.rwKeywordSet().asTable("FREQUENCIES");
323 Vector<Double> rstf;
324 t.keywordSet().get("RESTFREQS",rstf);
325 Bool canDo = True;
326 Unit u1("km/s");Unit u2("Hz");
327 if (Unit(un) == u1) {
328 Vector<Double> rstf;
329 t.keywordSet().get("RESTFREQS",rstf);
330 if (rstf.nelements() == 0) {
331 throw(AipsError("Can't set unit to km/s if no restfrequencies are specified"));
332 }
333 } else if (Unit(un) != u2 && un != "") {
334 throw(AipsError("Unit not conformant with Spectral Coordinates"));
335 }
336 t.rwKeywordSet().define("UNIT", un);
337//
338 MFrequency::Types mdr;
339 if (!MFrequency::getType(mdr, rfrm)) {
340
341 Int a,b;const uInt* c;
342 const String* valid = MFrequency::allMyTypes(a, b, c);
343 String pfix = "Please specify a legal frame type. Types are\n";
344 throw(AipsError(pfix+(*valid)));
345 } else {
346 t.rwKeywordSet().define("REFFRAME",rfrm);
347 }
348//
349 MDoppler::Types dtype;
350 dpl.upcase();
351 if (!MDoppler::getType(dtype, dpl)) {
352 throw(AipsError("Doppler type unknown"));
353 } else {
354 t.rwKeywordSet().define("DOPPLER",dpl);
355 }
356}
357
358std::vector<double> SDMemTable::getAbcissa(Int whichRow) const
359{
360 std::vector<double> absc(nChan());
361 Vector<Double> absc1(nChan());
362 indgen(absc1);
363 ROArrayColumn<uInt> fid(table_, "FREQID");
364 Vector<uInt> v;
365 fid.get(whichRow, v);
366 uInt specidx = v(IFSel_);
367 SpectralCoordinate spc = getCoordinate(specidx);
368 Table t = table_.keywordSet().asTable("FREQUENCIES");
369
370 MDirection direct = getDirection(whichRow);
371
372 ROScalarColumn<Double> tme(table_, "TIME");
373 Double obstime;
374 tme.get(whichRow,obstime);
375 MVEpoch tm2(Quantum<Double>(obstime, Unit(String("d"))));
376 MEpoch::Types met = getTimeReference();
377 MEpoch epoch(tm2, met);
378
379 Vector<Double> antpos;
380 table_.keywordSet().get("AntennaPosition", antpos);
381 MVPosition mvpos(antpos(0),antpos(1),antpos(2));
382 MPosition pos(mvpos);
383 String sunit;
384 t.keywordSet().get("UNIT",sunit);
385 if (sunit == "") sunit = "pixel";
386 Unit u(sunit);
387 String frm;
388 t.keywordSet().get("REFFRAME",frm);
389 if (frm == "") frm = "TOPO";
390 String dpl;
391 t.keywordSet().get("DOPPLER",dpl);
392 if (dpl == "") dpl = "RADIO";
393 MFrequency::Types mtype;
394 if (!MFrequency::getType(mtype, frm)) {
395 cout << "Frequency type unknown assuming TOPO" << endl; // SHould never happen
396 mtype = MFrequency::TOPO;
397 }
398 MDoppler::Types dtype;
399 if (!MDoppler::getType(dtype, dpl)) {
400 cout << "Doppler type unknown assuming RADIO" << endl; // SHould never happen
401 dtype = MDoppler::RADIO;
402 }
403
404 if (!spc.setReferenceConversion(mtype,epoch,pos,direct)) {
405 throw(AipsError("Couldn't convert frequency frame."));
406 }
407
408 if ( u == Unit("km/s") ) {
409 Vector<Double> rstf;
410 t.keywordSet().get("RESTFREQS",rstf);
411 if (rstf.nelements() > 0) {
412 if (rstf.nelements() >= nIF())
413 spc.selectRestFrequency(uInt(IFSel_));
414 spc.setVelocity(u.getName(),dtype);
415 Vector<Double> wrld;
416 spc.pixelToVelocity(wrld,absc1);
417 std::vector<double>::iterator it;
418 uInt i = 0;
419 for (it = absc.begin(); it != absc.end(); ++it) {
420 (*it) = wrld[i];
421 i++;
422 }
423 }
424 } else if (u == Unit("Hz")) {
425 Vector<String> wau(1); wau = u.getName();
426 spc.setWorldAxisUnits(wau);
427 std::vector<double>::iterator it;
428 Double tmp;
429 uInt i = 0;
430 for (it = absc.begin(); it != absc.end(); ++it) {
431 spc.toWorld(tmp,absc1[i]);
432 (*it) = tmp;
433 i++;
434 }
435
436 } else {
437 // assume channels/pixels
438 std::vector<double>::iterator it;
439 uInt i=0;
440 for (it = absc.begin(); it != absc.end(); ++it) {
441 (*it) = Double(i++);
442 }
443 }
444 return absc;
445}
446
447std::string SDMemTable::getAbcissaString(Int whichRow) const
448{
449 ROArrayColumn<uInt> fid(table_, "FREQID");
450 Table t = table_.keywordSet().asTable("FREQUENCIES");
451 String sunit;
452 t.keywordSet().get("UNIT",sunit);
453 if (sunit == "") sunit = "pixel";
454 Unit u(sunit);
455 Vector<uInt> v;
456 fid.get(whichRow, v);
457 uInt specidx = v(IFSel_);
458 SpectralCoordinate spc = getCoordinate(specidx);
459 String frm;
460 t.keywordSet().get("REFFRAME",frm);
461//
462 MFrequency::Types mtype;
463 if (!MFrequency::getType(mtype, frm)) {
464 cout << "Frequency type unknown assuming TOPO" << endl;
465 mtype = MFrequency::TOPO;
466 }
467 spc.setFrequencySystem(mtype);
468//
469 String dpl;
470 t.keywordSet().get("DOPPLER",dpl);
471 MDoppler::Types dtype;
472 MDoppler::getType(dtype, dpl); // Can't fail
473//
474 String s = "Channel";
475 if (u == Unit("km/s")) {
476 spc.setVelocity(u.getName(), dtype);
477 s = CoordinateUtil::axisLabel(spc,0,True,True,True);
478 } else if (u == Unit("Hz")) {
479 Vector<String> wau(1);wau = u.getName();
480 spc.setWorldAxisUnits(wau);
481 s = CoordinateUtil::axisLabel(spc);
482 }
483 return s;
484}
485
486void SDMemTable::setSpectrum(std::vector<float> spectrum, int whichRow)
487{
488 ArrayColumn<Float> spec(table_, "SPECTRA");
489 Array<Float> arr;
490 spec.get(whichRow, arr);
491 if (spectrum.size() != arr.shape()(3)) {
492 throw(AipsError("Attempting to set spectrum with incorrect length."));
493 }
494
495 ArrayAccessor<Float, Axis<asap::BeamAxis> > aa0(arr);
496 aa0.reset(aa0.begin(uInt(beamSel_)));//go to beam
497 ArrayAccessor<Float, Axis<asap::IFAxis> > aa1(aa0);
498 aa1.reset(aa1.begin(uInt(IFSel_)));// go to IF
499 ArrayAccessor<Float, Axis<asap::PolAxis> > aa2(aa1);
500 aa2.reset(aa2.begin(uInt(polSel_)));// go to pol
501
502 std::vector<float>::iterator it = spectrum.begin();
503 for (ArrayAccessor<Float, Axis<asap::ChanAxis> > i(aa2); i != i.end(); ++i) {
504 (*i) = Float(*it);
505 it++;
506 }
507 spec.put(whichRow, arr);
508}
509
510void SDMemTable::getSpectrum(Vector<Float>& spectrum, Int whichRow) const
511{
512 ROArrayColumn<Float> spec(table_, "SPECTRA");
513 Array<Float> arr;
514 spec.get(whichRow, arr);
515 spectrum.resize(arr.shape()(3));
516 ArrayAccessor<Float, Axis<asap::BeamAxis> > aa0(arr);
517 aa0.reset(aa0.begin(uInt(beamSel_)));//go to beam
518 ArrayAccessor<Float, Axis<asap::IFAxis> > aa1(aa0);
519 aa1.reset(aa1.begin(uInt(IFSel_)));// go to IF
520 ArrayAccessor<Float, Axis<asap::PolAxis> > aa2(aa1);
521 aa2.reset(aa2.begin(uInt(polSel_)));// go to pol
522
523 ArrayAccessor<Float, Axis<asap::BeamAxis> > va(spectrum);
524 for (ArrayAccessor<Float, Axis<asap::ChanAxis> > i(aa2); i != i.end(); ++i) {
525 (*va) = (*i);
526 va++;
527 }
528}
529/*
530void SDMemTable::getMask(Vector<Bool>& mask, Int whichRow) const {
531 ROArrayColumn<uChar> spec(table_, "FLAGTRA");
532 Array<uChar> arr;
533 spec.get(whichRow, arr);
534 mask.resize(arr.shape()(3));
535
536 ArrayAccessor<uChar, Axis<asap::BeamAxis> > aa0(arr);
537 aa0.reset(aa0.begin(uInt(beamSel_)));//go to beam
538 ArrayAccessor<uChar, Axis<asap::IFAxis> > aa1(aa0);
539 aa1.reset(aa1.begin(uInt(IFSel_)));// go to IF
540 ArrayAccessor<uChar, Axis<asap::PolAxis> > aa2(aa1);
541 aa2.reset(aa2.begin(uInt(polSel_)));// go to pol
542
543 Bool useUserMask = ( chanMask_.size() == arr.shape()(3) );
544
545 ArrayAccessor<Bool, Axis<asap::BeamAxis> > va(mask);
546 std::vector<bool> tmp;
547 tmp = chanMask_; // WHY the fxxx do I have to make a copy here. The
548 // iterator should work on chanMask_??
549 std::vector<bool>::iterator miter;
550 miter = tmp.begin();
551
552 for (ArrayAccessor<uChar, Axis<asap::ChanAxis> > i(aa2); i != i.end(); ++i) {
553 bool out =!static_cast<bool>(*i);
554 if (useUserMask) {
555 out = out && (*miter);
556 miter++;
557 }
558 (*va) = out;
559 va++;
560 }
561}
562*/
563MaskedArray<Float> SDMemTable::rowAsMaskedArray(uInt whichRow,
564 Bool useSelection) const
565{
566 ROArrayColumn<Float> spec(table_, "SPECTRA");
567 Array<Float> arr;
568 ROArrayColumn<uChar> flag(table_, "FLAGTRA");
569 Array<uChar> farr;
570 spec.get(whichRow, arr);
571 flag.get(whichRow, farr);
572 Array<Bool> barr(farr.shape());convertArray(barr, farr);
573 MaskedArray<Float> marr;
574 if (useSelection) {
575 ArrayAccessor<Float, Axis<asap::BeamAxis> > aa0(arr);
576 aa0.reset(aa0.begin(uInt(beamSel_)));//go to beam
577 ArrayAccessor<Float, Axis<asap::IFAxis> > aa1(aa0);
578 aa1.reset(aa1.begin(uInt(IFSel_)));// go to IF
579 ArrayAccessor<Float, Axis<asap::PolAxis> > aa2(aa1);
580 aa2.reset(aa2.begin(uInt(polSel_)));// go to pol
581
582 ArrayAccessor<Bool, Axis<asap::BeamAxis> > baa0(barr);
583 baa0.reset(baa0.begin(uInt(beamSel_)));//go to beam
584 ArrayAccessor<Bool, Axis<asap::IFAxis> > baa1(baa0);
585 baa1.reset(baa1.begin(uInt(IFSel_)));// go to IF
586 ArrayAccessor<Bool, Axis<asap::PolAxis> > baa2(baa1);
587 baa2.reset(baa2.begin(uInt(polSel_)));// go to pol
588
589 Vector<Float> a(arr.shape()(3));
590 Vector<Bool> b(barr.shape()(3));
591 ArrayAccessor<Float, Axis<asap::BeamAxis> > a0(a);
592 ArrayAccessor<Bool, Axis<asap::BeamAxis> > b0(b);
593
594 ArrayAccessor<Bool, Axis<asap::ChanAxis> > j(baa2);
595 for (ArrayAccessor<Float, Axis<asap::ChanAxis> > i(aa2);
596 i != i.end(); ++i) {
597 (*a0) = (*i);
598 (*b0) = !(*j);
599 j++;
600 a0++;
601 b0++;
602 }
603 marr.setData(a,b);
604 } else {
605 marr.setData(arr,!barr);
606 }
607 return marr;
608}
609
610Float SDMemTable::getTsys(Int whichRow) const
611{
612 ROArrayColumn<Float> ts(table_, "TSYS");
613 Array<Float> arr;
614 ts.get(whichRow, arr);
615 Float out;
616 IPosition ip(arr.shape());
617 ip(0) = beamSel_;ip(1) = IFSel_;ip(2) = polSel_;ip(3)=0;
618 out = arr(ip);
619 return out;
620}
621
622MDirection SDMemTable::getDirection(Int whichRow) const
623{
624 MDirection::Types mdr = getDirectionReference();
625 ROArrayColumn<Double> dir(table_, "DIRECTION");
626 Array<Double> posit;
627 dir.get(whichRow,posit);
628 Vector<Double> wpos(2);
629 wpos[0] = posit(IPosition(2,beamSel_,0));
630 wpos[1] = posit(IPosition(2,beamSel_,1));
631 Quantum<Double> lon(wpos[0],Unit(String("rad")));
632 Quantum<Double> lat(wpos[1],Unit(String("rad")));
633 MDirection direct(lon, lat, mdr);
634 return direct;
635}
636
637SpectralCoordinate SDMemTable::getCoordinate(uInt whichIdx) const
638{
639
640 Table t = table_.keywordSet().asTable("FREQUENCIES");
641 if (whichIdx > t.nrow() ) {
642 throw(AipsError("SDMemTable::getCoordinate - whichIdx out of range"));
643 }
644
645 Double rp,rv,inc;
646 String rf;
647 Vector<Double> vec;
648 ROScalarColumn<Double> rpc(t, "REFPIX");
649 ROScalarColumn<Double> rvc(t, "REFVAL");
650 ROScalarColumn<Double> incc(t, "INCREMENT");
651 t.keywordSet().get("RESTFREQS",vec);
652 t.keywordSet().get("BASEREFFRAME",rf);
653
654 MFrequency::Types mft;
655 if (!MFrequency::getType(mft, rf)) {
656 cerr << "Frequency type unknown assuming TOPO" << endl;
657 mft = MFrequency::TOPO;
658 }
659 rpc.get(whichIdx, rp);
660 rvc.get(whichIdx, rv);
661 incc.get(whichIdx, inc);
662 SpectralCoordinate spec(mft,rv,inc,rp);
663 if (vec.nelements() > 0)
664 spec.setRestFrequencies(vec);
665 return spec;
666}
667
668Bool SDMemTable::setCoordinate(const SpectralCoordinate& speccord,
669 uInt whichIdx) {
670 Table t = table_.rwKeywordSet().asTable("FREQUENCIES");
671 if (whichIdx > t.nrow() ) {
672 throw(AipsError("SDMemTable::setCoordinate - coord no out of range"));
673 }
674 ScalarColumn<Double> rpc(t, "REFPIX");
675 ScalarColumn<Double> rvc(t, "REFVAL");
676 ScalarColumn<Double> incc(t, "INCREMENT");
677
678 rpc.put(whichIdx, speccord.referencePixel()[0]);
679 rvc.put(whichIdx, speccord.referenceValue()[0]);
680 incc.put(whichIdx, speccord.increment()[0]);
681
682 return True;
683}
684
685Int SDMemTable::nCoordinates() const
686{
687 return table_.keywordSet().asTable("FREQUENCIES").nrow();
688}
689
690void SDMemTable::setRestFreqs(std::vector<double> freqs,
691 const std::string& theunit)
692{
693 Vector<Double> tvec(freqs);
694 Quantum<Vector<Double> > q(tvec, String(theunit));
695 tvec.resize();
696 tvec = q.getValue("Hz");
697 Table t = table_.keywordSet().asTable("FREQUENCIES");
698 t.rwKeywordSet().define("RESTFREQS",tvec);
699}
700
701std::vector<double> SDMemTable::getRestFreqs() const
702{
703 Table t = table_.keywordSet().asTable("FREQUENCIES");
704 Vector<Double> tvec;
705 t.keywordSet().get("RESTFREQS",tvec);
706 std::vector<double> stlout;
707 tvec.tovector(stlout);
708 return stlout;
709}
710
711bool SDMemTable::putSDFreqTable(const SDFrequencyTable& sdft)
712{
713 TableDesc td("", "1", TableDesc::Scratch);
714 td.addColumn(ScalarColumnDesc<Double>("REFPIX"));
715 td.addColumn(ScalarColumnDesc<Double>("REFVAL"));
716 td.addColumn(ScalarColumnDesc<Double>("INCREMENT"));
717 SetupNewTable aNewTab("freqs", td, Table::New);
718 Table aTable (aNewTab, Table::Memory, sdft.length());
719 ScalarColumn<Double> sc0(aTable, "REFPIX");
720 ScalarColumn<Double> sc1(aTable, "REFVAL");
721 ScalarColumn<Double> sc2(aTable, "INCREMENT");
722 for (uInt i=0; i < sdft.length(); ++i) {
723 sc0.put(i,sdft.referencePixel(i));
724 sc1.put(i,sdft.referenceValue(i));
725 sc2.put(i,sdft.increment(i));
726 }
727 String rf = sdft.refFrame();
728 if (rf.contains("TOPO")) rf = "TOPO";
729
730 aTable.rwKeywordSet().define("BASEREFFRAME", rf);
731 aTable.rwKeywordSet().define("REFFRAME", rf);
732 aTable.rwKeywordSet().define("EQUINOX", sdft.equinox());
733 aTable.rwKeywordSet().define("UNIT", String(""));
734 aTable.rwKeywordSet().define("DOPPLER", String("RADIO"));
735 Vector<Double> rfvec;
736 String rfunit;
737 sdft.restFrequencies(rfvec,rfunit);
738 Quantum<Vector<Double> > q(rfvec, rfunit);
739 rfvec.resize();
740 rfvec = q.getValue("Hz");
741 aTable.rwKeywordSet().define("RESTFREQS", rfvec);
742 table_.rwKeywordSet().defineTable ("FREQUENCIES", aTable);
743 return True;
744}
745
746SDFrequencyTable SDMemTable::getSDFreqTable() const
747{
748 // TODO !!!!! implement this properly USE with care
749 const Table& t = table_.keywordSet().asTable("FREQUENCIES");
750 SDFrequencyTable sdft;
751 sdft.setLength(t.nrow());
752 return sdft;
753}
754
755bool SDMemTable::putSDContainer(const SDContainer& sdc)
756{
757 ScalarColumn<Double> mjd(table_, "TIME");
758 ScalarColumn<String> srcn(table_, "SRCNAME");
759 ScalarColumn<String> fldn(table_, "FIELDNAME");
760 ArrayColumn<Float> spec(table_, "SPECTRA");
761 ArrayColumn<uChar> flags(table_, "FLAGTRA");
762 ArrayColumn<Float> ts(table_, "TSYS");
763 ScalarColumn<Int> scan(table_, "SCANID");
764 ScalarColumn<Double> integr(table_, "INTERVAL");
765 ArrayColumn<uInt> freqid(table_, "FREQID");
766 ArrayColumn<Double> dir(table_, "DIRECTION");
767 ScalarColumn<Int> rbeam(table_, "REFBEAM");
768 ArrayColumn<Float> tcal(table_, "TCAL");
769 ScalarColumn<String> tcalt(table_, "TCALTIME");
770 ScalarColumn<Float> az(table_, "AZIMUTH");
771 ScalarColumn<Float> el(table_, "ELEVATION");
772 ScalarColumn<Float> para(table_, "PARANGLE");
773 ArrayColumn<String> hist(table_, "HISTORY");
774
775 uInt rno = table_.nrow();
776 table_.addRow();
777
778 mjd.put(rno, sdc.timestamp);
779 srcn.put(rno, sdc.sourcename);
780 fldn.put(rno, sdc.fieldname);
781 spec.put(rno, sdc.getSpectrum());
782 flags.put(rno, sdc.getFlags());
783 ts.put(rno, sdc.getTsys());
784 scan.put(rno, sdc.scanid);
785 integr.put(rno, sdc.interval);
786 freqid.put(rno, sdc.getFreqMap());
787 dir.put(rno, sdc.getDirection());
788 rbeam.put(rno, sdc.refbeam);
789 tcal.put(rno, sdc.tcal);
790 tcalt.put(rno, sdc.tcaltime);
791 az.put(rno, sdc.azimuth);
792 el.put(rno, sdc.elevation);
793 para.put(rno, sdc.parangle);
794 hist.put(rno, sdc.getHistory());
795
796 return true;
797}
798
799SDContainer SDMemTable::getSDContainer(uInt whichRow) const
800{
801 ROScalarColumn<Double> mjd(table_, "TIME");
802 ROScalarColumn<String> srcn(table_, "SRCNAME");
803 ROScalarColumn<String> fldn(table_, "FIELDNAME");
804 ROArrayColumn<Float> spec(table_, "SPECTRA");
805 ROArrayColumn<uChar> flags(table_, "FLAGTRA");
806 ROArrayColumn<Float> ts(table_, "TSYS");
807 ROScalarColumn<Int> scan(table_, "SCANID");
808 ROScalarColumn<Double> integr(table_, "INTERVAL");
809 ROArrayColumn<uInt> freqid(table_, "FREQID");
810 ROArrayColumn<Double> dir(table_, "DIRECTION");
811 ROScalarColumn<Int> rbeam(table_, "REFBEAM");
812 ROArrayColumn<Float> tcal(table_, "TCAL");
813 ROScalarColumn<String> tcalt(table_, "TCALTIME");
814 ROScalarColumn<Float> az(table_, "AZIMUTH");
815 ROScalarColumn<Float> el(table_, "ELEVATION");
816 ROScalarColumn<Float> para(table_, "PARANGLE");
817 ROArrayColumn<String> hist(table_, "HISTORY");
818
819 SDContainer sdc(nBeam(),nIF(),nPol(),nChan());
820 mjd.get(whichRow, sdc.timestamp);
821 srcn.get(whichRow, sdc.sourcename);
822 integr.get(whichRow, sdc.interval);
823 scan.get(whichRow, sdc.scanid);
824 fldn.get(whichRow, sdc.fieldname);
825 rbeam.get(whichRow, sdc.refbeam);
826 az.get(whichRow, sdc.azimuth);
827 el.get(whichRow, sdc.elevation);
828 para.get(whichRow, sdc.parangle);
829 Vector<Float> tc;
830 tcal.get(whichRow, tc);
831 sdc.tcal[0] = tc[0];sdc.tcal[1] = tc[1];
832 tcalt.get(whichRow, sdc.tcaltime);
833 Array<Float> spectrum;
834 Array<Float> tsys;
835 Array<uChar> flagtrum;
836 Vector<uInt> fmap;
837 Array<Double> direction;
838 Vector<String> histo;
839 spec.get(whichRow, spectrum);
840 sdc.putSpectrum(spectrum);
841 flags.get(whichRow, flagtrum);
842 sdc.putFlags(flagtrum);
843 ts.get(whichRow, tsys);
844 sdc.putTsys(tsys);
845 freqid.get(whichRow, fmap);
846 sdc.putFreqMap(fmap);
847 dir.get(whichRow, direction);
848 sdc.putDirection(direction);
849 hist.get(whichRow, histo);
850 sdc.putHistory(histo);
851 return sdc;
852}
853
854bool SDMemTable::putSDHeader(const SDHeader& sdh)
855{
856 table_.rwKeywordSet().define("nIF", sdh.nif);
857 table_.rwKeywordSet().define("nBeam", sdh.nbeam);
858 table_.rwKeywordSet().define("nPol", sdh.npol);
859 table_.rwKeywordSet().define("nChan", sdh.nchan);
860 table_.rwKeywordSet().define("Observer", sdh.observer);
861 table_.rwKeywordSet().define("Project", sdh.project);
862 table_.rwKeywordSet().define("Obstype", sdh.obstype);
863 table_.rwKeywordSet().define("AntennaName", sdh.antennaname);
864 table_.rwKeywordSet().define("AntennaPosition", sdh.antennaposition);
865 table_.rwKeywordSet().define("Equinox", sdh.equinox);
866 table_.rwKeywordSet().define("FreqRefFrame", sdh.freqref);
867 table_.rwKeywordSet().define("FreqRefVal", sdh.reffreq);
868 table_.rwKeywordSet().define("Bandwidth", sdh.bandwidth);
869 table_.rwKeywordSet().define("UTC", sdh.utc);
870 table_.rwKeywordSet().define("FluxUnit", sdh.fluxunit);
871 table_.rwKeywordSet().define("Epoch", sdh.epoch);
872 return true;
873}
874
875SDHeader SDMemTable::getSDHeader() const
876{
877 SDHeader sdh;
878 table_.keywordSet().get("nBeam",sdh.nbeam);
879 table_.keywordSet().get("nIF",sdh.nif);
880 table_.keywordSet().get("nPol",sdh.npol);
881 table_.keywordSet().get("nChan",sdh.nchan);
882 table_.keywordSet().get("Observer", sdh.observer);
883 table_.keywordSet().get("Project", sdh.project);
884 table_.keywordSet().get("Obstype", sdh.obstype);
885 table_.keywordSet().get("AntennaName", sdh.antennaname);
886 table_.keywordSet().get("AntennaPosition", sdh.antennaposition);
887 table_.keywordSet().get("Equinox", sdh.equinox);
888 table_.keywordSet().get("FreqRefFrame", sdh.freqref);
889 table_.keywordSet().get("FreqRefVal", sdh.reffreq);
890 table_.keywordSet().get("Bandwidth", sdh.bandwidth);
891 table_.keywordSet().get("UTC", sdh.utc);
892 table_.keywordSet().get("FluxUnit", sdh.fluxunit);
893 table_.keywordSet().get("Epoch", sdh.epoch);
894 return sdh;
895}
896void SDMemTable::makePersistent(const std::string& filename)
897{
898 table_.deepCopy(filename,Table::New);
899}
900
901Int SDMemTable::nScan() const {
902 Int n = 0;
903 ROScalarColumn<Int> scans(table_, "SCANID");
904 Int previous = -1;Int current=0;
905 for (uInt i=0; i< scans.nrow();i++) {
906 scans.getScalar(i,current);
907 if (previous != current) {
908 previous = current;
909 n++;
910 }
911 }
912 return n;
913}
914
915String SDMemTable::formatSec(Double x) const
916{
917 Double xcop = x;
918 MVTime mvt(xcop/24./3600.); // make days
919 if (x < 59.95)
920 return String(" ") + mvt.string(MVTime::TIME_CLEAN_NO_HM, 7)+"s";
921 return mvt.string(MVTime::TIME_CLEAN_NO_H, 7)+" ";
922};
923
924std::string SDMemTable::getFluxUnit() const
925{
926 String tmp;
927 table_.keywordSet().get("FluxUnit", tmp);
928 return tmp;
929}
930
931void SDMemTable::setFluxUnit(const std::string& unit)
932{
933 String tmp(unit);
934 Unit tU(tmp);
935 if (tU==Unit("K") || tU==Unit("Jy")) {
936 table_.rwKeywordSet().define(String("FluxUnit"), tmp);
937 } else {
938 throw AipsError("Illegal unit - must be compatible with Jy or K");
939 }
940}
941
942
943void SDMemTable::setInstrument(const std::string& name)
944{
945 Bool throwIt = True;
946 Instrument ins = convertInstrument (name, throwIt);
947 String nameU(name);
948 nameU.upcase();
949 table_.rwKeywordSet().define(String("AntennaName"), nameU);
950}
951
952std::string SDMemTable::summary() const {
953 ROScalarColumn<Int> scans(table_, "SCANID");
954 ROScalarColumn<String> srcs(table_, "SRCNAME");
955 ostringstream oss;
956 oss << endl;
957 oss << "--------------------------------------------------" << endl;
958 oss << " Scan Table Summary" << endl;
959 oss << "--------------------------------------------------" << endl;
960 oss.flags(std::ios_base::left);
961 oss << setw(15) << "Beams:" << setw(4) << nBeam() << endl
962 << setw(15) << "IFs:" << setw(4) << nIF() << endl
963 << setw(15) << "Polarisations:" << setw(4) << nPol() << endl
964 << setw(15) << "Channels:" << setw(4) << nChan() << endl;
965 oss << endl;
966 String tmp;
967 table_.keywordSet().get("Observer", tmp);
968 oss << setw(15) << "Observer:" << tmp << endl;
969 table_.keywordSet().get("Project", tmp);
970 oss << setw(15) << "Project:" << tmp << endl;
971 table_.keywordSet().get("Obstype", tmp);
972 oss << setw(15) << "Obs. Type:" << tmp << endl;
973 table_.keywordSet().get("AntennaName", tmp);
974 oss << setw(15) << "Antenna Name:" << tmp << endl;
975 table_.keywordSet().get("FluxUnit", tmp);
976 oss << setw(15) << "Flux Unit:" << tmp << endl;
977 Table t = table_.keywordSet().asTable("FREQUENCIES");
978 Vector<Double> vec;
979 t.keywordSet().get("RESTFREQS",vec);
980 oss << setw(15) << "Rest Freqs:";
981 if (vec.nelements() > 0) {
982 oss << setprecision(0) << vec << " [Hz]" << endl;
983 } else {
984 oss << "None set" << endl;
985 }
986 oss << setw(15) << "Abcissa:" << getAbcissaString() << endl;
987 oss << setw(15) << "Cursor:" << "Beam[" << getBeam() << "] "
988 << "IF[" << getIF() << "] " << "Pol[" << getPol() << "]" << endl;
989 oss << endl;
990 uInt count = 0;
991 String name;
992 Int previous = -1;Int current=0;
993 Int integ = 0;
994 oss << setw(6) << "Scan"
995 << setw(12) << "Source"
996 << setw(21) << "Time"
997 << setw(11) << "Integration" << endl;
998 oss << "--------------------------------------------------" << endl;
999 for (uInt i=0; i< scans.nrow();i++) {
1000 scans.getScalar(i,current);
1001 if (previous != current) {
1002 srcs.getScalar(i,name);
1003 previous = current;
1004 String t = formatSec(Double(getInterval(i)));
1005 oss << setw(6) << count << setw(12) << name << setw(21) << getTime(i)
1006 << setw(2) << setprecision(1)
1007 << t << endl;
1008 count++;
1009 } else {
1010 integ++;
1011 }
1012 }
1013 oss << endl;
1014 oss << "Table contains " << table_.nrow() << " integration(s)." << endl;
1015 oss << "in " << count << " scan(s)." << endl;
1016 oss << "--------------------------------------------------";
1017 return String(oss);
1018}
1019
1020Int SDMemTable::nBeam() const
1021{
1022 Int n;
1023 table_.keywordSet().get("nBeam",n);
1024 return n;
1025}
1026Int SDMemTable::nIF() const {
1027 Int n;
1028 table_.keywordSet().get("nIF",n);
1029 return n;
1030}
1031Int SDMemTable::nPol() const {
1032 Int n;
1033 table_.keywordSet().get("nPol",n);
1034 return n;
1035}
1036Int SDMemTable::nChan() const {
1037 Int n;
1038 table_.keywordSet().get("nChan",n);
1039 return n;
1040}
1041bool SDMemTable::appendHistory(const std::string& hist, int whichRow)
1042{
1043 ArrayColumn<String> histo(table_, "HISTORY");
1044 Vector<String> history;
1045 histo.get(whichRow, history);
1046 history.resize(history.nelements()+1,True);
1047 history[history.nelements()-1] = hist;
1048 histo.put(whichRow, history);
1049}
1050
1051std::vector<std::string> SDMemTable::history(int whichRow) const
1052{
1053 ROArrayColumn<String> hist(table_, "HISTORY");
1054 Vector<String> history;
1055 hist.get(whichRow, history);
1056 std::vector<std::string> stlout;
1057 // there is no Array<String>.tovector(std::vector<std::string>), so
1058 // do it by hand
1059 for (uInt i=0; i<history.nelements(); ++i) {
1060 stlout.push_back(history[i]);
1061 }
1062 return stlout;
1063}
1064/*
1065void SDMemTable::maskChannels(const std::vector<Int>& whichChans ) {
1066
1067 std::vector<int>::iterator it;
1068 ArrayAccessor<uChar, Axis<asap::PolAxis> > j(flags_);
1069 for (it = whichChans.begin(); it != whichChans.end(); it++) {
1070 j.reset(j.begin(uInt(*it)));
1071 for (ArrayAccessor<uChar, Axis<asap::BeamAxis> > i(j); i != i.end(); ++i) {
1072 for (ArrayAccessor<uChar, Axis<asap::IFAxis> > ii(i); ii != ii.end(); ++ii) {
1073 for (ArrayAccessor<uChar, Axis<asap::ChanAxis> > iii(ii);
1074 iii != iii.end(); ++iii) {
1075 (*iii) =
1076 }
1077 }
1078 }
1079 }
1080
1081}
1082*/
1083void SDMemTable::flag(int whichRow)
1084{
1085 ArrayColumn<uChar> spec(table_, "FLAGTRA");
1086 Array<uChar> arr;
1087 spec.get(whichRow, arr);
1088
1089 ArrayAccessor<uChar, Axis<asap::BeamAxis> > aa0(arr);
1090 aa0.reset(aa0.begin(uInt(beamSel_)));//go to beam
1091 ArrayAccessor<uChar, Axis<asap::IFAxis> > aa1(aa0);
1092 aa1.reset(aa1.begin(uInt(IFSel_)));// go to IF
1093 ArrayAccessor<uChar, Axis<asap::PolAxis> > aa2(aa1);
1094 aa2.reset(aa2.begin(uInt(polSel_)));// go to pol
1095
1096 for (ArrayAccessor<uChar, Axis<asap::ChanAxis> > i(aa2); i != i.end(); ++i) {
1097 (*i) = uChar(True);
1098 }
1099
1100 spec.put(whichRow, arr);
1101}
1102
1103MDirection::Types SDMemTable::getDirectionReference () const
1104{
1105 Float eq;
1106 table_.keywordSet().get("Equinox",eq);
1107 std::map<float,string> mp;
1108 mp[2000.0] = "J2000";
1109 mp[1950.0] = "B1950";
1110 MDirection::Types mdr;
1111 if (!MDirection::getType(mdr, mp[eq])) {
1112 mdr = MDirection::J2000;
1113 cerr << "Unknown equinox using J2000" << endl;
1114 }
1115//
1116 return mdr;
1117}
1118
1119MEpoch::Types SDMemTable::getTimeReference () const
1120{
1121 MEpoch::Types met;
1122 String ep;
1123 table_.keywordSet().get("Epoch",ep);
1124 if (!MEpoch::getType(met, ep)) {
1125 cerr << "Epoch type uknown - using UTC" << endl;
1126 met = MEpoch::UTC;
1127 }
1128//
1129 return met;
1130}
1131
1132
1133Instrument SDMemTable::convertInstrument (const String& instrument,
1134 Bool throwIt)
1135{
1136 String t(instrument);
1137 t.upcase();
1138//
1139 Instrument inst = asap::UNKNOWN;
1140 if (t==String("TID")) {
1141 inst = TIDBINBILLA;
1142 } else if (t==String("ATPKSMB")) {
1143 inst = PKSMULTIBEAM;
1144 } else if (t==String("ATPKSSB")) {
1145 inst = PKSSINGLEBEAM;
1146 } else if (t==String("MOPRA")) {
1147 inst = MOPRA;
1148 } else {
1149 if (throwIt) {
1150 throw AipsError("Unrecognized instrument - use function scan.set_instrument to set");
1151 }
1152 }
1153 return inst;
1154}
1155
Note: See TracBrowser for help on using the repository browser.