source: branches/Release12/src/SDAsciiWriter.cc@ 1423

Last change on this file since 1423 was 717, checked in by mar637, 19 years ago

implemented use of SDLog

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 5.9 KB
Line 
1//#---------------------------------------------------------------------------
2//# SDAsciiWriter.cc: ASAP class to write out single dish spectra as FITS images
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: SDAsciiWriter.cc 717 2005-11-17 03:37:54Z mar637 $
30//#---------------------------------------------------------------------------
31
32#include <casa/aips.h>
33#include <casa/Arrays/Array.h>
34#include <casa/Arrays/Vector.h>
35#include <casa/Arrays/VectorIter.h>
36#include <casa/Utilities/CountedPtr.h>
37#include <casa/Quanta/Quantum.h>
38#include <casa/Quanta/MVAngle.h>
39#include <casa/Utilities/Assert.h>
40
41#include <casa/iostream.h>
42#include <casa/fstream.h>
43#include <casa/sstream.h>
44
45#include <coordinates/Coordinates/CoordinateUtil.h>
46#include <coordinates/Coordinates/SpectralCoordinate.h>
47#include <coordinates/Coordinates/DirectionCoordinate.h>
48#include <coordinates/Coordinates/StokesCoordinate.h>
49
50#include <measures/Measures/MEpoch.h>
51
52#include <tables/Tables/Table.h>
53#include <tables/Tables/ScalarColumn.h>
54#include <tables/Tables/ArrayColumn.h>
55
56#include "SDDefs.h"
57#include "SDContainer.h"
58#include "SDMemTable.h"
59#include "SDAsciiWriter.h"
60
61using namespace casa;
62using namespace asap;
63
64
65SDAsciiWriter::SDAsciiWriter()
66{;}
67
68SDAsciiWriter::~SDAsciiWriter()
69{;}
70
71
72Bool SDAsciiWriter::write(const SDMemTable& sdTable, const String& fileName, Bool toStokes)
73{
74
75// Get global Header from Table
76
77 SDHeader header = sdTable.getSDHeader();
78 MEpoch::Ref timeRef(MEpoch::UTC); // Should be in header
79 MDirection::Types dirRef(MDirection::J2000); // Should be in header
80
81// Column keywords
82
83 Table tab = sdTable.table();
84 ROArrayColumn<Double> dir(tab, String("DIRECTION"));
85 ROScalarColumn<Double> time(tab, "TIME");
86 ROArrayColumn<uInt> freqid(tab, "FREQID");
87 ROScalarColumn<String> src(tab, "SRCNAME");
88
89// Temps
90
91 Vector<Int> whichStokes(1,1);
92 Array<Double> whichDir;
93 Vector<Double> lonLat(2);
94 IPosition posDir(2,0);
95 const Unit RAD(String("rad"));
96
97// Open and write header file
98
99 String rootName(fileName);
100 if (rootName.length()==0) rootName = String("ascii");
101 {
102 String fName = String(rootName) + String("_header.txt");
103 pushLog("Writing header to "+fName);
104 ofstream of(fName.chars(), ios::trunc);
105 std::string summary = sdTable.summary(true);
106 of << summary;
107 of.close();
108 }
109
110// Open data file
111
112 String fName = rootName + String(".txt");
113 ofstream of(fName.chars(), ios::trunc);
114
115// Write header
116
117 of << "row beam IF pol source longitude latitude time nchan spectrum mask"
118 << endl;
119
120// Loop over rows
121
122 const uInt nRows = sdTable.nRow();
123 for (uInt iRow=0; iRow<nRows; iRow++) {
124
125// Get data
126
127 const MaskedArray<Float>& dataIn(sdTable.rowAsMaskedArray(iRow,toStokes));
128 const Array<Float>& values = dataIn.getArray();
129 const Array<Bool>& mask = dataIn.getMask();
130
131// Get abcissa
132
133 std::vector<double> abcissa = sdTable.getAbcissa(Int(iRow));
134 const uInt n = abcissa.size();
135
136// Iterate through data in this row by spectra
137
138 ReadOnlyVectorIterator<Float> itData(values, asap::ChanAxis);
139 ReadOnlyVectorIterator<Bool> itMask(mask, asap::ChanAxis);
140 while (!itData.pastEnd()) {
141 const IPosition& pos = itData.pos();
142 AlwaysAssert(itData.vector().nelements()==n,AipsError);
143
144// FreqID
145
146 Vector<uInt> iTmp;
147 freqid.get(iRow, iTmp);
148
149// Direction
150
151 dir.get(iRow, whichDir);
152 posDir(0) = pos(asap::BeamAxis);
153 posDir(1) = 0;
154 lonLat[0] = whichDir(posDir);
155//
156 posDir(0) = pos(asap::BeamAxis);
157 posDir(1) = 1;
158 lonLat[1] = whichDir(posDir);
159
160// Write preamble
161
162 of << iRow << " " << pos(asap::BeamAxis) << " " << pos(asap::IFAxis) << " " <<
163 pos(asap::PolAxis) << " " <<
164 src(iRow) << " " << formatDirection(lonLat) << " " <<
165 sdTable.getTime(iRow,True) << " " << n << " ";
166
167// Write abcissa
168
169 of.setf(ios::fixed, ios::floatfield);
170 of.precision(4);
171 for (uInt i=0; i<n; i++) {
172 of << abcissa[i] << " ";
173 }
174
175// Write data
176
177 const Vector<Float>& data = itData.vector();
178 const Vector<Bool>& mask = itMask.vector();
179 for (uInt i=0; i<n; i++) {
180 of << data[i] << " ";
181 }
182// Write mask
183
184 for (uInt i=0; i<n; i++) {
185 of << mask[i] << " ";
186 }
187 of << endl;
188
189// Next spectrum
190
191 itData.next();
192 itMask.next();
193 }
194 }
195
196 of.close();
197 ostringstream oss;
198 oss << "Wrote " << nRows << " rows into file " << fileName;
199 pushLog(String(oss));
200 return True;
201}
202
203
204
205String SDAsciiWriter::formatDirection(const Vector<Double>& lonLat)
206{
207 MVAngle x1(lonLat(0));
208 String s1 = x1.string(MVAngle::TIME, 12);
209
210 MVAngle x2(lonLat(1));
211 String s2 = x2.string(MVAngle::ANGLE, 12);
212
213 String ss = s1 + String(" ") + s2;
214 return ss;
215}
216
Note: See TracBrowser for help on using the repository browser.