source: trunk/src/STFITSImageWriter.cpp@ 1443

Last change on this file since 1443 was 1443, checked in by Malte Marquarding, 17 years ago

Fix ticket #127; still have to add class header hack
File size: 6.8 KB
Line 
1//#---------------------------------------------------------------------------
2//# STFITSImageWriter.cc: Class to write out single dish spectra as FITS images
3//#---------------------------------------------------------------------------
4//# Copyright (C) 2008
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: $
30//#---------------------------------------------------------------------------
31
32#include <images/Images/TempImage.h>
33
34#include <lattices/Lattices/ArrayLattice.h>
35
36#include <measures/Measures/MEpoch.h>
37#include <measures/Measures/Stokes.h>
38
39#include <tables/Tables/Table.h>
40#include <tables/Tables/TableIter.h>
41#include <tables/Tables/TableRecord.h>
42#include <casa/Containers/RecordField.h>
43#include <tables/Tables/TableRow.h>
44#include <tables/Tables/ScalarColumn.h>
45#include <tables/Tables/ArrayColumn.h>
46
47
48#include <coordinates/Coordinates/CoordinateUtil.h>
49#include <coordinates/Coordinates/CoordinateSystem.h>
50#include <coordinates/Coordinates/SpectralCoordinate.h>
51#include <coordinates/Coordinates/DirectionCoordinate.h>
52#include <coordinates/Coordinates/StokesCoordinate.h>
53#include <coordinates/Coordinates/ObsInfo.h>
54
55#include <images/Images/ImageFITSConverter.h>
56#include <images/Images/TempImage.h>
57
58#include <lattices/Lattices/ArrayLattice.h>
59
60
61//#include "STDefs.h"
62#include "STHeader.h"
63#include "Scantable.h"
64#include "STFITSImageWriter.h"
65
66using namespace casa;
67using namespace asap;
68
69
70STFITSImageWriter::STFITSImageWriter() : isClass_(False)
71{;}
72
73STFITSImageWriter::~STFITSImageWriter()
74{;}
75
76
77Bool STFITSImageWriter::write(const Scantable& stable,
78 const String& fileName)
79{
80
81// Get global Header from Table
82
83 STHeader hdr = stable.getHeader();
84
85// Column keywords
86
87 Table tab = stable.table();
88
89// Temps
90
91 const Unit RAD(String("rad"));
92
93// Open and write header file
94
95 String rootName(fileName);
96 if (rootName.length()==0) rootName = String("fits");
97
98 ObsInfo oi;
99 oi.setTelescope(String(hdr.antennaname));
100 oi.setObserver(String(hdr.observer));
101
102 uInt maxMem = 128;
103 Bool preferVelocity = False;
104 Bool opticalVelocity = False;
105 Int bitPix = -32; // FLoating point
106 Float minPix = 1.0;
107 Float maxPix = -1.0;
108 Bool overWrite = True;
109 Bool degLast = False;
110 Bool reallyVerbose = False;
111 String errMsg;
112
113
114 Block<String> cols(5);
115 cols[0] = String("SCANNO");
116 cols[1] = String("CYCLENO");
117 cols[2] = String("BEAMNO");
118 cols[3] = String("IFNO");
119 cols[4] = String("POLNO");
120 TableIterator iter(tab, cols);
121 // Open data file
122 while ( !iter.pastEnd() ) {
123 Table t = iter.table();
124 ROTableRow row(t);
125 const TableRecord& rec = row.get(0);
126 String dirtype = stable.getDirectionRefString();
127 ostringstream onstr;
128 onstr << "SCAN" << rec.asuInt("SCANNO")
129 << "_CYCLE" << rec.asuInt("CYCLENO")
130 << "_BEAM" << rec.asuInt("BEAMNO")
131 << "_IF" << rec.asuInt("IFNO")
132 << "_POL" << rec.asuInt("POLNO");
133 String fileName = rootName + String(onstr) + String(".fits");
134 int row0 = t.rowNumbers(tab)[0];
135
136 const MPosition& mp = stable.getAntennaPosition();
137 const MDirection& md = stable.getDirection(row0);
138 const MEpoch& me = stable.getEpoch(row0);
139 oi.setObsDate(me);
140
141 const Double& rf =
142 stable.molecules().getRestFrequency(rec.asuInt("MOLECULE_ID") );
143 SpectralCoordinate sC =
144 stable.frequencies().getSpectralCoordinate(md, mp, me, rf,
145 rec.asuInt("FREQ_ID"));
146
147 Int polno = rec.asuInt("POLNO");
148 Stokes::StokesTypes stokes =
149 Stokes::type(stable.getPolarizationLabel(polno, stable.getPolType()));
150 Vector<Int> whichStokes(1);
151 whichStokes(0) = Int(stokes);
152 StokesCoordinate stC(whichStokes);
153
154
155 CoordinateSystem cSys;
156 Vector<Double> lonlat = md.getAngle().getValue();
157 DirectionCoordinate dC =
158 getDirectionCoordinate(stable.getDirectionRefString(),
159 lonlat[0], lonlat[1]);
160 cSys.addCoordinate(sC);
161 cSys.addCoordinate(dC);
162 cSys.addCoordinate(stC);
163 cSys.setObsInfo(oi);
164
165 Vector<Float> spec = rec.asArrayFloat("SPECTRA");
166 Vector<uChar> flag = rec.asArrayuChar("FLAGTRA");
167 Vector<Bool> bflag(flag.shape());
168 convertArray(bflag, flag);
169
170 // Create casacore Image
171 IPosition shp(4, 1);
172 shp(0)= spec.nelements();
173 cout << shp << spec << endl;
174 TempImage<Float> tIm(shp, cSys);
175 tIm.put(spec);
176 ArrayLattice<Bool> latMask(shp);
177 IPosition where(4,0);
178 IPosition stride(4,1);
179 latMask.putSlice(!bflag, where, stride);
180 tIm.attachMask(latMask);
181 Bool ok = ImageFITSConverter::ImageToFITS(errMsg, tIm, fileName,
182 maxMem, preferVelocity,
183 opticalVelocity, bitPix,
184 minPix, maxPix, overWrite,
185 degLast, reallyVerbose);
186 if (!ok) {
187 throw(AipsError(errMsg));
188 }
189 if (isClass_) {
190 classHackHeader(fileName);
191 }
192 //pushLog(String(oss));
193 ++iter;
194 }
195 return True;
196}
197
198DirectionCoordinate
199STFITSImageWriter::getDirectionCoordinate(const String& reff,
200 Double lon, Double lat)
201{
202
203 Projection proj(Projection::SIN); // What should we use ?
204 Matrix<Double> xForm(2,2);
205 xForm = 0.0; xForm.diagonal() = 1.0;
206 Vector<Double> incLonLat(2,0.0);
207 Vector<Double> refPixLonLat(2,0.0);
208 MDirection::Types mdt;
209 if (!MDirection::getType(mdt, reff)) {
210 throw(AipsError("Illegal Direction frame."));
211 }
212 return DirectionCoordinate(mdt, proj, lon, lat,
213 incLonLat[0], incLonLat[1], xForm,
214 refPixLonLat[0], refPixLonLat[1]);
215
216}
217
218void STFITSImageWriter::classHackHeader(const String& filename) {
219 cout << "Not yet implemented" << endl;
220}
Note: See TracBrowser for help on using the repository browser.