source: branches/Release2.1.2/src/SDFITSImageWriter.cc@ 1757

Last change on this file since 1757 was 901, checked in by mar637, 20 years ago

SDContainer -> STHeader
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 8.0 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 901 2006-03-10 02:29:22Z 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#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 "STDefs.h"
61#include "STHeader.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 STHeader header = sdTable.getSTHeader();
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 ostringstream oss;
123 oss << "Created directory '" << dirName2 << "' for output files";
124
125 // Temps
126
127 Vector<Int> whichStokes(1,1);
128 Array<Double> whichDir;
129 Vector<Double> lonLat(2);
130 IPosition posDir(2,0);
131 Projection proj(Projection::SIN); // What should we use ?
132 Matrix<Double> xForm(2,2);
133 xForm = 0.0; xForm.diagonal() = 1.0;
134 Vector<Double> incLonLat(2,0.0);
135 Vector<Double> refPixLonLat(2,0.0);
136
137// Do we have linear or circular ? No way to know yet...
138
139 Bool linear = True;
140
141// Loop over rows
142
143 uInt maxMem = 128;
144 Bool preferVelocity = False;
145 Bool opticalVelocity = False;
146 Int bitPix = -32; // FLoating point
147 Float minPix = 1.0;
148 Float maxPix = -1.0;
149 Bool overWrite = True;
150 Bool degLast = False;
151 Bool reallyVerbose = False;
152//
153 String errMsg;
154 const uInt nRows = sdTable.nRow();
155 for (uInt iRow=0; iRow<nRows; iRow++) {
156
157// Get data converted to Stokes
158
159 const MaskedArray<Float>& dataIn(sdTable.rowAsMaskedArray(iRow,toStokes));
160 const Array<Float>& values = dataIn.getArray();
161 const Array<Bool>& mask = dataIn.getMask();
162 const IPosition& shapeIn = dataIn.shape();
163 shapeOut(0) = shapeIn(chanAxis);
164 TiledShape tShapeOut(shapeOut);
165
166// Update ObsInfo (time changes per integration)
167
168 Double dTmp;
169 timeCol.get(iRow, dTmp);
170 MVEpoch tmp2(Quantum<Double>(dTmp, Unit(String("d"))));
171 MEpoch epoch(tmp2, timeRef);
172 oi.setObsDate(epoch);
173
174// Iterate through data in this row by spectra
175
176 ReadOnlyVectorIterator<Float> itData(values, chanAxis);
177 ReadOnlyVectorIterator<Bool> itMask(mask, chanAxis);
178 while (!itData.pastEnd()) {
179 const IPosition& pos = itData.pos();
180
181// Form SpectralCoordinate
182
183 Vector<uInt> iTmp;
184 freqidCol.get(iRow, iTmp);
185 uInt freqID = iTmp(pos(ifAxis));
186//
187 restfreqidCol.get(iRow, iTmp);
188 uInt restFreqID = iTmp(pos(ifAxis));
189//
190 SpectralCoordinate sC = sdTable.getSpectralCoordinate(freqID, restFreqID, iRow);
191
192// Form DirectionCoordinate
193
194 dirCol.get(iRow, whichDir);
195 posDir(0) = pos(beamAxis);
196 posDir(1) = 0;
197 lonLat[0] = whichDir(posDir);
198//
199 posDir(0) = pos(beamAxis);
200 posDir(1) = 1;
201 lonLat[1] = whichDir(posDir);
202 DirectionCoordinate dC(dirRef, proj, lonLat[0], lonLat[1],
203 incLonLat[0], incLonLat[1], xForm,
204 refPixLonLat[0], refPixLonLat[1]);
205
206// Form Stokes Coordinate (Stokes info still sketchy);
207
208 Stokes::StokesTypes stokes = SDPolUtil::convertStokes(pos(stokesAxis), toStokes, linear);
209 String stokesName = Stokes::name(stokes);
210 whichStokes(0) = Int(stokes);
211 StokesCoordinate stC(whichStokes);
212
213// Create CoordinateSystem
214
215 CoordinateSystem cSys;
216 cSys.addCoordinate(sC);
217 cSys.addCoordinate(dC);
218 cSys.addCoordinate(stC);
219 cSys.setObsInfo(oi);
220
221// Reform data into correct shape for output
222
223 Array<Float> t1(itData.array().reform(shapeOut));
224 Array<Bool> m1(itMask.array().reform(shapeOut));
225
226// Create aips++ Image
227
228 TempImage<Float> tIm(tShapeOut, cSys);
229 tIm.put(t1);
230//
231 ArrayLattice<Bool> latMask(m1);
232 tIm.attachMask(latMask);
233
234// Write out as FITS Image file
235
236 ostringstream oss2;
237 oss2 << "row" << iRow << "_beam" << pos(0) << "_if"
238 << pos(1) << "_" << stokesName << "_" << srcCol(iRow) << ".fits";
239 String tS(oss2);
240 String fileName = dirName2 + String("/") + tS;
241 if (verbose) oss << endl << "Writing row " << iRow
242 << " into file '" << tS << "'";
243//
244 Bool ok = ImageFITSConverter::ImageToFITS(errMsg, tIm, fileName,
245 maxMem, preferVelocity,
246 opticalVelocity, bitPix,
247 minPix, maxPix, overWrite,
248 degLast, reallyVerbose);
249 if (!ok) {
250 oss << endl<< "Error writing fits - " << errMsg;
251 }
252
253// Next spectrum
254
255 itData.next();
256 itMask.next();
257 }
258 }
259//
260 if (!verbose) {
261 oss << endl<< "Wrote " << nRows << " into individual FITS files";
262 }
263 pushLog(String(oss));
264 return True;
265}
266
Note: See TracBrowser for help on using the repository browser.