source: trunk/src/STFITSImageWriter.cpp@ 2859

Last change on this file since 2859 was 2658, checked in by Malte Marquarding, 12 years ago

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