source: trunk/python/scantable.py@ 1113

Last change on this file since 1113 was 1099, checked in by mar637, 18 years ago

moved average_channel into average_time as an addition weight

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