source: trunk/python/scantable.py@ 2031

Last change on this file since 2031 was 2029, checked in by Takeshi Nakazato, 14 years ago

New Development: No

JIRA Issue: Yes CAS-2718

Ready for Test: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No

Module(s): Module Names change impacts.

Description: Describe your changes here...

Replaced STWriter with MSWriter for MS output.


  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 111.7 KB
Line 
1"""This module defines the scantable class."""
2
3import os
4import numpy
5try:
6 from functools import wraps as wraps_dec
7except ImportError:
8 from asap.compatibility import wraps as wraps_dec
9
10from asap.env import is_casapy
11from asap._asap import Scantable
12from asap._asap import filler, msfiller
13from asap.parameters import rcParams
14from asap.logging import asaplog, asaplog_post_dec
15from asap.selector import selector
16from asap.linecatalog import linecatalog
17from asap.coordinate import coordinate
18from asap.utils import _n_bools, mask_not, mask_and, mask_or, page
19from asap.asapfitter import fitter
20
21
22def preserve_selection(func):
23 @wraps_dec(func)
24 def wrap(obj, *args, **kw):
25 basesel = obj.get_selection()
26 try:
27 val = func(obj, *args, **kw)
28 finally:
29 obj.set_selection(basesel)
30 return val
31 return wrap
32
33def is_scantable(filename):
34 """Is the given file a scantable?
35
36 Parameters:
37
38 filename: the name of the file/directory to test
39
40 """
41 if ( os.path.isdir(filename)
42 and os.path.exists(filename+'/table.info')
43 and os.path.exists(filename+'/table.dat') ):
44 f=open(filename+'/table.info')
45 l=f.readline()
46 f.close()
47 #if ( l.find('Scantable') != -1 ):
48 if ( l.find('Measurement Set') == -1 ):
49 return True
50 else:
51 return False
52 else:
53 return False
54## return (os.path.isdir(filename)
55## and not os.path.exists(filename+'/table.f1')
56## and os.path.exists(filename+'/table.info'))
57
58def is_ms(filename):
59 """Is the given file a MeasurementSet?
60
61 Parameters:
62
63 filename: the name of the file/directory to test
64
65 """
66 if ( os.path.isdir(filename)
67 and os.path.exists(filename+'/table.info')
68 and os.path.exists(filename+'/table.dat') ):
69 f=open(filename+'/table.info')
70 l=f.readline()
71 f.close()
72 if ( l.find('Measurement Set') != -1 ):
73 return True
74 else:
75 return False
76 else:
77 return False
78
79class scantable(Scantable):
80 """\
81 The ASAP container for scans (single-dish data).
82 """
83
84 @asaplog_post_dec
85 #def __init__(self, filename, average=None, unit=None, getpt=None,
86 # antenna=None, parallactify=None):
87 def __init__(self, filename, average=None, unit=None, parallactify=None, **args):
88 """\
89 Create a scantable from a saved one or make a reference
90
91 Parameters:
92
93 filename: the name of an asap table on disk
94 or
95 the name of a rpfits/sdfits/ms file
96 (integrations within scans are auto averaged
97 and the whole file is read) or
98 [advanced] a reference to an existing scantable
99
100 average: average all integrations withinb a scan on read.
101 The default (True) is taken from .asaprc.
102
103 unit: brightness unit; must be consistent with K or Jy.
104 Over-rides the default selected by the filler
105 (input rpfits/sdfits/ms) or replaces the value
106 in existing scantables
107
108 getpt: for MeasurementSet input data only:
109 If True, all pointing data are filled.
110 The deafult is False, which makes time to load
111 the MS data faster in some cases.
112
113 antenna: for MeasurementSet input data only:
114 Antenna selection. integer (id) or string (name or id).
115
116 parallactify: Indicate that the data had been parallatified. Default
117 is taken from rc file.
118
119 """
120 if average is None:
121 average = rcParams['scantable.autoaverage']
122 #if getpt is None:
123 # getpt = True
124 #if antenna is not None:
125 # asaplog.push("Antenna selection currently unsupported."
126 # "Using ''")
127 # asaplog.post('WARN')
128 #if antenna is None:
129 # antenna = ''
130 #elif type(antenna) == int:
131 # antenna = '%s' % antenna
132 #elif type(antenna) == list:
133 # tmpstr = ''
134 # for i in range( len(antenna) ):
135 # if type(antenna[i]) == int:
136 # tmpstr = tmpstr + ('%s,'%(antenna[i]))
137 # elif type(antenna[i]) == str:
138 # tmpstr=tmpstr+antenna[i]+','
139 # else:
140 # raise TypeError('Bad antenna selection.')
141 # antenna = tmpstr.rstrip(',')
142 parallactify = parallactify or rcParams['scantable.parallactify']
143 varlist = vars()
144 from asap._asap import stmath
145 self._math = stmath( rcParams['insitu'] )
146 if isinstance(filename, Scantable):
147 Scantable.__init__(self, filename)
148 else:
149 if isinstance(filename, str):
150 filename = os.path.expandvars(filename)
151 filename = os.path.expanduser(filename)
152 if not os.path.exists(filename):
153 s = "File '%s' not found." % (filename)
154 raise IOError(s)
155 if is_scantable(filename):
156 ondisk = rcParams['scantable.storage'] == 'disk'
157 Scantable.__init__(self, filename, ondisk)
158 if unit is not None:
159 self.set_fluxunit(unit)
160 if average:
161 self._assign( self.average_time( scanav=True ) )
162 # do not reset to the default freqframe
163 #self.set_freqframe(rcParams['scantable.freqframe'])
164 #elif os.path.isdir(filename) \
165 # and not os.path.exists(filename+'/table.f1'):
166 elif is_ms(filename):
167 # Measurement Set
168 opts={'ms': {}}
169 mskeys=['getpt','antenna']
170 for key in mskeys:
171 if key in args.keys():
172 opts['ms'][key] = args[key]
173 #self._fill([filename], unit, average, getpt, antenna)
174 self._fill([filename], unit, average, opts)
175 elif os.path.isfile(filename):
176 #self._fill([filename], unit, average, getpt, antenna)
177 self._fill([filename], unit, average)
178 else:
179 msg = "The given file '%s'is not a valid " \
180 "asap table." % (filename)
181 raise IOError(msg)
182 elif (isinstance(filename, list) or isinstance(filename, tuple)) \
183 and isinstance(filename[-1], str):
184 #self._fill(filename, unit, average, getpt, antenna)
185 self._fill(filename, unit, average)
186 self.parallactify(parallactify)
187 self._add_history("scantable", varlist)
188
189 @asaplog_post_dec
190 def save(self, name=None, format=None, overwrite=False):
191 """\
192 Store the scantable on disk. This can be an asap (aips++) Table,
193 SDFITS or MS2 format.
194
195 Parameters:
196
197 name: the name of the outputfile. For format "ASCII"
198 this is the root file name (data in 'name'.txt
199 and header in 'name'_header.txt)
200
201 format: an optional file format. Default is ASAP.
202 Allowed are:
203
204 * 'ASAP' (save as ASAP [aips++] Table),
205 * 'SDFITS' (save as SDFITS file)
206 * 'ASCII' (saves as ascii text file)
207 * 'MS2' (saves as an casacore MeasurementSet V2)
208 * 'FITS' (save as image FITS - not readable by class)
209 * 'CLASS' (save as FITS readable by CLASS)
210
211 overwrite: If the file should be overwritten if it exists.
212 The default False is to return with warning
213 without writing the output. USE WITH CARE.
214
215 Example::
216
217 scan.save('myscan.asap')
218 scan.save('myscan.sdfits', 'SDFITS')
219
220 """
221 from os import path
222 format = format or rcParams['scantable.save']
223 suffix = '.'+format.lower()
224 if name is None or name == "":
225 name = 'scantable'+suffix
226 msg = "No filename given. Using default name %s..." % name
227 asaplog.push(msg)
228 name = path.expandvars(name)
229 if path.isfile(name) or path.isdir(name):
230 if not overwrite:
231 msg = "File %s exists." % name
232 raise IOError(msg)
233 format2 = format.upper()
234 if format2 == 'ASAP':
235 self._save(name)
236 elif format2 == 'MS2':
237 msopt = {'ms': {'overwrite': overwrite } }
238 from asap._asap import mswriter
239 writer = mswriter( self )
240 writer.write( name, msopt )
241 else:
242 from asap._asap import stwriter as stw
243 writer = stw(format2)
244 writer.write(self, name)
245 return
246
247 def copy(self):
248 """Return a copy of this scantable.
249
250 *Note*:
251
252 This makes a full (deep) copy. scan2 = scan1 makes a reference.
253
254 Example::
255
256 copiedscan = scan.copy()
257
258 """
259 sd = scantable(Scantable._copy(self))
260 return sd
261
262 def drop_scan(self, scanid=None):
263 """\
264 Return a new scantable where the specified scan number(s) has(have)
265 been dropped.
266
267 Parameters:
268
269 scanid: a (list of) scan number(s)
270
271 """
272 from asap import _is_sequence_or_number as _is_valid
273 from asap import _to_list
274 from asap import unique
275 if not _is_valid(scanid):
276 raise RuntimeError( 'Please specify a scanno to drop from the scantable' )
277 scanid = _to_list(scanid)
278 allscans = unique([ self.getscan(i) for i in range(self.nrow())])
279 for sid in scanid: allscans.remove(sid)
280 if len(allscans) == 0:
281 raise ValueError("Can't remove all scans")
282 sel = selector(scans=allscans)
283 return self._select_copy(sel)
284
285 def _select_copy(self, selection):
286 orig = self.get_selection()
287 self.set_selection(orig+selection)
288 cp = self.copy()
289 self.set_selection(orig)
290 return cp
291
292 def get_scan(self, scanid=None):
293 """\
294 Return a specific scan (by scanno) or collection of scans (by
295 source name) in a new scantable.
296
297 *Note*:
298
299 See scantable.drop_scan() for the inverse operation.
300
301 Parameters:
302
303 scanid: a (list of) scanno or a source name, unix-style
304 patterns are accepted for source name matching, e.g.
305 '*_R' gets all 'ref scans
306
307 Example::
308
309 # get all scans containing the source '323p459'
310 newscan = scan.get_scan('323p459')
311 # get all 'off' scans
312 refscans = scan.get_scan('*_R')
313 # get a susbset of scans by scanno (as listed in scan.summary())
314 newscan = scan.get_scan([0, 2, 7, 10])
315
316 """
317 if scanid is None:
318 raise RuntimeError( 'Please specify a scan no or name to '
319 'retrieve from the scantable' )
320 try:
321 bsel = self.get_selection()
322 sel = selector()
323 if type(scanid) is str:
324 sel.set_name(scanid)
325 return self._select_copy(sel)
326 elif type(scanid) is int:
327 sel.set_scans([scanid])
328 return self._select_copy(sel)
329 elif type(scanid) is list:
330 sel.set_scans(scanid)
331 return self._select_copy(sel)
332 else:
333 msg = "Illegal scanid type, use 'int' or 'list' if ints."
334 raise TypeError(msg)
335 except RuntimeError:
336 raise
337
338 def __str__(self):
339 return Scantable._summary(self, True)
340
341 def summary(self, filename=None):
342 """\
343 Print a summary of the contents of this scantable.
344
345 Parameters:
346
347 filename: the name of a file to write the putput to
348 Default - no file output
349
350 """
351 info = Scantable._summary(self, True)
352 if filename is not None:
353 if filename is "":
354 filename = 'scantable_summary.txt'
355 from os.path import expandvars, isdir
356 filename = expandvars(filename)
357 if not isdir(filename):
358 data = open(filename, 'w')
359 data.write(info)
360 data.close()
361 else:
362 msg = "Illegal file name '%s'." % (filename)
363 raise IOError(msg)
364 return page(info)
365
366 def get_spectrum(self, rowno):
367 """Return the spectrum for the current row in the scantable as a list.
368
369 Parameters:
370
371 rowno: the row number to retrieve the spectrum from
372
373 """
374 return self._getspectrum(rowno)
375
376 def get_mask(self, rowno):
377 """Return the mask for the current row in the scantable as a list.
378
379 Parameters:
380
381 rowno: the row number to retrieve the mask from
382
383 """
384 return self._getmask(rowno)
385
386 def set_spectrum(self, spec, rowno):
387 """Set the spectrum for the current row in the scantable.
388
389 Parameters:
390
391 spec: the new spectrum
392
393 rowno: the row number to set the spectrum for
394
395 """
396 assert(len(spec) == self.nchan())
397 return self._setspectrum(spec, rowno)
398
399 def get_coordinate(self, rowno):
400 """Return the (spectral) coordinate for a a given 'rowno'.
401
402 *Note*:
403
404 * This coordinate is only valid until a scantable method modifies
405 the frequency axis.
406 * This coordinate does contain the original frequency set-up
407 NOT the new frame. The conversions however are done using the user
408 specified frame (e.g. LSRK/TOPO). To get the 'real' coordinate,
409 use scantable.freq_align first. Without it there is no closure,
410 i.e.::
411
412 c = myscan.get_coordinate(0)
413 c.to_frequency(c.get_reference_pixel()) != c.get_reference_value()
414
415 Parameters:
416
417 rowno: the row number for the spectral coordinate
418
419 """
420 return coordinate(Scantable.get_coordinate(self, rowno))
421
422 def get_selection(self):
423 """\
424 Get the selection object currently set on this scantable.
425
426 Example::
427
428 sel = scan.get_selection()
429 sel.set_ifs(0) # select IF 0
430 scan.set_selection(sel) # apply modified selection
431
432 """
433 return selector(self._getselection())
434
435 def set_selection(self, selection=None, **kw):
436 """\
437 Select a subset of the data. All following operations on this scantable
438 are only applied to thi selection.
439
440 Parameters:
441
442 selection: a selector object (default unset the selection), or
443 any combination of "pols", "ifs", "beams", "scans",
444 "cycles", "name", "query"
445
446 Examples::
447
448 sel = selector() # create a selection object
449 self.set_scans([0, 3]) # select SCANNO 0 and 3
450 scan.set_selection(sel) # set the selection
451 scan.summary() # will only print summary of scanno 0 an 3
452 scan.set_selection() # unset the selection
453 # or the equivalent
454 scan.set_selection(scans=[0,3])
455 scan.summary() # will only print summary of scanno 0 an 3
456 scan.set_selection() # unset the selection
457
458 """
459 if selection is None:
460 # reset
461 if len(kw) == 0:
462 selection = selector()
463 else:
464 # try keywords
465 for k in kw:
466 if k not in selector.fields:
467 raise KeyError("Invalid selection key '%s', valid keys are %s" % (k, selector.fields))
468 selection = selector(**kw)
469 self._setselection(selection)
470
471 def get_row(self, row=0, insitu=None):
472 """\
473 Select a row in the scantable.
474 Return a scantable with single row.
475
476 Parameters:
477
478 row: row no of integration, default is 0.
479 insitu: if False a new scantable is returned. Otherwise, the
480 scaling is done in-situ. The default is taken from .asaprc
481 (False)
482
483 """
484 if insitu is None: insitu = rcParams['insitu']
485 if not insitu:
486 workscan = self.copy()
487 else:
488 workscan = self
489 # Select a row
490 sel=selector()
491 sel.set_rows([row])
492 #sel.set_scans([workscan.getscan(row)])
493 #sel.set_cycles([workscan.getcycle(row)])
494 #sel.set_beams([workscan.getbeam(row)])
495 #sel.set_ifs([workscan.getif(row)])
496 #sel.set_polarisations([workscan.getpol(row)])
497 #sel.set_name(workscan._getsourcename(row))
498 workscan.set_selection(sel)
499 if not workscan.nrow() == 1:
500 msg = "Cloud not identify single row. %d rows selected."%(workscan.nrow())
501 raise RuntimeError(msg)
502 del sel
503 if insitu:
504 self._assign(workscan)
505 else:
506 return workscan
507
508 @asaplog_post_dec
509 def stats(self, stat='stddev', mask=None, form='3.3f', row=None):
510 """\
511 Determine the specified statistic of the current beam/if/pol
512 Takes a 'mask' as an optional parameter to specify which
513 channels should be excluded.
514
515 Parameters:
516
517 stat: 'min', 'max', 'min_abc', 'max_abc', 'sumsq', 'sum',
518 'mean', 'var', 'stddev', 'avdev', 'rms', 'median'
519
520 mask: an optional mask specifying where the statistic
521 should be determined.
522
523 form: format string to print statistic values
524
525 row: row number of spectrum to process.
526 (default is None: for all rows)
527
528 Example:
529 scan.set_unit('channel')
530 msk = scan.create_mask([100, 200], [500, 600])
531 scan.stats(stat='mean', mask=m)
532
533 """
534 mask = mask or []
535 if not self._check_ifs():
536 raise ValueError("Cannot apply mask as the IFs have different "
537 "number of channels. Please use setselection() "
538 "to select individual IFs")
539 rtnabc = False
540 if stat.lower().endswith('_abc'): rtnabc = True
541 getchan = False
542 if stat.lower().startswith('min') or stat.lower().startswith('max'):
543 chan = self._math._minmaxchan(self, mask, stat)
544 getchan = True
545 statvals = []
546 if not rtnabc:
547 if row == None:
548 statvals = self._math._stats(self, mask, stat)
549 else:
550 statvals = self._math._statsrow(self, mask, stat, int(row))
551
552 #def cb(i):
553 # return statvals[i]
554
555 #return self._row_callback(cb, stat)
556
557 label=stat
558 #callback=cb
559 out = ""
560 #outvec = []
561 sep = '-'*50
562
563 if row == None:
564 rows = xrange(self.nrow())
565 elif isinstance(row, int):
566 rows = [ row ]
567
568 for i in rows:
569 refstr = ''
570 statunit= ''
571 if getchan:
572 qx, qy = self.chan2data(rowno=i, chan=chan[i])
573 if rtnabc:
574 statvals.append(qx['value'])
575 refstr = ('(value: %'+form) % (qy['value'])+' ['+qy['unit']+'])'
576 statunit= '['+qx['unit']+']'
577 else:
578 refstr = ('(@ %'+form) % (qx['value'])+' ['+qx['unit']+'])'
579
580 tm = self._gettime(i)
581 src = self._getsourcename(i)
582 out += 'Scan[%d] (%s) ' % (self.getscan(i), src)
583 out += 'Time[%s]:\n' % (tm)
584 if self.nbeam(-1) > 1: out += ' Beam[%d] ' % (self.getbeam(i))
585 if self.nif(-1) > 1: out += ' IF[%d] ' % (self.getif(i))
586 if self.npol(-1) > 1: out += ' Pol[%d] ' % (self.getpol(i))
587 #outvec.append(callback(i))
588 if len(rows) > 1:
589 # out += ('= %'+form) % (outvec[i]) +' '+refstr+'\n'
590 out += ('= %'+form) % (statvals[i]) +' '+refstr+'\n'
591 else:
592 # out += ('= %'+form) % (outvec[0]) +' '+refstr+'\n'
593 out += ('= %'+form) % (statvals[0]) +' '+refstr+'\n'
594 out += sep+"\n"
595
596 import os
597 if os.environ.has_key( 'USER' ):
598 usr = os.environ['USER']
599 else:
600 import commands
601 usr = commands.getoutput( 'whoami' )
602 tmpfile = '/tmp/tmp_'+usr+'_casapy_asap_scantable_stats'
603 f = open(tmpfile,'w')
604 print >> f, sep
605 print >> f, ' %s %s' % (label, statunit)
606 print >> f, sep
607 print >> f, out
608 f.close()
609 f = open(tmpfile,'r')
610 x = f.readlines()
611 f.close()
612 asaplog.push(''.join(x), False)
613
614 return statvals
615
616 def chan2data(self, rowno=0, chan=0):
617 """\
618 Returns channel/frequency/velocity and spectral value
619 at an arbitrary row and channel in the scantable.
620
621 Parameters:
622
623 rowno: a row number in the scantable. Default is the
624 first row, i.e. rowno=0
625
626 chan: a channel in the scantable. Default is the first
627 channel, i.e. pos=0
628
629 """
630 if isinstance(rowno, int) and isinstance(chan, int):
631 qx = {'unit': self.get_unit(),
632 'value': self._getabcissa(rowno)[chan]}
633 qy = {'unit': self.get_fluxunit(),
634 'value': self._getspectrum(rowno)[chan]}
635 return qx, qy
636
637 def stddev(self, mask=None):
638 """\
639 Determine the standard deviation of the current beam/if/pol
640 Takes a 'mask' as an optional parameter to specify which
641 channels should be excluded.
642
643 Parameters:
644
645 mask: an optional mask specifying where the standard
646 deviation should be determined.
647
648 Example::
649
650 scan.set_unit('channel')
651 msk = scan.create_mask([100, 200], [500, 600])
652 scan.stddev(mask=m)
653
654 """
655 return self.stats(stat='stddev', mask=mask);
656
657
658 def get_column_names(self):
659 """\
660 Return a list of column names, which can be used for selection.
661 """
662 return list(Scantable.get_column_names(self))
663
664 def get_tsys(self, row=-1):
665 """\
666 Return the System temperatures.
667
668 Parameters:
669
670 row: the rowno to get the information for. (default all rows)
671
672 Returns:
673
674 a list of Tsys values for the current selection
675
676 """
677 if row > -1:
678 return self._get_column(self._gettsys, row)
679 return self._row_callback(self._gettsys, "Tsys")
680
681
682 def get_weather(self, row=-1):
683 """\
684 Return the weather informations.
685
686 Parameters:
687
688 row: the rowno to get the information for. (default all rows)
689
690 Returns:
691
692 a dict or list of of dicts of values for the current selection
693
694 """
695
696 values = self._get_column(self._get_weather, row)
697 if row > -1:
698 return {'temperature': values[0],
699 'pressure': values[1], 'humidity' : values[2],
700 'windspeed' : values[3], 'windaz' : values[4]
701 }
702 else:
703 out = []
704 for r in values:
705
706 out.append({'temperature': r[0],
707 'pressure': r[1], 'humidity' : r[2],
708 'windspeed' : r[3], 'windaz' : r[4]
709 })
710 return out
711
712 def _row_callback(self, callback, label):
713 out = ""
714 outvec = []
715 sep = '-'*50
716 for i in range(self.nrow()):
717 tm = self._gettime(i)
718 src = self._getsourcename(i)
719 out += 'Scan[%d] (%s) ' % (self.getscan(i), src)
720 out += 'Time[%s]:\n' % (tm)
721 if self.nbeam(-1) > 1:
722 out += ' Beam[%d] ' % (self.getbeam(i))
723 if self.nif(-1) > 1: out += ' IF[%d] ' % (self.getif(i))
724 if self.npol(-1) > 1: out += ' Pol[%d] ' % (self.getpol(i))
725 outvec.append(callback(i))
726 out += '= %3.3f\n' % (outvec[i])
727 out += sep+'\n'
728
729 asaplog.push(sep)
730 asaplog.push(" %s" % (label))
731 asaplog.push(sep)
732 asaplog.push(out)
733 asaplog.post()
734 return outvec
735
736 def _get_column(self, callback, row=-1, *args):
737 """
738 """
739 if row == -1:
740 return [callback(i, *args) for i in range(self.nrow())]
741 else:
742 if 0 <= row < self.nrow():
743 return callback(row, *args)
744
745
746 def get_time(self, row=-1, asdatetime=False, prec=-1):
747 """\
748 Get a list of time stamps for the observations.
749 Return a datetime object or a string (default) for each
750 integration time stamp in the scantable.
751
752 Parameters:
753
754 row: row no of integration. Default -1 return all rows
755
756 asdatetime: return values as datetime objects rather than strings
757
758 prec: number of digits shown. Default -1 to automatic calculation.
759 Note this number is equals to the digits of MVTime,
760 i.e., 0<prec<3: dates with hh:: only,
761 <5: with hh:mm:, <7 or 0: with hh:mm:ss,
762 and 6> : with hh:mm:ss.tt... (prec-6 t's added)
763
764 """
765 from datetime import datetime
766 if prec < 0:
767 # automagically set necessary precision +1
768 prec = 7 - numpy.floor(numpy.log10(numpy.min(self.get_inttime(row))))
769 prec = max(6, int(prec))
770 else:
771 prec = max(0, prec)
772 if asdatetime:
773 #precision can be 1 millisecond at max
774 prec = min(12, prec)
775
776 times = self._get_column(self._gettime, row, prec)
777 if not asdatetime:
778 return times
779 format = "%Y/%m/%d/%H:%M:%S.%f"
780 if prec < 7:
781 nsub = 1 + (((6-prec)/2) % 3)
782 substr = [".%f","%S","%M"]
783 for i in range(nsub):
784 format = format.replace(substr[i],"")
785 if isinstance(times, list):
786 return [datetime.strptime(i, format) for i in times]
787 else:
788 return datetime.strptime(times, format)
789
790
791 def get_inttime(self, row=-1):
792 """\
793 Get a list of integration times for the observations.
794 Return a time in seconds for each integration in the scantable.
795
796 Parameters:
797
798 row: row no of integration. Default -1 return all rows.
799
800 """
801 return self._get_column(self._getinttime, row)
802
803
804 def get_sourcename(self, row=-1):
805 """\
806 Get a list source names for the observations.
807 Return a string for each integration in the scantable.
808 Parameters:
809
810 row: row no of integration. Default -1 return all rows.
811
812 """
813 return self._get_column(self._getsourcename, row)
814
815 def get_elevation(self, row=-1):
816 """\
817 Get a list of elevations for the observations.
818 Return a float for each integration in the scantable.
819
820 Parameters:
821
822 row: row no of integration. Default -1 return all rows.
823
824 """
825 return self._get_column(self._getelevation, row)
826
827 def get_azimuth(self, row=-1):
828 """\
829 Get a list of azimuths for the observations.
830 Return a float for each integration in the scantable.
831
832 Parameters:
833 row: row no of integration. Default -1 return all rows.
834
835 """
836 return self._get_column(self._getazimuth, row)
837
838 def get_parangle(self, row=-1):
839 """\
840 Get a list of parallactic angles for the observations.
841 Return a float for each integration in the scantable.
842
843 Parameters:
844
845 row: row no of integration. Default -1 return all rows.
846
847 """
848 return self._get_column(self._getparangle, row)
849
850 def get_direction(self, row=-1):
851 """
852 Get a list of Positions on the sky (direction) for the observations.
853 Return a string for each integration in the scantable.
854
855 Parameters:
856
857 row: row no of integration. Default -1 return all rows
858
859 """
860 return self._get_column(self._getdirection, row)
861
862 def get_directionval(self, row=-1):
863 """\
864 Get a list of Positions on the sky (direction) for the observations.
865 Return a float for each integration in the scantable.
866
867 Parameters:
868
869 row: row no of integration. Default -1 return all rows
870
871 """
872 return self._get_column(self._getdirectionvec, row)
873
874 @asaplog_post_dec
875 def set_unit(self, unit='channel'):
876 """\
877 Set the unit for all following operations on this scantable
878
879 Parameters:
880
881 unit: optional unit, default is 'channel'. Use one of '*Hz',
882 'km/s', 'channel' or equivalent ''
883
884 """
885 varlist = vars()
886 if unit in ['', 'pixel', 'channel']:
887 unit = ''
888 inf = list(self._getcoordinfo())
889 inf[0] = unit
890 self._setcoordinfo(inf)
891 self._add_history("set_unit", varlist)
892
893 @asaplog_post_dec
894 def set_instrument(self, instr):
895 """\
896 Set the instrument for subsequent processing.
897
898 Parameters:
899
900 instr: Select from 'ATPKSMB', 'ATPKSHOH', 'ATMOPRA',
901 'DSS-43' (Tid), 'CEDUNA', and 'HOBART'
902
903 """
904 self._setInstrument(instr)
905 self._add_history("set_instument", vars())
906
907 @asaplog_post_dec
908 def set_feedtype(self, feedtype):
909 """\
910 Overwrite the feed type, which might not be set correctly.
911
912 Parameters:
913
914 feedtype: 'linear' or 'circular'
915
916 """
917 self._setfeedtype(feedtype)
918 self._add_history("set_feedtype", vars())
919
920 @asaplog_post_dec
921 def set_doppler(self, doppler='RADIO'):
922 """\
923 Set the doppler for all following operations on this scantable.
924
925 Parameters:
926
927 doppler: One of 'RADIO', 'OPTICAL', 'Z', 'BETA', 'GAMMA'
928
929 """
930 varlist = vars()
931 inf = list(self._getcoordinfo())
932 inf[2] = doppler
933 self._setcoordinfo(inf)
934 self._add_history("set_doppler", vars())
935
936 @asaplog_post_dec
937 def set_freqframe(self, frame=None):
938 """\
939 Set the frame type of the Spectral Axis.
940
941 Parameters:
942
943 frame: an optional frame type, default 'LSRK'. Valid frames are:
944 'TOPO', 'LSRD', 'LSRK', 'BARY',
945 'GEO', 'GALACTO', 'LGROUP', 'CMB'
946
947 Example::
948
949 scan.set_freqframe('BARY')
950
951 """
952 frame = frame or rcParams['scantable.freqframe']
953 varlist = vars()
954 # "REST" is not implemented in casacore
955 #valid = ['REST', 'TOPO', 'LSRD', 'LSRK', 'BARY', \
956 # 'GEO', 'GALACTO', 'LGROUP', 'CMB']
957 valid = ['TOPO', 'LSRD', 'LSRK', 'BARY', \
958 'GEO', 'GALACTO', 'LGROUP', 'CMB']
959
960 if frame in valid:
961 inf = list(self._getcoordinfo())
962 inf[1] = frame
963 self._setcoordinfo(inf)
964 self._add_history("set_freqframe", varlist)
965 else:
966 msg = "Please specify a valid freq type. Valid types are:\n", valid
967 raise TypeError(msg)
968
969 @asaplog_post_dec
970 def set_dirframe(self, frame=""):
971 """\
972 Set the frame type of the Direction on the sky.
973
974 Parameters:
975
976 frame: an optional frame type, default ''. Valid frames are:
977 'J2000', 'B1950', 'GALACTIC'
978
979 Example:
980
981 scan.set_dirframe('GALACTIC')
982
983 """
984 varlist = vars()
985 Scantable.set_dirframe(self, frame)
986 self._add_history("set_dirframe", varlist)
987
988 def get_unit(self):
989 """\
990 Get the default unit set in this scantable
991
992 Returns:
993
994 A unit string
995
996 """
997 inf = self._getcoordinfo()
998 unit = inf[0]
999 if unit == '': unit = 'channel'
1000 return unit
1001
1002 @asaplog_post_dec
1003 def get_abcissa(self, rowno=0):
1004 """\
1005 Get the abcissa in the current coordinate setup for the currently
1006 selected Beam/IF/Pol
1007
1008 Parameters:
1009
1010 rowno: an optional row number in the scantable. Default is the
1011 first row, i.e. rowno=0
1012
1013 Returns:
1014
1015 The abcissa values and the format string (as a dictionary)
1016
1017 """
1018 abc = self._getabcissa(rowno)
1019 lbl = self._getabcissalabel(rowno)
1020 return abc, lbl
1021
1022 @asaplog_post_dec
1023 def flag(self, row=-1, mask=None, unflag=False):
1024 """\
1025 Flag the selected data using an optional channel mask.
1026
1027 Parameters:
1028
1029 row: an optional row number in the scantable.
1030 Default -1 flags all rows
1031
1032 mask: an optional channel mask, created with create_mask. Default
1033 (no mask) is all channels.
1034
1035 unflag: if True, unflag the data
1036
1037 """
1038 varlist = vars()
1039 mask = mask or []
1040 self._flag(row, mask, unflag)
1041 self._add_history("flag", varlist)
1042
1043 @asaplog_post_dec
1044 def flag_row(self, rows=[], unflag=False):
1045 """\
1046 Flag the selected data in row-based manner.
1047
1048 Parameters:
1049
1050 rows: list of row numbers to be flagged. Default is no row
1051 (must be explicitly specified to execute row-based flagging).
1052
1053 unflag: if True, unflag the data.
1054
1055 """
1056 varlist = vars()
1057 self._flag_row(rows, unflag)
1058 self._add_history("flag_row", varlist)
1059
1060 @asaplog_post_dec
1061 def clip(self, uthres=None, dthres=None, clipoutside=True, unflag=False):
1062 """\
1063 Flag the selected data outside a specified range (in channel-base)
1064
1065 Parameters:
1066
1067 uthres: upper threshold.
1068
1069 dthres: lower threshold
1070
1071 clipoutside: True for flagging data outside the range [dthres:uthres].
1072 False for flagging data inside the range.
1073
1074 unflag: if True, unflag the data.
1075
1076 """
1077 varlist = vars()
1078 self._clip(uthres, dthres, clipoutside, unflag)
1079 self._add_history("clip", varlist)
1080
1081 @asaplog_post_dec
1082 def lag_flag(self, start, end, unit="MHz", insitu=None):
1083 """\
1084 Flag the data in 'lag' space by providing a frequency to remove.
1085 Flagged data in the scantable gets interpolated over the region.
1086 No taper is applied.
1087
1088 Parameters:
1089
1090 start: the start frequency (really a period within the
1091 bandwidth) or period to remove
1092
1093 end: the end frequency or period to remove
1094
1095 unit: the frequency unit (default "MHz") or "" for
1096 explicit lag channels
1097
1098 *Notes*:
1099
1100 It is recommended to flag edges of the band or strong
1101 signals beforehand.
1102
1103 """
1104 if insitu is None: insitu = rcParams['insitu']
1105 self._math._setinsitu(insitu)
1106 varlist = vars()
1107 base = { "GHz": 1000000000., "MHz": 1000000., "kHz": 1000., "Hz": 1.}
1108 if not (unit == "" or base.has_key(unit)):
1109 raise ValueError("%s is not a valid unit." % unit)
1110 if unit == "":
1111 s = scantable(self._math._lag_flag(self, start, end, "lags"))
1112 else:
1113 s = scantable(self._math._lag_flag(self, start*base[unit],
1114 end*base[unit], "frequency"))
1115 s._add_history("lag_flag", varlist)
1116 if insitu:
1117 self._assign(s)
1118 else:
1119 return s
1120
1121 @asaplog_post_dec
1122 def create_mask(self, *args, **kwargs):
1123 """\
1124 Compute and return a mask based on [min, max] windows.
1125 The specified windows are to be INCLUDED, when the mask is
1126 applied.
1127
1128 Parameters:
1129
1130 [min, max], [min2, max2], ...
1131 Pairs of start/end points (inclusive)specifying the regions
1132 to be masked
1133
1134 invert: optional argument. If specified as True,
1135 return an inverted mask, i.e. the regions
1136 specified are EXCLUDED
1137
1138 row: create the mask using the specified row for
1139 unit conversions, default is row=0
1140 only necessary if frequency varies over rows.
1141
1142 Examples::
1143
1144 scan.set_unit('channel')
1145 # a)
1146 msk = scan.create_mask([400, 500], [800, 900])
1147 # masks everything outside 400 and 500
1148 # and 800 and 900 in the unit 'channel'
1149
1150 # b)
1151 msk = scan.create_mask([400, 500], [800, 900], invert=True)
1152 # masks the regions between 400 and 500
1153 # and 800 and 900 in the unit 'channel'
1154
1155 # c)
1156 #mask only channel 400
1157 msk = scan.create_mask([400])
1158
1159 """
1160 row = kwargs.get("row", 0)
1161 data = self._getabcissa(row)
1162 u = self._getcoordinfo()[0]
1163 if u == "":
1164 u = "channel"
1165 msg = "The current mask window unit is %s" % u
1166 i = self._check_ifs()
1167 if not i:
1168 msg += "\nThis mask is only valid for IF=%d" % (self.getif(i))
1169 asaplog.push(msg)
1170 n = self.nchan()
1171 msk = _n_bools(n, False)
1172 # test if args is a 'list' or a 'normal *args - UGLY!!!
1173
1174 ws = (isinstance(args[-1][-1], int) or isinstance(args[-1][-1], float)) \
1175 and args or args[0]
1176 for window in ws:
1177 if len(window) == 1:
1178 window = [window[0], window[0]]
1179 if len(window) == 0 or len(window) > 2:
1180 raise ValueError("A window needs to be defined as [start(, end)]")
1181 if window[0] > window[1]:
1182 tmp = window[0]
1183 window[0] = window[1]
1184 window[1] = tmp
1185 for i in range(n):
1186 if data[i] >= window[0] and data[i] <= window[1]:
1187 msk[i] = True
1188 if kwargs.has_key('invert'):
1189 if kwargs.get('invert'):
1190 msk = mask_not(msk)
1191 return msk
1192
1193 def get_masklist(self, mask=None, row=0, silent=False):
1194 """\
1195 Compute and return a list of mask windows, [min, max].
1196
1197 Parameters:
1198
1199 mask: channel mask, created with create_mask.
1200
1201 row: calcutate the masklist using the specified row
1202 for unit conversions, default is row=0
1203 only necessary if frequency varies over rows.
1204
1205 Returns:
1206
1207 [min, max], [min2, max2], ...
1208 Pairs of start/end points (inclusive)specifying
1209 the masked regions
1210
1211 """
1212 if not (isinstance(mask,list) or isinstance(mask, tuple)):
1213 raise TypeError("The mask should be list or tuple.")
1214 if len(mask) < 2:
1215 raise TypeError("The mask elements should be > 1")
1216 if self.nchan() != len(mask):
1217 msg = "Number of channels in scantable != number of mask elements"
1218 raise TypeError(msg)
1219 data = self._getabcissa(row)
1220 u = self._getcoordinfo()[0]
1221 if u == "":
1222 u = "channel"
1223 msg = "The current mask window unit is %s" % u
1224 i = self._check_ifs()
1225 if not i:
1226 msg += "\nThis mask is only valid for IF=%d" % (self.getif(i))
1227 if not silent:
1228 asaplog.push(msg)
1229 masklist=[]
1230 ist, ien = None, None
1231 ist, ien=self.get_mask_indices(mask)
1232 if ist is not None and ien is not None:
1233 for i in xrange(len(ist)):
1234 range=[data[ist[i]],data[ien[i]]]
1235 range.sort()
1236 masklist.append([range[0],range[1]])
1237 return masklist
1238
1239 def get_mask_indices(self, mask=None):
1240 """\
1241 Compute and Return lists of mask start indices and mask end indices.
1242
1243 Parameters:
1244
1245 mask: channel mask, created with create_mask.
1246
1247 Returns:
1248
1249 List of mask start indices and that of mask end indices,
1250 i.e., [istart1,istart2,....], [iend1,iend2,....].
1251
1252 """
1253 if not (isinstance(mask,list) or isinstance(mask, tuple)):
1254 raise TypeError("The mask should be list or tuple.")
1255 if len(mask) < 2:
1256 raise TypeError("The mask elements should be > 1")
1257 istart=[]
1258 iend=[]
1259 if mask[0]: istart.append(0)
1260 for i in range(len(mask)-1):
1261 if not mask[i] and mask[i+1]:
1262 istart.append(i+1)
1263 elif mask[i] and not mask[i+1]:
1264 iend.append(i)
1265 if mask[len(mask)-1]: iend.append(len(mask)-1)
1266 if len(istart) != len(iend):
1267 raise RuntimeError("Numbers of mask start != mask end.")
1268 for i in range(len(istart)):
1269 if istart[i] > iend[i]:
1270 raise RuntimeError("Mask start index > mask end index")
1271 break
1272 return istart,iend
1273
1274 @asaplog_post_dec
1275 def parse_maskexpr(self,maskstring):
1276 """
1277 Parse CASA type mask selection syntax (IF dependent).
1278
1279 Parameters:
1280 maskstring : A string mask selection expression.
1281 A comma separated selections mean different IF -
1282 channel combinations. IFs and channel selections
1283 are partitioned by a colon, ':'.
1284 examples:
1285 '' = all IFs (all channels)
1286 '<2,4~6,9' = IFs 0,1,4,5,6,9 (all channels)
1287 '3:3~45;60' = channels 3 to 45 and 60 in IF 3
1288 '0~1:2~6,8' = channels 2 to 6 in IFs 0,1, and
1289 all channels in IF8
1290 Returns:
1291 A dictionary of selected (valid) IF and masklist pairs,
1292 e.g. {'0': [[50,250],[350,462]], '2': [[100,400],[550,974]]}
1293 """
1294 if not isinstance(maskstring,str):
1295 asaplog.post()
1296 asaplog.push("Invalid mask expression")
1297 asaplog.post("ERROR")
1298
1299 valid_ifs = self.getifnos()
1300 frequnit = self.get_unit()
1301 seldict = {}
1302 if maskstring == "":
1303 maskstring = str(valid_ifs)[1:-1]
1304 ## split each selection
1305 sellist = maskstring.split(',')
1306 for currselstr in sellist:
1307 selset = currselstr.split(':')
1308 # spw and mask string (may include ~, < or >)
1309 spwmasklist = self._parse_selection(selset[0],typestr='integer',
1310 offset=1,minval=min(valid_ifs),
1311 maxval=max(valid_ifs))
1312 for spwlist in spwmasklist:
1313 selspws = []
1314 for ispw in range(spwlist[0],spwlist[1]+1):
1315 # Put into the list only if ispw exists
1316 if valid_ifs.count(ispw):
1317 selspws.append(ispw)
1318 del spwmasklist, spwlist
1319
1320 # parse frequency mask list
1321 if len(selset) > 1:
1322 freqmasklist = self._parse_selection(selset[1],typestr='float',
1323 offset=0.)
1324 else:
1325 # want to select the whole spectrum
1326 freqmasklist = [None]
1327
1328 ## define a dictionary of spw - masklist combination
1329 for ispw in selspws:
1330 #print "working on", ispw
1331 spwstr = str(ispw)
1332 if len(selspws) == 0:
1333 # empty spw
1334 continue
1335 else:
1336 ## want to get min and max of the spw and
1337 ## offset to set for '<' and '>'
1338 if frequnit == 'channel':
1339 minfreq = 0
1340 maxfreq = self.nchan(ifno=ispw)
1341 offset = 0.5
1342 else:
1343 ## This is ugly part. need improvement
1344 for ifrow in xrange(self.nrow()):
1345 if self.getif(ifrow) == ispw:
1346 #print "IF",ispw,"found in row =",ifrow
1347 break
1348 freqcoord = self.get_coordinate(ifrow)
1349 freqs = self._getabcissa(ifrow)
1350 minfreq = min(freqs)
1351 maxfreq = max(freqs)
1352 if len(freqs) == 1:
1353 offset = 0.5
1354 elif frequnit.find('Hz') > 0:
1355 offset = abs(freqcoord.to_frequency(1,unit=frequnit)
1356 -freqcoord.to_frequency(0,unit=frequnit))*0.5
1357 elif frequnit.find('m/s') > 0:
1358 offset = abs(freqcoord.to_velocity(1,unit=frequnit)
1359 -freqcoord.to_velocity(0,unit=frequnit))*0.5
1360 else:
1361 asaplog.post()
1362 asaplog.push("Invalid frequency unit")
1363 asaplog.post("ERROR")
1364 del freqs, freqcoord, ifrow
1365 for freq in freqmasklist:
1366 selmask = freq or [minfreq, maxfreq]
1367 if selmask[0] == None:
1368 ## selection was "<freq[1]".
1369 if selmask[1] < minfreq:
1370 ## avoid adding region selection
1371 selmask = None
1372 else:
1373 selmask = [minfreq,selmask[1]-offset]
1374 elif selmask[1] == None:
1375 ## selection was ">freq[0]"
1376 if selmask[0] > maxfreq:
1377 ## avoid adding region selection
1378 selmask = None
1379 else:
1380 selmask = [selmask[0]+offset,maxfreq]
1381 if selmask:
1382 if not seldict.has_key(spwstr):
1383 # new spw selection
1384 seldict[spwstr] = []
1385 seldict[spwstr] += [selmask]
1386 del minfreq,maxfreq,offset,freq,selmask
1387 del spwstr
1388 del freqmasklist
1389 del valid_ifs
1390 if len(seldict) == 0:
1391 asaplog.post()
1392 asaplog.push("No valid selection in the mask expression: "+maskstring)
1393 asaplog.post("WARN")
1394 return None
1395 msg = "Selected masklist:\n"
1396 for sif, lmask in seldict.iteritems():
1397 msg += " IF"+sif+" - "+str(lmask)+"\n"
1398 asaplog.push(msg)
1399 return seldict
1400
1401 def _parse_selection(self,selstr,typestr='float',offset=0.,minval=None,maxval=None):
1402 """
1403 Parameters:
1404 selstr : The Selection string, e.g., '<3;5~7;100~103;9'
1405 typestr : The type of the values in returned list
1406 ('integer' or 'float')
1407 offset : The offset value to subtract from or add to
1408 the boundary value if the selection string
1409 includes '<' or '>'
1410 minval, maxval : The minimum/maximum values to set if the
1411 selection string includes '<' or '>'.
1412 The list element is filled with None by default.
1413 Returns:
1414 A list of min/max pair of selections.
1415 Example:
1416 _parseSelection('<3;5~7;9',typestr='int',offset=1,minval=0)
1417 returns [[0,2],[5,7],[9,9]]
1418 """
1419 selgroups = selstr.split(';')
1420 sellists = []
1421 if typestr.lower().startswith('int'):
1422 formatfunc = int
1423 else:
1424 formatfunc = float
1425
1426 for currsel in selgroups:
1427 if currsel.find('~') > 0:
1428 minsel = formatfunc(currsel.split('~')[0].strip())
1429 maxsel = formatfunc(currsel.split('~')[1].strip())
1430 elif currsel.strip().startswith('<'):
1431 minsel = minval
1432 maxsel = formatfunc(currsel.split('<')[1].strip()) \
1433 - formatfunc(offset)
1434 elif currsel.strip().startswith('>'):
1435 minsel = formatfunc(currsel.split('>')[1].strip()) \
1436 + formatfunc(offset)
1437 maxsel = maxval
1438 else:
1439 minsel = formatfunc(currsel)
1440 maxsel = formatfunc(currsel)
1441 sellists.append([minsel,maxsel])
1442 return sellists
1443
1444# def get_restfreqs(self):
1445# """
1446# Get the restfrequency(s) stored in this scantable.
1447# The return value(s) are always of unit 'Hz'
1448# Parameters:
1449# none
1450# Returns:
1451# a list of doubles
1452# """
1453# return list(self._getrestfreqs())
1454
1455 def get_restfreqs(self, ids=None):
1456 """\
1457 Get the restfrequency(s) stored in this scantable.
1458 The return value(s) are always of unit 'Hz'
1459
1460 Parameters:
1461
1462 ids: (optional) a list of MOLECULE_ID for that restfrequency(s) to
1463 be retrieved
1464
1465 Returns:
1466
1467 dictionary containing ids and a list of doubles for each id
1468
1469 """
1470 if ids is None:
1471 rfreqs={}
1472 idlist = self.getmolnos()
1473 for i in idlist:
1474 rfreqs[i]=list(self._getrestfreqs(i))
1475 return rfreqs
1476 else:
1477 if type(ids)==list or type(ids)==tuple:
1478 rfreqs={}
1479 for i in ids:
1480 rfreqs[i]=list(self._getrestfreqs(i))
1481 return rfreqs
1482 else:
1483 return list(self._getrestfreqs(ids))
1484 #return list(self._getrestfreqs(ids))
1485
1486 def set_restfreqs(self, freqs=None, unit='Hz'):
1487 """\
1488 Set or replace the restfrequency specified and
1489 if the 'freqs' argument holds a scalar,
1490 then that rest frequency will be applied to all the selected
1491 data. If the 'freqs' argument holds
1492 a vector, then it MUST be of equal or smaller length than
1493 the number of IFs (and the available restfrequencies will be
1494 replaced by this vector). In this case, *all* data have
1495 the restfrequency set per IF according
1496 to the corresponding value you give in the 'freqs' vector.
1497 E.g. 'freqs=[1e9, 2e9]' would mean IF 0 gets restfreq 1e9 and
1498 IF 1 gets restfreq 2e9.
1499
1500 You can also specify the frequencies via a linecatalog.
1501
1502 Parameters:
1503
1504 freqs: list of rest frequency values or string idenitfiers
1505
1506 unit: unit for rest frequency (default 'Hz')
1507
1508
1509 Example::
1510
1511 # set the given restfrequency for the all currently selected IFs
1512 scan.set_restfreqs(freqs=1.4e9)
1513 # set restfrequencies for the n IFs (n > 1) in the order of the
1514 # list, i.e
1515 # IF0 -> 1.4e9, IF1 -> 1.41e9, IF3 -> 1.42e9
1516 # len(list_of_restfreqs) == nIF
1517 # for nIF == 1 the following will set multiple restfrequency for
1518 # that IF
1519 scan.set_restfreqs(freqs=[1.4e9, 1.41e9, 1.42e9])
1520 # set multiple restfrequencies per IF. as a list of lists where
1521 # the outer list has nIF elements, the inner s arbitrary
1522 scan.set_restfreqs(freqs=[[1.4e9, 1.41e9], [1.67e9]])
1523
1524 *Note*:
1525
1526 To do more sophisticate Restfrequency setting, e.g. on a
1527 source and IF basis, use scantable.set_selection() before using
1528 this function::
1529
1530 # provided your scantable is called scan
1531 selection = selector()
1532 selection.set_name("ORION*")
1533 selection.set_ifs([1])
1534 scan.set_selection(selection)
1535 scan.set_restfreqs(freqs=86.6e9)
1536
1537 """
1538 varlist = vars()
1539 from asap import linecatalog
1540 # simple value
1541 if isinstance(freqs, int) or isinstance(freqs, float):
1542 self._setrestfreqs([freqs], [""], unit)
1543 # list of values
1544 elif isinstance(freqs, list) or isinstance(freqs, tuple):
1545 # list values are scalars
1546 if isinstance(freqs[-1], int) or isinstance(freqs[-1], float):
1547 if len(freqs) == 1:
1548 self._setrestfreqs(freqs, [""], unit)
1549 else:
1550 # allow the 'old' mode of setting mulitple IFs
1551 sel = selector()
1552 savesel = self._getselection()
1553 iflist = self.getifnos()
1554 if len(freqs)>len(iflist):
1555 raise ValueError("number of elements in list of list "
1556 "exeeds the current IF selections")
1557 iflist = self.getifnos()
1558 for i, fval in enumerate(freqs):
1559 sel.set_ifs(iflist[i])
1560 self._setselection(sel)
1561 self._setrestfreqs([fval], [""], unit)
1562 self._setselection(savesel)
1563
1564 # list values are dict, {'value'=, 'name'=)
1565 elif isinstance(freqs[-1], dict):
1566 values = []
1567 names = []
1568 for d in freqs:
1569 values.append(d["value"])
1570 names.append(d["name"])
1571 self._setrestfreqs(values, names, unit)
1572 elif isinstance(freqs[-1], list) or isinstance(freqs[-1], tuple):
1573 sel = selector()
1574 savesel = self._getselection()
1575 iflist = self.getifnos()
1576 if len(freqs)>len(iflist):
1577 raise ValueError("number of elements in list of list exeeds"
1578 " the current IF selections")
1579 for i, fval in enumerate(freqs):
1580 sel.set_ifs(iflist[i])
1581 self._setselection(sel)
1582 self._setrestfreqs(fval, [""], unit)
1583 self._setselection(savesel)
1584 # freqs are to be taken from a linecatalog
1585 elif isinstance(freqs, linecatalog):
1586 sel = selector()
1587 savesel = self._getselection()
1588 for i in xrange(freqs.nrow()):
1589 sel.set_ifs(iflist[i])
1590 self._setselection(sel)
1591 self._setrestfreqs([freqs.get_frequency(i)],
1592 [freqs.get_name(i)], "MHz")
1593 # ensure that we are not iterating past nIF
1594 if i == self.nif()-1: break
1595 self._setselection(savesel)
1596 else:
1597 return
1598 self._add_history("set_restfreqs", varlist)
1599
1600 def shift_refpix(self, delta):
1601 """\
1602 Shift the reference pixel of the Spectra Coordinate by an
1603 integer amount.
1604
1605 Parameters:
1606
1607 delta: the amount to shift by
1608
1609 *Note*:
1610
1611 Be careful using this with broadband data.
1612
1613 """
1614 Scantable.shift_refpix(self, delta)
1615
1616 @asaplog_post_dec
1617 def history(self, filename=None):
1618 """\
1619 Print the history. Optionally to a file.
1620
1621 Parameters:
1622
1623 filename: The name of the file to save the history to.
1624
1625 """
1626 hist = list(self._gethistory())
1627 out = "-"*80
1628 for h in hist:
1629 if h.startswith("---"):
1630 out = "\n".join([out, h])
1631 else:
1632 items = h.split("##")
1633 date = items[0]
1634 func = items[1]
1635 items = items[2:]
1636 out += "\n"+date+"\n"
1637 out += "Function: %s\n Parameters:" % (func)
1638 for i in items:
1639 if i == '':
1640 continue
1641 s = i.split("=")
1642 out += "\n %s = %s" % (s[0], s[1])
1643 out = "\n".join([out, "-"*80])
1644 if filename is not None:
1645 if filename is "":
1646 filename = 'scantable_history.txt'
1647 import os
1648 filename = os.path.expandvars(os.path.expanduser(filename))
1649 if not os.path.isdir(filename):
1650 data = open(filename, 'w')
1651 data.write(out)
1652 data.close()
1653 else:
1654 msg = "Illegal file name '%s'." % (filename)
1655 raise IOError(msg)
1656 return page(out)
1657 #
1658 # Maths business
1659 #
1660 @asaplog_post_dec
1661 def average_time(self, mask=None, scanav=False, weight='tint', align=False):
1662 """\
1663 Return the (time) weighted average of a scan.
1664
1665 *Note*:
1666
1667 in channels only - align if necessary
1668
1669 Parameters:
1670
1671 mask: an optional mask (only used for 'var' and 'tsys'
1672 weighting)
1673
1674 scanav: True averages each scan separately
1675 False (default) averages all scans together,
1676
1677 weight: Weighting scheme.
1678 'none' (mean no weight)
1679 'var' (1/var(spec) weighted)
1680 'tsys' (1/Tsys**2 weighted)
1681 'tint' (integration time weighted)
1682 'tintsys' (Tint/Tsys**2)
1683 'median' ( median averaging)
1684 The default is 'tint'
1685
1686 align: align the spectra in velocity before averaging. It takes
1687 the time of the first spectrum as reference time.
1688
1689 Example::
1690
1691 # time average the scantable without using a mask
1692 newscan = scan.average_time()
1693
1694 """
1695 varlist = vars()
1696 weight = weight or 'TINT'
1697 mask = mask or ()
1698 scanav = (scanav and 'SCAN') or 'NONE'
1699 scan = (self, )
1700
1701 if align:
1702 scan = (self.freq_align(insitu=False), )
1703 s = None
1704 if weight.upper() == 'MEDIAN':
1705 s = scantable(self._math._averagechannel(scan[0], 'MEDIAN',
1706 scanav))
1707 else:
1708 s = scantable(self._math._average(scan, mask, weight.upper(),
1709 scanav))
1710 s._add_history("average_time", varlist)
1711 return s
1712
1713 @asaplog_post_dec
1714 def convert_flux(self, jyperk=None, eta=None, d=None, insitu=None):
1715 """\
1716 Return a scan where all spectra are converted to either
1717 Jansky or Kelvin depending upon the flux units of the scan table.
1718 By default the function tries to look the values up internally.
1719 If it can't find them (or if you want to over-ride), you must
1720 specify EITHER jyperk OR eta (and D which it will try to look up
1721 also if you don't set it). jyperk takes precedence if you set both.
1722
1723 Parameters:
1724
1725 jyperk: the Jy / K conversion factor
1726
1727 eta: the aperture efficiency
1728
1729 d: the geometric diameter (metres)
1730
1731 insitu: if False a new scantable is returned.
1732 Otherwise, the scaling is done in-situ
1733 The default is taken from .asaprc (False)
1734
1735 """
1736 if insitu is None: insitu = rcParams['insitu']
1737 self._math._setinsitu(insitu)
1738 varlist = vars()
1739 jyperk = jyperk or -1.0
1740 d = d or -1.0
1741 eta = eta or -1.0
1742 s = scantable(self._math._convertflux(self, d, eta, jyperk))
1743 s._add_history("convert_flux", varlist)
1744 if insitu: self._assign(s)
1745 else: return s
1746
1747 @asaplog_post_dec
1748 def gain_el(self, poly=None, filename="", method="linear", insitu=None):
1749 """\
1750 Return a scan after applying a gain-elevation correction.
1751 The correction can be made via either a polynomial or a
1752 table-based interpolation (and extrapolation if necessary).
1753 You specify polynomial coefficients, an ascii table or neither.
1754 If you specify neither, then a polynomial correction will be made
1755 with built in coefficients known for certain telescopes (an error
1756 will occur if the instrument is not known).
1757 The data and Tsys are *divided* by the scaling factors.
1758
1759 Parameters:
1760
1761 poly: Polynomial coefficients (default None) to compute a
1762 gain-elevation correction as a function of
1763 elevation (in degrees).
1764
1765 filename: The name of an ascii file holding correction factors.
1766 The first row of the ascii file must give the column
1767 names and these MUST include columns
1768 "ELEVATION" (degrees) and "FACTOR" (multiply data
1769 by this) somewhere.
1770 The second row must give the data type of the
1771 column. Use 'R' for Real and 'I' for Integer.
1772 An example file would be
1773 (actual factors are arbitrary) :
1774
1775 TIME ELEVATION FACTOR
1776 R R R
1777 0.1 0 0.8
1778 0.2 20 0.85
1779 0.3 40 0.9
1780 0.4 60 0.85
1781 0.5 80 0.8
1782 0.6 90 0.75
1783
1784 method: Interpolation method when correcting from a table.
1785 Values are "nearest", "linear" (default), "cubic"
1786 and "spline"
1787
1788 insitu: if False a new scantable is returned.
1789 Otherwise, the scaling is done in-situ
1790 The default is taken from .asaprc (False)
1791
1792 """
1793
1794 if insitu is None: insitu = rcParams['insitu']
1795 self._math._setinsitu(insitu)
1796 varlist = vars()
1797 poly = poly or ()
1798 from os.path import expandvars
1799 filename = expandvars(filename)
1800 s = scantable(self._math._gainel(self, poly, filename, method))
1801 s._add_history("gain_el", varlist)
1802 if insitu:
1803 self._assign(s)
1804 else:
1805 return s
1806
1807 @asaplog_post_dec
1808 def freq_align(self, reftime=None, method='cubic', insitu=None):
1809 """\
1810 Return a scan where all rows have been aligned in frequency/velocity.
1811 The alignment frequency frame (e.g. LSRK) is that set by function
1812 set_freqframe.
1813
1814 Parameters:
1815
1816 reftime: reference time to align at. By default, the time of
1817 the first row of data is used.
1818
1819 method: Interpolation method for regridding the spectra.
1820 Choose from "nearest", "linear", "cubic" (default)
1821 and "spline"
1822
1823 insitu: if False a new scantable is returned.
1824 Otherwise, the scaling is done in-situ
1825 The default is taken from .asaprc (False)
1826
1827 """
1828 if insitu is None: insitu = rcParams["insitu"]
1829 self._math._setinsitu(insitu)
1830 varlist = vars()
1831 reftime = reftime or ""
1832 s = scantable(self._math._freq_align(self, reftime, method))
1833 s._add_history("freq_align", varlist)
1834 if insitu: self._assign(s)
1835 else: return s
1836
1837 @asaplog_post_dec
1838 def opacity(self, tau=None, insitu=None):
1839 """\
1840 Apply an opacity correction. The data
1841 and Tsys are multiplied by the correction factor.
1842
1843 Parameters:
1844
1845 tau: (list of) opacity from which the correction factor is
1846 exp(tau*ZD)
1847 where ZD is the zenith-distance.
1848 If a list is provided, it has to be of length nIF,
1849 nIF*nPol or 1 and in order of IF/POL, e.g.
1850 [opif0pol0, opif0pol1, opif1pol0 ...]
1851 if tau is `None` the opacities are determined from a
1852 model.
1853
1854 insitu: if False a new scantable is returned.
1855 Otherwise, the scaling is done in-situ
1856 The default is taken from .asaprc (False)
1857
1858 """
1859 if insitu is None: insitu = rcParams['insitu']
1860 self._math._setinsitu(insitu)
1861 varlist = vars()
1862 if not hasattr(tau, "__len__"):
1863 tau = [tau]
1864 s = scantable(self._math._opacity(self, tau))
1865 s._add_history("opacity", varlist)
1866 if insitu: self._assign(s)
1867 else: return s
1868
1869 @asaplog_post_dec
1870 def bin(self, width=5, insitu=None):
1871 """\
1872 Return a scan where all spectra have been binned up.
1873
1874 Parameters:
1875
1876 width: The bin width (default=5) in pixels
1877
1878 insitu: if False a new scantable is returned.
1879 Otherwise, the scaling is done in-situ
1880 The default is taken from .asaprc (False)
1881
1882 """
1883 if insitu is None: insitu = rcParams['insitu']
1884 self._math._setinsitu(insitu)
1885 varlist = vars()
1886 s = scantable(self._math._bin(self, width))
1887 s._add_history("bin", varlist)
1888 if insitu:
1889 self._assign(s)
1890 else:
1891 return s
1892
1893 @asaplog_post_dec
1894 def resample(self, width=5, method='cubic', insitu=None):
1895 """\
1896 Return a scan where all spectra have been binned up.
1897
1898 Parameters:
1899
1900 width: The bin width (default=5) in pixels
1901
1902 method: Interpolation method when correcting from a table.
1903 Values are "nearest", "linear", "cubic" (default)
1904 and "spline"
1905
1906 insitu: if False a new scantable is returned.
1907 Otherwise, the scaling is done in-situ
1908 The default is taken from .asaprc (False)
1909
1910 """
1911 if insitu is None: insitu = rcParams['insitu']
1912 self._math._setinsitu(insitu)
1913 varlist = vars()
1914 s = scantable(self._math._resample(self, method, width))
1915 s._add_history("resample", varlist)
1916 if insitu: self._assign(s)
1917 else: return s
1918
1919 @asaplog_post_dec
1920 def average_pol(self, mask=None, weight='none'):
1921 """\
1922 Average the Polarisations together.
1923
1924 Parameters:
1925
1926 mask: An optional mask defining the region, where the
1927 averaging will be applied. The output will have all
1928 specified points masked.
1929
1930 weight: Weighting scheme. 'none' (default), 'var' (1/var(spec)
1931 weighted), or 'tsys' (1/Tsys**2 weighted)
1932
1933 """
1934 varlist = vars()
1935 mask = mask or ()
1936 s = scantable(self._math._averagepol(self, mask, weight.upper()))
1937 s._add_history("average_pol", varlist)
1938 return s
1939
1940 @asaplog_post_dec
1941 def average_beam(self, mask=None, weight='none'):
1942 """\
1943 Average the Beams together.
1944
1945 Parameters:
1946 mask: An optional mask defining the region, where the
1947 averaging will be applied. The output will have all
1948 specified points masked.
1949
1950 weight: Weighting scheme. 'none' (default), 'var' (1/var(spec)
1951 weighted), or 'tsys' (1/Tsys**2 weighted)
1952
1953 """
1954 varlist = vars()
1955 mask = mask or ()
1956 s = scantable(self._math._averagebeams(self, mask, weight.upper()))
1957 s._add_history("average_beam", varlist)
1958 return s
1959
1960 def parallactify(self, pflag):
1961 """\
1962 Set a flag to indicate whether this data should be treated as having
1963 been 'parallactified' (total phase == 0.0)
1964
1965 Parameters:
1966
1967 pflag: Bool indicating whether to turn this on (True) or
1968 off (False)
1969
1970 """
1971 varlist = vars()
1972 self._parallactify(pflag)
1973 self._add_history("parallactify", varlist)
1974
1975 @asaplog_post_dec
1976 def convert_pol(self, poltype=None):
1977 """\
1978 Convert the data to a different polarisation type.
1979 Note that you will need cross-polarisation terms for most conversions.
1980
1981 Parameters:
1982
1983 poltype: The new polarisation type. Valid types are:
1984 "linear", "circular", "stokes" and "linpol"
1985
1986 """
1987 varlist = vars()
1988 s = scantable(self._math._convertpol(self, poltype))
1989 s._add_history("convert_pol", varlist)
1990 return s
1991
1992 @asaplog_post_dec
1993 def smooth(self, kernel="hanning", width=5.0, order=2, plot=False, insitu=None):
1994 """\
1995 Smooth the spectrum by the specified kernel (conserving flux).
1996
1997 Parameters:
1998
1999 kernel: The type of smoothing kernel. Select from
2000 'hanning' (default), 'gaussian', 'boxcar', 'rmedian'
2001 or 'poly'
2002
2003 width: The width of the kernel in pixels. For hanning this is
2004 ignored otherwise it defauls to 5 pixels.
2005 For 'gaussian' it is the Full Width Half
2006 Maximum. For 'boxcar' it is the full width.
2007 For 'rmedian' and 'poly' it is the half width.
2008
2009 order: Optional parameter for 'poly' kernel (default is 2), to
2010 specify the order of the polnomial. Ignored by all other
2011 kernels.
2012
2013 plot: plot the original and the smoothed spectra.
2014 In this each indivual fit has to be approved, by
2015 typing 'y' or 'n'
2016
2017 insitu: if False a new scantable is returned.
2018 Otherwise, the scaling is done in-situ
2019 The default is taken from .asaprc (False)
2020
2021 """
2022 if insitu is None: insitu = rcParams['insitu']
2023 self._math._setinsitu(insitu)
2024 varlist = vars()
2025
2026 if plot: orgscan = self.copy()
2027
2028 s = scantable(self._math._smooth(self, kernel.lower(), width, order))
2029 s._add_history("smooth", varlist)
2030
2031 if plot:
2032 if rcParams['plotter.gui']:
2033 from asap.asaplotgui import asaplotgui as asaplot
2034 else:
2035 from asap.asaplot import asaplot
2036 self._p=asaplot()
2037 self._p.set_panels()
2038 ylab=s._get_ordinate_label()
2039 #self._p.palette(0,["#777777","red"])
2040 for r in xrange(s.nrow()):
2041 xsm=s._getabcissa(r)
2042 ysm=s._getspectrum(r)
2043 xorg=orgscan._getabcissa(r)
2044 yorg=orgscan._getspectrum(r)
2045 self._p.clear()
2046 self._p.hold()
2047 self._p.set_axes('ylabel',ylab)
2048 self._p.set_axes('xlabel',s._getabcissalabel(r))
2049 self._p.set_axes('title',s._getsourcename(r))
2050 self._p.set_line(label='Original',color="#777777")
2051 self._p.plot(xorg,yorg)
2052 self._p.set_line(label='Smoothed',color="red")
2053 self._p.plot(xsm,ysm)
2054 ### Ugly part for legend
2055 for i in [0,1]:
2056 self._p.subplots[0]['lines'].append([self._p.subplots[0]['axes'].lines[i]])
2057 self._p.release()
2058 ### Ugly part for legend
2059 self._p.subplots[0]['lines']=[]
2060 res = raw_input("Accept smoothing ([y]/n): ")
2061 if res.upper() == 'N':
2062 s._setspectrum(yorg, r)
2063 self._p.unmap()
2064 self._p = None
2065 del orgscan
2066
2067 if insitu: self._assign(s)
2068 else: return s
2069
2070
2071 @asaplog_post_dec
2072 def cspline_baseline(self, insitu=None, mask=None, npiece=None, clipthresh=None, clipniter=None, plot=None, outlog=None, blfile=None):
2073 """\
2074 Return a scan which has been baselined (all rows) by cubic spline function (piecewise cubic polynomial).
2075 Parameters:
2076 insitu: If False a new scantable is returned.
2077 Otherwise, the scaling is done in-situ
2078 The default is taken from .asaprc (False)
2079 mask: An optional mask
2080 npiece: Number of pieces. (default is 2)
2081 clipthresh: Clipping threshold. (default is 3.0, unit: sigma)
2082 clipniter: maximum number of iteration of 'clipthresh'-sigma clipping (default is 1)
2083 plot: *** CURRENTLY UNAVAILABLE, ALWAYS FALSE ***
2084 plot the fit and the residual. In this each
2085 indivual fit has to be approved, by typing 'y'
2086 or 'n'
2087 outlog: Output the coefficients of the best-fit
2088 function to logger (default is False)
2089 blfile: Name of a text file in which the best-fit
2090 parameter values to be written
2091 (default is "": no file/logger output)
2092
2093 Example:
2094 # return a scan baselined by a cubic spline consisting of 2 pieces (i.e., 1 internal knot),
2095 # also with 3-sigma clipping, iteration up to 4 times
2096 bscan = scan.cspline_baseline(npiece=2,clipthresh=3.0,clipniter=4)
2097 """
2098
2099 varlist = vars()
2100
2101 if insitu is None: insitu = rcParams["insitu"]
2102 if insitu:
2103 workscan = self
2104 else:
2105 workscan = self.copy()
2106
2107 nchan = workscan.nchan()
2108
2109 if mask is None: mask = [True for i in xrange(nchan)]
2110 if npiece is None: npiece = 2
2111 if clipthresh is None: clipthresh = 3.0
2112 if clipniter is None: clipniter = 1
2113 if plot is None: plot = False
2114 if outlog is None: outlog = False
2115 if blfile is None: blfile = ""
2116
2117 outblfile = (blfile != "") and os.path.exists(os.path.expanduser(os.path.expandvars(blfile)))
2118
2119 try:
2120 #CURRENTLY, PLOT=true UNAVAILABLE UNTIL cubic spline fitting is implemented as a fitter method.
2121 workscan._cspline_baseline(mask, npiece, clipthresh, clipniter, outlog, blfile)
2122
2123 workscan._add_history("cspline_baseline", varlist)
2124
2125 if insitu:
2126 self._assign(workscan)
2127 else:
2128 return workscan
2129
2130 except RuntimeError, e:
2131 msg = "The fit failed, possibly because it didn't converge."
2132 if rcParams["verbose"]:
2133 asaplog.push(str(e))
2134 asaplog.push(str(msg))
2135 return
2136 else:
2137 raise RuntimeError(str(e)+'\n'+msg)
2138
2139
2140 def auto_cspline_baseline(self, insitu=None, mask=None, npiece=None, clipthresh=None,
2141 clipniter=None, edge=None, threshold=None,
2142 chan_avg_limit=None, plot=None, outlog=None, blfile=None):
2143 """\
2144 Return a scan which has been baselined (all rows) by cubic spline
2145 function (piecewise cubic polynomial).
2146 Spectral lines are detected first using linefinder and masked out
2147 to avoid them affecting the baseline solution.
2148
2149 Parameters:
2150 insitu: if False a new scantable is returned.
2151 Otherwise, the scaling is done in-situ
2152 The default is taken from .asaprc (False)
2153 mask: an optional mask retreived from scantable
2154 npiece: Number of pieces. (default is 2)
2155 clipthresh: Clipping threshold. (default is 3.0, unit: sigma)
2156 clipniter: maximum number of iteration of 'clipthresh'-sigma clipping (default is 1)
2157 edge: an optional number of channel to drop at
2158 the edge of spectrum. If only one value is
2159 specified, the same number will be dropped
2160 from both sides of the spectrum. Default
2161 is to keep all channels. Nested tuples
2162 represent individual edge selection for
2163 different IFs (a number of spectral channels
2164 can be different)
2165 threshold: the threshold used by line finder. It is
2166 better to keep it large as only strong lines
2167 affect the baseline solution.
2168 chan_avg_limit:
2169 a maximum number of consequtive spectral
2170 channels to average during the search of
2171 weak and broad lines. The default is no
2172 averaging (and no search for weak lines).
2173 If such lines can affect the fitted baseline
2174 (e.g. a high order polynomial is fitted),
2175 increase this parameter (usually values up
2176 to 8 are reasonable). Most users of this
2177 method should find the default value sufficient.
2178 plot: *** CURRENTLY UNAVAILABLE, ALWAYS FALSE ***
2179 plot the fit and the residual. In this each
2180 indivual fit has to be approved, by typing 'y'
2181 or 'n'
2182 outlog: Output the coefficients of the best-fit
2183 function to logger (default is False)
2184 blfile: Name of a text file in which the best-fit
2185 parameter values to be written
2186 (default is "": no file/logger output)
2187
2188 Example:
2189 bscan = scan.auto_cspline_baseline(npiece=3, insitu=False)
2190 """
2191
2192 varlist = vars()
2193
2194 if insitu is None: insitu = rcParams['insitu']
2195 if insitu:
2196 workscan = self
2197 else:
2198 workscan = self.copy()
2199
2200 nchan = workscan.nchan()
2201
2202 if mask is None: mask = [True for i in xrange(nchan)]
2203 if npiece is None: npiece = 2
2204 if clipthresh is None: clipthresh = 3.0
2205 if clipniter is None: clipniter = 1
2206 if edge is None: edge = (0, 0)
2207 if threshold is None: threshold = 3
2208 if chan_avg_limit is None: chan_avg_limit = 1
2209 if plot is None: plot = False
2210 if outlog is None: outlog = False
2211 if blfile is None: blfile = ""
2212
2213 outblfile = (blfile != "") and os.path.exists(os.path.expanduser(os.path.expandvars(blfile)))
2214
2215 from asap.asaplinefind import linefinder
2216 from asap import _is_sequence_or_number as _is_valid
2217
2218 if not (isinstance(edge, list) or isinstance(edge, tuple)): edge = [ edge ]
2219 individualedge = False;
2220 if len(edge) > 1: individualedge = isinstance(edge[0], list) or isinstance(edge[0], tuple)
2221
2222 if individualedge:
2223 for edgepar in edge:
2224 if not _is_valid(edgepar, int):
2225 raise ValueError, "Each element of the 'edge' tuple has \
2226 to be a pair of integers or an integer."
2227 else:
2228 if not _is_valid(edge, int):
2229 raise ValueError, "Parameter 'edge' has to be an integer or a \
2230 pair of integers specified as a tuple. \
2231 Nested tuples are allowed \
2232 to make individual selection for different IFs."
2233
2234 if len(edge) > 1:
2235 curedge = edge
2236 else:
2237 curedge = edge + edge
2238
2239 try:
2240 #CURRENTLY, PLOT=true UNAVAILABLE UNTIL cubic spline fitting is implemented as a fitter method.
2241 if individualedge:
2242 curedge = []
2243 for i in xrange(len(edge)):
2244 curedge += edge[i]
2245
2246 workscan._auto_cspline_baseline(mask, npiece, clipthresh, clipniter, curedge, threshold, chan_avg_limit, outlog, blfile)
2247
2248 workscan._add_history("auto_cspline_baseline", varlist)
2249
2250 if insitu:
2251 self._assign(workscan)
2252 else:
2253 return workscan
2254
2255 except RuntimeError, e:
2256 msg = "The fit failed, possibly because it didn't converge."
2257 if rcParams["verbose"]:
2258 asaplog.push(str(e))
2259 asaplog.push(str(msg))
2260 return
2261 else:
2262 raise RuntimeError(str(e)+'\n'+msg)
2263
2264
2265 @asaplog_post_dec
2266 def poly_baseline(self, insitu=None, mask=None, order=None, plot=None, outlog=None, blfile=None):
2267 """\
2268 Return a scan which has been baselined (all rows) by a polynomial.
2269 Parameters:
2270 insitu: if False a new scantable is returned.
2271 Otherwise, the scaling is done in-situ
2272 The default is taken from .asaprc (False)
2273 mask: an optional mask
2274 order: the order of the polynomial (default is 0)
2275 plot: plot the fit and the residual. In this each
2276 indivual fit has to be approved, by typing 'y'
2277 or 'n'
2278 outlog: Output the coefficients of the best-fit
2279 function to logger (default is False)
2280 blfile: Name of a text file in which the best-fit
2281 parameter values to be written
2282 (default is "": no file/logger output)
2283
2284 Example:
2285 # return a scan baselined by a third order polynomial,
2286 # not using a mask
2287 bscan = scan.poly_baseline(order=3)
2288 """
2289
2290 varlist = vars()
2291
2292 if insitu is None: insitu = rcParams["insitu"]
2293 if insitu:
2294 workscan = self
2295 else:
2296 workscan = self.copy()
2297
2298 nchan = workscan.nchan()
2299
2300 if mask is None: mask = [True for i in xrange(nchan)]
2301 if order is None: order = 0
2302 if plot is None: plot = False
2303 if outlog is None: outlog = False
2304 if blfile is None: blfile = ""
2305
2306 outblfile = (blfile != "") and os.path.exists(os.path.expanduser(os.path.expandvars(blfile)))
2307
2308 try:
2309 rows = xrange(workscan.nrow())
2310
2311 #if len(rows) > 0: workscan._init_blinfo()
2312
2313 if plot:
2314 if outblfile: blf = open(blfile, "a")
2315
2316 f = fitter()
2317 f.set_function(lpoly=order)
2318 for r in rows:
2319 f.x = workscan._getabcissa(r)
2320 f.y = workscan._getspectrum(r)
2321 f.mask = mask_and(mask, workscan._getmask(r)) # (CAS-1434)
2322 f.data = None
2323 f.fit()
2324
2325 f.plot(residual=True)
2326 accept_fit = raw_input("Accept fit ( [y]/n ): ")
2327 if accept_fit.upper() == "N":
2328 #workscan._append_blinfo(None, None, None)
2329 continue
2330
2331 blpars = f.get_parameters()
2332 masklist = workscan.get_masklist(f.mask, row=r, silent=True)
2333 #workscan._append_blinfo(blpars, masklist, f.mask)
2334 workscan._setspectrum(f.fitter.getresidual(), r)
2335
2336 if outblfile:
2337 rms = workscan.get_rms(f.mask, r)
2338 dataout = workscan.format_blparams_row(blpars["params"], blpars["fixed"], rms, str(masklist), r, True)
2339 blf.write(dataout)
2340
2341 f._p.unmap()
2342 f._p = None
2343
2344 if outblfile: blf.close()
2345 else:
2346 workscan._poly_baseline(mask, order, outlog, blfile)
2347
2348 workscan._add_history("poly_baseline", varlist)
2349
2350 if insitu:
2351 self._assign(workscan)
2352 else:
2353 return workscan
2354
2355 except RuntimeError, e:
2356 msg = "The fit failed, possibly because it didn't converge."
2357 if rcParams["verbose"]:
2358 asaplog.push(str(e))
2359 asaplog.push(str(msg))
2360 return
2361 else:
2362 raise RuntimeError(str(e)+'\n'+msg)
2363
2364
2365 def auto_poly_baseline(self, insitu=None, mask=None, order=None, edge=None, threshold=None,
2366 chan_avg_limit=None, plot=None, outlog=None, blfile=None):
2367 """\
2368 Return a scan which has been baselined (all rows) by a polynomial.
2369 Spectral lines are detected first using linefinder and masked out
2370 to avoid them affecting the baseline solution.
2371
2372 Parameters:
2373 insitu: if False a new scantable is returned.
2374 Otherwise, the scaling is done in-situ
2375 The default is taken from .asaprc (False)
2376 mask: an optional mask retreived from scantable
2377 order: the order of the polynomial (default is 0)
2378 edge: an optional number of channel to drop at
2379 the edge of spectrum. If only one value is
2380 specified, the same number will be dropped
2381 from both sides of the spectrum. Default
2382 is to keep all channels. Nested tuples
2383 represent individual edge selection for
2384 different IFs (a number of spectral channels
2385 can be different)
2386 threshold: the threshold used by line finder. It is
2387 better to keep it large as only strong lines
2388 affect the baseline solution.
2389 chan_avg_limit:
2390 a maximum number of consequtive spectral
2391 channels to average during the search of
2392 weak and broad lines. The default is no
2393 averaging (and no search for weak lines).
2394 If such lines can affect the fitted baseline
2395 (e.g. a high order polynomial is fitted),
2396 increase this parameter (usually values up
2397 to 8 are reasonable). Most users of this
2398 method should find the default value sufficient.
2399 plot: plot the fit and the residual. In this each
2400 indivual fit has to be approved, by typing 'y'
2401 or 'n'
2402 outlog: Output the coefficients of the best-fit
2403 function to logger (default is False)
2404 blfile: Name of a text file in which the best-fit
2405 parameter values to be written
2406 (default is "": no file/logger output)
2407
2408 Example:
2409 bscan = scan.auto_poly_baseline(order=7, insitu=False)
2410 """
2411
2412 varlist = vars()
2413
2414 if insitu is None: insitu = rcParams['insitu']
2415 if insitu:
2416 workscan = self
2417 else:
2418 workscan = self.copy()
2419
2420 nchan = workscan.nchan()
2421
2422 if mask is None: mask = [True for i in xrange(nchan)]
2423 if order is None: order = 0
2424 if edge is None: edge = (0, 0)
2425 if threshold is None: threshold = 3
2426 if chan_avg_limit is None: chan_avg_limit = 1
2427 if plot is None: plot = False
2428 if outlog is None: outlog = False
2429 if blfile is None: blfile = ""
2430
2431 outblfile = (blfile != "") and os.path.exists(os.path.expanduser(os.path.expandvars(blfile)))
2432
2433 from asap.asaplinefind import linefinder
2434 from asap import _is_sequence_or_number as _is_valid
2435
2436 if not (isinstance(edge, list) or isinstance(edge, tuple)): edge = [ edge ]
2437 individualedge = False;
2438 if len(edge) > 1: individualedge = isinstance(edge[0], list) or isinstance(edge[0], tuple)
2439
2440 if individualedge:
2441 for edgepar in edge:
2442 if not _is_valid(edgepar, int):
2443 raise ValueError, "Each element of the 'edge' tuple has \
2444 to be a pair of integers or an integer."
2445 else:
2446 if not _is_valid(edge, int):
2447 raise ValueError, "Parameter 'edge' has to be an integer or a \
2448 pair of integers specified as a tuple. \
2449 Nested tuples are allowed \
2450 to make individual selection for different IFs."
2451
2452 if len(edge) > 1:
2453 curedge = edge
2454 else:
2455 curedge = edge + edge
2456
2457 try:
2458 rows = xrange(workscan.nrow())
2459
2460 #if len(rows) > 0: workscan._init_blinfo()
2461
2462 if plot:
2463 if outblfile: blf = open(blfile, "a")
2464
2465 fl = linefinder()
2466 fl.set_options(threshold=threshold,avg_limit=chan_avg_limit)
2467 fl.set_scan(workscan)
2468 f = fitter()
2469 f.set_function(lpoly=order)
2470
2471 for r in rows:
2472 if individualedge:
2473 if len(edge) <= workscan.getif(r):
2474 raise RuntimeError, "Number of edge elements appear to " \
2475 "be less than the number of IFs"
2476 else:
2477 curedge = edge[workscan.getif(r)]
2478
2479 fl.find_lines(r, mask_and(mask, workscan._getmask(r)), curedge) # (CAS-1434)
2480
2481 f.x = workscan._getabcissa(r)
2482 f.y = workscan._getspectrum(r)
2483 f.mask = fl.get_mask()
2484 f.data = None
2485 f.fit()
2486
2487 f.plot(residual=True)
2488 accept_fit = raw_input("Accept fit ( [y]/n ): ")
2489 if accept_fit.upper() == "N":
2490 #workscan._append_blinfo(None, None, None)
2491 continue
2492
2493 blpars = f.get_parameters()
2494 masklist = workscan.get_masklist(f.mask, row=r, silent=True)
2495 #workscan._append_blinfo(blpars, masklist, f.mask)
2496 workscan._setspectrum(f.fitter.getresidual(), r)
2497
2498 if outblfile:
2499 rms = workscan.get_rms(f.mask, r)
2500 dataout = workscan.format_blparams_row(blpars["params"], blpars["fixed"], rms, str(masklist), r, True)
2501 blf.write(dataout)
2502
2503 f._p.unmap()
2504 f._p = None
2505
2506 if outblfile: blf.close()
2507
2508 else:
2509 if individualedge:
2510 curedge = []
2511 for i in xrange(len(edge)):
2512 curedge += edge[i]
2513
2514 workscan._auto_poly_baseline(mask, order, curedge, threshold, chan_avg_limit, outlog, blfile)
2515
2516 workscan._add_history("auto_poly_baseline", varlist)
2517
2518 if insitu:
2519 self._assign(workscan)
2520 else:
2521 return workscan
2522
2523 except RuntimeError, e:
2524 msg = "The fit failed, possibly because it didn't converge."
2525 if rcParams["verbose"]:
2526 asaplog.push(str(e))
2527 asaplog.push(str(msg))
2528 return
2529 else:
2530 raise RuntimeError(str(e)+'\n'+msg)
2531
2532
2533 ### OBSOLETE ##################################################################
2534 @asaplog_post_dec
2535 def old_poly_baseline(self, mask=None, order=0, plot=False, uselin=False, insitu=None, rows=None):
2536 """
2537 Return a scan which has been baselined (all rows) by a polynomial.
2538
2539 Parameters:
2540
2541 mask: an optional mask
2542
2543 order: the order of the polynomial (default is 0)
2544
2545 plot: plot the fit and the residual. In this each
2546 indivual fit has to be approved, by typing 'y'
2547 or 'n'
2548
2549 uselin: use linear polynomial fit
2550
2551 insitu: if False a new scantable is returned.
2552 Otherwise, the scaling is done in-situ
2553 The default is taken from .asaprc (False)
2554
2555 rows: row numbers of spectra to be processed.
2556 (default is None: for all rows)
2557
2558 Example:
2559 # return a scan baselined by a third order polynomial,
2560 # not using a mask
2561 bscan = scan.poly_baseline(order=3)
2562
2563 """
2564 if insitu is None: insitu = rcParams['insitu']
2565 if not insitu:
2566 workscan = self.copy()
2567 else:
2568 workscan = self
2569 varlist = vars()
2570 if mask is None:
2571 mask = [True for i in xrange(self.nchan())]
2572
2573 try:
2574 f = fitter()
2575 if uselin:
2576 f.set_function(lpoly=order)
2577 else:
2578 f.set_function(poly=order)
2579
2580 if rows == None:
2581 rows = xrange(workscan.nrow())
2582 elif isinstance(rows, int):
2583 rows = [ rows ]
2584
2585 if len(rows) > 0:
2586 self.blpars = []
2587 self.masklists = []
2588 self.actualmask = []
2589
2590 for r in rows:
2591 f.x = workscan._getabcissa(r)
2592 f.y = workscan._getspectrum(r)
2593 f.mask = mask_and(mask, workscan._getmask(r)) # (CAS-1434)
2594 f.data = None
2595 f.fit()
2596 if plot:
2597 f.plot(residual=True)
2598 x = raw_input("Accept fit ( [y]/n ): ")
2599 if x.upper() == 'N':
2600 self.blpars.append(None)
2601 self.masklists.append(None)
2602 self.actualmask.append(None)
2603 continue
2604 workscan._setspectrum(f.fitter.getresidual(), r)
2605 self.blpars.append(f.get_parameters())
2606 self.masklists.append(workscan.get_masklist(f.mask, row=r, silent=True))
2607 self.actualmask.append(f.mask)
2608
2609 if plot:
2610 f._p.unmap()
2611 f._p = None
2612 workscan._add_history("poly_baseline", varlist)
2613 if insitu:
2614 self._assign(workscan)
2615 else:
2616 return workscan
2617 except RuntimeError:
2618 msg = "The fit failed, possibly because it didn't converge."
2619 raise RuntimeError(msg)
2620
2621 def _init_blinfo(self):
2622 """\
2623 Initialise the following three auxiliary members:
2624 blpars : parameters of the best-fit baseline,
2625 masklists : mask data (edge positions of masked channels) and
2626 actualmask : mask data (in boolean list),
2627 to keep for use later (including output to logger/text files).
2628 Used by poly_baseline() and auto_poly_baseline() in case of
2629 'plot=True'.
2630 """
2631 self.blpars = []
2632 self.masklists = []
2633 self.actualmask = []
2634 return
2635
2636 def _append_blinfo(self, data_blpars, data_masklists, data_actualmask):
2637 """\
2638 Append baseline-fitting related info to blpars, masklist and
2639 actualmask.
2640 """
2641 self.blpars.append(data_blpars)
2642 self.masklists.append(data_masklists)
2643 self.actualmask.append(data_actualmask)
2644 return
2645
2646 @asaplog_post_dec
2647 def rotate_linpolphase(self, angle):
2648 """\
2649 Rotate the phase of the complex polarization O=Q+iU correlation.
2650 This is always done in situ in the raw data. So if you call this
2651 function more than once then each call rotates the phase further.
2652
2653 Parameters:
2654
2655 angle: The angle (degrees) to rotate (add) by.
2656
2657 Example::
2658
2659 scan.rotate_linpolphase(2.3)
2660
2661 """
2662 varlist = vars()
2663 self._math._rotate_linpolphase(self, angle)
2664 self._add_history("rotate_linpolphase", varlist)
2665 return
2666
2667 @asaplog_post_dec
2668 def rotate_xyphase(self, angle):
2669 """\
2670 Rotate the phase of the XY correlation. This is always done in situ
2671 in the data. So if you call this function more than once
2672 then each call rotates the phase further.
2673
2674 Parameters:
2675
2676 angle: The angle (degrees) to rotate (add) by.
2677
2678 Example::
2679
2680 scan.rotate_xyphase(2.3)
2681
2682 """
2683 varlist = vars()
2684 self._math._rotate_xyphase(self, angle)
2685 self._add_history("rotate_xyphase", varlist)
2686 return
2687
2688 @asaplog_post_dec
2689 def swap_linears(self):
2690 """\
2691 Swap the linear polarisations XX and YY, or better the first two
2692 polarisations as this also works for ciculars.
2693 """
2694 varlist = vars()
2695 self._math._swap_linears(self)
2696 self._add_history("swap_linears", varlist)
2697 return
2698
2699 @asaplog_post_dec
2700 def invert_phase(self):
2701 """\
2702 Invert the phase of the complex polarisation
2703 """
2704 varlist = vars()
2705 self._math._invert_phase(self)
2706 self._add_history("invert_phase", varlist)
2707 return
2708
2709 @asaplog_post_dec
2710 def add(self, offset, insitu=None):
2711 """\
2712 Return a scan where all spectra have the offset added
2713
2714 Parameters:
2715
2716 offset: the offset
2717
2718 insitu: if False a new scantable is returned.
2719 Otherwise, the scaling is done in-situ
2720 The default is taken from .asaprc (False)
2721
2722 """
2723 if insitu is None: insitu = rcParams['insitu']
2724 self._math._setinsitu(insitu)
2725 varlist = vars()
2726 s = scantable(self._math._unaryop(self, offset, "ADD", False))
2727 s._add_history("add", varlist)
2728 if insitu:
2729 self._assign(s)
2730 else:
2731 return s
2732
2733 @asaplog_post_dec
2734 def scale(self, factor, tsys=True, insitu=None):
2735 """\
2736
2737 Return a scan where all spectra are scaled by the given 'factor'
2738
2739 Parameters:
2740
2741 factor: the scaling factor (float or 1D float list)
2742
2743 insitu: if False a new scantable is returned.
2744 Otherwise, the scaling is done in-situ
2745 The default is taken from .asaprc (False)
2746
2747 tsys: if True (default) then apply the operation to Tsys
2748 as well as the data
2749
2750 """
2751 if insitu is None: insitu = rcParams['insitu']
2752 self._math._setinsitu(insitu)
2753 varlist = vars()
2754 s = None
2755 import numpy
2756 if isinstance(factor, list) or isinstance(factor, numpy.ndarray):
2757 if isinstance(factor[0], list) or isinstance(factor[0], numpy.ndarray):
2758 from asapmath import _array2dOp
2759 s = _array2dOp( self.copy(), factor, "MUL", tsys )
2760 else:
2761 s = scantable( self._math._arrayop( self.copy(), factor, "MUL", tsys ) )
2762 else:
2763 s = scantable(self._math._unaryop(self.copy(), factor, "MUL", tsys))
2764 s._add_history("scale", varlist)
2765 if insitu:
2766 self._assign(s)
2767 else:
2768 return s
2769
2770 def set_sourcetype(self, match, matchtype="pattern",
2771 sourcetype="reference"):
2772 """\
2773 Set the type of the source to be an source or reference scan
2774 using the provided pattern.
2775
2776 Parameters:
2777
2778 match: a Unix style pattern, regular expression or selector
2779
2780 matchtype: 'pattern' (default) UNIX style pattern or
2781 'regex' regular expression
2782
2783 sourcetype: the type of the source to use (source/reference)
2784
2785 """
2786 varlist = vars()
2787 basesel = self.get_selection()
2788 stype = -1
2789 if sourcetype.lower().startswith("r"):
2790 stype = 1
2791 elif sourcetype.lower().startswith("s"):
2792 stype = 0
2793 else:
2794 raise ValueError("Illegal sourcetype use s(ource) or r(eference)")
2795 if matchtype.lower().startswith("p"):
2796 matchtype = "pattern"
2797 elif matchtype.lower().startswith("r"):
2798 matchtype = "regex"
2799 else:
2800 raise ValueError("Illegal matchtype, use p(attern) or r(egex)")
2801 sel = selector()
2802 if isinstance(match, selector):
2803 sel = match
2804 else:
2805 sel.set_query("SRCNAME == %s('%s')" % (matchtype, match))
2806 self.set_selection(basesel+sel)
2807 self._setsourcetype(stype)
2808 self.set_selection(basesel)
2809 self._add_history("set_sourcetype", varlist)
2810
2811 @asaplog_post_dec
2812 @preserve_selection
2813 def auto_quotient(self, preserve=True, mode='paired', verify=False):
2814 """\
2815 This function allows to build quotients automatically.
2816 It assumes the observation to have the same number of
2817 "ons" and "offs"
2818
2819 Parameters:
2820
2821 preserve: you can preserve (default) the continuum or
2822 remove it. The equations used are
2823
2824 preserve: Output = Toff * (on/off) - Toff
2825
2826 remove: Output = Toff * (on/off) - Ton
2827
2828 mode: the on/off detection mode
2829 'paired' (default)
2830 identifies 'off' scans by the
2831 trailing '_R' (Mopra/Parkes) or
2832 '_e'/'_w' (Tid) and matches
2833 on/off pairs from the observing pattern
2834 'time'
2835 finds the closest off in time
2836
2837 .. todo:: verify argument is not implemented
2838
2839 """
2840 varlist = vars()
2841 modes = ["time", "paired"]
2842 if not mode in modes:
2843 msg = "please provide valid mode. Valid modes are %s" % (modes)
2844 raise ValueError(msg)
2845 s = None
2846 if mode.lower() == "paired":
2847 sel = self.get_selection()
2848 sel.set_query("SRCTYPE==psoff")
2849 self.set_selection(sel)
2850 offs = self.copy()
2851 sel.set_query("SRCTYPE==pson")
2852 self.set_selection(sel)
2853 ons = self.copy()
2854 s = scantable(self._math._quotient(ons, offs, preserve))
2855 elif mode.lower() == "time":
2856 s = scantable(self._math._auto_quotient(self, mode, preserve))
2857 s._add_history("auto_quotient", varlist)
2858 return s
2859
2860 @asaplog_post_dec
2861 def mx_quotient(self, mask = None, weight='median', preserve=True):
2862 """\
2863 Form a quotient using "off" beams when observing in "MX" mode.
2864
2865 Parameters:
2866
2867 mask: an optional mask to be used when weight == 'stddev'
2868
2869 weight: How to average the off beams. Default is 'median'.
2870
2871 preserve: you can preserve (default) the continuum or
2872 remove it. The equations used are:
2873
2874 preserve: Output = Toff * (on/off) - Toff
2875
2876 remove: Output = Toff * (on/off) - Ton
2877
2878 """
2879 mask = mask or ()
2880 varlist = vars()
2881 on = scantable(self._math._mx_extract(self, 'on'))
2882 preoff = scantable(self._math._mx_extract(self, 'off'))
2883 off = preoff.average_time(mask=mask, weight=weight, scanav=False)
2884 from asapmath import quotient
2885 q = quotient(on, off, preserve)
2886 q._add_history("mx_quotient", varlist)
2887 return q
2888
2889 @asaplog_post_dec
2890 def freq_switch(self, insitu=None):
2891 """\
2892 Apply frequency switching to the data.
2893
2894 Parameters:
2895
2896 insitu: if False a new scantable is returned.
2897 Otherwise, the swictching is done in-situ
2898 The default is taken from .asaprc (False)
2899
2900 """
2901 if insitu is None: insitu = rcParams['insitu']
2902 self._math._setinsitu(insitu)
2903 varlist = vars()
2904 s = scantable(self._math._freqswitch(self))
2905 s._add_history("freq_switch", varlist)
2906 if insitu:
2907 self._assign(s)
2908 else:
2909 return s
2910
2911 @asaplog_post_dec
2912 def recalc_azel(self):
2913 """Recalculate the azimuth and elevation for each position."""
2914 varlist = vars()
2915 self._recalcazel()
2916 self._add_history("recalc_azel", varlist)
2917 return
2918
2919 @asaplog_post_dec
2920 def __add__(self, other):
2921 varlist = vars()
2922 s = None
2923 if isinstance(other, scantable):
2924 s = scantable(self._math._binaryop(self, other, "ADD"))
2925 elif isinstance(other, float):
2926 s = scantable(self._math._unaryop(self, other, "ADD", False))
2927 else:
2928 raise TypeError("Other input is not a scantable or float value")
2929 s._add_history("operator +", varlist)
2930 return s
2931
2932 @asaplog_post_dec
2933 def __sub__(self, other):
2934 """
2935 implicit on all axes and on Tsys
2936 """
2937 varlist = vars()
2938 s = None
2939 if isinstance(other, scantable):
2940 s = scantable(self._math._binaryop(self, other, "SUB"))
2941 elif isinstance(other, float):
2942 s = scantable(self._math._unaryop(self, other, "SUB", False))
2943 else:
2944 raise TypeError("Other input is not a scantable or float value")
2945 s._add_history("operator -", varlist)
2946 return s
2947
2948 @asaplog_post_dec
2949 def __mul__(self, other):
2950 """
2951 implicit on all axes and on Tsys
2952 """
2953 varlist = vars()
2954 s = None
2955 if isinstance(other, scantable):
2956 s = scantable(self._math._binaryop(self, other, "MUL"))
2957 elif isinstance(other, float):
2958 s = scantable(self._math._unaryop(self, other, "MUL", False))
2959 else:
2960 raise TypeError("Other input is not a scantable or float value")
2961 s._add_history("operator *", varlist)
2962 return s
2963
2964
2965 @asaplog_post_dec
2966 def __div__(self, other):
2967 """
2968 implicit on all axes and on Tsys
2969 """
2970 varlist = vars()
2971 s = None
2972 if isinstance(other, scantable):
2973 s = scantable(self._math._binaryop(self, other, "DIV"))
2974 elif isinstance(other, float):
2975 if other == 0.0:
2976 raise ZeroDivisionError("Dividing by zero is not recommended")
2977 s = scantable(self._math._unaryop(self, other, "DIV", False))
2978 else:
2979 raise TypeError("Other input is not a scantable or float value")
2980 s._add_history("operator /", varlist)
2981 return s
2982
2983 @asaplog_post_dec
2984 def get_fit(self, row=0):
2985 """\
2986 Print or return the stored fits for a row in the scantable
2987
2988 Parameters:
2989
2990 row: the row which the fit has been applied to.
2991
2992 """
2993 if row > self.nrow():
2994 return
2995 from asap.asapfit import asapfit
2996 fit = asapfit(self._getfit(row))
2997 asaplog.push( '%s' %(fit) )
2998 return fit.as_dict()
2999
3000 def flag_nans(self):
3001 """\
3002 Utility function to flag NaN values in the scantable.
3003 """
3004 import numpy
3005 basesel = self.get_selection()
3006 for i in range(self.nrow()):
3007 sel = self.get_row_selector(i)
3008 self.set_selection(basesel+sel)
3009 nans = numpy.isnan(self._getspectrum(0))
3010 if numpy.any(nans):
3011 bnans = [ bool(v) for v in nans]
3012 self.flag(bnans)
3013 self.set_selection(basesel)
3014
3015 def get_row_selector(self, rowno):
3016 #return selector(beams=self.getbeam(rowno),
3017 # ifs=self.getif(rowno),
3018 # pols=self.getpol(rowno),
3019 # scans=self.getscan(rowno),
3020 # cycles=self.getcycle(rowno))
3021 return selector(rows=[rowno])
3022
3023 def _add_history(self, funcname, parameters):
3024 if not rcParams['scantable.history']:
3025 return
3026 # create date
3027 sep = "##"
3028 from datetime import datetime
3029 dstr = datetime.now().strftime('%Y/%m/%d %H:%M:%S')
3030 hist = dstr+sep
3031 hist += funcname+sep#cdate+sep
3032 if parameters.has_key('self'): del parameters['self']
3033 for k, v in parameters.iteritems():
3034 if type(v) is dict:
3035 for k2, v2 in v.iteritems():
3036 hist += k2
3037 hist += "="
3038 if isinstance(v2, scantable):
3039 hist += 'scantable'
3040 elif k2 == 'mask':
3041 if isinstance(v2, list) or isinstance(v2, tuple):
3042 hist += str(self._zip_mask(v2))
3043 else:
3044 hist += str(v2)
3045 else:
3046 hist += str(v2)
3047 else:
3048 hist += k
3049 hist += "="
3050 if isinstance(v, scantable):
3051 hist += 'scantable'
3052 elif k == 'mask':
3053 if isinstance(v, list) or isinstance(v, tuple):
3054 hist += str(self._zip_mask(v))
3055 else:
3056 hist += str(v)
3057 else:
3058 hist += str(v)
3059 hist += sep
3060 hist = hist[:-2] # remove trailing '##'
3061 self._addhistory(hist)
3062
3063
3064 def _zip_mask(self, mask):
3065 mask = list(mask)
3066 i = 0
3067 segments = []
3068 while mask[i:].count(1):
3069 i += mask[i:].index(1)
3070 if mask[i:].count(0):
3071 j = i + mask[i:].index(0)
3072 else:
3073 j = len(mask)
3074 segments.append([i, j])
3075 i = j
3076 return segments
3077
3078 def _get_ordinate_label(self):
3079 fu = "("+self.get_fluxunit()+")"
3080 import re
3081 lbl = "Intensity"
3082 if re.match(".K.", fu):
3083 lbl = "Brightness Temperature "+ fu
3084 elif re.match(".Jy.", fu):
3085 lbl = "Flux density "+ fu
3086 return lbl
3087
3088 def _check_ifs(self):
3089 #nchans = [self.nchan(i) for i in range(self.nif(-1))]
3090 nchans = [self.nchan(i) for i in self.getifnos()]
3091 nchans = filter(lambda t: t > 0, nchans)
3092 return (sum(nchans)/len(nchans) == nchans[0])
3093
3094 @asaplog_post_dec
3095 #def _fill(self, names, unit, average, getpt, antenna):
3096 def _fill(self, names, unit, average, opts={}):
3097 first = True
3098 fullnames = []
3099 for name in names:
3100 name = os.path.expandvars(name)
3101 name = os.path.expanduser(name)
3102 if not os.path.exists(name):
3103 msg = "File '%s' does not exists" % (name)
3104 raise IOError(msg)
3105 fullnames.append(name)
3106 if average:
3107 asaplog.push('Auto averaging integrations')
3108 stype = int(rcParams['scantable.storage'].lower() == 'disk')
3109 for name in fullnames:
3110 tbl = Scantable(stype)
3111 if is_ms( name ):
3112 r = msfiller( tbl )
3113 else:
3114 r = filler( tbl )
3115 rx = rcParams['scantable.reference']
3116 r.setreferenceexpr(rx)
3117 #r = filler(tbl)
3118 #rx = rcParams['scantable.reference']
3119 #r.setreferenceexpr(rx)
3120 msg = "Importing %s..." % (name)
3121 asaplog.push(msg, False)
3122 #opts = {'ms': {'antenna' : antenna, 'getpt': getpt} }
3123 r.open(name, opts)# antenna, -1, -1, getpt)
3124 r.fill()
3125 if average:
3126 tbl = self._math._average((tbl, ), (), 'NONE', 'SCAN')
3127 if not first:
3128 tbl = self._math._merge([self, tbl])
3129 Scantable.__init__(self, tbl)
3130 r.close()
3131 del r, tbl
3132 first = False
3133 #flush log
3134 asaplog.post()
3135 if unit is not None:
3136 self.set_fluxunit(unit)
3137 if not is_casapy():
3138 self.set_freqframe(rcParams['scantable.freqframe'])
3139
3140
3141 def __getitem__(self, key):
3142 if key < 0:
3143 key += self.nrow()
3144 if key >= self.nrow():
3145 raise IndexError("Row index out of range.")
3146 return self._getspectrum(key)
3147
3148 def __setitem__(self, key, value):
3149 if key < 0:
3150 key += self.nrow()
3151 if key >= self.nrow():
3152 raise IndexError("Row index out of range.")
3153 if not hasattr(value, "__len__") or \
3154 len(value) > self.nchan(self.getif(key)):
3155 raise ValueError("Spectrum length doesn't match.")
3156 return self._setspectrum(value, key)
3157
3158 def __len__(self):
3159 return self.nrow()
3160
3161 def __iter__(self):
3162 for i in range(len(self)):
3163 yield self[i]
Note: See TracBrowser for help on using the repository browser.