source: trunk/src/STWriter.cpp@ 1446

Last change on this file since 1446 was 1443, checked in by Malte Marquarding, 16 years ago

Fix ticket #127; still have to add class header hack

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 11.5 KB
RevLine 
[470]1//#---------------------------------------------------------------------------
[822]2//# STWriter.cc: ASAP class to write out single dish spectra.
[28]3//#---------------------------------------------------------------------------
4//# Copyright (C) 2004
[125]5//# ATNF
[28]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 1443 2008-09-01 01:54:04Z MalteMarquarding $
30//#---------------------------------------------------------------------------
31
32#include <string>
33
[81]34#include <casa/aips.h>
[470]35#include <casa/Arrays/Array.h>
36#include <casa/Arrays/Vector.h>
[81]37#include <casa/BasicSL/Complex.h>
38#include <casa/Utilities/CountedPtr.h>
[470]39#include <casa/Utilities/Assert.h>
[28]40
41#include <atnf/PKSIO/PKSMS2writer.h>
42#include <atnf/PKSIO/PKSSDwriter.h>
43
[827]44#include <tables/Tables/Table.h>
45#include <tables/Tables/TableIter.h>
46#include <tables/Tables/TableRow.h>
[324]47#include <tables/Tables/ArrayColumn.h>
48
[1443]49#include "STFITSImageWriter.h"
[988]50#include "STAsciiWriter.h"
[901]51#include "STHeader.h"
[324]52
[822]53#include "STWriter.h"
[28]54
[125]55using namespace casa;
[827]56namespace asap {
[28]57
[822]58STWriter::STWriter(const std::string &format)
[28]59{
[988]60 format_ = format;
61 String t(format_);
[450]62 t.upcase();
[1443]63 if (t == "MS2") {
[988]64 writer_ = new PKSMS2writer();
[1443]65 } else if (t == "SDFITS") {
[988]66 writer_ = new PKSSDwriter();
[1443]67 } else if (t == "ASCII" || t == "FITS" || t == "CLASS") {
[988]68 writer_ = 0;
[450]69 } else {
[988]70 throw (AipsError("Unrecognized export format"));
[28]71 }
72}
73
[822]74STWriter::~STWriter()
[28]75{
[988]76 if (writer_) {
77 delete writer_;
[192]78 }
[28]79}
80
[822]81Int STWriter::setFormat(const std::string &format)
[28]82{
[988]83 if (format != format_) {
84 if (writer_) delete writer_;
[450]85 }
[827]86
[988]87 format_ = format;
88 String t(format_);
[450]89 t.upcase();
90 if (t== "MS2") {
[988]91 writer_ = new PKSMS2writer();
[450]92 } else if (t== "SDFITS") {
[988]93 writer_ = new PKSSDwriter();
[1443]94 } else if (t == "ASCII" || t == "FITS" || t == "CLASS") {
[988]95 writer_ = 0;
[450]96 } else {
97 throw (AipsError("Unrecognized Format"));
[28]98 }
99 return 0;
100}
101
[827]102Int STWriter::write(const CountedPtr<Scantable> in,
103 const std::string &filename)
[28]104{
[192]105
[1280]106 if (format_=="ASCII") {
107 STAsciiWriter iw;
108 if (iw.write(*in, filename)) {
109 return 0;
110 } else {
111 return 1;
112 }
[1443]113 } else if ( format_ == "FITS" || format_ == "CLASS") {
114 STFITSImageWriter iw;
115 if (format_ == "CLASS") {
116 iw.setClass(True);
117 }
118 if (iw.write(*in, filename)) {
119 return 0;
120 }
[199]121 }
[192]122
[827]123 // MS or SDFITS
[192]124
[28]125 // Extract the header from the table.
[1391]126 // this is a little different from what I have done
127 // before. Need to check with the Offline User Test data
[901]128 STHeader hdr = in->getHeader();
[1060]129 //const Int nPol = hdr.npol;
130 //const Int nChan = hdr.nchan;
[1295]131 std::vector<uint> ifs = in->getIFNos();
132 int nIF = in->nif();//ifs.size();
[1060]133 Vector<uInt> nPol(nIF),nChan(nIF);
134 Vector<Bool> havexpol(nIF);
[1391]135 String fluxUnit = hdr.fluxunit;
136
[1295]137 nPol = 0;nChan = 0; havexpol = False;
[1305]138 for (uint i=0;i<ifs.size();++i) {
[1295]139 nPol(ifs[i]) = in->npol();
140 nChan(ifs[i]) = in->nchan(ifs[i]);
141 havexpol(ifs[i]) = nPol(ifs[i]) > 2;
[1060]142 }
[28]143
[324]144 const Table table = in->table();
145
[28]146 // Create the output file and write static data.
147 Int status;
[1060]148 status = writer_->create(String(filename), hdr.observer, hdr.project,
[1410]149 hdr.antennaname, hdr.antennaposition,
150 hdr.obstype, hdr.fluxunit,
151 hdr.equinox, hdr.freqref,
152 nChan, nPol, havexpol, False);
[996]153 if ( status ) {
[717]154 throw(AipsError("Failed to create output file"));
[28]155 }
156
[999]157 Double srcVel;
[822]158
159 String fieldName, srcName, tcalTime;
160 Vector<Float> calFctr, sigma, tcal, tsys;
[999]161 Vector<Double> direction(2), scanRate(2), srcDir(2), srcPM(2);
[822]162 Matrix<Float> spectra;
163 Matrix<uChar> flagtra;
164 Complex xCalFctr;
[28]165 Int count = 0;
[822]166 Int scanno = 1;
167 // use spearate iterators to ensure renumbering of all numbers
168 TableIterator scanit(table, "SCANNO");
169 while (!scanit.pastEnd() ) {
170 Table stable = scanit.table();
[988]171 TableIterator beamit(stable, "BEAMNO");
[822]172 Int beamno = 1;
173 while (!beamit.pastEnd() ) {
174 Table btable = beamit.table();
[988]175 TableIterator cycit(btable, "CYCLENO");
[999]176 ROArrayColumn<Double> srateCol(btable, "SCANRATE");
177 srateCol.get(0, scanRate);
178 ROArrayColumn<Double> spmCol(btable, "SRCPROPERMOTION");
179 spmCol.get(0, srcPM);
180 ROArrayColumn <Double> sdirCol(btable, "SRCDIRECTION");
181 sdirCol.get(0, srcDir);
182 ROScalarColumn<Double> svelCol(btable, "SRCVELOCITY");
183 svelCol.get(0, srcVel);
[1295]184 ROScalarColumn<uInt> bCol(btable, "BEAMNO");
185 beamno = bCol(0)+1;
[988]186 Int cycno = 1;
187 while (!cycit.pastEnd() ) {
188 Table ctable = cycit.table();
189 TableIterator ifit(ctable, "IFNO");
[1390]190 MDirection::ScalarColumn dirCol(ctable, "DIRECTION");
191 Vector<Double> direction = dirCol(0).getAngle("rad").getValue();
[988]192 Int ifno = 1;
193 while (!ifit.pastEnd() ) {
194 Table itable = ifit.table();
195 TableRow row(itable);
[822]196 // use the first row to fill in all the "metadata"
[827]197 const TableRecord& rec = row.get(0);
[988]198 ROArrayColumn<Float> specCol(itable, "SPECTRA");
[1295]199 ifno = rec.asuInt("IFNO")+1;
[822]200 uInt nchan = specCol(0).nelements();
201 Double cdelt,crval,crpix, restfreq;
[827]202 Float focusAxi, focusTan, focusRot,
203 temperature, pressure, humidity, windSpeed, windAz;
[988]204 Float tmp0,tmp1,tmp2,tmp3,tmp4;
[827]205 Vector<Float> tcalval;
[988]206 String stmp0,stmp1, tcalt;
[822]207 in->frequencies().getEntry(crpix,crval,cdelt, rec.asuInt("FREQ_ID"));
[988]208 in->focus().getEntry(focusAxi, focusTan, focusRot,
209 tmp0,tmp1,tmp2,tmp3,tmp4,
210 rec.asuInt("FOCUS_ID"));
211 in->molecules().getEntry(restfreq,stmp0,stmp1,rec.asuInt("MOLECULE_ID"));
[827]212 in->tcal().getEntry(tcalt,tcalval,rec.asuInt("TCAL_ID"));
213 in->weather().getEntry(temperature, pressure, humidity,
214 windSpeed, windAz,
215 rec.asuInt("WEATHER_ID"));
[822]216 Double pixel = Double(nchan/2);
217 Double refFreqNew = (pixel-crpix)*cdelt + crval;
218 // ok, now we have nrows for the n polarizations in this table
219 Matrix<Float> specs;
220 Matrix<uChar> flags;
221 Vector<Complex> xpol;
[988]222 polConversion(specs, flags, xpol, itable);
223 Vector<Float> tsys = tsysFromTable(itable);
[822]224 // dummy data
[827]225 uInt npol = specs.ncolumn();
[988]226
[822]227 Matrix<Float> baseLin(npol,2, 0.0f);
228 Matrix<Float> baseSub(npol,9, 0.0f);
229 Complex xCalFctr;
230 Vector<Double> scanRate(2, 0.0);
231 Vector<Float> sigma(npol, 0.0f);
232 Vector<Float> calFctr(npol, 0.0f);
[996]233 status = writer_->write(scanno, cycno, rec.asDouble("TIME"),
[822]234 rec.asDouble("INTERVAL"),
235 rec.asString("FIELDNAME"),
236 rec.asString("SRCNAME"),
[1072]237 srcDir, srcPM, srcVel,hdr.obstype,
[822]238 ifno,
239 refFreqNew, nchan*abs(cdelt), cdelt,
240 restfreq,
[827]241 tcal,
242 tcalt,
[822]243 rec.asFloat("AZIMUTH"),
244 rec.asFloat("ELEVATION"),
245 rec.asFloat("PARANGLE"),
[988]246 focusAxi, focusTan, focusRot,
[827]247 temperature,
248 pressure, humidity, windSpeed, windAz,
[988]249 rec.asInt("REFBEAMNO")+1, beamno,
[827]250 direction,
[999]251 scanRate,
[827]252 tsys,
[822]253 sigma, calFctr,// not in scantable
254 baseLin, baseSub,// not in scantable
255 specs, flags,
256 xCalFctr,//
[996]257 xpol);
258 if ( status ) {
[988]259 writer_->close();
260 throw(AipsError("STWriter: Failed to export Scantable."));
[822]261 }
[999]262 ++count;
[1295]263 //++ifno;
[988]264 ++ifit;
[470]265 }
[988]266 ++cycno;
267 ++cycit;
[28]268 }
[1295]269 //++beamno;
[822]270 ++beamit;
[28]271 }
[822]272 ++scanno;
273 ++scanit;
[28]274 }
[717]275 ostringstream oss;
[999]276 oss << "STWriter: wrote " << count << " rows to " << filename;
[717]277 pushLog(String(oss));
[988]278 writer_->close();
[1391]279 //if MS2 delete POINTING table exists and copy the one in the keyword
280 if ( format_ == "MS2" ) {
281 replacePtTab(table, filename);
282 }
[28]283 return 0;
284}
[470]285
[827]286Vector<Float> STWriter::tsysFromTable(const Table& tab)
287{
288 ROArrayColumn<Float> tsysCol(tab, "TSYS");
289 Vector<Float> out(tab.nrow());
290 Vector<Float> tmp;
291 for (uInt i=0; i<tab.nrow(); ++i) {
[988]292 tmp.resize();
[827]293 tmp = tsysCol(i);
294 out[i] = tmp[0];
295 }
[988]296 return out;
[827]297}
298
299void STWriter::polConversion( Matrix< Float >& specs, Matrix< uChar >& flags,
[822]300 Vector< Complex > & xpol, const Table & tab )
301{
302 String poltype = tab.keywordSet().asString("POLTYPE");
[1259]303 if ( poltype == "stokes") {
[822]304 String msg = "poltype = " + poltype + " not yet supported in output.";
[988]305 throw(AipsError(msg));
306 }
[827]307 ROArrayColumn<Float> specCol(tab, "SPECTRA");
308 ROArrayColumn<uChar> flagCol(tab, "FLAGTRA");
[822]309 uInt nchan = specCol(0).nelements();
[988]310 uInt ncol = (tab.nrow()==1 ? 1: 2 );
311 specs.resize(nchan, ncol);
312 flags.resize(nchan, ncol);
[822]313 // the linears
[988]314 for (uInt i=0; i<ncol; ++i) {
[822]315 specs.column(i) = specCol(i);
316 flags.column(i) = flagCol(i);
317 }
318 // now the complex if exists
319 Bool hasxpol = False;
320 xpol.resize();
[827]321 if ( tab.nrow() == 4 ) {
[822]322 hasxpol = True;
323 xpol.resize(nchan);
324 Vector<Float> reals, imags;
325 reals = specCol(2); imags = specCol(3);
326 for (uInt k=0; k < nchan; ++k) {
[827]327 xpol[k] = Complex(reals[k], imags[k]);
[822]328 }
329 }
330}
[827]331
[1391]332// For writing MS data, if there is the reference to
333// original pointing table it replace it by it.
334void STWriter::replacePtTab (const Table& tab, const std::string& fname)
335{
336 String oldPtTabName = fname;
337 oldPtTabName.append("/POINTING");
338 if ( tab.keywordSet().isDefined("POINTING") ) {
339 String PointingTab = tab.keywordSet().asString("POINTING");
340 if ( Table::isReadable(PointingTab) ) {
341 Table newPtTab(PointingTab, Table::Old);
342 newPtTab.copy(oldPtTabName, Table::New);
343 ostringstream oss;
344 oss << "STWriter: copied " <<PointingTab << " to " << fname;
345 pushLog(String(oss));
346 }
347 }
348}
[827]349
350}
Note: See TracBrowser for help on using the repository browser.