source: branches/Release-1-fixes/src/SDFITSImageWriter.cc@ 1898

Last change on this file since 1898 was 448, checked in by kil064, 20 years ago

track change to SDmemTable::rowAsMaskedArray interface

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 7.9 KB
Line 
1//#---------------------------------------------------------------------------
2//# SDFITSImageWriter.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: SDFITSImageWriter.cc 448 2005-02-15 07:39:54Z kil064 $
30//#---------------------------------------------------------------------------
31
32#include "SDFITSImageWriter.h"
33
34#include <casa/aips.h>
35#include <casa/Arrays/Array.h>
36#include <casa/Arrays/Vector.h>
37#include <casa/Arrays/VectorIter.h>
38#include <casa/Utilities/CountedPtr.h>
39#include <casa/OS/Directory.h>
40
41#include <coordinates/Coordinates/CoordinateUtil.h>
42#include <coordinates/Coordinates/CoordinateSystem.h>
43#include <coordinates/Coordinates/SpectralCoordinate.h>
44#include <coordinates/Coordinates/DirectionCoordinate.h>
45#include <coordinates/Coordinates/StokesCoordinate.h>
46#include <coordinates/Coordinates/ObsInfo.h>
47
48#include <images/Images/ImageFITSConverter.h>
49#include <images/Images/TempImage.h>
50
51#include <lattices/Lattices/ArrayLattice.h>
52
53#include <measures/Measures/MEpoch.h>
54#include <measures/Measures/Stokes.h>
55
56#include <tables/Tables/Table.h>
57#include <tables/Tables/ScalarColumn.h>
58#include <tables/Tables/ArrayColumn.h>
59
60#include "SDDefs.h"
61#include "SDContainer.h"
62#include "SDMemTable.h"
63#include "SDPol.h"
64
65using namespace casa;
66using namespace asap;
67
68
69SDFITSImageWriter::SDFITSImageWriter()
70{;}
71
72SDFITSImageWriter::~SDFITSImageWriter()
73{;}
74
75
76Bool SDFITSImageWriter::write(const SDMemTable& sdTable,
77 const String& dirName, Bool toStokes,
78 Bool verbose)
79{
80
81// Get global Header from Table
82
83 SDHeader header = sdTable.getSDHeader();
84 MEpoch::Ref timeRef(sdTable.getTimeReference());
85 MDirection::Types dirRef = sdTable.getDirectionReference();
86
87// Prepare initial ObsInfo
88
89 ObsInfo oi;
90 oi.setTelescope(String(header.antennaname));
91 oi.setObserver(String(header.observer));
92
93// Column keywords
94
95 Table tab = sdTable.table();
96 ROArrayColumn<Double> dirCol(tab, String("DIRECTION"));
97 ROScalarColumn<Double> timeCol(tab, "TIME");
98 ROArrayColumn<uInt> freqidCol(tab, "FREQID");
99 ROArrayColumn<uInt> restfreqidCol(tab, "RESTFREQID");
100 ROScalarColumn<String> srcCol(tab, "SRCNAME");
101
102// Output Image Shape; spectral axis to be updated
103
104 IPosition shapeOut(4,1); // spectral, direction (2), stokes
105
106// Axes (should be in header)
107
108 const uInt beamAxis = asap::BeamAxis;
109 const uInt ifAxis = asap::IFAxis;
110 const uInt stokesAxis = asap::PolAxis;
111 const uInt chanAxis = asap::ChanAxis;
112 const Unit RAD(String("rad"));
113
114// Output Directory
115
116 String dirName2(dirName);
117 if (dirName.length()==0) {
118 dirName2 = String("asap_FITS");
119 }
120 Directory dir(dirName2);
121 dir.create(True);
122 cout << "Created directory '" << dirName2 << "' for output files" << endl;
123
124// Temps
125
126 Vector<Int> whichStokes(1,1);
127 Array<Double> whichDir;
128 Vector<Double> lonLat(2);
129 IPosition posDir(2,0);
130 Projection proj(Projection::SIN); // What should we use ?
131 Matrix<Double> xForm(2,2);
132 xForm = 0.0; xForm.diagonal() = 1.0;
133 Vector<Double> incLonLat(2,0.0);
134 Vector<Double> refPixLonLat(2,0.0);
135
136// Do we have linear or circular ? No way to know yet...
137
138 Bool linear = True;
139
140// Loop over rows
141
142 uInt maxMem = 128;
143 Bool preferVelocity = False;
144 Bool opticalVelocity = False;
145 Int bitPix = -32; // FLoating point
146 Float minPix = 1.0;
147 Float maxPix = -1.0;
148 Bool overWrite = True;
149 Bool degLast = False;
150 Bool reallyVerbose = False;
151//
152 String errMsg;
153 const uInt nRows = sdTable.nRow();
154 for (uInt iRow=0; iRow<nRows; iRow++) {
155
156// Get data converted to Stokes
157
158 const MaskedArray<Float>& dataIn(sdTable.rowAsMaskedArray(iRow,toStokes));
159 const Array<Float>& values = dataIn.getArray();
160 const Array<Bool>& mask = dataIn.getMask();
161 const IPosition& shapeIn = dataIn.shape();
162 shapeOut(0) = shapeIn(chanAxis);
163 TiledShape tShapeOut(shapeOut);
164
165// Update ObsInfo (time changes per integration)
166
167 Double dTmp;
168 timeCol.get(iRow, dTmp);
169 MVEpoch tmp2(Quantum<Double>(dTmp, Unit(String("d"))));
170 MEpoch epoch(tmp2, timeRef);
171 oi.setObsDate(epoch);
172
173// Iterate through data in this row by spectra
174
175 ReadOnlyVectorIterator<Float> itData(values, chanAxis);
176 ReadOnlyVectorIterator<Bool> itMask(mask, chanAxis);
177 while (!itData.pastEnd()) {
178 const IPosition& pos = itData.pos();
179
180// Form SpectralCoordinate
181
182 Vector<uInt> iTmp;
183 freqidCol.get(iRow, iTmp);
184 uInt freqID = iTmp(pos(ifAxis));
185//
186 restfreqidCol.get(iRow, iTmp);
187 uInt restFreqID = iTmp(pos(ifAxis));
188//
189 SpectralCoordinate sC = sdTable.getSpectralCoordinate(freqID, restFreqID, iRow);
190
191// Form DirectionCoordinate
192
193 dirCol.get(iRow, whichDir);
194 posDir(0) = pos(beamAxis);
195 posDir(1) = 0;
196 lonLat[0] = whichDir(posDir);
197//
198 posDir(0) = pos(beamAxis);
199 posDir(1) = 1;
200 lonLat[1] = whichDir(posDir);
201 DirectionCoordinate dC(dirRef, proj, lonLat[0], lonLat[1],
202 incLonLat[0], incLonLat[1], xForm,
203 refPixLonLat[0], refPixLonLat[1]);
204
205// Form Stokes Coordinate (Stokes info still sketchy);
206
207 Stokes::StokesTypes stokes = SDPolUtil::convertStokes(pos(stokesAxis), toStokes, linear);
208 String stokesName = Stokes::name(stokes);
209 whichStokes(0) = Int(stokes);
210 StokesCoordinate stC(whichStokes);
211
212// Create CoordinateSystem
213
214 CoordinateSystem cSys;
215 cSys.addCoordinate(sC);
216 cSys.addCoordinate(dC);
217 cSys.addCoordinate(stC);
218 cSys.setObsInfo(oi);
219
220// Reform data into correct shape for output
221
222 Array<Float> t1(itData.array().reform(shapeOut));
223 Array<Bool> m1(itMask.array().reform(shapeOut));
224
225// Create aips++ Image
226
227 TempImage<Float> tIm(tShapeOut, cSys);
228 tIm.put(t1);
229//
230 ArrayLattice<Bool> latMask(m1);
231 tIm.attachMask(latMask);
232
233// Write out as FITS Image file
234
235 ostringstream oss;
236 oss << "row" << iRow << "_beam" << pos(0) << "_if"
237 << pos(1) << "_" << stokesName << "_" << srcCol(iRow) << ".fits";
238 String tS(oss);
239 String fileName = dirName2 + String("/") + tS;
240 if (verbose) cout << "Writing row " << iRow
241 << " into file '" << tS << "'" << endl;
242//
243 Bool ok = ImageFITSConverter::ImageToFITS(errMsg, tIm, fileName,
244 maxMem, preferVelocity,
245 opticalVelocity, bitPix,
246 minPix, maxPix, overWrite,
247 degLast, reallyVerbose);
248 if (!ok) {
249 cerr << "Error writing fits - " << errMsg << endl;
250 }
251
252// Next spectrum
253
254 itData.next();
255 itMask.next();
256 }
257 }
258//
259 if (!verbose) {
260 cout << "Wrote " << nRows << " into individual FITS files" << endl;
261 }
262//
263 return True;
264}
265
266
Note: See TracBrowser for help on using the repository browser.