source: trunk/python/scantable.py@ 836

Last change on this file since 836 was 794, checked in by mar637, 19 years ago

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