source: trunk/src/STWriter.cpp @ 2004

Last change on this file since 2004 was 2004, checked in by Takeshi Nakazato, 13 years ago

New Development: No

JIRA Issue: Yes .CAS-2718

Ready for Test: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No?

Module(s): Module Names change impacts.

Description: Describe your changes here...

The MSFiller is called instead of PKSFiller when input data is MS.
I have tested all task regressions as well as sdsave unit test and passed.

A few modification was needed for STMath::dototalpower() and
STWriter::write().
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 14.8 KB
Line 
1//#---------------------------------------------------------------------------
2//# STWriter.cc: ASAP class to write out single dish spectra.
3//#---------------------------------------------------------------------------
4//# Copyright (C) 2004
5//# ATNF
6//#
7//# This program is free software; you can redistribute it and/or modify it
8//# under the terms of the GNU General Public License as published by the Free
9//# Software Foundation; either version 2 of the License, or (at your option)
10//# any later version.
11//#
12//# This program is distributed in the hope that it will be useful, but
13//# WITHOUT ANY WARRANTY; without even the implied warranty of
14//# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General
15//# Public License for more details.
16//#
17//# You should have received a copy of the GNU General Public License along
18//# with this program; if not, write to the Free Software Foundation, Inc.,
19//# 675 Massachusetts Ave, Cambridge, MA 02139, USA.
20//#
21//# Correspondence concerning this software should be addressed as follows:
22//#        Internet email: Malte.Marquarding@csiro.au
23//#        Postal address: Malte Marquarding,
24//#                        Australia Telescope National Facility,
25//#                        P.O. Box 76,
26//#                        Epping, NSW, 2121,
27//#                        AUSTRALIA
28//#
29//# $Id: STWriter.cpp 2004 2011-02-22 09:24:33Z TakeshiNakazato $
30//#---------------------------------------------------------------------------
31
32#include <string>
33
34#include <casa/aips.h>
35#include <casa/Arrays/Array.h>
36#include <casa/Arrays/Vector.h>
37#include <casa/BasicSL/Complex.h>
38#include <casa/Utilities/CountedPtr.h>
39#include <casa/Utilities/Assert.h>
40
41#include <atnf/PKSIO/PKSrecord.h>
42#include <atnf/PKSIO/PKSMS2writer.h>
43#include <atnf/PKSIO/PKSSDwriter.h>
44#include <atnf/PKSIO/SrcType.h>
45
46#include <tables/Tables/Table.h>
47#include <tables/Tables/TableIter.h>
48#include <tables/Tables/TableRow.h>
49#include <tables/Tables/ArrayColumn.h>
50
51#include "STFITSImageWriter.h"
52#include "STAsciiWriter.h"
53#include "STHeader.h"
54#include "STMath.h"
55
56
57#include "STWriter.h"
58
59using namespace casa;
60namespace asap {
61
62STWriter::STWriter(const std::string &format)
63{
64  format_ = format;
65  String t(format_);
66  t.upcase();
67  if (t == "MS2") {
68    writer_ = new PKSMS2writer();
69  } else if (t == "SDFITS") {
70    writer_ = new PKSSDwriter();
71  } else if (t == "ASCII" || t == "FITS" || t == "CLASS") {
72    writer_ = 0;
73  } else {
74    throw (AipsError("Unrecognized export format"));
75  }
76}
77
78STWriter::~STWriter()
79{
80   if (writer_) {
81     delete writer_;
82   }
83}
84
85Int STWriter::setFormat(const std::string &format)
86{
87  if (format != format_) {
88    if (writer_) delete writer_;
89  }
90
91  format_ = format;
92  String t(format_);
93  t.upcase();
94  if (t== "MS2") {
95    writer_ = new PKSMS2writer();
96  } else if (t== "SDFITS") {
97    writer_ = new PKSSDwriter();
98  } else if (t == "ASCII" || t == "FITS" || t == "CLASS") {
99    writer_ = 0;
100  } else {
101    throw (AipsError("Unrecognized Format"));
102  }
103  return 0;
104}
105
106Int STWriter::write(const CountedPtr<Scantable> in,
107                    const std::string &filename)
108{
109  // If we write out foreign formats we have to convert the frequency system
110  // into the output frame, as we do everything related to SPectarlCoordinates
111  // in asap on-the-fly.
112
113  CountedPtr<Scantable> inst = in;
114  if (in->frequencies().getFrame(true) != in->frequencies().getFrame(false)) {
115    STMath stm(false);
116    inst = stm.frequencyAlign(in);
117  }
118
119  if (format_=="ASCII") {
120    STAsciiWriter iw;
121    // ASCII calls SpectralCoordinate::toWorld so no freqAlign use 'in'
122    if (iw.write(*in, filename)) {
123      return 0;
124    } else {
125      return 1;
126    }
127  } else if ( format_ == "FITS" || format_ == "CLASS") {
128    STFITSImageWriter iw;
129    if (format_ == "CLASS") {
130      iw.setClass(True);
131    }
132    if (iw.write(*inst, filename)) {
133      return 0;
134    }
135  }
136
137  // MS or SDFITS
138
139  // Extract the header from the table.
140  // this is a little different from what I have done
141  // before. Need to check with the Offline User Test data
142  STHeader hdr = inst->getHeader();
143  //const Int nPol  = hdr.npol;
144  //const Int nChan = hdr.nchan;
145  std::vector<uint> ifs = inst->getIFNos();
146  int nIF = inst->nif();//ifs.size();
147  Vector<uInt> nPol(nIF),nChan(nIF);
148  Vector<Bool> havexpol(nIF);
149  String fluxUnit = hdr.fluxunit;
150  fluxUnit.upcase();
151  nPol = 0;nChan = 0; havexpol = False;
152  for (uint i=0;i<ifs.size();++i) {
153    nPol(ifs[i]) = inst->npol();
154    nChan(ifs[i]) = inst->nchan(ifs[i]);
155    havexpol(ifs[i]) = nPol(ifs[i]) > 2;
156  }
157//   Vector<String> obstypes(2);
158//   obstypes(0) = "TR";//on
159//   obstypes(1) = "RF TR";//off
160  const Table table = inst->table();
161
162
163  // Create the output file and write static data.
164  Int status;
165  status = writer_->create(String(filename), hdr.observer, hdr.project,
166                           hdr.antennaname, hdr.antennaposition,
167                           hdr.obstype, hdr.fluxunit,
168                           hdr.equinox, hdr.freqref,
169                           nChan, nPol, havexpol, False);
170  if ( status ) {
171    throw(AipsError("Failed to create output file"));
172  }
173
174  Int count = 0;
175  PKSrecord pksrec;
176  pksrec.scanNo = 1;
177  // use spearate iterators to ensure renumbering of all numbers
178  TableIterator scanit(table, "SCANNO");
179  while (!scanit.pastEnd() ) {
180    Table stable = scanit.table();
181    TableIterator beamit(stable, "BEAMNO");
182    pksrec.beamNo = 1;
183    while (!beamit.pastEnd() ) {
184      Table btable = beamit.table();
185      //MDirection::ScalarColumn dirCol(btable, "DIRECTION");
186      //pksrec.direction = dirCol(0).getAngle("rad").getValue();
187      TableIterator cycit(btable, "CYCLENO");
188      ROArrayColumn<Double> srateCol(btable, "SCANRATE");
189      Vector<Double> sratedbl;
190      srateCol.get(0, sratedbl);
191      Vector<Float> srateflt(sratedbl.nelements());
192      convertArray(srateflt, sratedbl);
193      //pksrec.scanRate = srateflt;
194      pksrec.scanRate = sratedbl;
195      ROArrayColumn<Double> spmCol(btable, "SRCPROPERMOTION");
196      spmCol.get(0, pksrec.srcPM);
197      ROArrayColumn <Double> sdirCol(btable, "SRCDIRECTION");
198      sdirCol.get(0, pksrec.srcDir);
199      ROScalarColumn<Double> svelCol(btable, "SRCVELOCITY");
200      svelCol.get(0, pksrec.srcVel);
201      ROScalarColumn<uInt> bCol(btable, "BEAMNO");
202      pksrec.beamNo = bCol(0)+1;
203      pksrec.cycleNo = 1;
204      while (!cycit.pastEnd() ) {
205        Table ctable = cycit.table();
206        TableIterator ifit(ctable, "IFNO", TableIterator::Ascending, TableIterator::HeapSort);
207        MDirection::ScalarColumn dirCol(ctable, "DIRECTION");
208        pksrec.direction = dirCol(0).getAngle("rad").getValue();
209        pksrec.IFno = 1;
210        while (!ifit.pastEnd() ) {
211          Table itable = ifit.table();
212          TableRow row(itable);
213          // use the first row to fill in all the "metadata"
214          const TableRecord& rec = row.get(0);
215          ROArrayColumn<Float> specCol(itable, "SPECTRA");
216          pksrec.IFno = rec.asuInt("IFNO")+1;
217          uInt nchan = specCol(0).nelements();
218          Double crval,crpix;
219          //Vector<Double> restfreq;
220          Float tmp0,tmp1,tmp2,tmp3,tmp4;
221          String tcalt;
222          Vector<String> stmp0, stmp1;
223          inst->frequencies().getEntry(crpix,crval, pksrec.freqInc,
224                                     rec.asuInt("FREQ_ID"));
225          inst->focus().getEntry(pksrec.parAngle, pksrec.focusAxi, pksrec.focusTan,
226                                 pksrec.focusRot, tmp0,tmp1,tmp2,tmp3,tmp4,
227                                 rec.asuInt("FOCUS_ID"));
228          inst->molecules().getEntry(pksrec.restFreq,stmp0,stmp1,
229                                   rec.asuInt("MOLECULE_ID"));
230//           inst->tcal().getEntry(pksrec.tcalTime, pksrec.tcal,
231//                               rec.asuInt("TCAL_ID"));
232          inst->weather().getEntry(pksrec.temperature, pksrec.pressure,
233                                 pksrec.humidity, pksrec.windSpeed,
234                                 pksrec.windAz, rec.asuInt("WEATHER_ID"));
235          Double pixel = Double(nchan/2);
236          pksrec.refFreq = (pixel-crpix)*pksrec.freqInc + crval;
237          // ok, now we have nrows for the n polarizations in this table
238          polConversion(pksrec.spectra, pksrec.flagged, pksrec.xPol, itable);
239          pksrec.tsys = tsysFromTable(itable);
240          // dummy data
241          uInt npol = pksrec.spectra.ncolumn();
242
243          // TCAL
244          inst->tcal().getEntry( pksrec.tcalTime, pksrec.tcal,
245                                 rec.asuInt("TCAL_ID") ) ;
246          if ( pksrec.tcal.nelements() == 1 ) {
247            ROScalarColumn<uInt> uintCol( itable, "TCAL_ID" ) ;
248            Vector<uInt> tcalids = uintCol.getColumn() ;
249            pksrec.tcal.resize( npol ) ;
250            Vector<Float> dummyA ;
251            String dummyS ;
252            for ( uInt ipol = 0 ; ipol < npol ; ipol++ ) {
253              inst->tcal().getEntry( dummyS, dummyA, tcalids[ipol] ) ;
254              pksrec.tcal[ipol] = dummyA[0] ;
255            }
256          }
257          else if ( pksrec.tcal.nelements() == nchan ) {
258            ROScalarColumn<uInt> uintCol( itable, "TCAL_ID" ) ;
259            Vector<uInt> tcalids = uintCol.getColumn() ;
260            pksrec.tcal.resize( npol ) ;
261            Vector<Float> dummyA ;
262            String dummyS ;
263            for ( uInt ipol = 0 ; ipol < npol ; ipol++ ) {
264              inst->tcal().getEntry( dummyS, dummyA, tcalids[ipol] ) ;
265              pksrec.tcal[ipol] = mean( dummyA ) ;
266            }
267          }
268          //LogIO os ;
269          //os << "npol = " << npol << " pksrec.tcal = " << pksrec.tcal << LogIO::POST ;
270
271          pksrec.mjd       = rec.asDouble("TIME");
272          pksrec.interval  = rec.asDouble("INTERVAL");
273          pksrec.fieldName = rec.asString("FIELDNAME");
274          pksrec.srcName   = rec.asString("SRCNAME");
275          //pksrec.obsType   = obstypes[rec.asInt("SRCTYPE")];
276          pksrec.obsType = getObsTypes( rec.asInt("SRCTYPE") ) ;
277          pksrec.bandwidth = nchan * abs(pksrec.freqInc);
278          pksrec.azimuth   = rec.asFloat("AZIMUTH");
279          pksrec.elevation = rec.asFloat("ELEVATION");
280          pksrec.refBeam   = rec.asInt("REFBEAMNO") + 1;
281          pksrec.sigma.resize(npol);
282          pksrec.sigma     = 0.0f;
283          pksrec.calFctr.resize(npol);
284          pksrec.calFctr   = 0.0f;
285          pksrec.baseLin.resize(npol,2);
286          pksrec.baseLin   = 0.0f;
287          pksrec.baseSub.resize(npol,9);
288          pksrec.baseSub   = 0.0f;
289          pksrec.xCalFctr  = 0.0;
290          pksrec.flagrow = rec.asuInt("FLAGROW");
291
292          status = writer_->write(pksrec);
293          if ( status ) {
294            writer_->close();
295            throw(AipsError("STWriter: Failed to export Scantable."));
296          }
297          ++count;
298          //++pksrec.IFno;
299          ++ifit;
300        }
301        ++pksrec.cycleNo;
302        ++cycit;
303      }
304      //++pksrec.beamNo;
305      ++beamit;
306    }
307    ++pksrec.scanNo;
308    ++scanit;
309  }
310  ostringstream oss;
311  oss << "STWriter: wrote " << count << " rows to " << filename;
312  pushLog(String(oss));
313  writer_->close();
314  //if MS2 delete POINTING table exists and copy the one in the keyword
315  if ( format_ == "MS2" ) {
316    replacePtTab(table, filename);
317  }
318  return 0;
319}
320
321Vector<Float> STWriter::tsysFromTable(const Table& tab)
322{
323  ROArrayColumn<Float> tsysCol(tab, "TSYS");
324  Vector<Float> out(tab.nrow());
325  Vector<Float> tmp;
326  for (uInt i=0; i<tab.nrow(); ++i) {
327    tmp.resize();
328    tmp = tsysCol(i);
329    out[i] = tmp[0];
330  }
331  return out;
332}
333
334void STWriter::polConversion( Matrix< Float >& specs, Matrix< uChar >& flags,
335                              Vector< Complex > & xpol, const Table & tab )
336{
337  String poltype = tab.keywordSet().asString("POLTYPE");
338// Full stokes is not supported. Just allow stokes I
339  if ( poltype == "stokes" && tab.nrow() != 1) {
340    String msg = "poltype = " + poltype + " not yet supported in output.";
341    throw(AipsError(msg));
342  }
343  ROArrayColumn<Float> specCol(tab, "SPECTRA");
344  ROArrayColumn<uChar> flagCol(tab, "FLAGTRA");
345  uInt nchan = specCol(0).nelements();
346  uInt ncol = (tab.nrow()==1 ? 1: 2 );
347  specs.resize(nchan, ncol);
348  flags.resize(nchan, ncol);
349  // the linears
350  for (uInt i=0; i<ncol; ++i) {
351    specs.column(i) = specCol(i);
352    flags.column(i) = flagCol(i);
353  }
354  // now the complex if exists
355  Bool hasxpol = False;
356  xpol.resize();
357  if ( tab.nrow() == 4 ) {
358    hasxpol = True;
359    xpol.resize(nchan);
360    Vector<Float> reals, imags;
361    reals = specCol(2); imags = specCol(3);
362    for (uInt k=0; k < nchan; ++k) {
363      xpol[k] = Complex(reals[k], imags[k]);
364    }
365  }
366}
367
368// For writing MS data, if there is the reference to
369// original pointing table it replace it by it.
370void STWriter::replacePtTab (const Table& tab, const std::string& fname)
371{
372  String oldPtTabName = fname;
373  oldPtTabName.append("/POINTING");
374  if ( tab.keywordSet().isDefined("POINTING") ) {
375    String PointingTab = tab.keywordSet().asString("POINTING");
376    if ( Table::isReadable(PointingTab) ) {
377      Table newPtTab(PointingTab, Table::Old);
378      newPtTab.copy(oldPtTabName, Table::New);
379      ostringstream oss;
380      oss << "STWriter: copied  " <<PointingTab  << " to " << fname;
381      pushLog(String(oss));
382    }
383  }
384}
385
386// get obsType string from SRCTYPE value
387String STWriter::getObsTypes( Int srctype )
388{
389  String obsType ;
390  switch( srctype ) {
391  case Int(SrcType::PSON):
392    obsType = "PSON" ;
393    break ;
394  case Int(SrcType::PSOFF):
395    obsType = "PSOFF" ;
396    break ;
397  case Int(SrcType::NOD):
398    obsType = "NOD" ;
399    break ;
400  case Int(SrcType::FSON):
401    obsType = "FSON" ;
402    break ;
403  case Int(SrcType::FSOFF):
404    obsType = "FSOFF" ;
405    break ;
406  case Int(SrcType::SKY):
407    obsType = "SKY" ;
408    break ;
409  case Int(SrcType::HOT):
410    obsType = "HOT" ;
411    break ;
412  case Int(SrcType::WARM):
413    obsType = "WARM" ;
414    break ;
415  case Int(SrcType::COLD):
416    obsType = "COLD" ;
417    break ;
418  case Int(SrcType::PONCAL):
419    obsType = "PSON:CALON" ;
420    break ;
421  case Int(SrcType::POFFCAL):
422    obsType = "PSOFF:CALON" ;
423    break ;
424  case Int(SrcType::NODCAL):
425    obsType = "NOD:CALON" ;
426    break ;
427  case Int(SrcType::FONCAL):
428    obsType = "FSON:CALON" ;
429    break ;
430  case Int(SrcType::FOFFCAL):
431    obsType = "FSOFF:CALOFF" ;
432    break ;
433  case Int(SrcType::FSLO):
434    obsType = "FSLO" ;
435    break ;
436  case Int(SrcType::FLOOFF):
437    obsType = "FS:LOWER:OFF" ;
438    break ;
439  case Int(SrcType::FLOSKY):
440    obsType = "FS:LOWER:SKY" ;
441    break ;
442  case Int(SrcType::FLOHOT):
443    obsType = "FS:LOWER:HOT" ;
444    break ;
445  case Int(SrcType::FLOWARM):
446    obsType = "FS:LOWER:WARM" ;
447    break ;
448  case Int(SrcType::FLOCOLD):
449    obsType = "FS:LOWER:COLD" ;
450    break ;
451  case Int(SrcType::FSHI):
452    obsType = "FSHI" ;
453    break ;
454  case Int(SrcType::FHIOFF):
455    obsType = "FS:HIGHER:OFF" ;
456    break ;
457  case Int(SrcType::FHISKY):
458    obsType = "FS:HIGHER:SKY" ;
459    break ;
460  case Int(SrcType::FHIHOT):
461    obsType = "FS:HIGHER:HOT" ;
462    break ;
463  case Int(SrcType::FHIWARM):
464    obsType = "FS:HIGHER:WARM" ;
465    break ;
466  case Int(SrcType::FHICOLD):
467    obsType = "FS:HIGHER:COLD" ;
468    break ;
469  default:
470    obsType = "NOTYPE" ;
471  }
472
473  return obsType ;
474}
475
476}
Note: See TracBrowser for help on using the repository browser.