source: branches/Release2.0/src/STWriter.cpp@ 1734

Last change on this file since 1734 was 999, checked in by mar637, 19 years ago

added "redundant" fields from reader. The aren't used in asap but should be passed on for consistency.

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