Context Navigation

source: trunk/src/SDMemTable.cc@ 163

Visit:

Last change on this file since 163 was 158, checked in by kil064, 20 years ago
Correct spelling 'abscissa' -> 'abcissa'
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 29.3 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: