source: branches/alma/src/STFITSImageWriter.cpp @ 1604

Last change on this file since 1604 was 1604, checked in by TakTsutsumi, 15 years ago

New Development: No (asap2.3.1 upgrade)

JIRA Issue: No

Ready to Release: Yes

Interface Changes: Yes

What Interface Changed:

Test Programs: List test programs

Put in Release Notes: No

Module(s): single dish

Description: new files added from asap2.3.1


File size: 9.5 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  while ( !iter.pastEnd() ) {
124    Table t = iter.table();
125    ROTableRow row(t);
126    const TableRecord& rec = row.get(0);
127    String dirtype = stable.getDirectionRefString();
128    ostringstream onstr;
129    onstr << "SCAN" << rec.asuInt("SCANNO")
130    << "_CYCLE" << rec.asuInt("CYCLENO")
131    << "_BEAM" << rec.asuInt("BEAMNO")
132    << "_IF" << rec.asuInt("IFNO")
133    << "_POL" << rec.asuInt("POLNO");
134    String fileName = rootName + String(onstr) + String(".fits");
135    int row0 = t.rowNumbers(tab)[0];
136
137    const MPosition& mp = stable.getAntennaPosition();
138    const MDirection& md = stable.getDirection(row0);
139    const MEpoch& me = stable.getEpoch(row0);
140    oi.setObsDate(me);
141   
142    //const Double& rf = 
143    //  stable.molecules().getRestFrequency(rec.asuInt("MOLECULE_ID") );
144    const std::vector<double>& rf = 
145      stable.molecules().getRestFrequency(rec.asuInt("MOLECULE_ID") );
146    SpectralCoordinate sC =
147      stable.frequencies().getSpectralCoordinate(md, mp, me, rf,
148                                                 rec.asuInt("FREQ_ID"));
149
150    Int polno = rec.asuInt("POLNO");
151    Stokes::StokesTypes stokes =
152      Stokes::type(stable.getPolarizationLabel(polno, stable.getPolType()));
153    Vector<Int> whichStokes(1);
154    whichStokes(0) = Int(stokes);
155    StokesCoordinate stC(whichStokes);
156
157
158    CoordinateSystem cSys;
159    Vector<Double> lonlat = md.getAngle().getValue();
160    DirectionCoordinate dC =
161      getDirectionCoordinate(stable.getDirectionRefString(),
162                             lonlat[0], lonlat[1]);
163    cSys.addCoordinate(sC);
164    cSys.addCoordinate(dC);
165    cSys.addCoordinate(stC);
166    cSys.setObsInfo(oi);
167   
168    Vector<Float> spec = rec.asArrayFloat("SPECTRA");
169    Vector<uChar> flag = rec.asArrayuChar("FLAGTRA");
170    Vector<Bool> bflag(flag.shape());
171    convertArray(bflag, flag);
172
173    // Create casacore Image
174    IPosition shp(4, 1);
175    shp(0)= spec.nelements();
176    TempImage<Float> tIm(shp, cSys);
177    tIm.put(spec);
178    ArrayLattice<Bool> latMask(shp);
179    IPosition where(4,0);
180    IPosition stride(4,1);
181    latMask.putSlice(!bflag, where, stride);
182    tIm.attachMask(latMask);
183    if (isClass_) {
184      preferVelocity = True;
185    }
186    Bool ok = ImageFITSConverter::ImageToFITS(errMsg, tIm, fileName,
187                                              maxMem, preferVelocity,
188                                              opticalVelocity, bitPix,
189                                              minPix, maxPix, overWrite,
190                                              degLast, reallyVerbose);
191    if (!ok) {
192      throw(AipsError(errMsg));
193    }
194   
195    int status = 0;
196    fitsfile *fptr;     
197    if( fits_open_file(&fptr, fileName.c_str(), READWRITE, &status) )
198      throw AipsError("Coudn't open fits file for modification");
199
200    if (isClass_) {
201      // Add CLASS specific information
202      // Apply hacks to fits file so it can be read by class
203      // modifications are
204      // CTYPE1 :  FREQ-XYZ -> FREQ
205      // CRVAL1 :  -> CRVAL1-RESTFREQ
206     
207      // Use preferVelocity which seems to work.
208      // When no restfrequency is given this will set up the correct
209      // frequency axis. Otherwise it uses a velocity axis
210      char ctyp[84];
211      fits_read_key(fptr, TSTRING, "CTYPE1", ctyp, NULL, &status);
212      if (String(ctyp).contains("FREQ", 0)) {
213        float restf,refval;
214        fits_read_key(fptr, TFLOAT, "CRVAL1",
215                      &refval, NULL, &status);
216        restf = 0.0;
217        if ( fits_update_key(fptr, TFLOAT, "CRVAL1", &restf,
218                             NULL, &status) ) {
219          throw AipsError("Couldn't modify CRVAL1");       
220        }
221        if ( fits_update_key(fptr, TFLOAT, "RESTFREQ", &refval,
222                             NULL, &status) ) {
223          throw AipsError("Couldn't modify RESTFREQ");       
224        }
225      } else {
226        float refval = sC.referenceValue()[0];
227        if ( fits_update_key(fptr, TFLOAT, "RESTFREQ", &refval,
228                             NULL, &status) ) {
229          throw AipsError("Couldn't modify RESTFREQ");       
230        }
231       
232      }
233    }
234
235    Float tsys = stable.getTsys(row0);
236    if ( fits_update_key(fptr, TFLOAT, "TSYS", &tsys,
237                         "System temperature", &status) ) {
238      throw AipsError("Couldn't modify TSYS");       
239    }
240    Float otime = Float(stable.getIntTime(row0));
241    if ( fits_update_key(fptr, TFLOAT, "OBSTIME", &otime,
242                         "Integration time", &status) ) {
243      throw AipsError("Couldn't modify OBSTIME");       
244    }
245   
246    const char* oname = stable.getSourceName(row0).c_str();
247    if ( fits_update_key(fptr, TSTRING, "OBJECT", (char *)oname,
248                         NULL, &status) ) {
249      throw AipsError("Couldn't modify OBJECT");       
250    }
251    Float elev = stable.getElevation(row0)/C::pi*180.0;
252    if ( fits_update_key(fptr, TFLOAT, "ELEVATIO", &elev,
253                         "Telescope elevation (degrees)", &status) ) {
254      throw AipsError("Couldn't modify ELEVATIO");       
255    }
256    Float azi = stable.getAzimuth(row0)/C::pi*180.0;
257    if ( fits_update_key(fptr, TFLOAT, "AZIMUTH", &azi,
258                         "Telescope azimuth (degrees)", &status) ) {
259      throw AipsError("Couldn't modify AZIMUTH");       
260    }
261    fits_close_file(fptr, &status);
262 
263    //pushLog(String(oss));
264    ++iter;
265  }
266  return True;
267}
268
269DirectionCoordinate
270STFITSImageWriter::getDirectionCoordinate(const String& reff,
271                                          Double lon, Double lat)
272{
273   
274   Projection proj(Projection::SIN);                   // What should we use ?
275   Matrix<Double> xForm(2,2);
276   xForm = 0.0; xForm.diagonal() = 1.0;
277   Vector<Double> incLonLat(2,0.0);
278   Vector<Double> refPixLonLat(2,0.0);
279   MDirection::Types mdt;
280   if (!MDirection::getType(mdt, reff)) {
281     throw(AipsError("Illegal Direction frame."));
282   }
283   return DirectionCoordinate(mdt, proj, lon, lat,
284                              incLonLat[0], incLonLat[1], xForm,
285                              refPixLonLat[0], refPixLonLat[1]);
286   
287}
288
Note: See TracBrowser for help on using the repository browser.