source: branches/Release-1-fixes/src/SDAsciiWriter.cc@ 1926

Last change on this file since 1926 was 496, checked in by kil064, 20 years ago

write header and data into seaprate files
write abcissa to data file

  • 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 496 2005-02-22 00:21:13Z kil064 $
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
44#include <coordinates/Coordinates/CoordinateUtil.h>
45#include <coordinates/Coordinates/SpectralCoordinate.h>
46#include <coordinates/Coordinates/DirectionCoordinate.h>
47#include <coordinates/Coordinates/StokesCoordinate.h>
48
49#include <measures/Measures/MEpoch.h>
50
51#include <tables/Tables/Table.h>
52#include <tables/Tables/ScalarColumn.h>
53#include <tables/Tables/ArrayColumn.h>
54
55#include "SDDefs.h"
56#include "SDContainer.h"
57#include "SDMemTable.h"
58#include "SDAsciiWriter.h"
59
60using namespace casa;
61using namespace asap;
62
63
64SDAsciiWriter::SDAsciiWriter()
65{;}
66
67SDAsciiWriter::~SDAsciiWriter()
68{;}
69
70
71Bool SDAsciiWriter::write(const SDMemTable& sdTable, const String& fileName, Bool toStokes)
72{
73
74// Get global Header from Table
75
76 SDHeader header = sdTable.getSDHeader();
77 MEpoch::Ref timeRef(MEpoch::UTC); // Should be in header
78 MDirection::Types dirRef(MDirection::J2000); // Should be in header
79
80// Column keywords
81
82 Table tab = sdTable.table();
83 ROArrayColumn<Double> dir(tab, String("DIRECTION"));
84 ROScalarColumn<Double> time(tab, "TIME");
85 ROArrayColumn<uInt> freqid(tab, "FREQID");
86 ROScalarColumn<String> src(tab, "SRCNAME");
87
88// Temps
89
90 Vector<Int> whichStokes(1,1);
91 Array<Double> whichDir;
92 Vector<Double> lonLat(2);
93 IPosition posDir(2,0);
94 const Unit RAD(String("rad"));
95
96// Open and write header file
97
98
99 String rootName(fileName);
100 if (rootName.length()==0) rootName = String("ascii");
101 {
102 String fName = String(rootName) + String("_header.txt");
103 cout << "Writing header to " << fName << endl;
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 cout << "Wrote " << nRows << " rows into file " << fileName << endl;
198//
199 return True;
200}
201
202
203
204String SDAsciiWriter::formatDirection(const Vector<Double>& lonLat)
205{
206 MVAngle x1(lonLat(0));
207 String s1 = x1.string(MVAngle::TIME, 12);
208//
209 MVAngle x2(lonLat(1));
210 String s2 = x2.string(MVAngle::ANGLE, 12);
211//
212 String ss = s1 + String(" ") + s2;
213 return ss;
214}
215
Note: See TracBrowser for help on using the repository browser.