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

Last change on this file since 799 was 798, checked in by phi196, 19 years ago

Added swap_pol & invert_phase
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 62.8 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 def invert_phase(self, allaxes=None):
1271 """
1272 Take the complex conjugate of the complex polarisation cross
1273 correlation. This is always done in situ in the data. If
1274 you call this function twice, you end up with the original phase
1275 Parameters:
1276 allaxes: If True apply to all spectra. Otherwise
1277 apply only to the selected (beam/pol/if)spectra only.
1278 The default is taken from .asaprc (True if none)
1279 Examples:
1280 scan.invert_phase()
1281 """
1282 if allaxes is None: allaxes = rcParams['scantable.allaxes']
1283 varlist = vars()
1284 from asap._asap import _invert_phase
1285 _invert_phase(self, allaxes)
1286 self._add_history("invert_phase", varlist)
1287 print_log()
1288 return
1289
1290 def swap_ppol(self, allaxes=None):
1291 """
1292 Swaps the order of parallel polarisation products (ie AA, BB
1293 becomes BB, AA) . This is always done in situ in the data.
1294 If you call this function twice, you end up with the original
1295 order
1296 Parameters:
1297 allaxes: If True apply to all spectra. Otherwise
1298 apply only to the selected (beam/pol/if)spectra only.
1299 The default is taken from .asaprc (True if none)
1300 Examples:
1301 scan.swap_pol()
1302 """
1303 if allaxes is None: allaxes = rcParams['scantable.allaxes']
1304 varlist = vars()
1305 from asap._asap import _swap_pol
1306 _swap_pol(self, allaxes)
1307 self._add_history("swap_phase", varlist)
1308 print_log()
1309 return
1310
1311
1312 def add(self, offset, insitu=None, allaxes=None):
1313 """
1314 Return a scan where all spectra have the offset added
1315 Parameters:
1316 offset: the offset
1317 insitu: if False a new scantable is returned.
1318 Otherwise, the scaling is done in-situ
1319 The default is taken from .asaprc (False)
1320 allaxes: if True apply to all spectra. Otherwise
1321 apply only to the selected (beam/pol/if)spectra only
1322 The default is taken from .asaprc (True if none)
1323 """
1324 if allaxes is None: allaxes = rcParams['scantable.allaxes']
1325 if insitu is None: insitu = rcParams['insitu']
1326 varlist = vars()
1327 if not insitu:
1328 from asap._asap import add as _add
1329 s = scantable(_add(self, offset, allaxes))
1330 s._add_history("add",varlist)
1331 print_log()
1332 return s
1333 else:
1334 from asap._asap import add_insitu as _add
1335 _add(self, offset, allaxes)
1336 self._add_history("add",varlist)
1337 print_log()
1338 return
1339
1340 def scale(self, factor, insitu=None, allaxes=None, tsys=True):
1341 """
1342 Return a scan where all spectra are scaled by the give 'factor'
1343 Parameters:
1344 factor: the scaling factor
1345 insitu: if False a new scantable is returned.
1346 Otherwise, the scaling is done in-situ
1347 The default is taken from .asaprc (False)
1348 allaxes: if True apply to all spectra. Otherwise
1349 apply only to the selected (beam/pol/if)spectra only.
1350 The default is taken from .asaprc (True if none)
1351 tsys: if True (default) then apply the operation to Tsys
1352 as well as the data
1353 """
1354 if allaxes is None: allaxes = rcParams['scantable.allaxes']
1355 if insitu is None: insitu = rcParams['insitu']
1356 varlist = vars()
1357 if not insitu:
1358 from asap._asap import scale as _scale
1359 s = scantable(_scale(self, factor, allaxes, tsys))
1360 s._add_history("scale",varlist)
1361 print_log()
1362 return s
1363 else:
1364 from asap._asap import scale_insitu as _scale
1365 _scale(self, factor, allaxes, tsys)
1366 self._add_history("scale",varlist)
1367 print_log()
1368 return
1369
1370 def auto_quotient(self, mode='suffix', preserve=True):
1371 """
1372 This function allows to build quotients automatically.
1373 It assumes the observation to have the same numer of
1374 "ons" and "offs"
1375 It will support "closest off in time" in the future
1376 Parameters:
1377 mode: the on/off detection mode; 'suffix' (default)
1378 'suffix' identifies 'off' scans by the
1379 trailing '_R' (Mopra/Parkes) or
1380 '_e'/'_w' (Tid)
1381 preserve: you can preserve (default) the continuum or
1382 remove it. The equations used are
1383 preserve: Output = Toff * (on/off) - Toff
1384 remove: Output = Tref * (on/off) - Ton
1385 """
1386 modes = ["suffix","time"]
1387 if not mode in modes:
1388 print "please provide valid mode. Valid modes are %s" % (modes)
1389 return None
1390 from asap._asap import quotient as _quot
1391 if mode == "suffix":
1392 srcs = self.get_scan("*[^_ewR]")
1393 refs = self.get_scan("*[_ewR]")
1394 if isinstance(srcs,scantable) and isinstance(refs,scantable):
1395 from asap import asaplog
1396 ns,nr = srcs.nrow(),refs.nrow()
1397 msg = "Found %i Off and %i On scans" % (ns,nr)
1398 asaplog.push(msg)
1399 if nr > ns:
1400 asaplog("Found more Off integrations than On scans - dropping excess Offs.")
1401 refs = refs.get_scan(range(ns))
1402 print_log()
1403 return scantable(_quot(srcs,refs, preserve))
1404 else:
1405 msg = "Couldn't find any on/off pairs"
1406 if rcParams['verbose']:
1407 print msg
1408 return
1409 else:
1410 raise RuntimeError()
1411 else:
1412 if rcParams['verbose']: print "not yet implemented"
1413 return None
1414
1415 def quotient(self, other, isreference=True, preserve=True):
1416 """
1417 Return the quotient of a 'source' (on) scan and a 'reference' (off)
1418 scan.
1419 The reference can have just one row, even if the signal has many.
1420 Otherwise they must have the same number of rows.
1421 The cursor of the output scan is set to 0
1422 Parameters:
1423 other: the 'other' scan
1424 isreference: if the 'other' scan is the reference (default)
1425 or source
1426 preserve: you can preserve (default) the continuum or
1427 remove it. The equations used are
1428 preserve: Output = Toff * (on/off) - Toff
1429 remove: Output = Tref * (on/off) - Ton
1430 Example:
1431 # src is a scantable for an 'on' scan, ref for an 'off' scantable
1432 q1 = src.quotient(ref)
1433 q2 = ref.quotient(src, isreference=False)
1434 # gives the same result
1435 """
1436 from asap._asap import quotient as _quot
1437 try:
1438 s = None
1439 if isreference:
1440 s = scantable(_quot(self, other, preserve))
1441 else:
1442 s = scantable(_quot(other, self, preserve))
1443 print_log()
1444 return s
1445 except RuntimeError,e:
1446 if rcParams['verbose']:
1447 print e
1448 return
1449 else: raise
1450
1451 def freq_switch(self, insitu=None):
1452 """
1453 Apply frequency switching to the data.
1454 Parameters:
1455 insitu: if False a new scantable is returned.
1456 Otherwise, the swictching is done in-situ
1457 The default is taken from .asaprc (False)
1458 Example:
1459 none
1460 """
1461 if insitu is None: insitu = rcParams['insitu']
1462 varlist = vars()
1463 try:
1464 if insitu:
1465 from asap._asap import _frequency_switch_insitu
1466 _frequency_switch_insitu(self)
1467 self._add_history("freq_switch", varlist)
1468 print_log()
1469 return
1470 else:
1471 from asap._asap import _frequency_switch
1472 sf = scantable(_frequency_switch(self))
1473 sf._add_history("freq_switch", varlist)
1474 print_log()
1475 return sf
1476 except RuntimeError,e:
1477 if rcParams['verbose']: print e
1478 else: raise
1479
1480 def recalc_azel(self):
1481 """
1482 Recalculate the azimuth and elevation for each position.
1483 Parameters:
1484 none
1485 Example:
1486 """
1487 varlist = vars()
1488 self._recalc_azel()
1489 self._add_history("recalc_azel", varlist)
1490 print_log()
1491 return
1492
1493 def __add__(self, other):
1494 varlist = vars()
1495 s = None
1496 if isinstance(other, scantable):
1497 from asap._asap import b_operate as _bop
1498 s = scantable(_bop(self, other, 'add', True))
1499 elif isinstance(other, float):
1500 from asap._asap import add as _add
1501 s = scantable(_add(self, other, True))
1502 else:
1503 raise TypeError("Other input is not a scantable or float value")
1504 s._add_history("operator +", varlist)
1505 print_log()
1506 return s
1507
1508 def __sub__(self, other):
1509 """
1510 implicit on all axes and on Tsys
1511 """
1512 varlist = vars()
1513 s = None
1514 if isinstance(other, scantable):
1515 from asap._asap import b_operate as _bop
1516 s = scantable(_bop(self, other, 'sub', True))
1517 elif isinstance(other, float):
1518 from asap._asap import add as _add
1519 s = scantable(_add(self, -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 __mul__(self, other):
1527 """
1528 implicit on all axes and on Tsys
1529 """
1530 varlist = vars()
1531 s = None
1532 if isinstance(other, scantable):
1533 from asap._asap import b_operate as _bop
1534 s = scantable(_bop(self, other, 'mul', True))
1535 elif isinstance(other, float):
1536 if other == 0.0:
1537 raise ZeroDivisionError("Dividing by zero is not recommended")
1538 from asap._asap import scale as _sca
1539 s = scantable(_sca(self, other, True))
1540 else:
1541 raise TypeError("Other input is not a scantable or float value")
1542 s._add_history("operator *", varlist)
1543 print_log()
1544 return s
1545
1546
1547 def __div__(self, other):
1548 """
1549 implicit on all axes and on Tsys
1550 """
1551 varlist = vars()
1552 s = None
1553 if isinstance(other, scantable):
1554 from asap._asap import b_operate as _bop
1555 s = scantable(_bop(self, other, 'div', True))
1556 elif isinstance(other, float):
1557 if other == 0.0:
1558 raise ZeroDivisionError("Dividing by zero is not recommended")
1559 from asap._asap import scale as _sca
1560 s = scantable(_sca(self, 1.0/other, True))
1561 else:
1562 raise TypeError("Other input is not a scantable or float value")
1563 s._add_history("operator /", varlist)
1564 print_log()
1565 return s
1566
1567 def get_fit(self, row=0):
1568 """
1569 Print or return the stored fits for a row in the scantable
1570 Parameters:
1571 row: the row which the fit has been applied to.
1572 """
1573 if row > self.nrow():
1574 return
1575 from asap import asapfit
1576 fit = asapfit(self._getfit(row))
1577 if rcParams['verbose']:
1578 print fit
1579 return
1580 else:
1581 return fit.as_dict()
1582
1583 def _add_history(self, funcname, parameters):
1584 # create date
1585 sep = "##"
1586 from datetime import datetime
1587 dstr = datetime.now().strftime('%Y/%m/%d %H:%M:%S')
1588 hist = dstr+sep
1589 hist += funcname+sep#cdate+sep
1590 if parameters.has_key('self'): del parameters['self']
1591 for k,v in parameters.iteritems():
1592 if type(v) is dict:
1593 for k2,v2 in v.iteritems():
1594 hist += k2
1595 hist += "="
1596 if isinstance(v2,scantable):
1597 hist += 'scantable'
1598 elif k2 == 'mask':
1599 if isinstance(v2,list) or isinstance(v2,tuple):
1600 hist += str(self._zip_mask(v2))
1601 else:
1602 hist += str(v2)
1603 else:
1604 hist += str(v2)
1605 else:
1606 hist += k
1607 hist += "="
1608 if isinstance(v,scantable):
1609 hist += 'scantable'
1610 elif k == 'mask':
1611 if isinstance(v,list) or isinstance(v,tuple):
1612 hist += str(self._zip_mask(v))
1613 else:
1614 hist += str(v)
1615 else:
1616 hist += str(v)
1617 hist += sep
1618 hist = hist[:-2] # remove trailing '##'
1619 self._addhistory(hist)
1620
1621
1622 def _zip_mask(self, mask):
1623 mask = list(mask)
1624 i = 0
1625 segments = []
1626 while mask[i:].count(1):
1627 i += mask[i:].index(1)
1628 if mask[i:].count(0):
1629 j = i + mask[i:].index(0)
1630 else:
1631 j = len(mask)
1632 segments.append([i,j])
1633 i = j
1634 return segments
1635
1636 def _get_ordinate_label(self):
1637 fu = "("+self.get_fluxunit()+")"
1638 import re
1639 lbl = "Intensity"
1640 if re.match(".K.",fu):
1641 lbl = "Brightness Temperature "+ fu
1642 elif re.match(".Jy.",fu):
1643 lbl = "Flux density "+ fu
1644 return lbl
1645
Note: See TracBrowser for help on using the repository browser.