source: trunk/python/scantable.py@ 522

Last change on this file since 522 was 516, checked in by vor010, 20 years ago

LineFinder/automatic baseline fitter: a bug related to multiple-row scantable handling has been corrected.
Help is changed to describe a new interface. Parameter checking + vector to tuple conversion for the edge parameter have been added

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 55.0 KB
Line 
1from asap._asap import sdtable
2from asap import rcParams
3from numarray import ones,zeros
4import sys
5
6class scantable(sdtable):
7 """
8 The ASAP container for scans
9 """
10
11 def __init__(self, filename, average=None, unit=None):
12 """
13 Create a scantable from a saved one or make a reference
14 Parameters:
15 filename: the name of an asap table on disk
16 or
17 the name of a rpfits/sdfits/ms file
18 (integrations within scans are auto averaged
19 and the whole file is read)
20 or
21 [advanced] a reference to an existing
22 scantable
23 average: average all integrations withinb a scan on read.
24 The default (True) is taken from .asaprc.
25 unit: brightness unit; must be consistent with K or Jy.
26 Over-rides the default selected by the reader
27 (input rpfits/sdfits/ms) or replaces the value
28 in existing scantables
29 """
30 if average is None or type(average) is not bool:
31 average = rcParams['scantable.autoaverage']
32
33 varlist = vars()
34 self._vb = rcParams['verbose']
35 self._p = None
36
37 if isinstance(filename,sdtable):
38 sdtable.__init__(self, filename)
39 if unit is not None:
40 self.set_fluxunit(unit)
41 else:
42 import os.path
43 if not os.path.exists(filename):
44 print "File '%s' not found." % (filename)
45 return
46 filename = os.path.expandvars(filename)
47 if os.path.isdir(filename):
48 # crude check if asap table
49 if os.path.exists(filename+'/table.info'):
50 sdtable.__init__(self, filename)
51 if unit is not None:
52 self.set_fluxunit(unit)
53 else:
54 print "The given file '%s'is not a valid asap table." % (filename)
55 return
56 else:
57 from asap._asap import sdreader
58 ifSel = -1
59 beamSel = -1
60 r = sdreader(filename,ifSel,beamSel)
61 print 'Importing data...'
62 r._read([-1])
63 tbl = r._getdata()
64 if unit is not None:
65 tbl.set_fluxunit(unit)
66 if average:
67 from asap._asap import average as _av
68 print 'Auto averaging integrations...'
69 tbl2 = _av((tbl,),(),True,'none')
70 sdtable.__init__(self,tbl2)
71 del tbl2
72 else:
73 sdtable.__init__(self,tbl)
74 del r,tbl
75 self._add_history("scantable", varlist)
76
77 def save(self, name=None, format=None, stokes=False, overwrite=False):
78 """
79 Store the scantable on disk. This can be an asap (aips++) Table, SDFITS,
80 Image FITS or MS2 format.
81 Parameters:
82 name: the name of the outputfile. For format="FITS" this
83 is the directory file name into which all the files
84 will be written (default is 'asap_FITS'). For format
85 "ASCII" this is the root file name (data in 'name'.txt
86 and header in 'name'_header.txt)
87 format: an optional file format. Default is ASAP.
88 Allowed are - 'ASAP' (save as ASAP [aips++] Table),
89 'SDFITS' (save as SDFITS file)
90 'FITS' (saves each row as a FITS Image)
91 'ASCII' (saves as ascii text file)
92 'MS2' (saves as an aips++
93 MeasurementSet V2)
94 stokes: Convert to Stokes parameters (only available
95 currently with FITS and ASCII formats.
96 Default is False.
97 overwrite: If the file should be overwritten if it exists.
98 The default False is to return with warning
99 without writing the output. USE WITH CARE.
100 Example:
101 scan.save('myscan.asap')
102 scan.save('myscan.sdfits','SDFITS')
103 """
104 from os import path
105 if format is None: format = rcParams['scantable.save']
106 suffix = '.'+format.lower()
107 if name is None or name =="":
108 name = 'scantable'+suffix
109 print "No filename given. Using default name %s..." % name
110 name = path.expandvars(name)
111 if path.isfile(name) or path.isdir(name):
112 if not overwrite:
113 print "File %s already exists." % name
114 return
115 format2 = format.upper()
116 if format2 == 'ASAP':
117 self._save(name)
118 else:
119 from asap._asap import sdwriter as _sw
120 w = _sw(format2)
121 w.write(self, name, stokes)
122 return
123
124 def copy(self):
125 """
126 Return a copy of this scantable.
127 Parameters:
128 none
129 Example:
130 copiedscan = scan.copy()
131 """
132 sd = scantable(sdtable._copy(self))
133 return sd
134
135 def get_scan(self, scanid=None):
136 """
137 Return a specific scan (by scanno) or collection of scans (by
138 source name) in a new scantable.
139 Parameters:
140 scanid: a (list of) scanno or a source name, unix-style
141 patterns are accepted for source name matching, e.g.
142 '*_R' gets all 'ref scans
143 Example:
144 # get all scans containing the source '323p459'
145 newscan = scan.get_scan('323p459')
146 # get all 'off' scans
147 refscans = scan.get_scan('*_R')
148 # get a susbset of scans by scanno (as listed in scan.summary())
149 newscan = scan.get_scan([0,2,7,10])
150 """
151 if scanid is None:
152 print "Please specify a scan no or name to retrieve from the scantable"
153 try:
154 if type(scanid) is str:
155 s = sdtable._getsource(self,scanid)
156 return scantable(s)
157 elif type(scanid) is int:
158 s = sdtable._getscan(self,[scanid])
159 return scantable(s)
160 elif type(scanid) is list:
161 s = sdtable._getscan(self,scanid)
162 return scantable(s)
163 else:
164 print "Illegal scanid type, use 'int' or 'list' if ints."
165 except RuntimeError:
166 print "Couldn't find any match."
167
168 def __str__(self):
169 return sdtable._summary(self,True)
170
171 def summary(self,filename=None, verbose=None):
172 """
173 Print a summary of the contents of this scantable.
174 Parameters:
175 filename: the name of a file to write the putput to
176 Default - no file output
177 verbose: print extra info such as the frequency table
178 The default (False) is taken from .asaprc
179 """
180 info = sdtable._summary(self, verbose)
181 if verbose is None: verbose = rcParams['scantable.verbosesummary']
182 if filename is not None:
183 if filename is "":
184 filename = 'scantable_summary.txt'
185 from os.path import expandvars, isdir
186 filename = expandvars(filename)
187 if not isdir(filename):
188 data = open(filename, 'w')
189 data.write(info)
190 data.close()
191 else:
192 print "Illegal file name '%s'." % (filename)
193 print info
194
195 def set_cursor(self, thebeam=0,theif=0,thepol=0):
196 """
197 Set the spectrum for individual operations.
198 Parameters:
199 thebeam,theif,thepol: a number
200 Example:
201 scan.set_cursor(0,0,1)
202 pol1sig = scan.stats(all=False) # returns std dev for beam=0
203 # if=0, pol=1
204 """
205 varlist = vars()
206 self.setbeam(thebeam)
207 self.setpol(thepol)
208 self.setif(theif)
209 self._add_history("set_cursor",varlist)
210 return
211
212 def get_cursor(self):
213 """
214 Return/print a the current 'cursor' into the Beam/IF/Pol cube.
215 Parameters:
216 none
217 Returns:
218 a list of values (currentBeam,currentIF,currentPol)
219 Example:
220 none
221 """
222 i = self.getbeam()
223 j = self.getif()
224 k = self.getpol()
225 if self._vb:
226 print "--------------------------------------------------"
227 print " Cursor position"
228 print "--------------------------------------------------"
229 out = 'Beam=%d IF=%d Pol=%d ' % (i,j,k)
230 print out
231 return i,j,k
232
233 def stats(self, stat='stddev', mask=None, allaxes=None):
234 """
235 Determine the specified statistic of the current beam/if/pol
236 Takes a 'mask' as an optional parameter to specify which
237 channels should be excluded.
238 Parameters:
239 stat: 'min', 'max', 'sumsq', 'sum', 'mean'
240 'var', 'stddev', 'avdev', 'rms', 'median'
241 mask: an optional mask specifying where the statistic
242 should be determined.
243 allaxes: if True apply to all spectra. Otherwise
244 apply only to the selected (beam/pol/if)spectra only.
245 The default is taken from .asaprc (True if none)
246 Example:
247 scan.set_unit('channel')
248 msk = scan.create_mask([100,200],[500,600])
249 scan.stats(stat='mean', mask=m)
250 """
251 if allaxes is None: allaxes = rcParams['scantable.allaxes']
252 from asap._asap import stats as _stats
253 from numarray import array,zeros,Float
254 if mask == None:
255 mask = ones(self.nchan())
256 axes = ['Beam','IF','Pol','Time']
257
258 beamSel,IFSel,polSel = (self.getbeam(),self.getif(),self.getpol())
259 if allaxes:
260 n = self.nbeam()*self.nif()*self.npol()*self.nrow()
261 shp = [self.nbeam(),self.nif(),self.npol(),self.nrow()]
262 arr = array(zeros(n),shape=shp,type=Float)
263
264 for i in range(self.nbeam()):
265 self.setbeam(i)
266 for j in range(self.nif()):
267 self.setif(j)
268 for k in range(self.npol()):
269 self.setpol(k)
270 arr[i,j,k,:] = _stats(self,mask,stat,-1)
271 retval = {'axes': axes, 'data': arr, 'cursor':None}
272 tm = [self._gettime(val) for val in range(self.nrow())]
273 if self._vb:
274 self._print_values(retval,stat,tm)
275 self.setbeam(beamSel)
276 self.setif(IFSel)
277 self.setpol(polSel)
278 return retval
279
280 else:
281 statval = _stats(self,mask,stat,-1)
282 out = ''
283 for l in range(self.nrow()):
284 tm = self._gettime(l)
285 out += 'Time[%s]:\n' % (tm)
286 if self.nbeam() > 1: out += ' Beam[%d] ' % (beamSel)
287 if self.nif() > 1: out += ' IF[%d] ' % (IFSel)
288 if self.npol() > 1: out += ' Pol[%d] ' % (polSel)
289 out += '= %3.3f\n' % (statval[l])
290 out += "--------------------------------------------------\n"
291
292 if self._vb:
293 print "--------------------------------------------------"
294 print " ",stat
295 print "--------------------------------------------------"
296 print out
297 retval = {'axes': axes, 'data': array(statval), 'cursor':(i,j,k)}
298 return retval
299
300 def stddev(self,mask=None, allaxes=None):
301 """
302 Determine the standard deviation of the current beam/if/pol
303 Takes a 'mask' as an optional parameter to specify which
304 channels should be excluded.
305 Parameters:
306 mask: an optional mask specifying where the standard
307 deviation should be determined.
308 allaxes: optional flag to show all or a cursor selected
309 spectrum of Beam/IF/Pol. Default is all or taken
310 from .asaprc
311
312 Example:
313 scan.set_unit('channel')
314 msk = scan.create_mask([100,200],[500,600])
315 scan.stddev(mask=m)
316 """
317 if allaxes is None: allaxes = rcParams['scantable.allaxes']
318 return self.stats(stat='stddev',mask=mask, allaxes=allaxes);
319
320 def get_tsys(self, allaxes=None):
321 """
322 Return the System temperatures.
323 Parameters:
324 allaxes: if True apply to all spectra. Otherwise
325 apply only to the selected (beam/pol/if)spectra only.
326 The default is taken from .asaprc (True if none)
327 Returns:
328 a list of Tsys values.
329 """
330 if allaxes is None: allaxes = rcParams['scantable.allaxes']
331 from numarray import array,zeros,Float
332 axes = ['Beam','IF','Pol','Time']
333
334 if allaxes:
335 n = self.nbeam()*self.nif()*self.npol()*self.nrow()
336 shp = [self.nbeam(),self.nif(),self.npol(),self.nrow()]
337 arr = array(zeros(n),shape=shp,type=Float)
338
339 for i in range(self.nbeam()):
340 self.setbeam(i)
341 for j in range(self.nif()):
342 self.setif(j)
343 for k in range(self.npol()):
344 self.setpol(k)
345 arr[i,j,k,:] = self._gettsys()
346 retval = {'axes': axes, 'data': arr, 'cursor':None}
347 tm = [self._gettime(val) for val in range(self.nrow())]
348 if self._vb:
349 self._print_values(retval,'Tsys',tm)
350 return retval
351
352 else:
353 i,j,k = (self.getbeam(),self.getif(),self.getpol())
354 statval = self._gettsys()
355 out = ''
356 for l in range(self.nrow()):
357 tm = self._gettime(l)
358 out += 'Time[%s]:\n' % (tm)
359 if self.nbeam() > 1: out += ' Beam[%d] ' % (i)
360 if self.nif() > 1: out += ' IF[%d] ' % (j)
361 if self.npol() > 1: out += ' Pol[%d] ' % (k)
362 out += '= %3.3f\n' % (statval[l])
363 out += "--------------------------------------------------\n"
364
365 if self._vb:
366 print "--------------------------------------------------"
367 print " TSys"
368 print "--------------------------------------------------"
369 print out
370 retval = {'axes': axes, 'data': array(statval), 'cursor':(i,j,k)}
371 return retval
372
373 def get_time(self, row=-1):
374 """
375 Get a list of time stamps for the observations.
376 Return a string for each integration in the scantable.
377 Parameters:
378 row: row no of integration. Default -1 return all rows
379 Example:
380 none
381 """
382 out = []
383 if row == -1:
384 for i in range(self.nrow()):
385 out.append(self._gettime(i))
386 return out
387 else:
388 if row < self.nrow():
389 return self._gettime(row)
390
391 def set_unit(self, unit='channel'):
392 """
393 Set the unit for all following operations on this scantable
394 Parameters:
395 unit: optional unit, default is 'channel'
396 one of '*Hz','km/s','channel', ''
397 """
398 varlist = vars()
399 if unit in ['','pixel', 'channel']:
400 unit = ''
401 inf = list(self._getcoordinfo())
402 inf[0] = unit
403 self._setcoordinfo(inf)
404 if self._p: self.plot()
405 self._add_history("set_unit",varlist)
406
407 def set_instrument(self, instr):
408 """
409 Set the instrument for subsequent processing
410 Parameters:
411 instr: Select from 'ATPKSMB', 'ATPKSHOH', 'ATMOPRA',
412 'DSS-43' (Tid), 'CEDUNA', and 'HOBART'
413 """
414 self._setInstrument(instr)
415 self._add_history("set_instument",vars())
416
417 def set_doppler(self, doppler='RADIO'):
418 """
419 Set the doppler for all following operations on this scantable.
420 Parameters:
421 doppler: One of 'RADIO', 'OPTICAL', 'Z', 'BETA', 'GAMMA'
422 """
423 varlist = vars()
424 inf = list(self._getcoordinfo())
425 inf[2] = doppler
426 self._setcoordinfo(inf)
427 if self._p: self.plot()
428 self._add_history("set_doppler",vars())
429
430 def set_freqframe(self, frame=None):
431 """
432 Set the frame type of the Spectral Axis.
433 Parameters:
434 frame: an optional frame type, default 'LSRK'.
435 Examples:
436 scan.set_freqframe('BARY')
437 """
438 if frame is None: frame = rcParams['scantable.freqframe']
439 varlist = vars()
440 valid = ['REST','TOPO','LSRD','LSRK','BARY', \
441 'GEO','GALACTO','LGROUP','CMB']
442 if frame in valid:
443 inf = list(self._getcoordinfo())
444 inf[1] = frame
445 self._setcoordinfo(inf)
446 self._add_history("set_freqframe",varlist)
447 else:
448 print "Please specify a valid freq type. Valid types are:\n",valid
449
450 def get_unit(self):
451 """
452 Get the default unit set in this scantable
453 Parameters:
454 Returns:
455 A unit string
456 """
457 inf = self._getcoordinfo()
458 unit = inf[0]
459 if unit == '': unit = 'channel'
460 return unit
461
462 def get_abcissa(self, rowno=0):
463 """
464 Get the abcissa in the current coordinate setup for the currently
465 selected Beam/IF/Pol
466 Parameters:
467 rowno: an optional row number in the scantable. Default is the
468 first row, i.e. rowno=0
469 Returns:
470 The abcissa values and it's format string (as a dictionary)
471 """
472 abc = self._getabcissa(rowno)
473 lbl = self._getabcissalabel(rowno)
474 return abc, lbl
475 #return {'abcissa':abc,'label':lbl}
476
477 def create_mask(self, *args, **kwargs):
478 """
479 Compute and return a mask based on [min,max] windows.
480 The specified windows are to be INCLUDED, when the mask is
481 applied.
482 Parameters:
483 [min,max],[min2,max2],...
484 Pairs of start/end points specifying the regions
485 to be masked
486 invert: optional argument. If specified as True,
487 return an inverted mask, i.e. the regions
488 specified are EXCLUDED
489 row: create the mask using the specified row for
490 unit conversions, default is row=0
491 only necessary if frequency varies over rows.
492 Example:
493 scan.set_unit('channel')
494
495 a)
496 msk = scan.set_mask([400,500],[800,900])
497 # masks everything outside 400 and 500
498 # and 800 and 900 in the unit 'channel'
499
500 b)
501 msk = scan.set_mask([400,500],[800,900], invert=True)
502 # masks the regions between 400 and 500
503 # and 800 and 900 in the unit 'channel'
504
505 """
506 row = 0
507 if kwargs.has_key("row"):
508 row = kwargs.get("row")
509 data = self._getabcissa(row)
510 u = self._getcoordinfo()[0]
511 if self._vb:
512 if u == "": u = "channel"
513 print "The current mask window unit is", u
514 n = self.nchan()
515 msk = zeros(n)
516 for window in args:
517 if (len(window) != 2 or window[0] > window[1] ):
518 print "A window needs to be defined as [min,max]"
519 return
520 for i in range(n):
521 if data[i] >= window[0] and data[i] < window[1]:
522 msk[i] = 1
523 if kwargs.has_key('invert'):
524 if kwargs.get('invert'):
525 from numarray import logical_not
526 msk = logical_not(msk)
527 return msk
528
529 def get_restfreqs(self):
530 """
531 Get the restfrequency(s) stored in this scantable.
532 The return value(s) are always of unit 'Hz'
533 Parameters:
534 none
535 Returns:
536 a list of doubles
537 """
538 return list(self._getrestfreqs())
539
540 def lines(self):
541 """
542 Print the list of known spectral lines
543 """
544 sdtable._lines(self)
545
546 def set_restfreqs(self, freqs=None, unit='Hz', lines=None, source=None,
547 theif=None):
548 """
549 Select the restfrequency for the specified source and IF OR
550 replace for all IFs. If the 'freqs' argument holds a scalar,
551 then that rest frequency will be applied to the selected
552 data (and added to the list of available rest frequencies).
553 In this way, you can set a rest frequency for each
554 source and IF combination. If the 'freqs' argument holds
555 a vector, then it MUST be of length the number of IFs
556 (and the available restfrequencies will be replaced by
557 this vector). In this case, *all* data ('source' and
558 'theif' are ignored) have the restfrequency set per IF according
559 to the corresponding value you give in the 'freqs' vector.
560 E.g. 'freqs=[1e9,2e9]' would mean IF 0 gets restfreq 1e9 and
561 IF 1 gets restfreq 2e9.
562
563 You can also specify the frequencies via known line names
564 in the argument 'lines'. Use 'freqs' or 'lines'. 'freqs'
565 takes precedence. See the list of known names via function
566 scantable.lines()
567 Parameters:
568 freqs: list of rest frequencies
569 unit: unit for rest frequency (default 'Hz')
570 lines: list of known spectral lines names (alternative to freqs).
571 See possible list via scantable.lines()
572 source: Source name (blank means all)
573 theif: IF (-1 means all)
574 Example:
575 scan.set_restfreqs(freqs=1.4e9, source='NGC253', theif=2)
576 scan.set_restfreqs(freqs=[1.4e9,1.67e9])
577 """
578 varlist = vars()
579 if source is None:
580 source = ""
581 if theif is None:
582 theif = -1
583 t = type(freqs)
584 if t is int or t is float:
585 freqs = [freqs]
586 if freqs is None:
587 freqs = []
588 t = type(lines)
589 if t is str:
590 lines = [lines]
591 if lines is None:
592 lines = []
593 sdtable._setrestfreqs(self, freqs, unit, lines, source, theif)
594 self._add_history("set_restfreqs", varlist)
595
596
597
598 def flag_spectrum(self, thebeam, theif, thepol):
599 """
600 This flags a selected spectrum in the scan 'for good'.
601 USE WITH CARE - not reversible.
602 Use masks for non-permanent exclusion of channels.
603 Parameters:
604 thebeam,theif,thepol: all have to be explicitly
605 specified
606 Example:
607 scan.flag_spectrum(0,0,1)
608 flags the spectrum for Beam=0, IF=0, Pol=1
609 """
610 if (thebeam < self.nbeam() and
611 theif < self.nif() and
612 thepol < self.npol()):
613 sdtable.setbeam(self, thebeam)
614 sdtable.setif(self, theif)
615 sdtable.setpol(self, thepol)
616 sdtable._flag(self)
617 self._add_history("flag_spectrum", vars())
618 else:
619 print "Please specify a valid (Beam/IF/Pol)"
620 return
621
622 def plot(self, what='spectrum',col='Pol', panel=None):
623 """
624 Plot the spectra contained in the scan. Alternatively you can also
625 Plot Tsys vs Time
626 Parameters:
627 what: a choice of 'spectrum' (default) or 'tsys'
628 col: which out of Beams/IFs/Pols should be colour stacked
629 panel: set up multiple panels, currently not working.
630 """
631 print "Warning! Not fully functional. Use plotter.plot() instead"
632
633 validcol = {'Beam':self.nbeam(),'IF':self.nif(),'Pol':self.npol()}
634
635 validyax = ['spectrum','tsys']
636 from asap.asaplot import ASAPlot
637 if not self._p:
638 self._p = ASAPlot()
639 #print "Plotting not enabled"
640 #return
641 if self._p.is_dead:
642 del self._p
643 self._p = ASAPlot()
644 npan = 1
645 x = None
646 if what == 'tsys':
647 n = self.nrow()
648 if n < 2:
649 print "Only one integration. Can't plot."
650 return
651 self._p.hold()
652 self._p.clear()
653 if panel == 'Time':
654 npan = self.nrow()
655 self._p.set_panels(rows=npan)
656 xlab,ylab,tlab = None,None,None
657 self._vb = False
658 sel = self.get_cursor()
659 for i in range(npan):
660 if npan > 1:
661 self._p.subplot(i)
662 for j in range(validcol[col]):
663 x = None
664 y = None
665 m = None
666 tlab = self._getsourcename(i)
667 import re
668 tlab = re.sub('_S','',tlab)
669 if col == 'Beam':
670 self.setbeam(j)
671 elif col == 'IF':
672 self.setif(j)
673 elif col == 'Pol':
674 self.setpol(j)
675 if what == 'tsys':
676 x = range(self.nrow())
677 xlab = 'Time [pixel]'
678 m = list(ones(len(x)))
679 y = []
680 ylab = r'$T_{sys}$'
681 for k in range(len(x)):
682 y.append(self._gettsys(k))
683 else:
684 x,xlab = self.get_abcissa(i)
685 y = self._getspectrum(i)
686 ylab = r'Flux'
687 m = self._getmask(i)
688 llab = col+' '+str(j)
689 self._p.set_line(label=llab)
690 self._p.plot(x,y,m)
691 self._p.set_axes('xlabel',xlab)
692 self._p.set_axes('ylabel',ylab)
693 self._p.set_axes('title',tlab)
694 self._p.release()
695 self.set_cursor(sel[0],sel[1],sel[2])
696 self._vb = rcParams['verbose']
697 return
698
699 print out
700
701 def _print_values(self, dat, label='', timestamps=[]):
702 d = dat['data']
703 a = dat['axes']
704 shp = d.getshape()
705 out = ''
706 for i in range(shp[3]):
707 out += '%s [%s]:\n' % (a[3],timestamps[i])
708 t = d[:,:,:,i]
709 for j in range(shp[0]):
710 if shp[0] > 1: out += ' %s[%d] ' % (a[0],j)
711 for k in range(shp[1]):
712 if shp[1] > 1: out += ' %s[%d] ' % (a[1],k)
713 for l in range(shp[2]):
714 if shp[2] > 1: out += ' %s[%d] ' % (a[2],l)
715 out += '= %3.3f\n' % (t[j,k,l])
716 out += "-"*80
717 out += "\n"
718 print "-"*80
719 print " ", label
720 print "-"*80
721 print out
722
723 def history(self):
724 hist = list(self._gethistory())
725 print "-"*80
726 for h in hist:
727 if h.startswith("---"):
728 print h
729 else:
730 items = h.split("##")
731 date = items[0]
732 func = items[1]
733 items = items[2:]
734 print date
735 print "Function: %s\n Parameters:" % (func)
736 for i in items:
737 s = i.split("=")
738 print " %s = %s" % (s[0],s[1])
739 print "-"*80
740 return
741
742 #
743 # Maths business
744 #
745
746 def average_time(self, mask=None, scanav=True, weight=None):
747 """
748 Return the (time) average of a scan, or apply it 'insitu'.
749 Note:
750 in channels only
751 The cursor of the output scan is set to 0.
752 Parameters:
753 one scan or comma separated scans
754 mask: an optional mask (only used for 'var' and 'tsys'
755 weighting)
756 scanav: True (default) averages each scan separately
757 False averages all scans together,
758 weight: Weighting scheme. 'none' (default), 'var' (variance
759 weighted), 'tsys'
760
761 Example:
762 # time average the scantable without using a mask
763 newscan = scan.average_time()
764 """
765 varlist = vars()
766 if weight is None: weight = 'none'
767 if mask is None: mask = ()
768 from asap._asap import average as _av
769 s = scantable(_av((self,), mask, scanav, weight))
770 s._add_history("average_time",varlist)
771 return s
772
773 def convert_flux(self, jyperk=None, eta=None, d=None, insitu=None,
774 allaxes=None):
775 """
776 Return a scan where all spectra are converted to either
777 Jansky or Kelvin depending upon the flux units of the scan table.
778 By default the function tries to look the values up internally.
779 If it can't find them (or if you want to over-ride), you must
780 specify EITHER jyperk OR eta (and D which it will try to look up
781 also if you don't set it). jyperk takes precedence if you set both.
782 Parameters:
783 jyperk: the Jy / K conversion factor
784 eta: the aperture efficiency
785 d: the geomtric diameter (metres)
786 insitu: if False a new scantable is returned.
787 Otherwise, the scaling is done in-situ
788 The default is taken from .asaprc (False)
789 allaxes: if True apply to all spectra. Otherwise
790 apply only to the selected (beam/pol/if)spectra only
791 The default is taken from .asaprc (True if none)
792 """
793 if allaxes is None: allaxes = rcParams['scantable.allaxes']
794 if insitu is None: insitu = rcParams['insitu']
795 varlist = vars()
796 if jyperk is None: jyperk = -1.0
797 if d is None: d = -1.0
798 if eta is None: eta = -1.0
799 if not insitu:
800 from asap._asap import convertflux as _convert
801 s = scantable(_convert(self, d, eta, jyperk, allaxes))
802 s._add_history("convert_flux", varlist)
803 return s
804 else:
805 from asap._asap import convertflux_insitu as _convert
806 _convert(self, d, eta, jyperk, allaxes)
807 self._add_history("convert_flux", varlist)
808 return
809
810 def gain_el(self, poly=None, filename="", method="linear",
811 insitu=None, allaxes=None):
812 """
813 Return a scan after applying a gain-elevation correction.
814 The correction can be made via either a polynomial or a
815 table-based interpolation (and extrapolation if necessary).
816 You specify polynomial coefficients, an ascii table or neither.
817 If you specify neither, then a polynomial correction will be made
818 with built in coefficients known for certain telescopes (an error
819 will occur if the instrument is not known).
820 The data and Tsys are *divided* by the scaling factors.
821 Parameters:
822 poly: Polynomial coefficients (default None) to compute a
823 gain-elevation correction as a function of
824 elevation (in degrees).
825 filename: The name of an ascii file holding correction factors.
826 The first row of the ascii file must give the column
827 names and these MUST include columns
828 "ELEVATION" (degrees) and "FACTOR" (multiply data
829 by this) somewhere.
830 The second row must give the data type of the
831 column. Use 'R' for Real and 'I' for Integer.
832 An example file would be
833 (actual factors are arbitrary) :
834
835 TIME ELEVATION FACTOR
836 R R R
837 0.1 0 0.8
838 0.2 20 0.85
839 0.3 40 0.9
840 0.4 60 0.85
841 0.5 80 0.8
842 0.6 90 0.75
843 method: Interpolation method when correcting from a table.
844 Values are "nearest", "linear" (default), "cubic"
845 and "spline"
846 insitu: if False a new scantable is returned.
847 Otherwise, the scaling is done in-situ
848 The default is taken from .asaprc (False)
849 allaxes: If True apply to all spectra. Otherwise
850 apply only to the selected (beam/pol/if) spectra only
851 The default is taken from .asaprc (True if none)
852 """
853
854 if allaxes is None: allaxes = rcParams['scantable.allaxes']
855 if insitu is None: insitu = rcParams['insitu']
856 varlist = vars()
857 if poly is None:
858 poly = ()
859 from os.path import expandvars
860 filename = expandvars(filename)
861 if not insitu:
862 from asap._asap import gainel as _gainEl
863 s = scantable(_gainEl(self, poly, filename, method, allaxes))
864 s._add_history("gain_el", varlist)
865 return s
866 else:
867 from asap._asap import gainel_insitu as _gainEl
868 _gainEl(self, poly, filename, method, allaxes)
869 self._add_history("gain_el", varlist)
870 return
871
872 def freq_align(self, reftime=None, method='cubic', perif=False,
873 insitu=None):
874 """
875 Return a scan where all rows have been aligned in frequency/velocity.
876 The alignment frequency frame (e.g. LSRK) is that set by function
877 set_freqframe.
878 Parameters:
879 reftime: reference time to align at. By default, the time of
880 the first row of data is used.
881 method: Interpolation method for regridding the spectra.
882 Choose from "nearest", "linear", "cubic" (default)
883 and "spline"
884 perif: Generate aligners per freqID (no doppler tracking) or
885 per IF (scan-based doppler tracking)
886 insitu: if False a new scantable is returned.
887 Otherwise, the scaling is done in-situ
888 The default is taken from .asaprc (False)
889 """
890 if insitu is None: insitu = rcParams['insitu']
891 varlist = vars()
892 if reftime is None: reftime = ''
893 perfreqid = not perif
894 if not insitu:
895 from asap._asap import freq_align as _align
896 s = scantable(_align(self, reftime, method, perfreqid))
897 s._add_history("freq_align", varlist)
898 return s
899 else:
900 from asap._asap import freq_align_insitu as _align
901 _align(self, reftime, method, perfreqid)
902 self._add_history("freq_align", varlist)
903 return
904
905 def opacity(self, tau, insitu=None, allaxes=None):
906 """
907 Apply an opacity correction. The data
908 and Tsys are multiplied by the correction factor.
909 Parameters:
910 tau: Opacity from which the correction factor is
911 exp(tau*ZD)
912 where ZD is the zenith-distance
913 insitu: if False a new scantable is returned.
914 Otherwise, the scaling is done in-situ
915 The default is taken from .asaprc (False)
916 allaxes: if True apply to all spectra. Otherwise
917 apply only to the selected (beam/pol/if)spectra only
918 The default is taken from .asaprc (True if none)
919 """
920 if allaxes is None: allaxes = rcParams['scantable.allaxes']
921 if insitu is None: insitu = rcParams['insitu']
922 varlist = vars()
923 if not insitu:
924 from asap._asap import opacity as _opacity
925 s = scantable(_opacity(self, tau, allaxes))
926 s._add_history("opacity", varlist)
927 return s
928 else:
929 from asap._asap import opacity_insitu as _opacity
930 _opacity(self, tau, allaxes)
931 self._add_history("opacity", varlist)
932 return
933
934 def bin(self, width=5, insitu=None):
935 """
936 Return a scan where all spectra have been binned up.
937 width: The bin width (default=5) in pixels
938 insitu: if False a new scantable is returned.
939 Otherwise, the scaling is done in-situ
940 The default is taken from .asaprc (False)
941 """
942 if insitu is None: insitu = rcParams['insitu']
943 varlist = vars()
944 if not insitu:
945 from asap._asap import bin as _bin
946 s = scantable(_bin(self, width))
947 s._add_history("bin",varlist)
948 return s
949 else:
950 from asap._asap import bin_insitu as _bin
951 _bin(self, width)
952 self._add_history("bin",varlist)
953 return
954
955
956 def resample(self, width=5, method='cubic', insitu=None):
957 """
958 Return a scan where all spectra have been binned up
959 width: The bin width (default=5) in pixels
960 method: Interpolation method when correcting from a table.
961 Values are "nearest", "linear", "cubic" (default)
962 and "spline"
963 insitu: if False a new scantable is returned.
964 Otherwise, the scaling is done in-situ
965 The default is taken from .asaprc (False)
966 """
967 if insitu is None: insitu = rcParams['insitu']
968 varlist = vars()
969 if not insitu:
970 from asap._asap import resample as _resample
971 s = scantable(_resample(self, method, width))
972 s._add_history("resample",varlist)
973 return s
974 else:
975 from asap._asap import resample_insitu as _resample
976 _resample(self, method, width)
977 self._add_history("resample",varlist)
978 return
979
980 def average_pol(self, mask=None, weight='none', insitu=None):
981 """
982 Average the Polarisations together.
983 The polarisation cursor of the output scan is set to 0
984 Parameters:
985 scan: The scantable
986 mask: An optional mask defining the region, where the
987 averaging will be applied. The output will have all
988 specified points masked.
989 weight: Weighting scheme. 'none' (default), or
990 'var' (variance weighted)
991 insitu: if False a new scantable is returned.
992 Otherwise, the scaling is done in-situ
993 The default is taken from .asaprc (False)
994 Example:
995 polav = average_pols(myscan)
996 """
997 if insitu is None: insitu = rcParams['insitu']
998 varlist = vars()
999
1000 if mask is None:
1001 mask = ()
1002 if not insitu:
1003 from asap._asap import averagepol as _avpol
1004 s = scantable(_avpol(self, mask, weight))
1005 s._add_history("average_pol",varlist)
1006 return s
1007 else:
1008 from asap._asap import averagepol_insitu as _avpol
1009 _avpol(self, mask, weight)
1010 self._add_history("average_pol",varlist)
1011 return
1012
1013 def smooth(self, kernel="hanning", width=5.0, insitu=None, allaxes=None):
1014 """
1015 Smooth the spectrum by the specified kernel (conserving flux).
1016 Parameters:
1017 scan: The input scan
1018 kernel: The type of smoothing kernel. Select from
1019 'hanning' (default), 'gaussian' and 'boxcar'.
1020 The first three characters are sufficient.
1021 width: The width of the kernel in pixels. For hanning this is
1022 ignored otherwise it defauls to 5 pixels.
1023 For 'gaussian' it is the Full Width Half
1024 Maximum. For 'boxcar' it is the full width.
1025 insitu: if False a new scantable is returned.
1026 Otherwise, the scaling is done in-situ
1027 The default is taken from .asaprc (False)
1028 allaxes: If True (default) apply to all spectra. Otherwise
1029 apply only to the selected (beam/pol/if)spectra only
1030 The default is taken from .asaprc (True if none)
1031 Example:
1032 none
1033 """
1034 if allaxes is None: allaxes = rcParams['scantable.allaxes']
1035 if insitu is None: insitu = rcParams['insitu']
1036 varlist = vars()
1037 if not insitu:
1038 from asap._asap import smooth as _smooth
1039 s = scantable(_smooth(self,kernel,width,allaxes))
1040 s._add_history("smooth", varlist)
1041 return s
1042 else:
1043 from asap._asap import smooth_insitu as _smooth
1044 _smooth(self,kernel,width,allaxes)
1045 self._add_history("smooth", varlist)
1046 return
1047
1048 def poly_baseline(self, mask=None, order=0, insitu=None):
1049 """
1050 Return a scan which has been baselined (all rows) by a polynomial.
1051 Parameters:
1052 scan: a scantable
1053 mask: an optional mask
1054 order: the order of the polynomial (default is 0)
1055 insitu: if False a new scantable is returned.
1056 Otherwise, the scaling is done in-situ
1057 The default is taken from .asaprc (False)
1058 Example:
1059 # return a scan baselined by a third order polynomial,
1060 # not using a mask
1061 bscan = scan.poly_baseline(order=3)
1062 """
1063 if insitu is None: insitu = rcParams['insitu']
1064 varlist = vars()
1065 if mask is None:
1066 from numarray import ones
1067 mask = list(ones(scan.nchan()))
1068 from asap.asapfitter import fitter
1069 f = fitter()
1070 f._verbose(True)
1071 f.set_scan(self, mask)
1072 f.set_function(poly=order)
1073 sf = f.auto_fit(insitu)
1074 if insitu:
1075 self._add_history("poly_baseline", varlist)
1076 return
1077 else:
1078 sf._add_history("poly_baseline", varlist)
1079 return sf
1080
1081 def auto_poly_baseline(self, mask=None, edge=(0,0), order=0,
1082 threshold=3,insitu=None):
1083 """
1084 Return a scan which has been baselined (all rows) by a polynomial.
1085 Spectral lines are detected first using linefinder and masked out
1086 to avoid them affecting the baseline solution.
1087
1088 Parameters:
1089 scan: a scantable
1090 mask: an optional mask retreived from scantable
1091 edge: an optional number of channel to drop at
1092 the edge of spectrum. If only one value is
1093 specified, the same number will be dropped from
1094 both sides of the spectrum. Default is to keep
1095 all channels
1096 order: the order of the polynomial (default is 0)
1097 threshold: the threshold used by line finder. It is better to
1098 keep it large as only strong lines affect the
1099 baseline solution.
1100 insitu: if False a new scantable is returned.
1101 Otherwise, the scaling is done in-situ
1102 The default is taken from .asaprc (False)
1103
1104 Example:
1105 scan2=scan.auto_poly_baseline(order=7)
1106 """
1107 if insitu is None: insitu = rcParams['insitu']
1108 varlist = vars()
1109 from asap.asapfitter import fitter
1110 from asap.asaplinefind import linefinder
1111 from asap import _is_sequence_or_number as _is_valid
1112
1113 if not _is_valid(edge, int):
1114 raise RuntimeError, "Parameter 'edge' has to be an integer or a \
1115 pair of integers specified as a tuple"
1116
1117 # setup fitter
1118 f = fitter()
1119 f._verbose(True)
1120 f.set_function(poly=order)
1121
1122 # setup line finder
1123 fl=linefinder()
1124 fl.set_options(threshold=threshold)
1125
1126 if not insitu:
1127 workscan=self.copy()
1128 else:
1129 workscan=self
1130
1131 vb=workscan._vb
1132 # remember the verbose parameter and selection
1133 workscan._vb=False
1134 sel=workscan.get_cursor()
1135 rows=range(workscan.nrow())
1136 for i in range(workscan.nbeam()):
1137 workscan.setbeam(i)
1138 for j in range(workscan.nif()):
1139 workscan.setif(j)
1140 for k in range(workscan.npol()):
1141 workscan.setpol(k)
1142 if f._vb:
1143 print "Processing:"
1144 print 'Beam[%d], IF[%d], Pol[%d]' % (i,j,k)
1145 for iRow in rows:
1146 fl.set_scan(workscan,mask,edge)
1147 fl.find_lines(iRow)
1148 f.set_scan(workscan, fl.get_mask())
1149 f.x=workscan._getabcissa(iRow)
1150 f.y=workscan._getspectrum(iRow)
1151 f.data=None
1152 f.fit()
1153 x=f.get_parameters()
1154 workscan._setspectrum(f.fitter.getresidual(),iRow)
1155 workscan.set_cursor(sel[0],sel[1],sel[2])
1156 workscan._vb = vb
1157 if not insitu:
1158 return workscan
1159
1160 def rotate_linpolphase(self, angle, allaxes=None):
1161 """
1162 Rotate the phase of the complex polarization O=Q+iU correlation.
1163 This is always done in situ in the raw data. So if you call this
1164 function more than once then each call rotates the phase further.
1165 Parameters:
1166 angle: The angle (degrees) to rotate (add) by.
1167 allaxes: If True apply to all spectra. Otherwise
1168 apply only to the selected (beam/pol/if)spectra only.
1169 The default is taken from .asaprc (True if none)
1170 Examples:
1171 scan.rotate_linpolphase(2.3)
1172 """
1173 if allaxes is None: allaxes = rcParams['scantable.allaxes']
1174 varlist = vars()
1175 from asap._asap import _rotate_linpolphase as _rotate
1176 _rotate(self, angle, allaxes)
1177 self._add_history("rotate_linpolphase", varlist)
1178 return
1179
1180
1181 def rotate_xyphase(self, angle, allaxes=None):
1182 """
1183 Rotate the phase of the XY correlation. This is always done in situ
1184 in the data. So if you call this function more than once
1185 then each call rotates the phase further.
1186 Parameters:
1187 angle: The angle (degrees) to rotate (add) by.
1188 allaxes: If True apply to all spectra. Otherwise
1189 apply only to the selected (beam/pol/if)spectra only.
1190 The default is taken from .asaprc (True if none)
1191 Examples:
1192 scan.rotate_xyphase(2.3)
1193 """
1194 if allaxes is None: allaxes = rcParams['scantable.allaxes']
1195 varlist = vars()
1196 from asap._asap import rotate_xyphase as _rotate_xyphase
1197 _rotate_xyphase(self, angle, allaxes)
1198 self._add_history("rotate_xyphase", varlist)
1199 return
1200
1201
1202 def add(self, offset, insitu=None, allaxes=None):
1203 """
1204 Return a scan where all spectra have the offset added
1205 Parameters:
1206 offset: the offset
1207 insitu: if False a new scantable is returned.
1208 Otherwise, the scaling is done in-situ
1209 The default is taken from .asaprc (False)
1210 allaxes: if True apply to all spectra. Otherwise
1211 apply only to the selected (beam/pol/if)spectra only
1212 The default is taken from .asaprc (True if none)
1213 """
1214 if allaxes is None: allaxes = rcParams['scantable.allaxes']
1215 if insitu is None: insitu = rcParams['insitu']
1216 varlist = vars()
1217 if not insitu:
1218 from asap._asap import add as _add
1219 s = scantable(_add(self, offset, allaxes))
1220 s._add_history("add",varlist)
1221 return s
1222 else:
1223 from asap._asap import add_insitu as _add
1224 _add(self, offset, allaxes)
1225 self._add_history("add",varlist)
1226 return
1227
1228 def scale(self, factor, insitu=None, allaxes=None, tsys=True):
1229 """
1230 Return a scan where all spectra are scaled by the give 'factor'
1231 Parameters:
1232 factor: the scaling factor
1233 insitu: if False a new scantable is returned.
1234 Otherwise, the scaling is done in-situ
1235 The default is taken from .asaprc (False)
1236 allaxes: if True apply to all spectra. Otherwise
1237 apply only to the selected (beam/pol/if)spectra only.
1238 The default is taken from .asaprc (True if none)
1239 tsys: if True (default) then apply the operation to Tsys
1240 as well as the data
1241 """
1242 if allaxes is None: allaxes = rcParams['scantable.allaxes']
1243 if insitu is None: insitu = rcParams['insitu']
1244 varlist = vars()
1245 if not insitu:
1246 from asap._asap import scale as _scale
1247 s = scantable(_scale(self, factor, allaxes, tsys))
1248 s._add_history("scale",varlist)
1249 return s
1250 else:
1251 from asap._asap import scale_insitu as _scale
1252 _scale(self, factor, allaxes, tsys)
1253 self._add_history("scale",varlist)
1254 return
1255
1256
1257 def quotient(self, other, isreference=True, preserve=True):
1258 """
1259 Return the quotient of a 'source' (signal) scan and a 'reference' scan.
1260 The reference can have just one row, even if the signal has many.
1261 Otherwise they must have the same number of rows.
1262 The cursor of the output scan is set to 0
1263 Parameters:
1264 source: the 'other' scan
1265 isreference: if the 'other' scan is the reference (default)
1266 or source
1267 preserve: you can preserve (default) the continuum or
1268 remove it. The equations used are
1269 preserve - Output = Tref * (sig/ref) - Tref
1270 remove - Output = Tref * (sig/ref) - Tsig
1271 Example:
1272 # src is a scantable for an 'on' scan, ref for an 'off' scantable
1273 q1 = src.quotient(ref)
1274 q2 = ref.quotient(src, isreference=False)
1275 # gives the same result
1276 """
1277 from asap._asap import quotient as _quot
1278 if isreference:
1279 return scantable(_quot(self, other, preserve))
1280 else:
1281 return scantable(_quot(other, self, preserve))
1282
1283 def __add__(self, other):
1284 varlist = vars()
1285 s = None
1286 if isinstance(other, scantable):
1287 from asap._asap import b_operate as _bop
1288 s = scantable(_bop(self, other, 'add', True))
1289 elif isinstance(other, float):
1290 from asap._asap import add as _add
1291 s = scantable(_add(self, other, True))
1292 else:
1293 print "Other input is not a scantable or float value"
1294 return
1295 s._add_history("operator +", varlist)
1296 return s
1297
1298
1299 def __sub__(self, other):
1300 """
1301 implicit on all axes and on Tsys
1302 """
1303 varlist = vars()
1304 s = None
1305 if isinstance(other, scantable):
1306 from asap._asap import b_operate as _bop
1307 s = scantable(_bop(self, other, 'sub', True))
1308 elif isinstance(other, float):
1309 from asap._asap import add as _add
1310 s = scantable(_add(self, -other, True))
1311 else:
1312 print "Other input is not a scantable or float value"
1313 return
1314 s._add_history("operator -", varlist)
1315 return s
1316
1317 def __mul__(self, other):
1318 """
1319 implicit on all axes and on Tsys
1320 """
1321 varlist = vars()
1322 s = None
1323 if isinstance(other, scantable):
1324 from asap._asap import b_operate as _bop
1325 s = scantable(_bop(self, other, 'mul', True))
1326 elif isinstance(other, float):
1327 if other == 0.0:
1328 print "Multiplying by zero is not recommended."
1329 return
1330 from asap._asap import scale as _sca
1331 s = scantable(_sca(self, other, True))
1332 else:
1333 print "Other input is not a scantable or float value"
1334 return
1335 s._add_history("operator *", varlist)
1336 return s
1337
1338
1339 def __div__(self, other):
1340 """
1341 implicit on all axes and on Tsys
1342 """
1343 varlist = vars()
1344 s = None
1345 if isinstance(other, scantable):
1346 from asap._asap import b_operate as _bop
1347 s = scantable(_bop(self, other, 'div', True))
1348 elif isinstance(other, float):
1349 if other == 0.0:
1350 print "Dividing by zero is not recommended"
1351 return
1352 from asap._asap import scale as _sca
1353 s = scantable(_sca(self, 1.0/other, True))
1354 else:
1355 print "Other input is not a scantable or float value"
1356 return
1357 s._add_history("operator /", varlist)
1358 return s
1359
1360 def _add_history(self, funcname, parameters):
1361 # create date
1362 sep = "##"
1363 from datetime import datetime
1364 dstr = datetime.now().strftime('%Y/%m/%d %H:%M:%S')
1365 hist = dstr+sep
1366 hist += funcname+sep#cdate+sep
1367 if parameters.has_key('self'): del parameters['self']
1368 for k,v in parameters.iteritems():
1369 if type(v) is dict:
1370 for k2,v2 in v.iteritems():
1371 hist += k2
1372 hist += "="
1373 if isinstance(v2,scantable):
1374 hist += 'scantable'
1375 elif k2 == 'mask':
1376 if isinstance(v2,list) or isinstance(v2,tuple):
1377 hist += str(self._zip_mask(v2))
1378 else:
1379 hist += str(v2)
1380 else:
1381 hist += str(v2)
1382 else:
1383 hist += k
1384 hist += "="
1385 if isinstance(v,scantable):
1386 hist += 'scantable'
1387 elif k == 'mask':
1388 if isinstance(v,list) or isinstance(v,tuple):
1389 hist += str(self._zip_mask(v))
1390 else:
1391 hist += str(v)
1392 else:
1393 hist += str(v)
1394 hist += sep
1395 hist = hist[:-2] # remove trailing '##'
1396 self._addhistory(hist)
1397
1398
1399 def _zip_mask(self, mask):
1400 mask = list(mask)
1401 i = 0
1402 segments = []
1403 while mask[i:].count(1):
1404 i += mask[i:].index(1)
1405 if mask[i:].count(0):
1406 j = i + mask[i:].index(0)
1407 else:
1408 j = len(mask)
1409 segments.append([i,j])
1410 i = j
1411 return segments
Note: See TracBrowser for help on using the repository browser.