source: branches/alma/src/STWriter.cpp@ 2793

Last change on this file since 2793 was 1757, checked in by Kana Sugimoto, 14 years ago

New Development: Yes

JIRA Issue: Yes (CAS-2211)

Ready for Test: Yes

Interface Changes: Yes

What Interface Changed: ASAP 3.0.0 interface changes

Test Programs:

Put in Release Notes: Yes

Module(s): all the CASA sd tools and tasks are affected.

Description: Merged ATNF-ASAP 3.0.0 developments to CASA (alma) branch.

Note you also need to update casa/code/atnf.


  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 13.6 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 1757 2010-06-09 09:03:06Z KanaSugimoto $
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
[1603]41#include <atnf/PKSIO/PKSrecord.h>
[28]42#include <atnf/PKSIO/PKSMS2writer.h>
43#include <atnf/PKSIO/PKSSDwriter.h>
[1757]44#include <atnf/PKSIO/SrcType.h>
[28]45
[827]46#include <tables/Tables/Table.h>
47#include <tables/Tables/TableIter.h>
48#include <tables/Tables/TableRow.h>
[324]49#include <tables/Tables/ArrayColumn.h>
50
[1603]51#include "STFITSImageWriter.h"
[988]52#include "STAsciiWriter.h"
[901]53#include "STHeader.h"
[1757]54#include "STMath.h"
[324]55
[1757]56
[822]57#include "STWriter.h"
[28]58
[125]59using namespace casa;
[827]60namespace asap {
[28]61
[822]62STWriter::STWriter(const std::string &format)
[28]63{
[988]64 format_ = format;
65 String t(format_);
[450]66 t.upcase();
[1603]67 if (t == "MS2") {
[988]68 writer_ = new PKSMS2writer();
[1603]69 } else if (t == "SDFITS") {
[988]70 writer_ = new PKSSDwriter();
[1603]71 } else if (t == "ASCII" || t == "FITS" || t == "CLASS") {
[988]72 writer_ = 0;
[450]73 } else {
[988]74 throw (AipsError("Unrecognized export format"));
[28]75 }
76}
77
[822]78STWriter::~STWriter()
[28]79{
[988]80 if (writer_) {
81 delete writer_;
[192]82 }
[28]83}
84
[822]85Int STWriter::setFormat(const std::string &format)
[28]86{
[988]87 if (format != format_) {
88 if (writer_) delete writer_;
[450]89 }
[827]90
[988]91 format_ = format;
92 String t(format_);
[450]93 t.upcase();
94 if (t== "MS2") {
[988]95 writer_ = new PKSMS2writer();
[450]96 } else if (t== "SDFITS") {
[988]97 writer_ = new PKSSDwriter();
[1603]98 } else if (t == "ASCII" || t == "FITS" || t == "CLASS") {
[988]99 writer_ = 0;
[450]100 } else {
101 throw (AipsError("Unrecognized Format"));
[28]102 }
103 return 0;
104}
105
[827]106Int STWriter::write(const CountedPtr<Scantable> in,
107 const std::string &filename)
[28]108{
[1757]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.
[192]112
[1757]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
[1280]119 if (format_=="ASCII") {
120 STAsciiWriter iw;
[1757]121 // ASCII calls SpectralCoordinate::toWorld so no freqAlign use 'in'
[1280]122 if (iw.write(*in, filename)) {
123 return 0;
124 } else {
125 return 1;
126 }
[1603]127 } else if ( format_ == "FITS" || format_ == "CLASS") {
128 STFITSImageWriter iw;
129 if (format_ == "CLASS") {
130 iw.setClass(True);
131 }
[1757]132 if (iw.write(*inst, filename)) {
[1603]133 return 0;
134 }
[199]135 }
[192]136
[827]137 // MS or SDFITS
[192]138
[28]139 // Extract the header from the table.
[1387]140 // this is a little different from what I have done
141 // before. Need to check with the Offline User Test data
[1757]142 STHeader hdr = inst->getHeader();
[1060]143 //const Int nPol = hdr.npol;
144 //const Int nChan = hdr.nchan;
[1757]145 std::vector<uint> ifs = inst->getIFNos();
146 int nIF = inst->nif();//ifs.size();
[1060]147 Vector<uInt> nPol(nIF),nChan(nIF);
148 Vector<Bool> havexpol(nIF);
[1387]149 String fluxUnit = hdr.fluxunit;
[1757]150 fluxUnit.upcase();
[1295]151 nPol = 0;nChan = 0; havexpol = False;
[1305]152 for (uint i=0;i<ifs.size();++i) {
[1757]153 nPol(ifs[i]) = inst->npol();
154 nChan(ifs[i]) = inst->nchan(ifs[i]);
[1295]155 havexpol(ifs[i]) = nPol(ifs[i]) > 2;
[1060]156 }
[1757]157// Vector<String> obstypes(2);
158// obstypes(0) = "TR";//on
159// obstypes(1) = "RF TR";//off
160 const Table table = inst->table();
[28]161
[324]162
[28]163 // Create the output file and write static data.
164 Int status;
[1060]165 status = writer_->create(String(filename), hdr.observer, hdr.project,
[1603]166 hdr.antennaname, hdr.antennaposition,
167 hdr.obstype, hdr.fluxunit,
168 hdr.equinox, hdr.freqref,
169 nChan, nPol, havexpol, False);
[996]170 if ( status ) {
[717]171 throw(AipsError("Failed to create output file"));
[28]172 }
173
174 Int count = 0;
[1603]175 PKSrecord pksrec;
176 pksrec.scanNo = 1;
[822]177 // use spearate iterators to ensure renumbering of all numbers
178 TableIterator scanit(table, "SCANNO");
179 while (!scanit.pastEnd() ) {
180 Table stable = scanit.table();
[988]181 TableIterator beamit(stable, "BEAMNO");
[1603]182 pksrec.beamNo = 1;
[822]183 while (!beamit.pastEnd() ) {
184 Table btable = beamit.table();
[1672]185 //MDirection::ScalarColumn dirCol(btable, "DIRECTION");
186 //pksrec.direction = dirCol(0).getAngle("rad").getValue();
[988]187 TableIterator cycit(btable, "CYCLENO");
[999]188 ROArrayColumn<Double> srateCol(btable, "SCANRATE");
[1603]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;
[999]195 ROArrayColumn<Double> spmCol(btable, "SRCPROPERMOTION");
[1603]196 spmCol.get(0, pksrec.srcPM);
[999]197 ROArrayColumn <Double> sdirCol(btable, "SRCDIRECTION");
[1603]198 sdirCol.get(0, pksrec.srcDir);
[999]199 ROScalarColumn<Double> svelCol(btable, "SRCVELOCITY");
[1603]200 svelCol.get(0, pksrec.srcVel);
[1295]201 ROScalarColumn<uInt> bCol(btable, "BEAMNO");
[1603]202 pksrec.beamNo = bCol(0)+1;
203 pksrec.cycleNo = 1;
[988]204 while (!cycit.pastEnd() ) {
205 Table ctable = cycit.table();
[1683]206 TableIterator ifit(ctable, "IFNO", TableIterator::Ascending, TableIterator::HeapSort);
[1672]207 MDirection::ScalarColumn dirCol(ctable, "DIRECTION");
208 pksrec.direction = dirCol(0).getAngle("rad").getValue();
[1603]209 pksrec.IFno = 1;
[988]210 while (!ifit.pastEnd() ) {
211 Table itable = ifit.table();
212 TableRow row(itable);
[822]213 // use the first row to fill in all the "metadata"
[827]214 const TableRecord& rec = row.get(0);
[988]215 ROArrayColumn<Float> specCol(itable, "SPECTRA");
[1603]216 pksrec.IFno = rec.asuInt("IFNO")+1;
[822]217 uInt nchan = specCol(0).nelements();
[1603]218 Double crval,crpix;
219 //Vector<Double> restfreq;
[988]220 Float tmp0,tmp1,tmp2,tmp3,tmp4;
[1446]221 String tcalt;
222 Vector<String> stmp0, stmp1;
[1757]223 inst->frequencies().getEntry(crpix,crval, pksrec.freqInc,
[1603]224 rec.asuInt("FREQ_ID"));
[1757]225 inst->focus().getEntry(pksrec.parAngle, pksrec.focusAxi, pksrec.focusTan,
226 pksrec.focusRot, tmp0,tmp1,tmp2,tmp3,tmp4,
227 rec.asuInt("FOCUS_ID"));
228 inst->molecules().getEntry(pksrec.restFreq,stmp0,stmp1,
[1603]229 rec.asuInt("MOLECULE_ID"));
[1757]230 inst->tcal().getEntry(pksrec.tcalTime, pksrec.tcal,
[1603]231 rec.asuInt("TCAL_ID"));
[1757]232 inst->weather().getEntry(pksrec.temperature, pksrec.pressure,
[1603]233 pksrec.humidity, pksrec.windSpeed,
234 pksrec.windAz, rec.asuInt("WEATHER_ID"));
[822]235 Double pixel = Double(nchan/2);
[1603]236 pksrec.refFreq = (pixel-crpix)*pksrec.freqInc + crval;
[822]237 // ok, now we have nrows for the n polarizations in this table
[1603]238 polConversion(pksrec.spectra, pksrec.flagged, pksrec.xPol, itable);
239 pksrec.tsys = tsysFromTable(itable);
[822]240 // dummy data
[1603]241 uInt npol = pksrec.spectra.ncolumn();
[988]242
[1603]243 pksrec.mjd = rec.asDouble("TIME");
244 pksrec.interval = rec.asDouble("INTERVAL");
245 pksrec.fieldName = rec.asString("FIELDNAME");
246 pksrec.srcName = rec.asString("SRCNAME");
[1757]247 //pksrec.obsType = obstypes[rec.asInt("SRCTYPE")];
248 pksrec.obsType = getObsTypes( rec.asInt("SRCTYPE") ) ;
[1603]249 pksrec.bandwidth = nchan * abs(pksrec.freqInc);
250 pksrec.azimuth = rec.asFloat("AZIMUTH");
251 pksrec.elevation = rec.asFloat("ELEVATION");
252 pksrec.refBeam = rec.asInt("REFBEAMNO") + 1;
253 pksrec.sigma.resize(npol);
254 pksrec.sigma = 0.0f;
255 pksrec.calFctr.resize(npol);
256 pksrec.calFctr = 0.0f;
257 pksrec.baseLin.resize(npol,2);
258 pksrec.baseLin = 0.0f;
259 pksrec.baseSub.resize(npol,9);
260 pksrec.baseSub = 0.0f;
261 pksrec.xCalFctr = 0.0;
[1661]262 pksrec.flagrow = rec.asuInt("FLAGROW");
[1603]263
264 status = writer_->write(pksrec);
[996]265 if ( status ) {
[988]266 writer_->close();
267 throw(AipsError("STWriter: Failed to export Scantable."));
[822]268 }
[999]269 ++count;
[1603]270 //++pksrec.IFno;
[988]271 ++ifit;
[470]272 }
[1603]273 ++pksrec.cycleNo;
[988]274 ++cycit;
[28]275 }
[1603]276 //++pksrec.beamNo;
[822]277 ++beamit;
[28]278 }
[1603]279 ++pksrec.scanNo;
[822]280 ++scanit;
[28]281 }
[717]282 ostringstream oss;
[999]283 oss << "STWriter: wrote " << count << " rows to " << filename;
[717]284 pushLog(String(oss));
[988]285 writer_->close();
[1387]286 //if MS2 delete POINTING table exists and copy the one in the keyword
287 if ( format_ == "MS2" ) {
288 replacePtTab(table, filename);
[1603]289 }
[28]290 return 0;
291}
[470]292
[827]293Vector<Float> STWriter::tsysFromTable(const Table& tab)
294{
295 ROArrayColumn<Float> tsysCol(tab, "TSYS");
296 Vector<Float> out(tab.nrow());
297 Vector<Float> tmp;
298 for (uInt i=0; i<tab.nrow(); ++i) {
[988]299 tmp.resize();
[827]300 tmp = tsysCol(i);
301 out[i] = tmp[0];
302 }
[988]303 return out;
[827]304}
305
306void STWriter::polConversion( Matrix< Float >& specs, Matrix< uChar >& flags,
[822]307 Vector< Complex > & xpol, const Table & tab )
308{
309 String poltype = tab.keywordSet().asString("POLTYPE");
[1757]310// Full stokes is not supported. Just allow stokes I
311 if ( poltype == "stokes" && tab.nrow() != 1) {
[822]312 String msg = "poltype = " + poltype + " not yet supported in output.";
[988]313 throw(AipsError(msg));
314 }
[827]315 ROArrayColumn<Float> specCol(tab, "SPECTRA");
316 ROArrayColumn<uChar> flagCol(tab, "FLAGTRA");
[822]317 uInt nchan = specCol(0).nelements();
[988]318 uInt ncol = (tab.nrow()==1 ? 1: 2 );
319 specs.resize(nchan, ncol);
320 flags.resize(nchan, ncol);
[822]321 // the linears
[988]322 for (uInt i=0; i<ncol; ++i) {
[822]323 specs.column(i) = specCol(i);
324 flags.column(i) = flagCol(i);
325 }
326 // now the complex if exists
327 Bool hasxpol = False;
328 xpol.resize();
[827]329 if ( tab.nrow() == 4 ) {
[822]330 hasxpol = True;
331 xpol.resize(nchan);
332 Vector<Float> reals, imags;
333 reals = specCol(2); imags = specCol(3);
334 for (uInt k=0; k < nchan; ++k) {
[827]335 xpol[k] = Complex(reals[k], imags[k]);
[822]336 }
337 }
338}
[827]339
[1387]340// For writing MS data, if there is the reference to
341// original pointing table it replace it by it.
342void STWriter::replacePtTab (const Table& tab, const std::string& fname)
343{
344 String oldPtTabName = fname;
345 oldPtTabName.append("/POINTING");
346 if ( tab.keywordSet().isDefined("POINTING") ) {
347 String PointingTab = tab.keywordSet().asString("POINTING");
348 if ( Table::isReadable(PointingTab) ) {
349 Table newPtTab(PointingTab, Table::Old);
350 newPtTab.copy(oldPtTabName, Table::New);
351 ostringstream oss;
352 oss << "STWriter: copied " <<PointingTab << " to " << fname;
353 pushLog(String(oss));
354 }
355 }
356}
[827]357
[1757]358// get obsType string from SRCTYPE value
359String STWriter::getObsTypes( Int srctype )
360{
361 String obsType ;
362 switch( srctype ) {
363 case Int(SrcType::PSON):
364 obsType = "PSON" ;
365 break ;
366 case Int(SrcType::PSOFF):
367 obsType = "PSOFF" ;
368 break ;
369 case Int(SrcType::NOD):
370 obsType = "NOD" ;
371 break ;
372 case Int(SrcType::FSON):
373 obsType = "FSON" ;
374 break ;
375 case Int(SrcType::FSOFF):
376 obsType = "FSOFF" ;
377 break ;
378 case Int(SrcType::SKY):
379 obsType = "SKY" ;
380 break ;
381 case Int(SrcType::HOT):
382 obsType = "HOT" ;
383 break ;
384 case Int(SrcType::WARM):
385 obsType = "WARM" ;
386 break ;
387 case Int(SrcType::COLD):
388 obsType = "COLD" ;
389 break ;
390 case Int(SrcType::PONCAL):
391 obsType = "PSON:CALON" ;
392 break ;
393 case Int(SrcType::POFFCAL):
394 obsType = "PSOFF:CALON" ;
395 break ;
396 case Int(SrcType::NODCAL):
397 obsType = "NOD:CALON" ;
398 break ;
399 case Int(SrcType::FONCAL):
400 obsType = "FSON:CALON" ;
401 break ;
402 case Int(SrcType::FOFFCAL):
403 obsType = "FSOFF:CALOFF" ;
404 break ;
405 case Int(SrcType::FSLO):
406 obsType = "FSLO" ;
407 break ;
408 case Int(SrcType::FLOOFF):
409 obsType = "FS:LOWER:OFF" ;
410 break ;
411 case Int(SrcType::FLOSKY):
412 obsType = "FS:LOWER:SKY" ;
413 break ;
414 case Int(SrcType::FLOHOT):
415 obsType = "FS:LOWER:HOT" ;
416 break ;
417 case Int(SrcType::FLOWARM):
418 obsType = "FS:LOWER:WARM" ;
419 break ;
420 case Int(SrcType::FLOCOLD):
421 obsType = "FS:LOWER:COLD" ;
422 break ;
423 case Int(SrcType::FSHI):
424 obsType = "FSHI" ;
425 break ;
426 case Int(SrcType::FHIOFF):
427 obsType = "FS:HIGHER:OFF" ;
428 break ;
429 case Int(SrcType::FHISKY):
430 obsType = "FS:HIGHER:SKY" ;
431 break ;
432 case Int(SrcType::FHIHOT):
433 obsType = "FS:HIGHER:HOT" ;
434 break ;
435 case Int(SrcType::FHIWARM):
436 obsType = "FS:HIGHER:WARM" ;
437 break ;
438 case Int(SrcType::FHICOLD):
439 obsType = "FS:HIGHER:COLD" ;
440 break ;
441 default:
442 obsType = "NOTYPE" ;
443 }
444
445 return obsType ;
[827]446}
[1757]447
448}
Note: See TracBrowser for help on using the repository browser.