source: trunk/src/SDFITSImageWriter.cc@ 212

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