source: trunk/src/STWriter.cpp@ 2603

Last change on this file since 2603 was 2577, checked in by Takeshi Nakazato, 12 years ago

New Development: No

JIRA Issue: No

Ready for Test: Yes

Interface Changes: 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...

Defined PKSMS option in cmake.
if PKSMS=ON, PKSMSreader/writer will be built while if PKSMS=OFF,
those classes will be ignored. Default is OFF.


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