source: trunk/src/STMath.cpp@ 977

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

Removed align option from average as it is buggy. The user has to provide a vector of aligned scantables.
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 36.7 KB
Line 
1//
2// C++ Implementation: STMath
3//
4// Description:
5//
6//
7// Author: Malte Marquarding <asap@atnf.csiro.au>, (C) 2006
8//
9// Copyright: See COPYING file that comes with this distribution
10//
11//
12
13#include <casa/iomanip.h>
14#include <casa/Exceptions/Error.h>
15#include <casa/Containers/Block.h>
16#include <casa/BasicSL/String.h>
17#include <tables/Tables/TableIter.h>
18#include <tables/Tables/TableCopy.h>
19#include <casa/Arrays/MaskArrLogi.h>
20#include <casa/Arrays/MaskArrMath.h>
21#include <casa/Arrays/ArrayLogical.h>
22#include <casa/Arrays/ArrayMath.h>
23#include <casa/Containers/RecordField.h>
24#include <tables/Tables/TableRow.h>
25#include <tables/Tables/TableVector.h>
26#include <tables/Tables/TabVecMath.h>
27#include <tables/Tables/ExprNode.h>
28#include <tables/Tables/TableRecord.h>
29#include <tables/Tables/ReadAsciiTable.h>
30
31#include <lattices/Lattices/LatticeUtilities.h>
32
33#include <coordinates/Coordinates/SpectralCoordinate.h>
34#include <coordinates/Coordinates/CoordinateSystem.h>
35#include <coordinates/Coordinates/CoordinateUtil.h>
36#include <coordinates/Coordinates/FrequencyAligner.h>
37
38#include <scimath/Mathematics/VectorKernel.h>
39#include <scimath/Mathematics/Convolver.h>
40#include <scimath/Functionals/Polynomial.h>
41
42#include "MathUtils.h"
43#include "RowAccumulator.h"
44#include "STAttr.h"
45#include "STMath.h"
46
47using namespace casa;
48
49using namespace asap;
50
51STMath::STMath(bool insitu) :
52 insitu_(insitu)
53{
54}
55
56
57STMath::~STMath()
58{
59}
60
61CountedPtr<Scantable>
62STMath::average( const std::vector<CountedPtr<Scantable> >& in,
63 const std::vector<bool>& mask,
64 const std::string& weight,
65 const std::string& avmode)
66{
67 if ( avmode == "SCAN" && in.size() != 1 )
68 throw(AipsError("Can't perform 'SCAN' averaging on multiple tables"));
69 WeightType wtype = stringToWeight(weight);
70
71 // output
72 // clone as this is non insitu
73 bool insitu = insitu_;
74 setInsitu(false);
75 CountedPtr< Scantable > out = getScantable(in[0], true);
76 setInsitu(insitu);
77 std::vector<CountedPtr<Scantable> >::const_iterator stit = in.begin();
78 ++stit;
79 while ( stit != in.end() ) {
80 out->appendToHistoryTable((*stit)->history());
81 ++stit;
82 }
83
84 Table& tout = out->table();
85
86 /// @todo check if all scantables are conformant
87
88 ArrayColumn<Float> specColOut(tout,"SPECTRA");
89 ArrayColumn<uChar> flagColOut(tout,"FLAGTRA");
90 ArrayColumn<Float> tsysColOut(tout,"TSYS");
91 ScalarColumn<Double> mjdColOut(tout,"TIME");
92 ScalarColumn<Double> intColOut(tout,"INTERVAL");
93
94 // set up the output table rows. These are based on the structure of the
95 // FIRST scantable in the vector
96 const Table& baset = in[0]->table();
97
98 Block<String> cols(3);
99 cols[0] = String("BEAMNO");
100 cols[1] = String("IFNO");
101 cols[2] = String("POLNO");
102 if ( avmode == "SOURCE" ) {
103 cols.resize(4);
104 cols[3] = String("SRCNAME");
105 }
106 if ( avmode == "SCAN" && in.size() == 1) {
107 cols.resize(4);
108 cols[3] = String("SCANNO");
109 }
110 uInt outrowCount = 0;
111 TableIterator iter(baset, cols);
112 while (!iter.pastEnd()) {
113 Table subt = iter.table();
114 // copy the first row of this selection into the new table
115 tout.addRow();
116 TableCopy::copyRows(tout, subt, outrowCount, 0, 1);
117 ++outrowCount;
118 ++iter;
119 }
120 RowAccumulator acc(wtype);
121 Vector<Bool> cmask(mask);
122 acc.setUserMask(cmask);
123 ROTableRow row(tout);
124 ROArrayColumn<Float> specCol, tsysCol;
125 ROArrayColumn<uChar> flagCol;
126 ROScalarColumn<Double> mjdCol, intCol;
127 ROScalarColumn<Int> scanIDCol;
128
129 for (uInt i=0; i < tout.nrow(); ++i) {
130 for ( int j=0; j < in.size(); ++j ) {
131 const Table& tin = in[j]->table();
132 const TableRecord& rec = row.get(i);
133 ROScalarColumn<Double> tmp(tin, "TIME");
134 Double td;tmp.get(0,td);
135 Table basesubt = tin(tin.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
136 && tin.col("IFNO") == Int(rec.asuInt("IFNO"))
137 && tin.col("POLNO") == Int(rec.asuInt("POLNO")) );
138 Table subt;
139 if ( avmode == "SOURCE") {
140 subt = basesubt( basesubt.col("SRCNAME") == rec.asString("SRCNAME") );
141 } else if (avmode == "SCAN") {
142 subt = basesubt( basesubt.col("SCANNO") == Int(rec.asuInt("SCANNO")) );
143 } else {
144 subt = basesubt;
145 }
146 specCol.attach(subt,"SPECTRA");
147 flagCol.attach(subt,"FLAGTRA");
148 tsysCol.attach(subt,"TSYS");
149 intCol.attach(subt,"INTERVAL");
150 mjdCol.attach(subt,"TIME");
151 Vector<Float> spec,tsys;
152 Vector<uChar> flag;
153 Double inter,time;
154 for (uInt k = 0; k < subt.nrow(); ++k ) {
155 flagCol.get(k, flag);
156 Vector<Bool> bflag(flag.shape());
157 convertArray(bflag, flag);
158 if ( allEQ(bflag, True) ) {
159 continue;//don't accumulate
160 }
161 specCol.get(k, spec);
162 tsysCol.get(k, tsys);
163 intCol.get(k, inter);
164 mjdCol.get(k, time);
165 // spectrum has to be added last to enable weighting by the other values
166 acc.add(spec, !bflag, tsys, inter, time);
167 }
168 }
169 //write out
170 specColOut.put(i, acc.getSpectrum());
171 const Vector<Bool>& msk = acc.getMask();
172 Vector<uChar> flg(msk.shape());
173 convertArray(flg, !msk);
174 flagColOut.put(i, flg);
175 tsysColOut.put(i, acc.getTsys());
176 intColOut.put(i, acc.getInterval());
177 mjdColOut.put(i, acc.getTime());
178 acc.reset();
179 }
180 return out;
181}
182
183CountedPtr< Scantable > STMath::getScantable(const CountedPtr< Scantable >& in,
184 bool droprows)
185{
186 if (insitu_) return in;
187 else {
188 // clone
189 Scantable* tabp = new Scantable(*in, Bool(droprows));
190 return CountedPtr<Scantable>(tabp);
191 }
192}
193
194CountedPtr< Scantable > STMath::unaryOperate( const CountedPtr< Scantable >& in,
195 float val,
196 const std::string& mode,
197 bool tsys )
198{
199 // modes are "ADD" and "MUL"
200 CountedPtr< Scantable > out = getScantable(in, false);
201 Table& tab = out->table();
202 ArrayColumn<Float> specCol(tab,"SPECTRA");
203 ArrayColumn<Float> tsysCol(tab,"TSYS");
204 for (uInt i=0; i<tab.nrow(); ++i) {
205 Vector<Float> spec;
206 Vector<Float> ts;
207 specCol.get(i, spec);
208 tsysCol.get(i, ts);
209 if (mode == "MUL") {
210 spec *= val;
211 specCol.put(i, spec);
212 if ( tsys ) {
213 ts *= val;
214 tsysCol.put(i, ts);
215 }
216 } else if ( mode == "ADD" ) {
217 spec += val;
218 specCol.put(i, spec);
219 if ( tsys ) {
220 ts += val;
221 tsysCol.put(i, ts);
222 }
223 }
224 }
225 return out;
226}
227
228MaskedArray<Float> STMath::maskedArray( const Vector<Float>& s,
229 const Vector<uChar>& f)
230{
231 Vector<Bool> mask;
232 mask.resize(f.shape());
233 convertArray(mask, f);
234 return MaskedArray<Float>(s,!mask);
235}
236
237Vector<uChar> STMath::flagsFromMA(const MaskedArray<Float>& ma)
238{
239 const Vector<Bool>& m = ma.getMask();
240 Vector<uChar> flags(m.shape());
241 convertArray(flags, !m);
242 return flags;
243}
244
245CountedPtr< Scantable > STMath::quotient( const CountedPtr< Scantable >& in,
246 const std::string & mode,
247 bool preserve )
248{
249 /// @todo make other modes available
250 /// modes should be "nearest", "pair"
251 // make this operation non insitu
252 const Table& tin = in->table();
253 Table ons = tin(tin.col("SRCTYPE") == Int(0));
254 Table offs = tin(tin.col("SRCTYPE") == Int(1));
255 if ( offs.nrow() == 0 )
256 throw(AipsError("No 'off' scans present."));
257 // put all "on" scans into output table
258
259 bool insitu = insitu_;
260 setInsitu(false);
261 CountedPtr< Scantable > out = getScantable(in, true);
262 setInsitu(insitu);
263 Table& tout = out->table();
264
265 TableCopy::copyRows(tout, ons);
266 TableRow row(tout);
267 ROScalarColumn<Double> offtimeCol(offs, "TIME");
268
269 ArrayColumn<Float> outspecCol(tout, "SPECTRA");
270 ROArrayColumn<Float> outtsysCol(tout, "TSYS");
271 ArrayColumn<uChar> outflagCol(tout, "FLAGTRA");
272 for (uInt i=0; i < tout.nrow(); ++i) {
273 const TableRecord& rec = row.get(i);
274 Double ontime = rec.asDouble("TIME");
275 ROScalarColumn<Double> offtimeCol(offs, "TIME");
276 Double mindeltat = min(abs(offtimeCol.getColumn() - ontime));
277 Table sel = offs( abs(offs.col("TIME")-ontime) <= mindeltat
278 && offs.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
279 && offs.col("IFNO") == Int(rec.asuInt("IFNO"))
280 && offs.col("POLNO") == Int(rec.asuInt("POLNO")) );
281
282 TableRow offrow(sel);
283 const TableRecord& offrec = offrow.get(0);//should only be one row
284 RORecordFieldPtr< Array<Float> > specoff(offrec, "SPECTRA");
285 RORecordFieldPtr< Array<Float> > tsysoff(offrec, "TSYS");
286 RORecordFieldPtr< Array<uChar> > flagoff(offrec, "FLAGTRA");
287 /// @fixme this assumes tsys is a scalar not vector
288 Float tsysoffscalar = (*tsysoff)(IPosition(1,0));
289 Vector<Float> specon, tsyson;
290 outtsysCol.get(i, tsyson);
291 outspecCol.get(i, specon);
292 Vector<uChar> flagon;
293 outflagCol.get(i, flagon);
294 MaskedArray<Float> mon = maskedArray(specon, flagon);
295 MaskedArray<Float> moff = maskedArray(*specoff, *flagoff);
296 MaskedArray<Float> quot = (tsysoffscalar * mon / moff);
297 if (preserve) {
298 quot -= tsysoffscalar;
299 } else {
300 quot -= tsyson[0];
301 }
302 outspecCol.put(i, quot.getArray());
303 outflagCol.put(i, flagsFromMA(quot));
304 }
305 // renumber scanno
306 TableIterator it(tout, "SCANNO");
307 uInt i = 0;
308 while ( !it.pastEnd() ) {
309 Table t = it.table();
310 TableVector<uInt> vec(t, "SCANNO");
311 vec = i;
312 ++i;
313 ++it;
314 }
315 return out;
316}
317
318CountedPtr< Scantable > STMath::freqSwitch( const CountedPtr< Scantable >& in )
319{
320 // make copy or reference
321 CountedPtr< Scantable > out = getScantable(in, false);
322 Table& tout = out->table();
323 Block<String> cols(3);
324 cols[0] = String("SCANNO");
325 cols[1] = String("BEAMNO");
326 cols[2] = String("POLNO");
327 TableIterator iter(tout, cols);
328 while (!iter.pastEnd()) {
329 Table subt = iter.table();
330 // this should leave us with two rows for the two IFs....if not ignore
331 if (subt.nrow() != 2 ) {
332 continue;
333 }
334 ArrayColumn<Float> specCol(tout, "SPECTRA");
335 ArrayColumn<Float> tsysCol(tout, "TSYS");
336 ArrayColumn<uChar> flagCol(tout, "FLAGTRA");
337 Vector<Float> onspec,offspec, ontsys, offtsys;
338 Vector<uChar> onflag, offflag;
339 tsysCol.get(0, ontsys); tsysCol.get(1, offtsys);
340 specCol.get(0, onspec); specCol.get(1, offspec);
341 flagCol.get(0, onflag); flagCol.get(1, offflag);
342 MaskedArray<Float> on = maskedArray(onspec, onflag);
343 MaskedArray<Float> off = maskedArray(offspec, offflag);
344 MaskedArray<Float> oncopy = on.copy();
345
346 on /= off; on -= 1.0f;
347 on *= ontsys[0];
348 off /= oncopy; off -= 1.0f;
349 off *= offtsys[0];
350 specCol.put(0, on.getArray());
351 const Vector<Bool>& m0 = on.getMask();
352 Vector<uChar> flags0(m0.shape());
353 convertArray(flags0, !m0);
354 flagCol.put(0, flags0);
355
356 specCol.put(1, off.getArray());
357 const Vector<Bool>& m1 = off.getMask();
358 Vector<uChar> flags1(m1.shape());
359 convertArray(flags1, !m1);
360 flagCol.put(1, flags1);
361 ++iter;
362 }
363
364 return out;
365}
366
367std::vector< float > STMath::statistic( const CountedPtr< Scantable > & in,
368 const std::vector< bool > & mask,
369 const std::string& which )
370{
371
372 Vector<Bool> m(mask);
373 const Table& tab = in->table();
374 ROArrayColumn<Float> specCol(tab, "SPECTRA");
375 ROArrayColumn<uChar> flagCol(tab, "FLAGTRA");
376 std::vector<float> out;
377 for (uInt i=0; i < tab.nrow(); ++i ) {
378 Vector<Float> spec; specCol.get(i, spec);
379 Vector<uChar> flag; flagCol.get(i, flag);
380 MaskedArray<Float> ma = maskedArray(spec, flag);
381 float outstat = 0.0;
382 if ( spec.nelements() == m.nelements() ) {
383 outstat = mathutil::statistics(which, ma(m));
384 } else {
385 outstat = mathutil::statistics(which, ma);
386 }
387 out.push_back(outstat);
388 }
389 return out;
390}
391
392CountedPtr< Scantable > STMath::bin( const CountedPtr< Scantable > & in,
393 int width )
394{
395 if ( !in->getSelection().empty() ) throw(AipsError("Can't bin subset of the data."));
396 CountedPtr< Scantable > out = getScantable(in, false);
397 Table& tout = out->table();
398 out->frequencies().rescale(width, "BIN");
399 ArrayColumn<Float> specCol(tout, "SPECTRA");
400 ArrayColumn<uChar> flagCol(tout, "FLAGTRA");
401 for (uInt i=0; i < tout.nrow(); ++i ) {
402 MaskedArray<Float> main = maskedArray(specCol(i), flagCol(i));
403 MaskedArray<Float> maout;
404 LatticeUtilities::bin(maout, main, 0, Int(width));
405 /// @todo implement channel based tsys binning
406 specCol.put(i, maout.getArray());
407 flagCol.put(i, flagsFromMA(maout));
408 // take only the first binned spectrum's length for the deprecated
409 // global header item nChan
410 if (i==0) tout.rwKeywordSet().define(String("nChan"),
411 Int(maout.getArray().nelements()));
412 }
413 return out;
414}
415
416CountedPtr< Scantable > STMath::resample( const CountedPtr< Scantable >& in,
417 const std::string& method,
418 float width )
419//
420// Should add the possibility of width being specified in km/s. This means
421// that for each freqID (SpectralCoordinate) we will need to convert to an
422// average channel width (say at the reference pixel). Then we would need
423// to be careful to make sure each spectrum (of different freqID)
424// is the same length.
425//
426{
427 InterpolateArray1D<Double,Float>::InterpolationMethod interp;
428 Int interpMethod(stringToIMethod(method));
429
430 CountedPtr< Scantable > out = getScantable(in, false);
431 Table& tout = out->table();
432
433// Resample SpectralCoordinates (one per freqID)
434 out->frequencies().rescale(width, "RESAMPLE");
435 TableIterator iter(tout, "IFNO");
436 TableRow row(tout);
437 while ( !iter.pastEnd() ) {
438 Table tab = iter.table();
439 ArrayColumn<Float> specCol(tab, "SPECTRA");
440 //ArrayColumn<Float> tsysCol(tout, "TSYS");
441 ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
442 Vector<Float> spec;
443 Vector<uChar> flag;
444 specCol.get(0,spec); // the number of channels should be constant per IF
445 uInt nChanIn = spec.nelements();
446 Vector<Float> xIn(nChanIn); indgen(xIn);
447 Int fac = Int(nChanIn/width);
448 Vector<Float> xOut(fac+10); // 10 to be safe - resize later
449 uInt k = 0;
450 Float x = 0.0;
451 while (x < Float(nChanIn) ) {
452 xOut(k) = x;
453 k++;
454 x += width;
455 }
456 uInt nChanOut = k;
457 xOut.resize(nChanOut, True);
458 // process all rows for this IFNO
459 Vector<Float> specOut;
460 Vector<Bool> maskOut;
461 Vector<uChar> flagOut;
462 for (uInt i=0; i < tab.nrow(); ++i) {
463 specCol.get(i, spec);
464 flagCol.get(i, flag);
465 Vector<Bool> mask(flag.nelements());
466 convertArray(mask, flag);
467
468 IPosition shapeIn(spec.shape());
469 //sh.nchan = nChanOut;
470 InterpolateArray1D<Float,Float>::interpolate(specOut, maskOut, xOut,
471 xIn, spec, mask,
472 interpMethod, True, True);
473 /// @todo do the same for channel based Tsys
474 flagOut.resize(maskOut.nelements());
475 convertArray(flagOut, maskOut);
476 specCol.put(i, specOut);
477 flagCol.put(i, flagOut);
478 }
479 ++iter;
480 }
481
482 return out;
483}
484
485STMath::imethod STMath::stringToIMethod(const std::string& in)
486{
487 static STMath::imap lookup;
488
489 // initialize the lookup table if necessary
490 if ( lookup.empty() ) {
491 lookup["nearest"] = InterpolateArray1D<Double,Float>::nearestNeighbour;
492 lookup["linear"] = InterpolateArray1D<Double,Float>::linear;
493 lookup["cubic"] = InterpolateArray1D<Double,Float>::cubic;
494 lookup["spline"] = InterpolateArray1D<Double,Float>::spline;
495 }
496
497 STMath::imap::const_iterator iter = lookup.find(in);
498
499 if ( lookup.end() == iter ) {
500 std::string message = in;
501 message += " is not a valid interpolation mode";
502 throw(AipsError(message));
503 }
504 return iter->second;
505}
506
507WeightType STMath::stringToWeight(const std::string& in)
508{
509 static std::map<std::string, WeightType> lookup;
510
511 // initialize the lookup table if necessary
512 if ( lookup.empty() ) {
513 lookup["NONE"] = asap::NONE;
514 lookup["TINT"] = asap::TINT;
515 lookup["TINTSYS"] = asap::TINTSYS;
516 lookup["TSYS"] = asap::TSYS;
517 lookup["VAR"] = asap::VAR;
518 }
519
520 std::map<std::string, WeightType>::const_iterator iter = lookup.find(in);
521
522 if ( lookup.end() == iter ) {
523 std::string message = in;
524 message += " is not a valid weighting mode";
525 throw(AipsError(message));
526 }
527 return iter->second;
528}
529
530CountedPtr< Scantable > STMath::gainElevation( const CountedPtr< Scantable >& in,
531 const vector< float > & coeff,
532 const std::string & filename,
533 const std::string& method)
534{
535 // Get elevation data from Scantable and convert to degrees
536 CountedPtr< Scantable > out = getScantable(in, false);
537 Table& tab = out->table();
538 ROScalarColumn<Float> elev(tab, "ELEVATION");
539 Vector<Float> x = elev.getColumn();
540 x *= Float(180 / C::pi); // Degrees
541
542 Vector<Float> coeffs(coeff);
543 const uInt nc = coeffs.nelements();
544 if ( filename.length() > 0 && nc > 0 ) {
545 throw(AipsError("You must choose either polynomial coefficients or an ascii file, not both"));
546 }
547
548 // Correct
549 if ( nc > 0 || filename.length() == 0 ) {
550 // Find instrument
551 Bool throwit = True;
552 Instrument inst =
553 STAttr::convertInstrument(tab.keywordSet().asString("AntennaName"),
554 throwit);
555
556 // Set polynomial
557 Polynomial<Float>* ppoly = 0;
558 Vector<Float> coeff;
559 String msg;
560 if ( nc > 0 ) {
561 ppoly = new Polynomial<Float>(nc);
562 coeff = coeffs;
563 msg = String("user");
564 } else {
565 STAttr sdAttr;
566 coeff = sdAttr.gainElevationPoly(inst);
567 ppoly = new Polynomial<Float>(3);
568 msg = String("built in");
569 }
570
571 if ( coeff.nelements() > 0 ) {
572 ppoly->setCoefficients(coeff);
573 } else {
574 delete ppoly;
575 throw(AipsError("There is no known gain-elevation polynomial known for this instrument"));
576 }
577 ostringstream oss;
578 oss << "Making polynomial correction with " << msg << " coefficients:" << endl;
579 oss << " " << coeff;
580 pushLog(String(oss));
581 const uInt nrow = tab.nrow();
582 Vector<Float> factor(nrow);
583 for ( uInt i=0; i < nrow; ++i ) {
584 factor[i] = 1.0 / (*ppoly)(x[i]);
585 }
586 delete ppoly;
587 scaleByVector(tab, factor, true);
588
589 } else {
590 // Read and correct
591 pushLog("Making correction from ascii Table");
592 scaleFromAsciiTable(tab, filename, method, x, true);
593 }
594 return out;
595}
596
597void STMath::scaleFromAsciiTable(Table& in, const std::string& filename,
598 const std::string& method,
599 const Vector<Float>& xout, bool dotsys)
600{
601
602// Read gain-elevation ascii file data into a Table.
603
604 String formatString;
605 Table tbl = readAsciiTable(formatString, Table::Memory, filename, "", "", False);
606 scaleFromTable(in, tbl, method, xout, dotsys);
607}
608
609void STMath::scaleFromTable(Table& in,
610 const Table& table,
611 const std::string& method,
612 const Vector<Float>& xout, bool dotsys)
613{
614
615 ROScalarColumn<Float> geElCol(table, "ELEVATION");
616 ROScalarColumn<Float> geFacCol(table, "FACTOR");
617 Vector<Float> xin = geElCol.getColumn();
618 Vector<Float> yin = geFacCol.getColumn();
619 Vector<Bool> maskin(xin.nelements(),True);
620
621 // Interpolate (and extrapolate) with desired method
622
623 //InterpolateArray1D<Double,Float>::InterpolationMethod method;
624 Int intmethod(stringToIMethod(method));
625
626 Vector<Float> yout;
627 Vector<Bool> maskout;
628 InterpolateArray1D<Float,Float>::interpolate(yout, maskout, xout,
629 xin, yin, maskin, intmethod,
630 True, True);
631
632 scaleByVector(in, Float(1.0)/yout, dotsys);
633}
634
635void STMath::scaleByVector( Table& in,
636 const Vector< Float >& factor,
637 bool dotsys )
638{
639 uInt nrow = in.nrow();
640 if ( factor.nelements() != nrow ) {
641 throw(AipsError("factors.nelements() != table.nelements()"));
642 }
643 ArrayColumn<Float> specCol(in, "SPECTRA");
644 ArrayColumn<uChar> flagCol(in, "FLAGTRA");
645 ArrayColumn<Float> tsysCol(in, "TSYS");
646 for (uInt i=0; i < nrow; ++i) {
647 MaskedArray<Float> ma = maskedArray(specCol(i), flagCol(i));
648 ma *= factor[i];
649 specCol.put(i, ma.getArray());
650 flagCol.put(i, flagsFromMA(ma));
651 if ( dotsys ) {
652 Vector<Float> tsys = tsysCol(i);
653 tsys *= factor[i];
654 tsysCol.put(i,tsys);
655 }
656 }
657}
658
659CountedPtr< Scantable > STMath::convertFlux( const CountedPtr< Scantable >& in,
660 float d, float etaap,
661 float jyperk )
662{
663 CountedPtr< Scantable > out = getScantable(in, false);
664 Table& tab = in->table();
665 Unit fluxUnit(tab.keywordSet().asString("FluxUnit"));
666 Unit K(String("K"));
667 Unit JY(String("Jy"));
668
669 bool tokelvin = true;
670 Double cfac = 1.0;
671
672 if ( fluxUnit == JY ) {
673 pushLog("Converting to K");
674 Quantum<Double> t(1.0,fluxUnit);
675 Quantum<Double> t2 = t.get(JY);
676 cfac = (t2 / t).getValue(); // value to Jy
677
678 tokelvin = true;
679 out->setFluxUnit("K");
680 } else if ( fluxUnit == K ) {
681 pushLog("Converting to Jy");
682 Quantum<Double> t(1.0,fluxUnit);
683 Quantum<Double> t2 = t.get(K);
684 cfac = (t2 / t).getValue(); // value to K
685
686 tokelvin = false;
687 out->setFluxUnit("Jy");
688 } else {
689 throw(AipsError("Unrecognized brightness units in Table - must be consistent with Jy or K"));
690 }
691 // Make sure input values are converted to either Jy or K first...
692 Float factor = cfac;
693
694 // Select method
695 if (jyperk > 0.0) {
696 factor *= jyperk;
697 if ( tokelvin ) factor = 1.0 / jyperk;
698 ostringstream oss;
699 oss << "Jy/K = " << jyperk;
700 pushLog(String(oss));
701 Vector<Float> factors(tab.nrow(), factor);
702 scaleByVector(tab,factors, false);
703 } else if ( etaap > 0.0) {
704 Instrument inst =
705 STAttr::convertInstrument(tab.keywordSet().asString("AntennaName"), True);
706 STAttr sda;
707 if (d < 0) d = sda.diameter(inst);
708 Float jyPerk = STAttr::findJyPerK(etaap, d);
709 ostringstream oss;
710 oss << "Jy/K = " << jyperk;
711 pushLog(String(oss));
712 factor *= jyperk;
713 if ( tokelvin ) {
714 factor = 1.0 / factor;
715 }
716 Vector<Float> factors(tab.nrow(), factor);
717 scaleByVector(tab, factors, False);
718 } else {
719
720 // OK now we must deal with automatic look up of values.
721 // We must also deal with the fact that the factors need
722 // to be computed per IF and may be different and may
723 // change per integration.
724
725 pushLog("Looking up conversion factors");
726 convertBrightnessUnits(out, tokelvin, cfac);
727 }
728
729 return out;
730}
731
732void STMath::convertBrightnessUnits( CountedPtr<Scantable>& in,
733 bool tokelvin, float cfac )
734{
735 Table& table = in->table();
736 Instrument inst =
737 STAttr::convertInstrument(table.keywordSet().asString("AntennaName"), True);
738 TableIterator iter(table, "FREQ_ID");
739 STFrequencies stfreqs = in->frequencies();
740 STAttr sdAtt;
741 while (!iter.pastEnd()) {
742 Table tab = iter.table();
743 ArrayColumn<Float> specCol(tab, "SPECTRA");
744 ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
745 ROScalarColumn<uInt> freqidCol(tab, "FREQ_ID");
746 MEpoch::ROScalarColumn timeCol(tab, "TIME");
747
748 uInt freqid; freqidCol.get(0, freqid);
749 Vector<Float> tmpspec; specCol.get(0, tmpspec);
750 // STAttr.JyPerK has a Vector interface... change sometime.
751 Vector<Float> freqs(1,stfreqs.getRefFreq(freqid, tmpspec.nelements()));
752 for ( uInt i=0; i<tab.nrow(); ++i) {
753 Float jyperk = (sdAtt.JyPerK(inst, timeCol(i), freqs))[0];
754 Float factor = cfac * jyperk;
755 if ( tokelvin ) factor = Float(1.0) / factor;
756 MaskedArray<Float> ma = maskedArray(specCol(i), flagCol(i));
757 ma *= factor;
758 specCol.put(i, ma.getArray());
759 flagCol.put(i, flagsFromMA(ma));
760 }
761 ++iter;
762 }
763}
764
765CountedPtr< Scantable > STMath::opacity( const CountedPtr< Scantable > & in,
766 float tau )
767{
768 CountedPtr< Scantable > out = getScantable(in, false);
769
770 Table tab = out->table();
771 ROScalarColumn<Float> elev(tab, "ELEVATION");
772 ArrayColumn<Float> specCol(tab, "SPECTRA");
773 ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
774 for ( uInt i=0; i<tab.nrow(); ++i) {
775 Float zdist = Float(C::pi_2) - elev(i);
776 Float factor = exp(tau)/cos(zdist);
777 MaskedArray<Float> ma = maskedArray(specCol(i), flagCol(i));
778 ma *= factor;
779 specCol.put(i, ma.getArray());
780 flagCol.put(i, flagsFromMA(ma));
781 }
782 return out;
783}
784
785CountedPtr< Scantable > STMath::smooth( const CountedPtr< Scantable >& in,
786 const std::string& kernel, float width )
787{
788 CountedPtr< Scantable > out = getScantable(in, false);
789 Table& table = in->table();
790 VectorKernel::KernelTypes type = VectorKernel::toKernelType(kernel);
791 // same IFNO should have same no of channels
792 // this saves overhead
793 TableIterator iter(table, "IFNO");
794 while (!iter.pastEnd()) {
795 Table tab = iter.table();
796 ArrayColumn<Float> specCol(tab, "SPECTRA");
797 ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
798 Vector<Float> tmpspec; specCol.get(0, tmpspec);
799 uInt nchan = tmpspec.nelements();
800 Vector<Float> kvec = VectorKernel::make(type, width, nchan, True, False);
801 Convolver<Float> conv(kvec, IPosition(1,nchan));
802 Vector<Float> spec;
803 Vector<uChar> flag;
804 for ( uInt i=0; i<tab.nrow(); ++i) {
805 specCol.get(i, spec);
806 flagCol.get(i, flag);
807 Vector<Bool> mask(flag.nelements());
808 convertArray(mask, flag);
809 Vector<Float> specout;
810 if ( type == VectorKernel::HANNING ) {
811 Vector<Bool> maskout;
812 mathutil::hanning(specout, maskout, spec , mask);
813 convertArray(flag, maskout);
814 flagCol.put(i, flag);
815 } else {
816 mathutil::replaceMaskByZero(specout, mask);
817 conv.linearConv(specout, spec);
818 }
819 specCol.put(i, specout);
820 }
821 ++iter;
822 }
823 return out;
824}
825
826CountedPtr< Scantable >
827 STMath::merge( const std::vector< CountedPtr < Scantable > >& in )
828{
829 if ( in.size() < 2 ) {
830 throw(AipsError("Need at least two scantables to perform a merge."));
831 }
832 std::vector<CountedPtr < Scantable > >::const_iterator it = in.begin();
833 bool insitu = insitu_;
834 setInsitu(false);
835 CountedPtr< Scantable > out = getScantable(*it, false);
836 setInsitu(insitu);
837 Table& tout = out->table();
838 ScalarColumn<uInt> freqidcol(tout,"FREQ_ID"), molidcol(tout, "MOLECULE_ID");
839 ScalarColumn<uInt> scannocol(tout,"SCANNO"), focusidcol(tout,"FOCUS_ID");
840 // Renumber SCANNO to be 0-based
841 Vector<uInt> scannos = scannocol.getColumn();
842 uInt offset = min(scannos);
843 scannos -= offset;
844 scannocol.putColumn(scannos);
845 uInt newscanno = max(scannos)+1;
846 ++it;
847 while ( it != in.end() ){
848 if ( ! (*it)->conformant(*out) ) {
849 // log message: "ignoring scantable i, as it isn't
850 // conformant with the other(s)"
851 cerr << "oh oh" << endl;
852 ++it;
853 continue;
854 }
855 out->appendToHistoryTable((*it)->history());
856 const Table& tab = (*it)->table();
857 TableIterator scanit(tab, "SCANNO");
858 while (!scanit.pastEnd()) {
859 TableIterator freqit(scanit.table(), "FREQ_ID");
860 while ( !freqit.pastEnd() ) {
861 Table thetab = freqit.table();
862 uInt nrow = tout.nrow();
863 //tout.addRow(thetab.nrow());
864 TableCopy::copyRows(tout, thetab, nrow, 0, thetab.nrow());
865 ROTableRow row(thetab);
866 for ( uInt i=0; i<thetab.nrow(); ++i) {
867 uInt k = nrow+i;
868 scannocol.put(k, newscanno);
869 const TableRecord& rec = row.get(i);
870 Double rv,rp,inc;
871 (*it)->frequencies().getEntry(rp, rv, inc, rec.asuInt("FREQ_ID"));
872 uInt id;
873 id = out->frequencies().addEntry(rp, rv, inc);
874 freqidcol.put(k,id);
875 String name,fname;Double rf;
876 (*it)->molecules().getEntry(rf, name, fname, rec.asuInt("MOLECULE_ID"));
877 id = out->molecules().addEntry(rf, name, fname);
878 molidcol.put(k, id);
879 Float frot,fax,ftan,fhand,fmount,fuser, fxy, fxyp;
880 (*it)->focus().getEntry(fax, ftan, frot, fhand,
881 fmount,fuser, fxy, fxyp,
882 rec.asuInt("FOCUS_ID"));
883 id = out->focus().addEntry(fax, ftan, frot, fhand,
884 fmount,fuser, fxy, fxyp);
885 focusidcol.put(k, id);
886 }
887 ++freqit;
888 }
889 ++newscanno;
890 ++scanit;
891 }
892 ++it;
893 }
894 return out;
895}
896
897CountedPtr< Scantable >
898 STMath::invertPhase( const CountedPtr < Scantable >& in )
899{
900 applyToPol(in, &STPol::invertPhase, Float(0.0));
901}
902
903CountedPtr< Scantable >
904 STMath::rotateXYPhase( const CountedPtr < Scantable >& in, float phase )
905{
906 return applyToPol(in, &STPol::rotatePhase, Float(phase));
907}
908
909CountedPtr< Scantable >
910 STMath::rotateLinPolPhase( const CountedPtr < Scantable >& in, float phase )
911{
912 return applyToPol(in, &STPol::rotateLinPolPhase, Float(phase));
913}
914
915CountedPtr< Scantable > STMath::applyToPol( const CountedPtr<Scantable>& in,
916 STPol::polOperation fptr,
917 Float phase )
918{
919 CountedPtr< Scantable > out = getScantable(in, false);
920 Table& tout = out->table();
921 Block<String> cols(4);
922 cols[0] = String("SCANNO");
923 cols[1] = String("BEAMNO");
924 cols[2] = String("IFNO");
925 cols[3] = String("CYCLENO");
926 TableIterator iter(tout, cols);
927 STPol* stpol = NULL;
928 stpol =STPol::getPolClass(out->factories_, out->getPolType() );
929 while (!iter.pastEnd()) {
930 Table t = iter.table();
931 ArrayColumn<Float> speccol(t, "SPECTRA");
932 Matrix<Float> pols = speccol.getColumn();
933 try {
934 stpol->setSpectra(pols);
935 (stpol->*fptr)(phase);
936 speccol.putColumn(stpol->getSpectra());
937 } catch (AipsError& e) {
938 delete stpol;stpol=0;
939 throw(e);
940 }
941 ++iter;
942 }
943 delete stpol;stpol=0;
944 return out;
945}
946
947CountedPtr< Scantable >
948 STMath::swapPolarisations( const CountedPtr< Scantable > & in )
949{
950 CountedPtr< Scantable > out = getScantable(in, false);
951 Table& tout = out->table();
952 Table t0 = tout(tout.col("POLNO") == 0);
953 Table t1 = tout(tout.col("POLNO") == 1);
954 if ( t0.nrow() != t1.nrow() )
955 throw(AipsError("Inconsistent number of polarisations"));
956 ArrayColumn<Float> speccol0(t0, "SPECTRA");
957 ArrayColumn<uChar> flagcol0(t0, "FLAGTRA");
958 ArrayColumn<Float> speccol1(t1, "SPECTRA");
959 ArrayColumn<uChar> flagcol1(t1, "FLAGTRA");
960 Matrix<Float> s0 = speccol0.getColumn();
961 Matrix<uChar> f0 = flagcol0.getColumn();
962 speccol0.putColumn(speccol1.getColumn());
963 flagcol0.putColumn(flagcol1.getColumn());
964 speccol1.putColumn(s0);
965 flagcol1.putColumn(f0);
966 return out;
967}
968
969CountedPtr< Scantable >
970 STMath::averagePolarisations( const CountedPtr< Scantable > & in,
971 const std::vector<bool>& mask,
972 const std::string& weight )
973{
974 if (in->getPolType() != "linear" || in->npol() != 2 )
975 throw(AipsError("averagePolarisations can only be applied to two linear polarisations."));
976 CountedPtr<Scantable> pol0( new Scantable(*in), false);
977 CountedPtr<Scantable> pol1( new Scantable(*in), false);
978 Table& tpol0 = pol0->table();
979 Table& tpol1 = pol1->table();
980 Vector<uInt> pol0rows = tpol0(tpol0.col("POLNO") == 0).rowNumbers();
981 Vector<uInt> pol1rows = tpol1(tpol1.col("POLNO") == 1).rowNumbers();
982 tpol0.removeRow(pol1rows);
983 tpol1.removeRow(pol0rows);
984 // give both tables the same POLNO
985 TableVector<uInt> vec(tpol1,"POLNO");
986 vec = 0;
987 std::vector<CountedPtr<Scantable> > pols;
988 pols.push_back(pol0);
989 pols.push_back(pol1);
990 CountedPtr< Scantable > out = average(pols, mask, weight, "NONE");
991 out->table_.rwKeywordSet().define("nPol",Int(1));
992 return out;
993}
994
995
996CountedPtr< Scantable >
997 asap::STMath::frequencyAlign( const CountedPtr< Scantable > & in,
998 const std::string & refTime,
999 const std::string & method)
1000{
1001 // clone as this is not working insitu
1002 bool insitu = insitu_;
1003 setInsitu(false);
1004 CountedPtr< Scantable > out = getScantable(in, false);
1005 setInsitu(insitu);
1006 Table& tout = out->table();
1007 // clear ouput frequency table
1008 //Table ftable = out->frequencies().table();
1009 //ftable.removeRow(ftable.rowNumbers());
1010 // Get reference Epoch to time of first row or given String
1011 Unit DAY(String("d"));
1012 MEpoch::Ref epochRef(in->getTimeReference());
1013 MEpoch refEpoch;
1014 if (refTime.length()>0) {
1015 Quantum<Double> qt;
1016 if (MVTime::read(qt,refTime)) {
1017 MVEpoch mv(qt);
1018 refEpoch = MEpoch(mv, epochRef);
1019 } else {
1020 throw(AipsError("Invalid format for Epoch string"));
1021 }
1022 } else {
1023 refEpoch = in->timeCol_(0);
1024 }
1025 MPosition refPos = in->getAntennaPosition();
1026
1027 InterpolateArray1D<Double,Float>::InterpolationMethod interp;
1028 Int interpMethod(stringToIMethod(method));
1029 // test if user frame is different to base frame
1030 if ( in->frequencies().getFrameString(true)
1031 == in->frequencies().getFrameString(false) ) {
1032 throw(AipsError("You have not set a frequency frame different from the initial - use function set_freqframe"));
1033 }
1034 MFrequency::Types system = in->frequencies().getFrame();
1035 MVTime mvt(refEpoch.getValue());
1036 String epochout = mvt.string(MVTime::YMD) + String(" (") + refEpoch.getRefString() + String(")");
1037 ostringstream oss;
1038 oss << "Aligned at reference Epoch " << epochout
1039 << " in frame " << MFrequency::showType(system);
1040 pushLog(String(oss));
1041 // set up the iterator
1042 Block<String> cols(4);
1043 // select by constant direction
1044 cols[0] = String("SRCNAME");
1045 cols[1] = String("BEAMNO");
1046 // select by IF ( no of channels varies over this )
1047 cols[2] = String("IFNO");
1048 // select by restfrequency
1049 cols[3] = String("MOLECULE_ID");
1050 TableIterator iter(tout, cols);
1051 while ( !iter.pastEnd() ) {
1052 Table t = iter.table();
1053 MDirection::ROScalarColumn dirCol(t, "DIRECTION");
1054 TableIterator fiter(t, "FREQ_ID");
1055 // determine nchan from the first row. This should work as
1056 // we are iterating over BEAMNO and IFNO // we should have constant direction
1057
1058 ROArrayColumn<Float> sCol(t, "SPECTRA");
1059 MDirection direction = dirCol(0);
1060 uInt nchan = sCol(0).nelements();
1061 while ( !fiter.pastEnd() ) {
1062 Table ftab = fiter.table();
1063 ScalarColumn<uInt> freqidCol(ftab, "FREQ_ID");
1064 // get the SpectralCoordinate for the freqid, which we are iterating over
1065 SpectralCoordinate sC = in->frequencies().getSpectralCoordinate(freqidCol(0));
1066 FrequencyAligner<Float> fa( sC, nchan, refEpoch,
1067 direction, refPos, system );
1068 // realign the SpectralCoordinate and put into the output Scantable
1069 Vector<String> units(1);
1070 units = String("Hz");
1071 Bool linear=True;
1072 SpectralCoordinate sc2 = fa.alignedSpectralCoordinate(linear);
1073 sc2.setWorldAxisUnits(units);
1074 uInt id = out->frequencies().addEntry(sc2.referencePixel()[0],
1075 sc2.referenceValue()[0],
1076 sc2.increment()[0]);
1077 TableVector<uInt> tvec(ftab, "FREQ_ID");
1078 tvec = id;
1079 // create the "global" abcissa for alignment with same FREQ_ID
1080 Vector<Double> abc(nchan);
1081 Double w;
1082 for (uInt i=0; i<nchan; i++) {
1083 sC.toWorld(w,Double(i));
1084 abc[i] = w;
1085 }
1086 // cache abcissa for same time stamps, so iterate over those
1087 TableIterator timeiter(ftab, "TIME");
1088 while ( !timeiter.pastEnd() ) {
1089 Table tab = timeiter.table();
1090 ArrayColumn<Float> specCol(tab, "SPECTRA");
1091 ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
1092 MEpoch::ROScalarColumn timeCol(tab, "TIME");
1093 // use align abcissa cache after the first row
1094 bool first = true;
1095 // these rows should be just be POLNO
1096 for (int i=0; i<tab.nrow(); ++i) {
1097 // input values
1098 Vector<uChar> flag = flagCol(i);
1099 Vector<Bool> mask(flag.shape());
1100 Vector<Float> specOut, spec;
1101 spec = specCol(i);
1102 Vector<Bool> maskOut;Vector<uChar> flagOut;
1103 convertArray(mask, flag);
1104 // alignment
1105 Bool ok = fa.align(specOut, maskOut, abc, spec,
1106 mask, timeCol(i), !first,
1107 interp, False);
1108 // back into scantable
1109 flagOut.resize(maskOut.nelements());
1110 convertArray(flagOut, maskOut);
1111 flagCol.put(i, flagOut);
1112 specCol.put(i, specOut);
1113 // start abcissa caching
1114 first = false;
1115 }
1116 // next timestamp
1117 ++timeiter;
1118 }
1119 // next FREQ_ID
1120 ++fiter;
1121 }
1122 // next aligner
1123 ++iter;
1124 }
1125 // set this afterwards to ensure we are doing insitu correctly.
1126 out->frequencies().setFrame(system, true);
1127 return out;
1128}
Note: See TracBrowser for help on using the repository browser.