source: trunk/python/scantable.py@ 1574

Last change on this file since 1574 was 1574, checked in by Malte Marquarding, 15 years ago

Ticket #167: python part of running polynomial smoothing

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