source: trunk/src/STMath.cpp@ 1004

Last change on this file since 1004 was 1000, checked in by mar637, 19 years ago

fixed bug due to incorrect use of Table::rowNumbers (RTFM!).
added flagging
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 38.6 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 < int(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 interp = stringToIMethod(method);
624
625 Vector<Float> yout;
626 Vector<Bool> maskout;
627 InterpolateArray1D<Float,Float>::interpolate(yout, maskout, xout,
628 xin, yin, maskin, interp,
629 True, True);
630
631 scaleByVector(in, Float(1.0)/yout, dotsys);
632}
633
634void STMath::scaleByVector( Table& in,
635 const Vector< Float >& factor,
636 bool dotsys )
637{
638 uInt nrow = in.nrow();
639 if ( factor.nelements() != nrow ) {
640 throw(AipsError("factors.nelements() != table.nelements()"));
641 }
642 ArrayColumn<Float> specCol(in, "SPECTRA");
643 ArrayColumn<uChar> flagCol(in, "FLAGTRA");
644 ArrayColumn<Float> tsysCol(in, "TSYS");
645 for (uInt i=0; i < nrow; ++i) {
646 MaskedArray<Float> ma = maskedArray(specCol(i), flagCol(i));
647 ma *= factor[i];
648 specCol.put(i, ma.getArray());
649 flagCol.put(i, flagsFromMA(ma));
650 if ( dotsys ) {
651 Vector<Float> tsys = tsysCol(i);
652 tsys *= factor[i];
653 tsysCol.put(i,tsys);
654 }
655 }
656}
657
658CountedPtr< Scantable > STMath::convertFlux( const CountedPtr< Scantable >& in,
659 float d, float etaap,
660 float jyperk )
661{
662 CountedPtr< Scantable > out = getScantable(in, false);
663 Table& tab = in->table();
664 Unit fluxUnit(tab.keywordSet().asString("FluxUnit"));
665 Unit K(String("K"));
666 Unit JY(String("Jy"));
667
668 bool tokelvin = true;
669 Double cfac = 1.0;
670
671 if ( fluxUnit == JY ) {
672 pushLog("Converting to K");
673 Quantum<Double> t(1.0,fluxUnit);
674 Quantum<Double> t2 = t.get(JY);
675 cfac = (t2 / t).getValue(); // value to Jy
676
677 tokelvin = true;
678 out->setFluxUnit("K");
679 } else if ( fluxUnit == K ) {
680 pushLog("Converting to Jy");
681 Quantum<Double> t(1.0,fluxUnit);
682 Quantum<Double> t2 = t.get(K);
683 cfac = (t2 / t).getValue(); // value to K
684
685 tokelvin = false;
686 out->setFluxUnit("Jy");
687 } else {
688 throw(AipsError("Unrecognized brightness units in Table - must be consistent with Jy or K"));
689 }
690 // Make sure input values are converted to either Jy or K first...
691 Float factor = cfac;
692
693 // Select method
694 if (jyperk > 0.0) {
695 factor *= jyperk;
696 if ( tokelvin ) factor = 1.0 / jyperk;
697 ostringstream oss;
698 oss << "Jy/K = " << jyperk;
699 pushLog(String(oss));
700 Vector<Float> factors(tab.nrow(), factor);
701 scaleByVector(tab,factors, false);
702 } else if ( etaap > 0.0) {
703 Instrument inst =
704 STAttr::convertInstrument(tab.keywordSet().asString("AntennaName"), True);
705 STAttr sda;
706 if (d < 0) d = sda.diameter(inst);
707 jyperk = STAttr::findJyPerK(etaap, d);
708 ostringstream oss;
709 oss << "Jy/K = " << jyperk;
710 pushLog(String(oss));
711 factor *= jyperk;
712 if ( tokelvin ) {
713 factor = 1.0 / factor;
714 }
715 Vector<Float> factors(tab.nrow(), factor);
716 scaleByVector(tab, factors, False);
717 } else {
718
719 // OK now we must deal with automatic look up of values.
720 // We must also deal with the fact that the factors need
721 // to be computed per IF and may be different and may
722 // change per integration.
723
724 pushLog("Looking up conversion factors");
725 convertBrightnessUnits(out, tokelvin, cfac);
726 }
727
728 return out;
729}
730
731void STMath::convertBrightnessUnits( CountedPtr<Scantable>& in,
732 bool tokelvin, float cfac )
733{
734 Table& table = in->table();
735 Instrument inst =
736 STAttr::convertInstrument(table.keywordSet().asString("AntennaName"), True);
737 TableIterator iter(table, "FREQ_ID");
738 STFrequencies stfreqs = in->frequencies();
739 STAttr sdAtt;
740 while (!iter.pastEnd()) {
741 Table tab = iter.table();
742 ArrayColumn<Float> specCol(tab, "SPECTRA");
743 ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
744 ROScalarColumn<uInt> freqidCol(tab, "FREQ_ID");
745 MEpoch::ROScalarColumn timeCol(tab, "TIME");
746
747 uInt freqid; freqidCol.get(0, freqid);
748 Vector<Float> tmpspec; specCol.get(0, tmpspec);
749 // STAttr.JyPerK has a Vector interface... change sometime.
750 Vector<Float> freqs(1,stfreqs.getRefFreq(freqid, tmpspec.nelements()));
751 for ( uInt i=0; i<tab.nrow(); ++i) {
752 Float jyperk = (sdAtt.JyPerK(inst, timeCol(i), freqs))[0];
753 Float factor = cfac * jyperk;
754 if ( tokelvin ) factor = Float(1.0) / factor;
755 MaskedArray<Float> ma = maskedArray(specCol(i), flagCol(i));
756 ma *= factor;
757 specCol.put(i, ma.getArray());
758 flagCol.put(i, flagsFromMA(ma));
759 }
760 ++iter;
761 }
762}
763
764CountedPtr< Scantable > STMath::opacity( const CountedPtr< Scantable > & in,
765 float tau )
766{
767 CountedPtr< Scantable > out = getScantable(in, false);
768
769 Table tab = out->table();
770 ROScalarColumn<Float> elev(tab, "ELEVATION");
771 ArrayColumn<Float> specCol(tab, "SPECTRA");
772 ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
773 for ( uInt i=0; i<tab.nrow(); ++i) {
774 Float zdist = Float(C::pi_2) - elev(i);
775 Float factor = exp(tau)/cos(zdist);
776 MaskedArray<Float> ma = maskedArray(specCol(i), flagCol(i));
777 ma *= factor;
778 specCol.put(i, ma.getArray());
779 flagCol.put(i, flagsFromMA(ma));
780 }
781 return out;
782}
783
784CountedPtr< Scantable > STMath::smooth( const CountedPtr< Scantable >& in,
785 const std::string& kernel, float width )
786{
787 CountedPtr< Scantable > out = getScantable(in, false);
788 Table& table = in->table();
789 VectorKernel::KernelTypes type = VectorKernel::toKernelType(kernel);
790 // same IFNO should have same no of channels
791 // this saves overhead
792 TableIterator iter(table, "IFNO");
793 while (!iter.pastEnd()) {
794 Table tab = iter.table();
795 ArrayColumn<Float> specCol(tab, "SPECTRA");
796 ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
797 Vector<Float> tmpspec; specCol.get(0, tmpspec);
798 uInt nchan = tmpspec.nelements();
799 Vector<Float> kvec = VectorKernel::make(type, width, nchan, True, False);
800 Convolver<Float> conv(kvec, IPosition(1,nchan));
801 Vector<Float> spec;
802 Vector<uChar> flag;
803 for ( uInt i=0; i<tab.nrow(); ++i) {
804 specCol.get(i, spec);
805 flagCol.get(i, flag);
806 Vector<Bool> mask(flag.nelements());
807 convertArray(mask, flag);
808 Vector<Float> specout;
809 if ( type == VectorKernel::HANNING ) {
810 Vector<Bool> maskout;
811 mathutil::hanning(specout, maskout, spec , mask);
812 convertArray(flag, maskout);
813 flagCol.put(i, flag);
814 } else {
815 mathutil::replaceMaskByZero(specout, mask);
816 conv.linearConv(specout, spec);
817 }
818 specCol.put(i, specout);
819 }
820 ++iter;
821 }
822 return out;
823}
824
825CountedPtr< Scantable >
826 STMath::merge( const std::vector< CountedPtr < Scantable > >& in )
827{
828 if ( in.size() < 2 ) {
829 throw(AipsError("Need at least two scantables to perform a merge."));
830 }
831 std::vector<CountedPtr < Scantable > >::const_iterator it = in.begin();
832 bool insitu = insitu_;
833 setInsitu(false);
834 CountedPtr< Scantable > out = getScantable(*it, false);
835 setInsitu(insitu);
836 Table& tout = out->table();
837 ScalarColumn<uInt> freqidcol(tout,"FREQ_ID"), molidcol(tout, "MOLECULE_ID");
838 ScalarColumn<uInt> scannocol(tout,"SCANNO"), focusidcol(tout,"FOCUS_ID");
839 // Renumber SCANNO to be 0-based
840 Vector<uInt> scannos = scannocol.getColumn();
841 uInt offset = min(scannos);
842 scannos -= offset;
843 scannocol.putColumn(scannos);
844 uInt newscanno = max(scannos)+1;
845 ++it;
846 while ( it != in.end() ){
847 if ( ! (*it)->conformant(*out) ) {
848 // log message: "ignoring scantable i, as it isn't
849 // conformant with the other(s)"
850 cerr << "oh oh" << endl;
851 ++it;
852 continue;
853 }
854 out->appendToHistoryTable((*it)->history());
855 const Table& tab = (*it)->table();
856 TableIterator scanit(tab, "SCANNO");
857 while (!scanit.pastEnd()) {
858 TableIterator freqit(scanit.table(), "FREQ_ID");
859 while ( !freqit.pastEnd() ) {
860 Table thetab = freqit.table();
861 uInt nrow = tout.nrow();
862 //tout.addRow(thetab.nrow());
863 TableCopy::copyRows(tout, thetab, nrow, 0, thetab.nrow());
864 ROTableRow row(thetab);
865 for ( uInt i=0; i<thetab.nrow(); ++i) {
866 uInt k = nrow+i;
867 scannocol.put(k, newscanno);
868 const TableRecord& rec = row.get(i);
869 Double rv,rp,inc;
870 (*it)->frequencies().getEntry(rp, rv, inc, rec.asuInt("FREQ_ID"));
871 uInt id;
872 id = out->frequencies().addEntry(rp, rv, inc);
873 freqidcol.put(k,id);
874 String name,fname;Double rf;
875 (*it)->molecules().getEntry(rf, name, fname, rec.asuInt("MOLECULE_ID"));
876 id = out->molecules().addEntry(rf, name, fname);
877 molidcol.put(k, id);
878 Float frot,fax,ftan,fhand,fmount,fuser, fxy, fxyp;
879 (*it)->focus().getEntry(fax, ftan, frot, fhand,
880 fmount,fuser, fxy, fxyp,
881 rec.asuInt("FOCUS_ID"));
882 id = out->focus().addEntry(fax, ftan, frot, fhand,
883 fmount,fuser, fxy, fxyp);
884 focusidcol.put(k, id);
885 }
886 ++freqit;
887 }
888 ++newscanno;
889 ++scanit;
890 }
891 ++it;
892 }
893 return out;
894}
895
896CountedPtr< Scantable >
897 STMath::invertPhase( const CountedPtr < Scantable >& in )
898{
899 return applyToPol(in, &STPol::invertPhase, Float(0.0));
900}
901
902CountedPtr< Scantable >
903 STMath::rotateXYPhase( const CountedPtr < Scantable >& in, float phase )
904{
905 return applyToPol(in, &STPol::rotatePhase, Float(phase));
906}
907
908CountedPtr< Scantable >
909 STMath::rotateLinPolPhase( const CountedPtr < Scantable >& in, float phase )
910{
911 return applyToPol(in, &STPol::rotateLinPolPhase, Float(phase));
912}
913
914CountedPtr< Scantable > STMath::applyToPol( const CountedPtr<Scantable>& in,
915 STPol::polOperation fptr,
916 Float phase )
917{
918 CountedPtr< Scantable > out = getScantable(in, false);
919 Table& tout = out->table();
920 Block<String> cols(4);
921 cols[0] = String("SCANNO");
922 cols[1] = String("BEAMNO");
923 cols[2] = String("IFNO");
924 cols[3] = String("CYCLENO");
925 TableIterator iter(tout, cols);
926 STPol* stpol = NULL;
927 stpol =STPol::getPolClass(out->factories_, out->getPolType() );
928 while (!iter.pastEnd()) {
929 Table t = iter.table();
930 ArrayColumn<Float> speccol(t, "SPECTRA");
931 Matrix<Float> pols = speccol.getColumn();
932 try {
933 stpol->setSpectra(pols);
934 (stpol->*fptr)(phase);
935 speccol.putColumn(stpol->getSpectra());
936 } catch (AipsError& e) {
937 delete stpol;stpol=0;
938 throw(e);
939 }
940 ++iter;
941 }
942 delete stpol;stpol=0;
943 return out;
944}
945
946CountedPtr< Scantable >
947 STMath::swapPolarisations( const CountedPtr< Scantable > & in )
948{
949 CountedPtr< Scantable > out = getScantable(in, false);
950 Table& tout = out->table();
951 Table t0 = tout(tout.col("POLNO") == 0);
952 Table t1 = tout(tout.col("POLNO") == 1);
953 if ( t0.nrow() != t1.nrow() )
954 throw(AipsError("Inconsistent number of polarisations"));
955 ArrayColumn<Float> speccol0(t0, "SPECTRA");
956 ArrayColumn<uChar> flagcol0(t0, "FLAGTRA");
957 ArrayColumn<Float> speccol1(t1, "SPECTRA");
958 ArrayColumn<uChar> flagcol1(t1, "FLAGTRA");
959 Matrix<Float> s0 = speccol0.getColumn();
960 Matrix<uChar> f0 = flagcol0.getColumn();
961 speccol0.putColumn(speccol1.getColumn());
962 flagcol0.putColumn(flagcol1.getColumn());
963 speccol1.putColumn(s0);
964 flagcol1.putColumn(f0);
965 return out;
966}
967
968CountedPtr< Scantable >
969 STMath::averagePolarisations( const CountedPtr< Scantable > & in,
970 const std::vector<bool>& mask,
971 const std::string& weight )
972{
973 if (in->getPolType() != "linear" || in->npol() != 2 )
974 throw(AipsError("averagePolarisations can only be applied to two linear polarisations."));
975 CountedPtr<Scantable> pol0( new Scantable(*in), false);
976 CountedPtr<Scantable> pol1( new Scantable(*in), false);
977 Table& tpol0 = pol0->table();
978 Table& tpol1 = pol1->table();
979 Vector<uInt> pol0rows = tpol0(tpol0.col("POLNO") == 0).rowNumbers();
980 Vector<uInt> pol1rows = tpol1(tpol1.col("POLNO") == 1).rowNumbers();
981 tpol0.removeRow(pol1rows);
982 tpol1.removeRow(pol0rows);
983 // give both tables the same POLNO
984 TableVector<uInt> vec(tpol1,"POLNO");
985 vec = 0;
986 std::vector<CountedPtr<Scantable> > pols;
987 pols.push_back(pol0);
988 pols.push_back(pol1);
989 CountedPtr< Scantable > out = average(pols, mask, weight, "NONE");
990 out->table_.rwKeywordSet().define("nPol",Int(1));
991 return out;
992}
993
994
995CountedPtr< Scantable >
996 asap::STMath::frequencyAlign( const CountedPtr< Scantable > & in,
997 const std::string & refTime,
998 const std::string & method)
999{
1000 // clone as this is not working insitu
1001 bool insitu = insitu_;
1002 setInsitu(false);
1003 CountedPtr< Scantable > out = getScantable(in, false);
1004 setInsitu(insitu);
1005 Table& tout = out->table();
1006 // Get reference Epoch to time of first row or given String
1007 Unit DAY(String("d"));
1008 MEpoch::Ref epochRef(in->getTimeReference());
1009 MEpoch refEpoch;
1010 if (refTime.length()>0) {
1011 Quantum<Double> qt;
1012 if (MVTime::read(qt,refTime)) {
1013 MVEpoch mv(qt);
1014 refEpoch = MEpoch(mv, epochRef);
1015 } else {
1016 throw(AipsError("Invalid format for Epoch string"));
1017 }
1018 } else {
1019 refEpoch = in->timeCol_(0);
1020 }
1021 MPosition refPos = in->getAntennaPosition();
1022
1023 InterpolateArray1D<Double,Float>::InterpolationMethod interp = stringToIMethod(method);
1024 // test if user frame is different to base frame
1025 if ( in->frequencies().getFrameString(true)
1026 == in->frequencies().getFrameString(false) ) {
1027 throw(AipsError("Can't convert as no output frame has been set"
1028 " (use set_freqframe) or it is aligned already."));
1029 }
1030 MFrequency::Types system = in->frequencies().getFrame();
1031 MVTime mvt(refEpoch.getValue());
1032 String epochout = mvt.string(MVTime::YMD) + String(" (") + refEpoch.getRefString() + String(")");
1033 ostringstream oss;
1034 oss << "Aligned at reference Epoch " << epochout
1035 << " in frame " << MFrequency::showType(system);
1036 pushLog(String(oss));
1037 // set up the iterator
1038 Block<String> cols(4);
1039 // select by constant direction
1040 cols[0] = String("SRCNAME");
1041 cols[1] = String("BEAMNO");
1042 // select by IF ( no of channels varies over this )
1043 cols[2] = String("IFNO");
1044 // select by restfrequency
1045 cols[3] = String("MOLECULE_ID");
1046 TableIterator iter(tout, cols);
1047 while ( !iter.pastEnd() ) {
1048 Table t = iter.table();
1049 MDirection::ROScalarColumn dirCol(t, "DIRECTION");
1050 TableIterator fiter(t, "FREQ_ID");
1051 // determine nchan from the first row. This should work as
1052 // we are iterating over BEAMNO and IFNO // we should have constant direction
1053
1054 ROArrayColumn<Float> sCol(t, "SPECTRA");
1055 MDirection direction = dirCol(0);
1056 uInt nchan = sCol(0).nelements();
1057 while ( !fiter.pastEnd() ) {
1058 Table ftab = fiter.table();
1059 ScalarColumn<uInt> freqidCol(ftab, "FREQ_ID");
1060 // get the SpectralCoordinate for the freqid, which we are iterating over
1061 SpectralCoordinate sC = in->frequencies().getSpectralCoordinate(freqidCol(0));
1062 FrequencyAligner<Float> fa( sC, nchan, refEpoch,
1063 direction, refPos, system );
1064 // realign the SpectralCoordinate and put into the output Scantable
1065 Vector<String> units(1);
1066 units = String("Hz");
1067 Bool linear=True;
1068 SpectralCoordinate sc2 = fa.alignedSpectralCoordinate(linear);
1069 sc2.setWorldAxisUnits(units);
1070 uInt id = out->frequencies().addEntry(sc2.referencePixel()[0],
1071 sc2.referenceValue()[0],
1072 sc2.increment()[0]);
1073 TableVector<uInt> tvec(ftab, "FREQ_ID");
1074 tvec = id;
1075 // create the "global" abcissa for alignment with same FREQ_ID
1076 Vector<Double> abc(nchan);
1077 Double w;
1078 for (uInt i=0; i<nchan; i++) {
1079 sC.toWorld(w,Double(i));
1080 abc[i] = w;
1081 }
1082 // cache abcissa for same time stamps, so iterate over those
1083 TableIterator timeiter(ftab, "TIME");
1084 while ( !timeiter.pastEnd() ) {
1085 Table tab = timeiter.table();
1086 ArrayColumn<Float> specCol(tab, "SPECTRA");
1087 ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
1088 MEpoch::ROScalarColumn timeCol(tab, "TIME");
1089 // use align abcissa cache after the first row
1090 bool first = true;
1091 // these rows should be just be POLNO
1092 for (int i=0; i<int(tab.nrow()); ++i) {
1093 // input values
1094 Vector<uChar> flag = flagCol(i);
1095 Vector<Bool> mask(flag.shape());
1096 Vector<Float> specOut, spec;
1097 spec = specCol(i);
1098 Vector<Bool> maskOut;Vector<uChar> flagOut;
1099 convertArray(mask, flag);
1100 // alignment
1101 Bool ok = fa.align(specOut, maskOut, abc, spec,
1102 mask, timeCol(i), !first,
1103 interp, False);
1104 // back into scantable
1105 flagOut.resize(maskOut.nelements());
1106 convertArray(flagOut, maskOut);
1107 flagCol.put(i, flagOut);
1108 specCol.put(i, specOut);
1109 // start abcissa caching
1110 first = false;
1111 }
1112 // next timestamp
1113 ++timeiter;
1114 }
1115 // next FREQ_ID
1116 ++fiter;
1117 }
1118 // next aligner
1119 ++iter;
1120 }
1121 // set this afterwards to ensure we are doing insitu correctly.
1122 out->frequencies().setFrame(system, true);
1123 return out;
1124}
1125
1126CountedPtr<Scantable>
1127 asap::STMath::convertPolarisation( const CountedPtr<Scantable>& in,
1128 const std::string & newtype )
1129{
1130 if (in->npol() != 2 && in->npol() != 4)
1131 throw(AipsError("Can only convert two or four polarisations."));
1132 if ( in->getPolType() == newtype )
1133 throw(AipsError("No need to convert."));
1134 if ( ! in->selector_.empty() )
1135 throw(AipsError("Can only convert whole scantable. Unset the selection."));
1136 bool insitu = insitu_;
1137 setInsitu(false);
1138 CountedPtr< Scantable > out = getScantable(in, true);
1139 setInsitu(insitu);
1140 Table& tout = out->table();
1141 tout.rwKeywordSet().define("POLTYPE", String(newtype));
1142
1143 Block<String> cols(4);
1144 cols[0] = "SCANNO";
1145 cols[1] = "CYCLENO";
1146 cols[2] = "BEAMNO";
1147 cols[3] = "IFNO";
1148 TableIterator it(in->originalTable_, cols);
1149 String basetype = in->getPolType();
1150 STPol* stpol = STPol::getPolClass(in->factories_, basetype);
1151 try {
1152 while ( !it.pastEnd() ) {
1153 Table tab = it.table();
1154 uInt row = tab.rowNumbers()[0];
1155 stpol->setSpectra(in->getPolMatrix(row));
1156 Float fang,fhand,parang;
1157 fang = in->focusTable_.getTotalFeedAngle(in->mfocusidCol_(row));
1158 fhand = in->focusTable_.getFeedHand(in->mfocusidCol_(row));
1159 parang = in->paraCol_(row);
1160 /// @todo re-enable this
1161 // disable total feed angle to support paralactifying Caswell style
1162 stpol->setPhaseCorrections(parang, -parang, fhand);
1163 Int npolout = 0;
1164 for (uInt i=0; i<tab.nrow(); ++i) {
1165 Vector<Float> outvec = stpol->getSpectrum(i, newtype);
1166 if ( outvec.nelements() > 0 ) {
1167 tout.addRow();
1168 TableCopy::copyRows(tout, tab, tout.nrow()-1, 0, 1);
1169 ArrayColumn<Float> sCol(tout,"SPECTRA");
1170 ScalarColumn<uInt> pCol(tout,"POLNO");
1171 sCol.put(tout.nrow()-1 ,outvec);
1172 pCol.put(tout.nrow()-1 ,uInt(npolout));
1173 npolout++;
1174 }
1175 }
1176 tout.rwKeywordSet().define("nPol", npolout);
1177 ++it;
1178 }
1179 } catch (AipsError& e) {
1180 delete stpol;
1181 throw(e);
1182 }
1183 delete stpol;
1184 return out;
1185}
Note: See TracBrowser for help on using the repository browser.