source: trunk/src/STWriter.cpp@ 3122

Last change on this file since 3122 was 3106, checked in by Takeshi Nakazato, 8 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...


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