source: trunk/src/STWriter.cpp@ 1508

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

fixed unititialised vector length vor convertArray

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