source: branches/Release12/python/scantable.py@ 789

Last change on this file since 789 was 788, checked in by mar637, 19 years ago

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