source: trunk/src/STWriter.cpp@ 881

Last change on this file since 881 was 844, checked in by mar637, 19 years ago

changed getSDHeader -> getHeader to reflect interface changes in Scantable
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 9.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 844 2006-02-23 00:23:57Z mar637 $
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/PKSMS2writer.h>
42#include <atnf/PKSIO/PKSSDwriter.h>
43
44#include <tables/Tables/Table.h>
45#include <tables/Tables/TableIter.h>
46#include <tables/Tables/TableRow.h>
47#include <tables/Tables/ArrayColumn.h>
48
49//#include "SDFITSImageWriter.h"
50//#include "STAsciiWriter.h"
51//#include "SDPol.h"
52#include "SDContainer.h"
53
54#include "STWriter.h"
55
56using namespace casa;
57namespace asap {
58
59STWriter::STWriter(const std::string &format)
60{
61 cFormat = format;
62 String t(cFormat);
63 t.upcase();
64 if (t== "MS2") {
65 cWriter = new PKSMS2writer();
66 } else if (t== "SDFITS") {
67 cWriter = new PKSSDwriter();
68 } else if (t== "FITS") {
69 cWriter = 0;
70 } else if (t== "ASCII") {
71 cWriter = 0;
72 } else {
73 throw (AipsError("Unrecognized Format"));
74 }
75}
76
77STWriter::~STWriter()
78{
79 if (cWriter) {
80 delete cWriter;
81 }
82}
83
84Int STWriter::setFormat(const std::string &format)
85{
86 if (format != cFormat) {
87 if (cWriter) delete cWriter;
88 }
89
90 cFormat = format;
91 String t(cFormat);
92 t.upcase();
93 if (t== "MS2") {
94 cWriter = new PKSMS2writer();
95 } else if (t== "SDFITS") {
96 cWriter = new PKSSDwriter();
97 } else if (t== "FITS") {
98 cWriter = 0;
99 } else if (t== "ASCII") {
100 cWriter = 0;
101 } else {
102 throw (AipsError("Unrecognized Format"));
103 }
104 return 0;
105}
106
107Int STWriter::write(const CountedPtr<Scantable> in,
108 const std::string &filename)
109{
110
111// Image FITS
112
113 if (cFormat=="FITS") {
114// Bool verbose = True;
115// SDFITSImageWriter iw;
116// if (iw.write(*in, filename, verbose)) {
117// return 0;
118// } else {
119// return 1;
120// }
121 } else if (cFormat=="ASCII") {
122 /*SDAsciiWriter iw;
123 if (iw.write(*in, filename)) {
124 return 0;
125 } else {
126 return 1;
127 }*/
128 }
129
130 // MS or SDFITS
131
132 // Extract the header from the table.
133 SDHeader hdr = in->getHeader();
134 const Int nPol = hdr.npol;
135 const Int nChan = hdr.nchan;
136
137 const Table table = in->table();
138// ROArrayColumn<uInt> freqIDCol(table, "FREQ_ID");
139// Vector<uInt> freqIDs;
140
141 // Create the output file and write static data.
142 Int status;
143 if (status = cWriter->create(filename, hdr.observer, hdr.project,
144 hdr.antennaname, hdr.antennaposition,
145 hdr.obstype, hdr.equinox, hdr.freqref,
146 nChan, nPol, False, False)) {
147 throw(AipsError("Failed to create output file"));
148 }
149
150 Double srcVel = 0.0;
151
152 String fieldName, srcName, tcalTime;
153 Vector<Float> calFctr, sigma, tcal, tsys;
154 Vector<Double> direction(2), scanRate(2), srcDir(2), srcPM(2,0.0);
155 Matrix<Float> spectra;
156 Matrix<uChar> flagtra;
157 Complex xCalFctr;
158 Int count = 0;
159 Int scanno = 1;
160 // use spearate iterators to ensure renumbering of all numbers
161 TableIterator scanit(table, "SCANNO");
162 while (!scanit.pastEnd() ) {
163 Table stable = scanit.table();
164 TableIterator beamit(table, "BEAMNO");
165 Int beamno = 1;
166 while (!beamit.pastEnd() ) {
167 Table btable = beamit.table();
168 // position only varies by beam
169 MDirection::ScalarColumn dirCol(btable, "DIRECTION");
170 Vector<Double> direction = dirCol(0).getAngle("rad").getValue();
171 TableIterator ifit(btable, "IFNO");
172 Int ifno = 1;
173 while (!ifit.pastEnd() ) {
174 Table itable = ifit.table();
175 TableIterator cycit(itable, "CYCLENO");
176 Int cycno = 1;
177 while (!cycit.pastEnd() ) {
178 Table ctable = cycit.table();
179 TableRow row(ctable);
180 // use the first row to fill in all the "metadata"
181 const TableRecord& rec = row.get(0);
182 ROArrayColumn<Float> specCol(ctable, "SPECTRA");
183 uInt nchan = specCol(0).nelements();
184 Double cdelt,crval,crpix, restfreq;
185 Float focusAxi, focusTan, focusRot,
186 temperature, pressure, humidity, windSpeed, windAz;
187 Vector<Float> tcalval;
188 String tmp,tmp2, tcalt;
189 in->frequencies().getEntry(crpix,crval,cdelt, rec.asuInt("FREQ_ID"));
190 in->molecules().getEntry(restfreq,tmp,tmp2,rec.asuInt("RESTFREQ_ID"));
191 in->tcal().getEntry(tcalt,tcalval,rec.asuInt("TCAL_ID"));
192 in->weather().getEntry(temperature, pressure, humidity,
193 windSpeed, windAz,
194 rec.asuInt("WEATHER_ID"));
195 Double pixel = Double(nchan/2);
196 Double refFreqNew = (pixel-crpix)*cdelt + crval;
197 // ok, now we have nrows for the n polarizations in this table
198 Matrix<Float> specs;
199 Matrix<uChar> flags;
200 Vector<Complex> xpol;
201 polConversion(specs, flags, xpol, ctable);
202 Vector<Float> tsys = tsysFromTable(ctable);
203 // dummy data
204 uInt npol = specs.ncolumn();
205 Matrix<Float> baseLin(npol,2, 0.0f);
206 Matrix<Float> baseSub(npol,9, 0.0f);
207 Complex xCalFctr;
208 Vector<Double> scanRate(2, 0.0);
209 Vector<Float> sigma(npol, 0.0f);
210 Vector<Float> calFctr(npol, 0.0f);
211
212
213 if (status = cWriter->write(scanno, cycno, rec.asDouble("TIME"),
214 rec.asDouble("INTERVAL"),
215 rec.asString("FIELDNAME"),
216 rec.asString("SRCNAME"),
217 direction,
218 srcPM, srcVel, // not in scantable yet
219 ifno,
220 refFreqNew, nchan*abs(cdelt), cdelt,
221 restfreq,
222 tcal,
223 tcalt,
224 rec.asFloat("AZIMUTH"),
225 rec.asFloat("ELEVATION"),
226 rec.asFloat("PARANGLE"),
227 focusAxi, focusTan, focusRot,//
228 temperature,
229 pressure, humidity, windSpeed, windAz,
230 rec.asInt("REFBEAM"), beamno,
231 direction,
232 scanRate,// not in scantable
233 tsys,
234 sigma, calFctr,// not in scantable
235 baseLin, baseSub,// not in scantable
236 specs, flags,
237 xCalFctr,//
238 xpol)
239 ) {
240 cerr << "Error writing output file." << endl;
241 return 1;
242 }
243
244 ++cycno;
245 ++cycit;
246 }
247 ++ifno;
248 ++ifit;
249 }
250 ++beamno;
251 ++beamit;
252 }
253 ++scanno;
254 ++scanit;
255 }
256 ostringstream oss;
257 oss << "STWriter: wrote " << count << " rows to " << filename << endl;
258 pushLog(String(oss));
259 cWriter->close();
260
261 return 0;
262}
263
264Vector<Float> STWriter::tsysFromTable(const Table& tab)
265{
266 ROArrayColumn<Float> tsysCol(tab, "TSYS");
267 Vector<Float> out(tab.nrow());
268 Vector<Float> tmp;
269 for (uInt i=0; i<tab.nrow(); ++i) {
270 tmp = tsysCol(i);
271 out[i] = tmp[0];
272 }
273}
274
275void STWriter::polConversion( Matrix< Float >& specs, Matrix< uChar >& flags,
276 Vector< Complex > & xpol, const Table & tab )
277{
278 TableRow row(tab);
279 String poltype = tab.keywordSet().asString("POLTYPE");
280 if ( poltype != "linear")
281 String msg = "poltype = " + poltype + " not yet supported in output.";
282 throw(AipsError("msg"));
283 // use the first row to fill in all the "metadata"
284 const TableRecord& rec = row.get(0);
285 ROArrayColumn<Float> specCol(tab, "SPECTRA");
286 ROArrayColumn<uChar> flagCol(tab, "FLAGTRA");
287 uInt nchan = specCol(0).nelements();
288 uInt ncols = (tab.nrow()==1 ? 1: 2 );
289 specs.resize(nchan, ncols);
290 flags.resize(nchan, ncols);
291 // the linears
292 for (uInt i=0; i<ncols; ++i) {
293 specs.column(i) = specCol(i);
294 flags.column(i) = flagCol(i);
295 }
296 // now the complex if exists
297 Bool hasxpol = False;
298 xpol.resize();
299 if ( tab.nrow() == 4 ) {
300 hasxpol = True;
301 xpol.resize(nchan);
302 Vector<Float> reals, imags;
303 reals = specCol(2); imags = specCol(3);
304 for (uInt k=0; k < nchan; ++k) {
305 xpol[k] = Complex(reals[k], imags[k]);
306 }
307 }
308}
309
310
311}
Note: See TracBrowser for help on using the repository browser.