source: trunk/python/scantable.py@ 120

Last change on this file since 120 was 116, checked in by mar637, 20 years ago

Added SDFITS writing.
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 13.8 KB
Line 
1from asap._asap import sdtable
2from numarray import ones
3import sys
4
5class scantable(sdtable):
6 """
7 The ASAP container for scans
8 """
9
10 def __init__(self, filename):
11 """
12 Create a scantable from a saved one or make a reference
13 Parameters:
14 filename: the name of an asap table on disk, or
15 [advanced] a refernce to an existing
16 scantable
17 """
18 self._vb = True
19 self._p = None
20 sdtable.__init__(self, filename)
21
22 def save(self, name, format='ASAP'):
23 """
24 Store the scantable on disk. This can be a asap file or SDFITS/MS2.
25 Parameters:
26 name: the name of the outputfile
27 format: an optional file format. Default is ASAP.
28 Alllowed are 'ASAP', 'SDFITS' and 'MS2'
29 Example:
30 scan.save('myscan.asap')
31 scan.save('myscan.sdfits','SDFITS')
32 """
33 if format == 'ASAP':
34 self._save(name)
35 else:
36 from asap._asap import sdwriter as _sw
37 w = _sw()
38 if format == 'SDFITS':
39 w.write(self, name)
40 return
41
42 def _verbose(self, *args):
43 """
44 Set the verbose level True or False, to indicate if output
45 should be printed as well as returned.
46 """
47 if type(args[0]) is bool:
48 self._vb = args[0]
49 return
50 elif len(args) == 0:
51 return self._vb
52
53
54 def copy(self):
55 """
56 Return a copy of this scantable.
57 Parameters:
58 none
59 Example:
60 copiedscan = scan.copy()
61 """
62 sd = scantable(sdtable._copy(self))
63 return sd
64
65 def get_scan(self, scanid=None):
66 """
67 Return a specific scan (by scanno) or collection of scans (by
68 source name) in a new scantable.
69 Parameters:
70 scanid: a scanno or a source name
71 Example:
72 scan.get_scan('323p459')
73 # gets all scans containing the source '323p459'
74 """
75 if scanid is None:
76 print "Please specify a scan no or name to retrieve from the scantable"
77 try:
78 if type(scanid) is str:
79 s = sdtable._getsource(self,scanid)
80 return scantable(s)
81 elif type(scanid) is int:
82 s = sdtable._getscan(self,scanid)
83 return scantable(s)
84 except RuntimeError:
85 print "Couldn't find any match."
86
87 def __str__(self):
88 return sdtable.summary(self)
89
90 def summary(self,filename=None):
91 """
92 Print a summary of the contents of this scantable.
93 Parameters:
94 filename: the name of a file to write the putput to
95 Default - no file output
96 """
97 info = sdtable.summary(self)
98 if filename is not None:
99 data = open(filename, 'w')
100 data.write(info)
101 data.close()
102 print info
103
104 def set_selection(self, thebeam=0,theif=0,thepol=0):
105 """
106 Set the spectrum for individual operations.
107 Parameters:
108 thebeam,theif,thepol: a number
109 Example:
110 scan.set_selection(0,0,1)
111 pol1rms = scan.rms(all=False) # returns rms for beam=0
112 # if=0, pol=1
113 """
114 self.setbeam(thebeam)
115 self.setpol(thepol)
116 self.setif(theif)
117 return
118
119 def get_selection(self):
120 """
121 Return/print a the current 'cursor' into the Beam/IF/Pol cube.
122 Parameters:
123 none
124 Returns:
125 a list of values (currentBeam,currentIF,currentPol)
126 Example:
127 none
128 """
129 i = self.getbeam()
130 j = self.getif()
131 k = self.getpol()
132 if self._vb:
133 out = 'Beam=%d IF=%d Pol=%d '% (i,j,k)
134 print out
135 return i,j,k
136
137 def rms(self,mask=None, all=True):
138 """
139 Determine the root mean square of the current beam/if/pol
140 Takes a 'mask' as an optional parameter to specify which
141 channels should be excluded.
142 Parameters:
143 mask: an optional mask specifying where the rms
144 should be determined.
145 all: optional flag to show all or a selected
146 spectrum of Beam/IF/Pol
147
148 Example:
149 scan.set_unit('channel')
150 msk = scan.create_mask([100,200],[500,600])
151 scan.rms(mask=m)
152 """
153 from asap._asap import rms as _rms
154 if mask == None:
155 mask = ones(self.nchan())
156 if all:
157 out = ''
158 tmp = []
159 for i in range(self.nbeam()):
160 self.setbeam(i)
161 for j in range(self.nif()):
162 self.setif(j)
163 for k in range(self.npol()):
164 self.setpol(k)
165 rmsval = _rms(self,mask)
166 tmp.append(rmsval)
167 out += 'Beam[%d], IF[%d], Pol[%d] = %3.3f\n' % (i,j,k,rmsval)
168 if self._vb:
169 print out
170 return tmp
171
172 else:
173 i = self.getbeam()
174 j = self.getif()
175 k = self.getpol()
176 rmsval = _rms(self,mask)
177 out = 'Beam[%d], IF[%d], Pol[%d] = %3.3f' % (i,j,k,rmsval)
178 if self._vb:
179 print out
180 return rmsval
181
182 def get_tsys(self, all=True):
183 """
184 Return the System temperatures.
185 Parameters:
186 all: optional parameter to get the Tsys values for all
187 Beams/IFs/Pols (default) or just the one selected
188 with scantable.set_selection()
189 [True or False]
190 Returns:
191 a list of Tsys values.
192 """
193 if all:
194 tmp = []
195 out = ''
196 for i in range(self.nbeam()):
197 self.setbeam(i)
198 for j in range(self.nif()):
199 self.setif(j)
200 for k in range(self.npol()):
201 self.setpol(k)
202 ts = self._gettsys()
203 tmp.append(ts)
204 out += 'TSys: Beam[%d], IF[%d], Pol[%d] = %3.3f\n' % (i,j,k,ts)
205 if self._vb:
206 print out
207 return tmp
208 else:
209 i = self.getbeam()
210 j = self.getif()
211 k = self.getpol()
212 ts = self._gettsys()
213 out = 'TSys: Beam[%d], IF[%d], Pol[%d] = %3.3f' % (i,j,k,ts)
214 if self._vb:
215 print out
216 return ts
217
218 def get_time(self):
219 """
220 Get a list of time stamps for the observations.
221 Return a string for each intergration in the scantable.
222 Parameters:
223 none
224 Example:
225 none
226 """
227 out = []
228 for i in range(self.nrow()):
229 out.append(self._gettime(i))
230 return out
231
232 def set_unit(self, unit='channel'):
233 """
234 Set the unit for all following operations on this scantable
235 Parameters:
236 unit: optional unit, default is 'channel'
237 one of '*Hz','km/s','channel', ''
238 """
239
240 if unit in ['','pixel', 'channel']:
241 unit = ''
242 inf = list(self._getcoordinfo())
243 inf[0] = unit
244 self._setcoordinfo(inf)
245 if self._p: self.plot()
246
247 def set_freqframe(self, frame='LSRK'):
248 """
249 Set the frame type of the Spectral Axis.
250 Parameters:
251 frame: an optional frame type, default 'LSRK'.
252 Examples:
253 scan.set_freqframe('BARY')
254 """
255 valid = ['REST','TOPO','LSRD','LSRK','BARY', \
256 'GEO','GALACTO','LGROUP','CMB']
257 if frame in valid:
258 inf = list(self._getcoordinfo())
259 inf[1] = frame
260 self._setcoordinfo(inf)
261 else:
262 print "Please specify a valid freq type. Valid types are:\n",valid
263
264 def get_unit(self):
265 """
266 Get the default unit set in this scantable
267 Parameters:
268 Returns:
269 A unit string
270 """
271 inf = self._getcoordinfo()
272 unit = inf[0]
273 if unit == '': unit = 'channel'
274 return unit
275
276 def get_abscissa(self, rowno=0):
277 """
278 Get the abscissa in the current coordinate setup for the currently
279 selected Beam/IF/Pol
280 Parameters:
281 none
282 Returns:
283 The abscissa values and it's format string.
284 """
285 absc = self.getabscissa(rowno)
286 lbl = self.getabscissalabel(rowno)
287 return absc, lbl
288
289 def create_mask(self, *args, **kwargs):
290 """
291 Compute and return a mask based on [min,max] windows.
292 The specified windows are to be EXCLUDED, when the mask is
293 applied.
294 Parameters:
295 [min,max],[min2,max2],...
296 Pairs of start/end points specifying the regions
297 to be masked
298 invert: return an inverted mask, i.e. the regions
299 specified are not masked (INCLUDED)
300 Example:
301 scan.set_unit('channel')
302
303 a)
304 msk = scan.set_mask([400,500],[800,900])
305 # masks the regions between 400 and 500
306 # and 800 and 900 in the unit 'channel'
307
308 b)
309 msk = scan.set_mask([400,500],[800,900], invert=True)
310 # masks the regions outside 400 and 500
311 # and 800 and 900 in the unit 'channel'
312
313 """
314 u = self._getcoordinfo()[0]
315 if self._vb:
316 if u == "": u = "channel"
317 print "The current mask window unit is", u
318 n = self.nchan()
319 data = self.getabscissa()
320 msk = ones(n)
321 for window in args:
322 if (len(window) != 2 or window[0] > window[1] ):
323 print "A window needs to be defined as [min,max]"
324 return
325 for i in range(n):
326 if data[i] >= window[0] and data[i] < window[1]:
327 msk[i] = 0
328 if kwargs.has_key('invert'):
329 if kwargs.get('invert'):
330 from numarray import logical_not
331 msk = logical_not(msk)
332 return msk
333
334 def set_restfreqs(self, freqs, unit='Hz'):
335 """
336 Set the restfrequency(s) for this scantable.
337 Parameters:
338 freqs: one or more frequencies
339 unit: optional 'unit', default 'Hz'
340 Example:
341 scan.set_restfreqs([1000000000.0])
342 """
343 if type(freqs) is float or int:
344 freqs = (freqs)
345 sdtable._setrestfreqs(self,freqs, unit)
346 return
347
348 def flag_spectrum(self, thebeam, theif, thepol):
349 """
350 This flags a selected spectrum in the scan 'for good'.
351 USE WITH CARE - not reversible.
352 Use masks for non-permanent exclusion of channels.
353 Parameters:
354 thebeam,theif,thepol: all have to be explicitly
355 specified
356 Example:
357 scan.flag_spectrum(0,0,1)
358 flags the spectrum for Beam=0, IF=0, Pol=1
359 """
360 if (thebeam < self.nbeam() and theif < self.nif() and thepol < self.npol()):
361 stable.setbeam(thebeam)
362 stable.setif(theif)
363 stable.setpol(thepol)
364 stable._flag(self)
365 else:
366 print "Please specify a valid (Beam/IF/Pol)"
367 return
368
369 def plot(self, what='spectrum',col='Pol', panel=None):
370 """
371 Plot the spectra contained in the scan. Alternatively you can also
372 Plot Tsys vs Time
373 Parameters:
374 what: a choice of 'spectrum' (default) or 'tsys'
375 col: which out of Beams/IFs/Pols should be colour stacked
376 panel: set up multiple panels, currently not working.
377 """
378 validcol = {'Beam':self.nbeam(),'IF':self.nif(),'Pol':self.npol()}
379 validyax = ['spectrum','tsys']
380 if not self._p:
381 from asap.asaplot import ASAPlot
382 self._p = ASAPlot()
383# print "Plotting not enabled"
384# return
385 npan = 1
386 x = None
387 if what == 'tsys':
388 n = self.nrow()
389 if n < 2:
390 print "Only one integration. Can't plot."
391 return
392 if panel == 'Time':
393 npan = self.nrow()
394 self._p.hold()
395 self._p.clear()
396 xlab,ylab,tlab = None,None,None
397 vb = self._verbose
398 self._verbose(False)
399 sel = self.get_selection()
400 for i in range(npan):
401 for j in range(validcol[col]):
402 x = None
403 y = None
404 m = None
405 tlab = self._getsourcename(i)
406 if col == 'Beam':
407 self.setbeam(j)
408 elif col == 'IF':
409 self.setif(j)
410 elif col == 'Pol':
411 self.setpol(j)
412 if what == 'tsys':
413 x = range(self.nrow())
414 xlab = 'Time [pixel]'
415 m = list(ones(len(x)))
416 y = []
417 ylab = r'$T_{sys}$'
418 for k in range(len(x)):
419 y.append(self._gettsys(k))
420 else:
421 x,xlab = self.get_abscissa(i)
422 y = self.getspectrum(i)
423 ylab = r'Flux'
424 m = self.getmask(i)
425 llab = col+' '+str(j)
426 self._p.set_line(label=llab)
427 self._p.plot(x,y,m)
428 self._p.set_axes('xlabel',xlab)
429 self._p.set_axes('ylabel',ylab)
430 self._p.set_axes('title',tlab)
431 self._p.release()
432 self.set_selection(sel[0],sel[1],sel[2])
433 self._verbose(vb)
434 return
Note: See TracBrowser for help on using the repository browser.