source: trunk/src/STMath.cpp @ 1033

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

Fix for Ticket #33; was confusing mask and flag again. == flag svn diff ../python/scantable.py!

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