source: trunk/python/scantable.py@ 713

Last change on this file since 713 was 710, checked in by mar637, 19 years ago

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