source: trunk/python/scantable.py@ 1172

Last change on this file since 1172 was 1157, checked in by mar637, 18 years ago

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