source: trunk/src/STWriter.cpp@ 2657

Last change on this file since 2657 was 2657, checked in by Malte Marquarding, 13 years ago

Ticket #280: Use FITSSpectralUtil::specsysFromFrame to write correct FITS output frame

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 13.6 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 2657 2012-10-10 01:17:12Z 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/PKSSDwriter.h>
43#include <atnf/PKSIO/SrcType.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 <fits/FITS/FITSSpectralUtil.h>
51
52#include "STFITSImageWriter.h"
53#include "STAsciiWriter.h"
54#include "STHeader.h"
55#include "STMath.h"
56
57
58#include "STWriter.h"
59
60using namespace casa;
61namespace asap {
62
63STWriter::STWriter(const std::string &format)
64{
65 format_ = format;
66 String t(format_);
67 t.upcase();
68 if (t == "MS2") {
69 throw (AipsError("MS2 OUTPUT FORMAT IS NO LONGER SUPPORTED"));
70 } else if (t == "SDFITS") {
71 writer_ = new PKSSDwriter();
72 } else if (t == "ASCII" || t == "FITS" || t == "CLASS") {
73 writer_ = 0;
74 } else {
75 throw (AipsError("Unrecognized export format"));
76 }
77}
78
79STWriter::~STWriter()
80{
81 if (writer_) {
82 delete writer_;
83 }
84}
85
86Int STWriter::setFormat(const std::string &format)
87{
88 if (format != format_) {
89 if (writer_) delete writer_;
90 }
91
92 format_ = format;
93 String t(format_);
94 t.upcase();
95 if (t== "MS2") {
96 throw (AipsError("MS2 OUTPUT FORMAT IS NO LONGER SUPPORTED"));
97 } else if (t== "SDFITS") {
98 writer_ = new PKSSDwriter();
99 } else if (t == "ASCII" || t == "FITS" || t == "CLASS") {
100 writer_ = 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 // If we write out foreign formats we have to convert the frequency system
111 // into the output frame, as we do everything related to SPectarlCoordinates
112 // in asap on-the-fly.
113
114 String freqframe;
115 FITSSpectralUtil::specsysFromFrame(freqframe,
116 in->frequencies().getFrame(true));
117 CountedPtr<Scantable> inst = in;
118 if (in->frequencies().getFrame(true) != in->frequencies().getFrame(false)) {
119 FITSSpectralUtil::specsysFromFrame(freqframe,
120 in->frequencies().getFrame(false));
121
122 STMath stm(false);
123 inst = stm.frequencyAlign(in);
124 }
125
126 if (format_=="ASCII") {
127 STAsciiWriter iw;
128 // ASCII calls SpectralCoordinate::toWorld so no freqAlign use 'in'
129 if (iw.write(*in, filename)) {
130 return 0;
131 } else {
132 return 1;
133 }
134 } else if ( format_ == "FITS" || format_ == "CLASS") {
135 STFITSImageWriter iw;
136 if (format_ == "CLASS") {
137 iw.setClass(True);
138 }
139 if (iw.write(*inst, filename)) {
140 return 0;
141 }
142 }
143
144 // MS or SDFITS
145
146 // Extract the header from the table.
147 // this is a little different from what I have done
148 // before. Need to check with the Offline User Test data
149 STHeader hdr = inst->getHeader();
150 std::vector<uint> ifs = inst->getIFNos();
151 int nIF = inst->nif();//ifs.size();
152 Vector<uInt> nPol(nIF),nChan(nIF);
153 Vector<Bool> havexpol(nIF);
154 String fluxUnit = hdr.fluxunit;
155 fluxUnit.upcase();
156 nPol = 0;nChan = 0; havexpol = False;
157 for (uint i=0;i<ifs.size();++i) {
158 if ( inst->npol() > 2 ) {
159 havexpol(ifs[i]) = True;
160 nPol(ifs[i]) = 2;
161 } else {
162 nPol(ifs[i]) = inst->npol();
163 }
164 nChan(ifs[i]) = inst->nchan(ifs[i]);
165 }
166// Vector<String> obstypes(2);
167// obstypes(0) = "TR";//on
168// obstypes(1) = "RF TR";//off
169 const Table table = inst->table();
170
171
172 // Create the output file and write static data.
173 Int status;
174 status = writer_->create(String(filename), hdr.observer, hdr.project,
175 inst->getAntennaName(), hdr.antennaposition,
176 hdr.obstype, hdr.fluxunit,
177 hdr.equinox, freqframe,
178 nChan, nPol, havexpol, False);
179 if ( status ) {
180 throw(AipsError("Failed to create output file"));
181 }
182
183 Int count = 0;
184 PKSrecord pksrec;
185 pksrec.scanNo = 1;
186 // use spearate iterators to ensure renumbering of all numbers
187 TableIterator scanit(table, "SCANNO");
188 while (!scanit.pastEnd() ) {
189 Table stable = scanit.table();
190 TableIterator beamit(stable, "BEAMNO");
191 pksrec.beamNo = 1;
192 while (!beamit.pastEnd() ) {
193 Table btable = beamit.table();
194 //MDirection::ScalarColumn dirCol(btable, "DIRECTION");
195 //pksrec.direction = dirCol(0).getAngle("rad").getValue();
196 TableIterator cycit(btable, "CYCLENO");
197 ROArrayColumn<Double> srateCol(btable, "SCANRATE");
198 Vector<Double> sratedbl;
199 srateCol.get(0, sratedbl);
200 Vector<Float> srateflt(sratedbl.nelements());
201 convertArray(srateflt, sratedbl);
202 //pksrec.scanRate = srateflt;
203 pksrec.scanRate = sratedbl;
204 ROArrayColumn<Double> spmCol(btable, "SRCPROPERMOTION");
205 spmCol.get(0, pksrec.srcPM);
206 ROArrayColumn <Double> sdirCol(btable, "SRCDIRECTION");
207 sdirCol.get(0, pksrec.srcDir);
208 ROScalarColumn<Double> svelCol(btable, "SRCVELOCITY");
209 svelCol.get(0, pksrec.srcVel);
210 ROScalarColumn<uInt> bCol(btable, "BEAMNO");
211 pksrec.beamNo = bCol(0)+1;
212 pksrec.cycleNo = 1;
213 while (!cycit.pastEnd() ) {
214 Table ctable = cycit.table();
215 TableIterator typeit( ctable, "SRCTYPE" ) ;
216 while(!typeit.pastEnd() ) {
217 Table ttable = typeit.table() ;
218 TableIterator ifit(ttable, "IFNO",
219 TableIterator::Ascending, TableIterator::HeapSort);
220 MDirection::ScalarColumn dirCol(ctable, "DIRECTION");
221 pksrec.direction = dirCol(0).getAngle("rad").getValue();
222 pksrec.IFno = 1;
223 while (!ifit.pastEnd() ) {
224 Table itable = ifit.table();
225 TableRow row(itable);
226 // use the first row to fill in all the "metadata"
227 const TableRecord& rec = row.get(0);
228 ROArrayColumn<Float> specCol(itable, "SPECTRA");
229 pksrec.IFno = rec.asuInt("IFNO")+1;
230 uInt nchan = specCol(0).nelements();
231 Double crval,crpix;
232 //Vector<Double> restfreq;
233 Float tmp0,tmp1,tmp2,tmp3,tmp4;
234 String tcalt;
235 Vector<String> stmp0, stmp1;
236 inst->frequencies().getEntry(crpix,crval, pksrec.freqInc,
237 rec.asuInt("FREQ_ID"));
238 inst->focus().getEntry(pksrec.parAngle, pksrec.focusAxi,
239 pksrec.focusTan,
240 pksrec.focusRot, tmp0,tmp1,tmp2,tmp3,tmp4,
241 rec.asuInt("FOCUS_ID"));
242 inst->molecules().getEntry(pksrec.restFreq,stmp0,stmp1,
243 rec.asuInt("MOLECULE_ID"));
244// inst->tcal().getEntry(pksrec.tcalTime, pksrec.tcal,
245// rec.asuInt("TCAL_ID"));
246 inst->weather().getEntry(pksrec.temperature, pksrec.pressure,
247 pksrec.humidity, pksrec.windSpeed,
248 pksrec.windAz, rec.asuInt("WEATHER_ID"));
249 Double pixel = Double(nchan/2);
250 pksrec.refFreq = (pixel-crpix)*pksrec.freqInc + crval;
251 // ok, now we have nrows for the n polarizations in this table
252 polConversion(pksrec.spectra, pksrec.flagged, pksrec.xPol, itable);
253 pksrec.tsys = tsysFromTable(itable);
254 // dummy data
255 uInt npol = pksrec.spectra.ncolumn();
256
257 // TCAL
258 inst->tcal().getEntry( pksrec.tcalTime, pksrec.tcal,
259 rec.asuInt("TCAL_ID") ) ;
260 if ( pksrec.tcal.nelements() == 1 ) {
261 ROScalarColumn<uInt> uintCol( itable, "TCAL_ID" ) ;
262 Vector<uInt> tcalids = uintCol.getColumn() ;
263 pksrec.tcal.resize( npol ) ;
264 Vector<Float> dummyA ;
265 String dummyS ;
266 for ( uInt ipol = 0 ; ipol < npol ; ipol++ ) {
267 inst->tcal().getEntry( dummyS, dummyA, tcalids[ipol] ) ;
268 pksrec.tcal[ipol] = dummyA[0] ;
269 }
270 }
271 else if ( pksrec.tcal.nelements() == nchan ) {
272 ROScalarColumn<uInt> uintCol( itable, "TCAL_ID" ) ;
273 Vector<uInt> tcalids = uintCol.getColumn() ;
274 pksrec.tcal.resize( npol ) ;
275 Vector<Float> dummyA ;
276 String dummyS ;
277 for ( uInt ipol = 0 ; ipol < npol ; ipol++ ) {
278 inst->tcal().getEntry( dummyS, dummyA, tcalids[ipol] ) ;
279 pksrec.tcal[ipol] = mean( dummyA ) ;
280 }
281 }
282 pksrec.mjd = rec.asDouble("TIME");
283 pksrec.interval = rec.asDouble("INTERVAL");
284 pksrec.fieldName = rec.asString("FIELDNAME");
285 pksrec.srcName = rec.asString("SRCNAME");
286 //pksrec.obsType = obstypes[rec.asInt("SRCTYPE")];
287 pksrec.obsType = getObsTypes( rec.asInt("SRCTYPE") ) ;
288 pksrec.bandwidth = nchan * abs(pksrec.freqInc);
289 pksrec.azimuth = rec.asFloat("AZIMUTH");
290 pksrec.elevation = rec.asFloat("ELEVATION");
291 pksrec.refBeam = rec.asInt("REFBEAMNO") + 1;
292 pksrec.sigma.resize(npol);
293 pksrec.sigma = 0.0f;
294 pksrec.calFctr.resize(npol);
295 pksrec.calFctr = 0.0f;
296 pksrec.baseLin.resize(npol,2);
297 pksrec.baseLin = 0.0f;
298 pksrec.baseSub.resize(npol,9);
299 pksrec.baseSub = 0.0f;
300 pksrec.xCalFctr = 0.0;
301 pksrec.flagrow = rec.asuInt("FLAGROW");
302
303 status = writer_->write(pksrec);
304 if ( status ) {
305 writer_->close();
306 throw(AipsError("STWriter: Failed to export Scantable."));
307 }
308 ++count;
309 //++pksrec.IFno;
310 ++ifit;
311 }
312 ++typeit ;
313 }
314 ++pksrec.cycleNo;
315 ++cycit;
316 }
317 //++pksrec.beamNo;
318 ++beamit;
319 }
320 ++pksrec.scanNo;
321 ++scanit;
322 }
323 ostringstream oss;
324 oss << "STWriter: wrote " << count << " rows to " << filename;
325 pushLog(String(oss));
326
327 writer_->close();
328 //if MS2 delete POINTING table exists and copy the one in the keyword
329 if ( format_ == "MS2" ) {
330 replacePtTab(table, filename);
331 }
332 return 0;
333}
334
335Vector<Float> STWriter::tsysFromTable(const Table& tab)
336{
337 ROArrayColumn<Float> tsysCol(tab, "TSYS");
338 Vector<Float> out(tab.nrow());
339 Vector<Float> tmp;
340 for (uInt i=0; i<tab.nrow(); ++i) {
341 tmp.resize();
342 tmp = tsysCol(i);
343 out[i] = tmp[0];
344 }
345 return out;
346}
347
348void STWriter::polConversion( Matrix< Float >& specs, Matrix< uChar >& flags,
349 Vector< Complex > & xpol, const Table & tab )
350{
351 String poltype = tab.keywordSet().asString("POLTYPE");
352// Full stokes is not supported. Just allow stokes I
353 if ( poltype == "stokes" && tab.nrow() != 1) {
354 String msg = "poltype = " + poltype + " not yet supported in output.";
355 throw(AipsError(msg));
356 }
357 ROArrayColumn<Float> specCol(tab, "SPECTRA");
358 ROArrayColumn<uChar> flagCol(tab, "FLAGTRA");
359 uInt nchan = specCol(0).nelements();
360 uInt ncol = (tab.nrow()==1 ? 1: 2 );
361 specs.resize(nchan, ncol);
362 flags.resize(nchan, ncol);
363 // the linears
364 for (uInt i=0; i<ncol; ++i) {
365 specs.column(i) = specCol(i);
366 flags.column(i) = flagCol(i);
367 }
368 // now the complex if exists
369 Bool hasxpol = False;
370 xpol.resize();
371 if ( tab.nrow() == 4 ) {
372 hasxpol = True;
373 xpol.resize(nchan);
374 Vector<Float> reals, imags;
375 reals = specCol(2); imags = specCol(3);
376 for (uInt k=0; k < nchan; ++k) {
377 xpol[k] = Complex(reals[k], imags[k]);
378 }
379 }
380}
381
382// For writing MS data, if there is the reference to
383// original pointing table it replace it by it.
384void STWriter::replacePtTab (const Table& tab, const std::string& fname)
385{
386 String oldPtTabName = fname;
387 oldPtTabName.append("/POINTING");
388 if ( tab.keywordSet().isDefined("POINTING") ) {
389 String PointingTab = tab.keywordSet().asString("POINTING");
390 if ( Table::isReadable(PointingTab) ) {
391 Table newPtTab(PointingTab, Table::Old);
392 newPtTab.copy(oldPtTabName, Table::New);
393 ostringstream oss;
394 oss << "STWriter: copied " <<PointingTab << " to " << fname;
395 pushLog(String(oss));
396 }
397 }
398}
399
400// get obsType string from SRCTYPE value
401String STWriter::getObsTypes( Int srctype )
402{
403 return SrcType::getName(srctype) ;
404}
405
406}
Note: See TracBrowser for help on using the repository browser.