source: trunk/src/STWriter.cpp@ 2011

Last change on this file since 2011 was 2007, checked in by Takeshi Nakazato, 14 years ago

New Development: No

JIRA Issue: No

Ready for Test: Yes/No

Interface Changes: Yes/No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No

Module(s): Module Names change impacts.

Description: Describe your changes here...

Bug fix on OBS_MODE.


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