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 988 2006-04-05 04:31:47Z 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 |
|
---|
55 | using namespace casa;
|
---|
56 | namespace asap {
|
---|
57 |
|
---|
58 | STWriter::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 |
|
---|
76 | STWriter::~STWriter()
|
---|
77 | {
|
---|
78 | if (writer_) {
|
---|
79 | delete writer_;
|
---|
80 | }
|
---|
81 | }
|
---|
82 |
|
---|
83 | Int 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 |
|
---|
106 | Int 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 | if (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 | 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(stable, "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 cycit(btable, "CYCLENO");
|
---|
172 | Int cycno = 1;
|
---|
173 | while (!cycit.pastEnd() ) {
|
---|
174 | Table ctable = cycit.table();
|
---|
175 | TableIterator ifit(ctable, "IFNO");
|
---|
176 | Int ifno = 1;
|
---|
177 | while (!ifit.pastEnd() ) {
|
---|
178 | Table itable = ifit.table();
|
---|
179 | TableRow row(itable);
|
---|
180 | // use the first row to fill in all the "metadata"
|
---|
181 | const TableRecord& rec = row.get(0);
|
---|
182 | ROArrayColumn<Float> specCol(itable, "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 | Float tmp0,tmp1,tmp2,tmp3,tmp4;
|
---|
188 | Vector<Float> tcalval;
|
---|
189 | String stmp0,stmp1, tcalt;
|
---|
190 | in->frequencies().getEntry(crpix,crval,cdelt, rec.asuInt("FREQ_ID"));
|
---|
191 | in->focus().getEntry(focusAxi, focusTan, focusRot,
|
---|
192 | tmp0,tmp1,tmp2,tmp3,tmp4,
|
---|
193 | rec.asuInt("FOCUS_ID"));
|
---|
194 | in->molecules().getEntry(restfreq,stmp0,stmp1,rec.asuInt("MOLECULE_ID"));
|
---|
195 | in->tcal().getEntry(tcalt,tcalval,rec.asuInt("TCAL_ID"));
|
---|
196 | in->weather().getEntry(temperature, pressure, humidity,
|
---|
197 | windSpeed, windAz,
|
---|
198 | rec.asuInt("WEATHER_ID"));
|
---|
199 | Double pixel = Double(nchan/2);
|
---|
200 | Double refFreqNew = (pixel-crpix)*cdelt + crval;
|
---|
201 | // ok, now we have nrows for the n polarizations in this table
|
---|
202 | Matrix<Float> specs;
|
---|
203 | Matrix<uChar> flags;
|
---|
204 | Vector<Complex> xpol;
|
---|
205 | polConversion(specs, flags, xpol, itable);
|
---|
206 | Vector<Float> tsys = tsysFromTable(itable);
|
---|
207 | // dummy data
|
---|
208 | uInt npol = specs.ncolumn();
|
---|
209 |
|
---|
210 | Matrix<Float> baseLin(npol,2, 0.0f);
|
---|
211 | Matrix<Float> baseSub(npol,9, 0.0f);
|
---|
212 | Complex xCalFctr;
|
---|
213 | Vector<Double> scanRate(2, 0.0);
|
---|
214 | Vector<Float> sigma(npol, 0.0f);
|
---|
215 | Vector<Float> calFctr(npol, 0.0f);
|
---|
216 | if (status = writer_->write(scanno, cycno, rec.asDouble("TIME"),
|
---|
217 | rec.asDouble("INTERVAL"),
|
---|
218 | rec.asString("FIELDNAME"),
|
---|
219 | rec.asString("SRCNAME"),
|
---|
220 | direction,
|
---|
221 | srcPM, srcVel, // not in scantable yet
|
---|
222 | ifno,
|
---|
223 | refFreqNew, nchan*abs(cdelt), cdelt,
|
---|
224 | restfreq,
|
---|
225 | tcal,
|
---|
226 | tcalt,
|
---|
227 | rec.asFloat("AZIMUTH"),
|
---|
228 | rec.asFloat("ELEVATION"),
|
---|
229 | rec.asFloat("PARANGLE"),
|
---|
230 | focusAxi, focusTan, focusRot,
|
---|
231 | temperature,
|
---|
232 | pressure, humidity, windSpeed, windAz,
|
---|
233 | rec.asInt("REFBEAMNO")+1, beamno,
|
---|
234 | direction,
|
---|
235 | scanRate,// not in scantable
|
---|
236 | tsys,
|
---|
237 | sigma, calFctr,// not in scantable
|
---|
238 | baseLin, baseSub,// not in scantable
|
---|
239 | specs, flags,
|
---|
240 | xCalFctr,//
|
---|
241 | xpol)
|
---|
242 | ) {
|
---|
243 | writer_->close();
|
---|
244 | throw(AipsError("STWriter: Failed to export Scantable."));
|
---|
245 | }
|
---|
246 |
|
---|
247 | ++ifno;
|
---|
248 | ++ifit;
|
---|
249 | }
|
---|
250 | ++cycno;
|
---|
251 | ++cycit;
|
---|
252 | }
|
---|
253 | ++beamno;
|
---|
254 | ++beamit;
|
---|
255 | }
|
---|
256 | ++scanno;
|
---|
257 | ++scanit;
|
---|
258 | }
|
---|
259 | ostringstream oss;
|
---|
260 | oss << "STWriter: wrote " << count << " rows to " << filename << endl;
|
---|
261 | pushLog(String(oss));
|
---|
262 | writer_->close();
|
---|
263 |
|
---|
264 | return 0;
|
---|
265 | }
|
---|
266 |
|
---|
267 | Vector<Float> STWriter::tsysFromTable(const Table& tab)
|
---|
268 | {
|
---|
269 | ROArrayColumn<Float> tsysCol(tab, "TSYS");
|
---|
270 | Vector<Float> out(tab.nrow());
|
---|
271 | Vector<Float> tmp;
|
---|
272 | for (uInt i=0; i<tab.nrow(); ++i) {
|
---|
273 | tmp.resize();
|
---|
274 | tmp = tsysCol(i);
|
---|
275 | out[i] = tmp[0];
|
---|
276 | }
|
---|
277 | return out;
|
---|
278 | }
|
---|
279 |
|
---|
280 | void STWriter::polConversion( Matrix< Float >& specs, Matrix< uChar >& flags,
|
---|
281 | Vector< Complex > & xpol, const Table & tab )
|
---|
282 | {
|
---|
283 | TableRow row(tab);
|
---|
284 | String poltype = tab.keywordSet().asString("POLTYPE");
|
---|
285 | if ( poltype != "linear") {
|
---|
286 | String msg = "poltype = " + poltype + " not yet supported in output.";
|
---|
287 | throw(AipsError(msg));
|
---|
288 | }
|
---|
289 | // use the first row to fill in all the "metadata"
|
---|
290 | const TableRecord& rec = row.get(0);
|
---|
291 | ROArrayColumn<Float> specCol(tab, "SPECTRA");
|
---|
292 | ROArrayColumn<uChar> flagCol(tab, "FLAGTRA");
|
---|
293 | uInt nchan = specCol(0).nelements();
|
---|
294 | uInt ncol = (tab.nrow()==1 ? 1: 2 );
|
---|
295 | specs.resize(nchan, ncol);
|
---|
296 | flags.resize(nchan, ncol);
|
---|
297 | // the linears
|
---|
298 | for (uInt i=0; i<ncol; ++i) {
|
---|
299 | specs.column(i) = specCol(i);
|
---|
300 | flags.column(i) = flagCol(i);
|
---|
301 | }
|
---|
302 | // now the complex if exists
|
---|
303 | Bool hasxpol = False;
|
---|
304 | xpol.resize();
|
---|
305 | if ( tab.nrow() == 4 ) {
|
---|
306 | hasxpol = True;
|
---|
307 | xpol.resize(nchan);
|
---|
308 | Vector<Float> reals, imags;
|
---|
309 | reals = specCol(2); imags = specCol(3);
|
---|
310 | for (uInt k=0; k < nchan; ++k) {
|
---|
311 | xpol[k] = Complex(reals[k], imags[k]);
|
---|
312 | }
|
---|
313 | }
|
---|
314 | }
|
---|
315 |
|
---|
316 |
|
---|
317 | }
|
---|