source: branches/newfiller/src/STFITSImageWriter.cpp@ 2380

Last change on this file since 2380 was 1757, checked in by Kana Sugimoto, 14 years ago

New Development: Yes

JIRA Issue: Yes (CAS-2211)

Ready for Test: Yes

Interface Changes: Yes

What Interface Changed: ASAP 3.0.0 interface changes

Test Programs:

Put in Release Notes: Yes

Module(s): all the CASA sd tools and tasks are affected.

Description: Merged ATNF-ASAP 3.0.0 developments to CASA (alma) branch.

Note you also need to update casa/code/atnf.
File size: 9.9 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 <fitsio.h>
33#include <images/Images/TempImage.h>
34
35#include <lattices/Lattices/ArrayLattice.h>
36
37#include <measures/Measures/MEpoch.h>
38#include <measures/Measures/Stokes.h>
39
40#include <tables/Tables/Table.h>
41#include <tables/Tables/TableIter.h>
42#include <tables/Tables/TableRecord.h>
43#include <casa/Containers/RecordField.h>
44#include <tables/Tables/TableRow.h>
45#include <tables/Tables/ScalarColumn.h>
46#include <tables/Tables/ArrayColumn.h>
47
48
49#include <coordinates/Coordinates/CoordinateUtil.h>
50#include <coordinates/Coordinates/CoordinateSystem.h>
51#include <coordinates/Coordinates/SpectralCoordinate.h>
52#include <coordinates/Coordinates/DirectionCoordinate.h>
53#include <coordinates/Coordinates/StokesCoordinate.h>
54#include <coordinates/Coordinates/ObsInfo.h>
55
56#include <images/Images/ImageFITSConverter.h>
57#include <images/Images/TempImage.h>
58
59#include <lattices/Lattices/ArrayLattice.h>
60
61
62//#include "STDefs.h"
63#include "STHeader.h"
64#include "Scantable.h"
65#include "STFITSImageWriter.h"
66
67using namespace casa;
68using namespace asap;
69
70
71STFITSImageWriter::STFITSImageWriter() : isClass_(False)
72{;}
73
74STFITSImageWriter::~STFITSImageWriter()
75{;}
76
77
78Bool STFITSImageWriter::write(const Scantable& stable,
79 const String& fileName)
80{
81
82// Get global Header from Table
83
84 STHeader hdr = stable.getHeader();
85
86// Column keywords
87
88 Table tab = stable.table();
89
90// Temps
91
92 const Unit RAD(String("rad"));
93
94// Open and write header file
95
96 String rootName(fileName);
97 if (rootName.length()==0) rootName = String("fits");
98
99 ObsInfo oi;
100 oi.setTelescope(String(hdr.antennaname));
101 oi.setObserver(String(hdr.observer));
102
103 uInt maxMem = 128;
104 Bool preferVelocity = False;
105 Bool opticalVelocity = False;
106 Int bitPix = -32; // FLoating point
107 Float minPix = 1.0;
108 Float maxPix = -1.0;
109 Bool overWrite = True;
110 Bool degLast = False;
111 Bool reallyVerbose = False;
112 String errMsg;
113
114
115 Block<String> cols(5);
116 cols[0] = String("SCANNO");
117 cols[1] = String("CYCLENO");
118 cols[2] = String("BEAMNO");
119 cols[3] = String("IFNO");
120 cols[4] = String("POLNO");
121 TableIterator iter(tab, cols);
122 // Open data file
123
124 while ( !iter.pastEnd() ) {
125 Table t = iter.table();
126 ROTableRow row(t);
127 const TableRecord& rec = row.get(0);
128 String dirtype = stable.getDirectionRefString();
129 ostringstream onstr;
130 if (rootName.length()==0) {
131 rootName = "fits";
132 }
133 if (tab.nrow() > 1) {
134 if (stable.nscan() > 1)
135 onstr << "_SCAN" << rec.asuInt("SCANNO");
136 if (stable.ncycle(rec.asuInt("SCANNO")) > 1)
137 onstr << "_CYCLE" << rec.asuInt("CYCLENO");
138 if (stable.nbeam(rec.asuInt("SCANNO")) > 1)
139 onstr << "_BEAM" << rec.asuInt("BEAMNO");
140 if (stable.nif(rec.asuInt("SCANNO")) > 1)
141 onstr << "_IF" << rec.asuInt("IFNO");
142 if (stable.npol(rec.asuInt("SCANNO")) > 1)
143 onstr << "_POL" << rec.asuInt("POLNO");
144 }
145 String fileName = rootName + String(onstr) + String(".fits");
146 int row0 = t.rowNumbers(tab)[0];
147
148 const MPosition& mp = stable.getAntennaPosition();
149 const MDirection& md = stable.getDirection(row0);
150 const MEpoch& me = stable.getEpoch(row0);
151 oi.setObsDate(me);
152
153 //const Double& rf =
154 // stable.molecules().getRestFrequency(rec.asuInt("MOLECULE_ID") );
155 const std::vector<double>& rf =
156 stable.molecules().getRestFrequency(rec.asuInt("MOLECULE_ID") );
157 SpectralCoordinate sC =
158 stable.frequencies().getSpectralCoordinate(md, mp, me, rf,
159 rec.asuInt("FREQ_ID"));
160
161 Int polno = rec.asuInt("POLNO");
162 Stokes::StokesTypes stokes =
163 Stokes::type(stable.getPolarizationLabel(polno, stable.getPolType()));
164 Vector<Int> whichStokes(1);
165 whichStokes(0) = Int(stokes);
166 StokesCoordinate stC(whichStokes);
167
168
169 CoordinateSystem cSys;
170 Vector<Double> lonlat = md.getAngle().getValue();
171 DirectionCoordinate dC =
172 getDirectionCoordinate(stable.getDirectionRefString(),
173 lonlat[0], lonlat[1]);
174 cSys.addCoordinate(sC);
175 cSys.addCoordinate(dC);
176 cSys.addCoordinate(stC);
177 cSys.setObsInfo(oi);
178
179 Vector<Float> spec = rec.asArrayFloat("SPECTRA");
180 Vector<uChar> flag = rec.asArrayuChar("FLAGTRA");
181 Vector<Bool> bflag(flag.shape());
182 convertArray(bflag, flag);
183
184 // Create casacore Image
185 IPosition shp(4, 1);
186 shp(0)= spec.nelements();
187 TempImage<Float> tIm(shp, cSys);
188 tIm.put(spec);
189 ArrayLattice<Bool> latMask(shp);
190 IPosition where(4,0);
191 IPosition stride(4,1);
192 latMask.putSlice(!bflag, where, stride);
193 tIm.attachMask(latMask);
194 if (isClass_) {
195 preferVelocity = True;
196 }
197 Bool ok = ImageFITSConverter::ImageToFITS(errMsg, tIm, fileName,
198 maxMem, preferVelocity,
199 opticalVelocity, bitPix,
200 minPix, maxPix, overWrite,
201 degLast, reallyVerbose);
202 if (!ok) {
203 throw(AipsError(errMsg));
204 }
205
206 int status = 0;
207 fitsfile *fptr;
208 if( fits_open_file(&fptr, fileName.c_str(), READWRITE, &status) )
209 throw AipsError("Coudn't open fits file for modification");
210
211 if (isClass_) {
212 // Add CLASS specific information
213 // Apply hacks to fits file so it can be read by class
214 // modifications are
215 // CTYPE1 : FREQ-XYZ -> FREQ
216 // CRVAL1 : -> CRVAL1-RESTFREQ
217
218 // Use preferVelocity which seems to work.
219 // When no restfrequency is given this will set up the correct
220 // frequency axis. Otherwise it uses a velocity axis
221 char ctyp[84];
222 fits_read_key(fptr, TSTRING, "CTYPE1", ctyp, NULL, &status);
223 if (String(ctyp).contains("FREQ", 0)) {
224 float restf,refval;
225 fits_read_key(fptr, TFLOAT, "CRVAL1",
226 &refval, NULL, &status);
227 restf = 0.0;
228 if ( fits_update_key(fptr, TFLOAT, "CRVAL1", &restf,
229 NULL, &status) ) {
230 throw AipsError("Couldn't modify CRVAL1");
231 }
232 if ( fits_update_key(fptr, TFLOAT, "RESTFREQ", &refval,
233 NULL, &status) ) {
234 throw AipsError("Couldn't modify RESTFREQ");
235 }
236 } else {
237 float refval = sC.referenceValue()[0];
238 if ( fits_update_key(fptr, TFLOAT, "RESTFREQ", &refval,
239 NULL, &status) ) {
240 throw AipsError("Couldn't modify RESTFREQ");
241 }
242
243 }
244 }
245
246 Float tsys = stable.getTsys(row0);
247 if ( fits_update_key(fptr, TFLOAT, "TSYS", &tsys,
248 "System temperature", &status) ) {
249 throw AipsError("Couldn't modify TSYS");
250 }
251 Float otime = Float(stable.getIntTime(row0));
252 if ( fits_update_key(fptr, TFLOAT, "OBSTIME", &otime,
253 "Integration time", &status) ) {
254 throw AipsError("Couldn't modify OBSTIME");
255 }
256
257 const char* oname = stable.getSourceName(row0).c_str();
258 if ( fits_update_key(fptr, TSTRING, "OBJECT", (char *)oname,
259 NULL, &status) ) {
260 throw AipsError("Couldn't modify OBJECT");
261 }
262 Float elev = stable.getElevation(row0)/C::pi*180.0;
263 if ( fits_update_key(fptr, TFLOAT, "ELEVATIO", &elev,
264 "Telescope elevation (degrees)", &status) ) {
265 throw AipsError("Couldn't modify ELEVATIO");
266 }
267 Float azi = stable.getAzimuth(row0)/C::pi*180.0;
268 if ( fits_update_key(fptr, TFLOAT, "AZIMUTH", &azi,
269 "Telescope azimuth (degrees)", &status) ) {
270 throw AipsError("Couldn't modify AZIMUTH");
271 }
272 fits_close_file(fptr, &status);
273 ostringstream oss;
274 oss << "Wrote " << fileName;
275 pushLog(String(oss));
276 //pushLog(String(oss));
277 ++iter;
278 }
279 return True;
280}
281
282DirectionCoordinate
283STFITSImageWriter::getDirectionCoordinate(const String& reff,
284 Double lon, Double lat)
285{
286
287 Projection proj(Projection::SIN); // What should we use ?
288 Matrix<Double> xForm(2,2);
289 xForm = 0.0; xForm.diagonal() = 1.0;
290 Vector<Double> incLonLat(2,0.0);
291 Vector<Double> refPixLonLat(2,0.0);
292 MDirection::Types mdt;
293 if (!MDirection::getType(mdt, reff)) {
294 throw(AipsError("Illegal Direction frame."));
295 }
296 return DirectionCoordinate(mdt, proj, lon, lat,
297 incLonLat[0], incLonLat[1], xForm,
298 refPixLonLat[0], refPixLonLat[1]);
299
300}
301
Note: See TracBrowser for help on using the repository browser.