source: trunk/python/scantable.py@ 1436

Last change on this file since 1436 was 1435, checked in by Malte Marquarding, 16 years ago

optionally suppress history

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