source: trunk/src/STWriter.cpp @ 2652

Last change on this file since 2652 was 2652, checked in by Malte Marquarding, 12 years ago

Fix for SDFITS part of issue #279: always limit to npol=2 if npol>2

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 13.9 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 2652 2012-09-28 01:06:32Z MalteMarquarding $
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#ifndef NOPKSMS
43#include <atnf/PKSIO/PKSMS2writer.h>
44#endif
45#include <atnf/PKSIO/PKSSDwriter.h>
46#include <atnf/PKSIO/SrcType.h>
47
48#include <tables/Tables/Table.h>
49#include <tables/Tables/TableIter.h>
50#include <tables/Tables/TableRow.h>
51#include <tables/Tables/ArrayColumn.h>
52
53#include "STFITSImageWriter.h"
54#include "STAsciiWriter.h"
55#include "STHeader.h"
56#include "STMath.h"
57
58
59#include "STWriter.h"
60
61using namespace casa;
62namespace asap {
63
64STWriter::STWriter(const std::string &format)
65{
66  format_ = format;
67  String t(format_);
68  t.upcase();
69  if (t == "MS2") {
70    #ifdef NOPKSMS
71    throw (AipsError("MS2 OUTPUT FORMAT IS NO LONGER SUPPORTED"));
72    #else
73    writer_ = new PKSMS2writer();
74    #endif
75  } else if (t == "SDFITS") {
76    writer_ = new PKSSDwriter();
77  } else if (t == "ASCII" || t == "FITS" || t == "CLASS") {
78    writer_ = 0;
79  } else {
80    throw (AipsError("Unrecognized export format"));
81  }
82}
83
84STWriter::~STWriter()
85{
86   if (writer_) {
87     delete writer_;
88   }
89}
90
91Int STWriter::setFormat(const std::string &format)
92{
93  if (format != format_) {
94    if (writer_) delete writer_;
95  }
96
97  format_ = format;
98  String t(format_);
99  t.upcase();
100  if (t== "MS2") {
101    #ifdef NOPKSMS
102    throw (AipsError("MS2 OUTPUT FORMAT IS NO LONGER SUPPORTED"));
103    #else
104    writer_ = new PKSMS2writer();
105    #endif
106  } else if (t== "SDFITS") {
107    writer_ = new PKSSDwriter();
108  } else if (t == "ASCII" || t == "FITS" || t == "CLASS") {
109    writer_ = 0;
110  } else {
111    throw (AipsError("Unrecognized Format"));
112  }
113  return 0;
114}
115
116Int STWriter::write(const CountedPtr<Scantable> in,
117                    const std::string &filename)
118{
119  // If we write out foreign formats we have to convert the frequency system
120  // into the output frame, as we do everything related to SPectarlCoordinates
121  // in asap on-the-fly.
122
123  CountedPtr<Scantable> inst = in;
124  if (in->frequencies().getFrame(true) != in->frequencies().getFrame(false)) {
125    STMath stm(false);
126    inst = stm.frequencyAlign(in);
127  }
128
129  if (format_=="ASCII") {
130    STAsciiWriter iw;
131    // ASCII calls SpectralCoordinate::toWorld so no freqAlign use 'in'
132    if (iw.write(*in, filename)) {
133      return 0;
134    } else {
135      return 1;
136    }
137  } else if ( format_ == "FITS" || format_ == "CLASS") {
138    STFITSImageWriter iw;
139    if (format_ == "CLASS") {
140      iw.setClass(True);
141    }
142    if (iw.write(*inst, filename)) {
143      return 0;
144    }
145  }
146
147  // MS or SDFITS
148
149  // Extract the header from the table.
150  // this is a little different from what I have done
151  // before. Need to check with the Offline User Test data
152  STHeader hdr = inst->getHeader();
153  //const Int nPol  = hdr.npol;
154  //const Int nChan = hdr.nchan;
155  std::vector<uint> ifs = inst->getIFNos();
156  int nIF = inst->nif();//ifs.size();
157  Vector<uInt> nPol(nIF),nChan(nIF);
158  Vector<Bool> havexpol(nIF);
159  String fluxUnit = hdr.fluxunit;
160  fluxUnit.upcase();
161  nPol = 0;nChan = 0; havexpol = False;
162  for (uint i=0;i<ifs.size();++i) {
163    if ( inst->npol() > 2 ) {
164      havexpol(ifs[i]) = True;     
165      nPol(ifs[i]) = 2;
166    } else {
167      nPol(ifs[i]) = inst->npol();
168    }
169    nChan(ifs[i]) = inst->nchan(ifs[i]);
170  }
171//   Vector<String> obstypes(2);
172//   obstypes(0) = "TR";//on
173//   obstypes(1) = "RF TR";//off
174  const Table table = inst->table();
175
176
177  // Create the output file and write static data.
178  Int status;
179//   status = writer_->create(String(filename), hdr.observer, hdr.project,
180//                            hdr.antennaname, hdr.antennaposition,
181//                            hdr.obstype, hdr.fluxunit,
182//                            hdr.equinox, hdr.freqref,
183//                            nChan, nPol, havexpol, False);
184  status = writer_->create(String(filename), hdr.observer, hdr.project,
185                           inst->getAntennaName(), hdr.antennaposition,
186                           hdr.obstype, hdr.fluxunit,
187                           hdr.equinox, hdr.freqref,
188                           nChan, nPol, havexpol, False);
189  if ( status ) {
190    throw(AipsError("Failed to create output file"));
191  }
192
193  Int count = 0;
194  PKSrecord pksrec;
195  pksrec.scanNo = 1;
196  // use spearate iterators to ensure renumbering of all numbers
197  TableIterator scanit(table, "SCANNO");
198  while (!scanit.pastEnd() ) {
199    Table stable = scanit.table();
200    TableIterator beamit(stable, "BEAMNO");
201    pksrec.beamNo = 1;
202    while (!beamit.pastEnd() ) {
203      Table btable = beamit.table();
204      //MDirection::ScalarColumn dirCol(btable, "DIRECTION");
205      //pksrec.direction = dirCol(0).getAngle("rad").getValue();
206      TableIterator cycit(btable, "CYCLENO");
207      ROArrayColumn<Double> srateCol(btable, "SCANRATE");
208      Vector<Double> sratedbl;
209      srateCol.get(0, sratedbl);
210      Vector<Float> srateflt(sratedbl.nelements());
211      convertArray(srateflt, sratedbl);
212      //pksrec.scanRate = srateflt;
213      pksrec.scanRate = sratedbl;
214      ROArrayColumn<Double> spmCol(btable, "SRCPROPERMOTION");
215      spmCol.get(0, pksrec.srcPM);
216      ROArrayColumn <Double> sdirCol(btable, "SRCDIRECTION");
217      sdirCol.get(0, pksrec.srcDir);
218      ROScalarColumn<Double> svelCol(btable, "SRCVELOCITY");
219      svelCol.get(0, pksrec.srcVel);
220      ROScalarColumn<uInt> bCol(btable, "BEAMNO");
221      pksrec.beamNo = bCol(0)+1;
222      pksrec.cycleNo = 1;
223      while (!cycit.pastEnd() ) {
224        Table ctable = cycit.table();
225        TableIterator typeit( ctable, "SRCTYPE" ) ;
226        while(!typeit.pastEnd() ) {
227          Table ttable = typeit.table() ;
228          TableIterator ifit(ttable, "IFNO",
229                             TableIterator::Ascending, TableIterator::HeapSort);
230          MDirection::ScalarColumn dirCol(ctable, "DIRECTION");
231          pksrec.direction = dirCol(0).getAngle("rad").getValue();
232          pksrec.IFno = 1;
233          while (!ifit.pastEnd() ) {
234            Table itable = ifit.table();
235            TableRow row(itable);
236            // use the first row to fill in all the "metadata"
237            const TableRecord& rec = row.get(0);
238            ROArrayColumn<Float> specCol(itable, "SPECTRA");
239            pksrec.IFno = rec.asuInt("IFNO")+1;
240            uInt nchan = specCol(0).nelements();
241            Double crval,crpix;
242            //Vector<Double> restfreq;
243            Float tmp0,tmp1,tmp2,tmp3,tmp4;
244            String tcalt;
245            Vector<String> stmp0, stmp1;
246            inst->frequencies().getEntry(crpix,crval, pksrec.freqInc,
247                                         rec.asuInt("FREQ_ID"));
248            inst->focus().getEntry(pksrec.parAngle, pksrec.focusAxi,
249                                   pksrec.focusTan,
250                                   pksrec.focusRot, tmp0,tmp1,tmp2,tmp3,tmp4,
251                                   rec.asuInt("FOCUS_ID"));
252            inst->molecules().getEntry(pksrec.restFreq,stmp0,stmp1,
253                                       rec.asuInt("MOLECULE_ID"));
254//             inst->tcal().getEntry(pksrec.tcalTime, pksrec.tcal,
255//                                   rec.asuInt("TCAL_ID"));
256            inst->weather().getEntry(pksrec.temperature, pksrec.pressure,
257                                     pksrec.humidity, pksrec.windSpeed,
258                                     pksrec.windAz, rec.asuInt("WEATHER_ID"));
259            Double pixel = Double(nchan/2);
260            pksrec.refFreq = (pixel-crpix)*pksrec.freqInc + crval;
261            // ok, now we have nrows for the n polarizations in this table
262            polConversion(pksrec.spectra, pksrec.flagged, pksrec.xPol, itable);
263            pksrec.tsys = tsysFromTable(itable);
264            // dummy data
265            uInt npol = pksrec.spectra.ncolumn();
266           
267            // TCAL
268            inst->tcal().getEntry( pksrec.tcalTime, pksrec.tcal,
269                                   rec.asuInt("TCAL_ID") ) ;
270            if ( pksrec.tcal.nelements() == 1 ) {
271              ROScalarColumn<uInt> uintCol( itable, "TCAL_ID" ) ;
272              Vector<uInt> tcalids = uintCol.getColumn() ;
273              pksrec.tcal.resize( npol ) ;
274              Vector<Float> dummyA ;
275              String dummyS ;
276              for ( uInt ipol = 0 ; ipol < npol ; ipol++ ) {
277                inst->tcal().getEntry( dummyS, dummyA, tcalids[ipol] ) ;
278                pksrec.tcal[ipol] = dummyA[0] ;
279              }
280            }
281            else if ( pksrec.tcal.nelements() == nchan ) {
282              ROScalarColumn<uInt> uintCol( itable, "TCAL_ID" ) ;
283              Vector<uInt> tcalids = uintCol.getColumn() ;
284              pksrec.tcal.resize( npol ) ;
285              Vector<Float> dummyA ;
286              String dummyS ;
287              for ( uInt ipol = 0 ; ipol < npol ; ipol++ ) {
288                inst->tcal().getEntry( dummyS, dummyA, tcalids[ipol] ) ;
289                pksrec.tcal[ipol] = mean( dummyA ) ;
290              }
291            }
292            pksrec.mjd       = rec.asDouble("TIME");
293            pksrec.interval  = rec.asDouble("INTERVAL");
294            pksrec.fieldName = rec.asString("FIELDNAME");
295            pksrec.srcName   = rec.asString("SRCNAME");
296            //pksrec.obsType   = obstypes[rec.asInt("SRCTYPE")];
297            pksrec.obsType = getObsTypes( rec.asInt("SRCTYPE") ) ;
298            pksrec.bandwidth = nchan * abs(pksrec.freqInc);
299            pksrec.azimuth   = rec.asFloat("AZIMUTH");
300            pksrec.elevation = rec.asFloat("ELEVATION");
301            pksrec.refBeam   = rec.asInt("REFBEAMNO") + 1;
302            pksrec.sigma.resize(npol);
303            pksrec.sigma     = 0.0f;
304            pksrec.calFctr.resize(npol);
305            pksrec.calFctr   = 0.0f;
306            pksrec.baseLin.resize(npol,2);
307            pksrec.baseLin   = 0.0f;
308            pksrec.baseSub.resize(npol,9);
309            pksrec.baseSub   = 0.0f;
310            pksrec.xCalFctr  = 0.0;
311            pksrec.flagrow = rec.asuInt("FLAGROW");
312           
313            status = writer_->write(pksrec);
314            if ( status ) {
315              writer_->close();
316              throw(AipsError("STWriter: Failed to export Scantable."));
317            }
318            ++count;
319            //++pksrec.IFno;
320            ++ifit;
321          }
322          ++typeit ;
323        }
324        ++pksrec.cycleNo;
325        ++cycit;
326      }
327      //++pksrec.beamNo;
328      ++beamit;
329    }
330    ++pksrec.scanNo;
331    ++scanit;
332  }
333  ostringstream oss;
334  oss << "STWriter: wrote " << count << " rows to " << filename;
335  pushLog(String(oss));
336 
337  writer_->close();
338  //if MS2 delete POINTING table exists and copy the one in the keyword
339  if ( format_ == "MS2" ) {
340    replacePtTab(table, filename);
341  }
342  return 0;
343}
344
345Vector<Float> STWriter::tsysFromTable(const Table& tab)
346{
347  ROArrayColumn<Float> tsysCol(tab, "TSYS");
348  Vector<Float> out(tab.nrow());
349  Vector<Float> tmp;
350  for (uInt i=0; i<tab.nrow(); ++i) {
351    tmp.resize();
352    tmp = tsysCol(i);
353    out[i] = tmp[0];
354  }
355  return out;
356}
357
358void STWriter::polConversion( Matrix< Float >& specs, Matrix< uChar >& flags,
359                              Vector< Complex > & xpol, const Table & tab )
360{
361  String poltype = tab.keywordSet().asString("POLTYPE");
362// Full stokes is not supported. Just allow stokes I
363  if ( poltype == "stokes" && tab.nrow() != 1) {
364    String msg = "poltype = " + poltype + " not yet supported in output.";
365    throw(AipsError(msg));
366  }
367  ROArrayColumn<Float> specCol(tab, "SPECTRA");
368  ROArrayColumn<uChar> flagCol(tab, "FLAGTRA");
369  uInt nchan = specCol(0).nelements();
370  uInt ncol = (tab.nrow()==1 ? 1: 2 );
371  specs.resize(nchan, ncol);
372  flags.resize(nchan, ncol);
373  // the linears
374  for (uInt i=0; i<ncol; ++i) {
375    specs.column(i) = specCol(i);
376    flags.column(i) = flagCol(i);
377  }
378  // now the complex if exists
379  Bool hasxpol = False;
380  xpol.resize();
381  if ( tab.nrow() == 4 ) {
382    hasxpol = True;
383    xpol.resize(nchan);
384    Vector<Float> reals, imags;
385    reals = specCol(2); imags = specCol(3);
386    for (uInt k=0; k < nchan; ++k) {
387      xpol[k] = Complex(reals[k], imags[k]);
388    }
389  }
390}
391
392// For writing MS data, if there is the reference to
393// original pointing table it replace it by it.
394void STWriter::replacePtTab (const Table& tab, const std::string& fname)
395{
396  String oldPtTabName = fname;
397  oldPtTabName.append("/POINTING");
398  if ( tab.keywordSet().isDefined("POINTING") ) {
399    String PointingTab = tab.keywordSet().asString("POINTING");
400    if ( Table::isReadable(PointingTab) ) {
401      Table newPtTab(PointingTab, Table::Old);
402      newPtTab.copy(oldPtTabName, Table::New);
403      ostringstream oss;
404      oss << "STWriter: copied  " <<PointingTab  << " to " << fname;
405      pushLog(String(oss));
406    }
407  }
408}
409
410// get obsType string from SRCTYPE value
411String STWriter::getObsTypes( Int srctype )
412{
413  return SrcType::getName(srctype) ;
414}
415
416}
Note: See TracBrowser for help on using the repository browser.