source: trunk/src/STWriter.cpp@ 1575

Last change on this file since 1575 was 1575, checked in by Malte Marquarding, 15 years ago

Ticket #168; align data to convert frequency frame before export

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