source: branches/parallel/python/scantable.py@ 2041

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

New Development: No

JIRA Issue: No

Ready for Test: Yes/No

Interface Changes: Yes/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...

Bug fix.
The average parameter of scantable constructor makes effective
for Scantable format.


  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 94.0 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# def get_restfreqs(self):
1275# """
1276# Get the restfrequency(s) stored in this scantable.
1277# The return value(s) are always of unit 'Hz'
1278# Parameters:
1279# none
1280# Returns:
1281# a list of doubles
1282# """
1283# return list(self._getrestfreqs())
1284
1285 def get_restfreqs(self, ids=None):
1286 """\
1287 Get the restfrequency(s) stored in this scantable.
1288 The return value(s) are always of unit 'Hz'
1289
1290 Parameters:
1291
1292 ids: (optional) a list of MOLECULE_ID for that restfrequency(s) to
1293 be retrieved
1294
1295 Returns:
1296
1297 dictionary containing ids and a list of doubles for each id
1298
1299 """
1300 if ids is None:
1301 rfreqs={}
1302 idlist = self.getmolnos()
1303 for i in idlist:
1304 rfreqs[i]=list(self._getrestfreqs(i))
1305 return rfreqs
1306 else:
1307 if type(ids)==list or type(ids)==tuple:
1308 rfreqs={}
1309 for i in ids:
1310 rfreqs[i]=list(self._getrestfreqs(i))
1311 return rfreqs
1312 else:
1313 return list(self._getrestfreqs(ids))
1314 #return list(self._getrestfreqs(ids))
1315
1316 def set_restfreqs(self, freqs=None, unit='Hz'):
1317 """\
1318 Set or replace the restfrequency specified and
1319 if the 'freqs' argument holds a scalar,
1320 then that rest frequency will be applied to all the selected
1321 data. If the 'freqs' argument holds
1322 a vector, then it MUST be of equal or smaller length than
1323 the number of IFs (and the available restfrequencies will be
1324 replaced by this vector). In this case, *all* data have
1325 the restfrequency set per IF according
1326 to the corresponding value you give in the 'freqs' vector.
1327 E.g. 'freqs=[1e9, 2e9]' would mean IF 0 gets restfreq 1e9 and
1328 IF 1 gets restfreq 2e9.
1329
1330 You can also specify the frequencies via a linecatalog.
1331
1332 Parameters:
1333
1334 freqs: list of rest frequency values or string idenitfiers
1335
1336 unit: unit for rest frequency (default 'Hz')
1337
1338
1339 Example::
1340
1341 # set the given restfrequency for the all currently selected IFs
1342 scan.set_restfreqs(freqs=1.4e9)
1343 # set restfrequencies for the n IFs (n > 1) in the order of the
1344 # list, i.e
1345 # IF0 -> 1.4e9, IF1 -> 1.41e9, IF3 -> 1.42e9
1346 # len(list_of_restfreqs) == nIF
1347 # for nIF == 1 the following will set multiple restfrequency for
1348 # that IF
1349 scan.set_restfreqs(freqs=[1.4e9, 1.41e9, 1.42e9])
1350 # set multiple restfrequencies per IF. as a list of lists where
1351 # the outer list has nIF elements, the inner s arbitrary
1352 scan.set_restfreqs(freqs=[[1.4e9, 1.41e9], [1.67e9]])
1353
1354 *Note*:
1355
1356 To do more sophisticate Restfrequency setting, e.g. on a
1357 source and IF basis, use scantable.set_selection() before using
1358 this function::
1359
1360 # provided your scantable is called scan
1361 selection = selector()
1362 selection.set_name("ORION*")
1363 selection.set_ifs([1])
1364 scan.set_selection(selection)
1365 scan.set_restfreqs(freqs=86.6e9)
1366
1367 """
1368 varlist = vars()
1369 from asap import linecatalog
1370 # simple value
1371 if isinstance(freqs, int) or isinstance(freqs, float):
1372 self._setrestfreqs([freqs], [""], unit)
1373 # list of values
1374 elif isinstance(freqs, list) or isinstance(freqs, tuple):
1375 # list values are scalars
1376 if isinstance(freqs[-1], int) or isinstance(freqs[-1], float):
1377 if len(freqs) == 1:
1378 self._setrestfreqs(freqs, [""], unit)
1379 else:
1380 # allow the 'old' mode of setting mulitple IFs
1381 sel = selector()
1382 savesel = self._getselection()
1383 iflist = self.getifnos()
1384 if len(freqs)>len(iflist):
1385 raise ValueError("number of elements in list of list "
1386 "exeeds the current IF selections")
1387 iflist = self.getifnos()
1388 for i, fval in enumerate(freqs):
1389 sel.set_ifs(iflist[i])
1390 self._setselection(sel)
1391 self._setrestfreqs([fval], [""], unit)
1392 self._setselection(savesel)
1393
1394 # list values are dict, {'value'=, 'name'=)
1395 elif isinstance(freqs[-1], dict):
1396 values = []
1397 names = []
1398 for d in freqs:
1399 values.append(d["value"])
1400 names.append(d["name"])
1401 self._setrestfreqs(values, names, unit)
1402 elif isinstance(freqs[-1], list) or isinstance(freqs[-1], tuple):
1403 sel = selector()
1404 savesel = self._getselection()
1405 iflist = self.getifnos()
1406 if len(freqs)>len(iflist):
1407 raise ValueError("number of elements in list of list exeeds"
1408 " the current IF selections")
1409 for i, fval in enumerate(freqs):
1410 sel.set_ifs(iflist[i])
1411 self._setselection(sel)
1412 self._setrestfreqs(fval, [""], unit)
1413 self._setselection(savesel)
1414 # freqs are to be taken from a linecatalog
1415 elif isinstance(freqs, linecatalog):
1416 sel = selector()
1417 savesel = self._getselection()
1418 for i in xrange(freqs.nrow()):
1419 sel.set_ifs(iflist[i])
1420 self._setselection(sel)
1421 self._setrestfreqs([freqs.get_frequency(i)],
1422 [freqs.get_name(i)], "MHz")
1423 # ensure that we are not iterating past nIF
1424 if i == self.nif()-1: break
1425 self._setselection(savesel)
1426 else:
1427 return
1428 self._add_history("set_restfreqs", varlist)
1429
1430 def shift_refpix(self, delta):
1431 """\
1432 Shift the reference pixel of the Spectra Coordinate by an
1433 integer amount.
1434
1435 Parameters:
1436
1437 delta: the amount to shift by
1438
1439 *Note*:
1440
1441 Be careful using this with broadband data.
1442
1443 """
1444 Scantable.shift_refpix(self, delta)
1445
1446 @asaplog_post_dec
1447 def history(self, filename=None):
1448 """\
1449 Print the history. Optionally to a file.
1450
1451 Parameters:
1452
1453 filename: The name of the file to save the history to.
1454
1455 """
1456 hist = list(self._gethistory())
1457 out = "-"*80
1458 for h in hist:
1459 if h.startswith("---"):
1460 out = "\n".join([out, h])
1461 else:
1462 items = h.split("##")
1463 date = items[0]
1464 func = items[1]
1465 items = items[2:]
1466 out += "\n"+date+"\n"
1467 out += "Function: %s\n Parameters:" % (func)
1468 for i in items:
1469 if i == '':
1470 continue
1471 s = i.split("=")
1472 out += "\n %s = %s" % (s[0], s[1])
1473 out = "\n".join([out, "-"*80])
1474 if filename is not None:
1475 if filename is "":
1476 filename = 'scantable_history.txt'
1477 import os
1478 filename = os.path.expandvars(os.path.expanduser(filename))
1479 if not os.path.isdir(filename):
1480 data = open(filename, 'w')
1481 data.write(out)
1482 data.close()
1483 else:
1484 msg = "Illegal file name '%s'." % (filename)
1485 raise IOError(msg)
1486 return page(out)
1487 #
1488 # Maths business
1489 #
1490 @asaplog_post_dec
1491 def average_time(self, mask=None, scanav=False, weight='tint', align=False):
1492 """\
1493 Return the (time) weighted average of a scan.
1494
1495 *Note*:
1496
1497 in channels only - align if necessary
1498
1499 Parameters:
1500
1501 mask: an optional mask (only used for 'var' and 'tsys'
1502 weighting)
1503
1504 scanav: True averages each scan separately
1505 False (default) averages all scans together,
1506
1507 weight: Weighting scheme.
1508 'none' (mean no weight)
1509 'var' (1/var(spec) weighted)
1510 'tsys' (1/Tsys**2 weighted)
1511 'tint' (integration time weighted)
1512 'tintsys' (Tint/Tsys**2)
1513 'median' ( median averaging)
1514 The default is 'tint'
1515
1516 align: align the spectra in velocity before averaging. It takes
1517 the time of the first spectrum as reference time.
1518
1519 Example::
1520
1521 # time average the scantable without using a mask
1522 newscan = scan.average_time()
1523
1524 """
1525 varlist = vars()
1526 weight = weight or 'TINT'
1527 mask = mask or ()
1528 scanav = (scanav and 'SCAN') or 'NONE'
1529 scan = (self, )
1530
1531 if align:
1532 scan = (self.freq_align(insitu=False), )
1533 s = None
1534 if weight.upper() == 'MEDIAN':
1535 s = scantable(self._math._averagechannel(scan[0], 'MEDIAN',
1536 scanav))
1537 else:
1538 s = scantable(self._math._average(scan, mask, weight.upper(),
1539 scanav))
1540 s._add_history("average_time", varlist)
1541 return s
1542
1543 @asaplog_post_dec
1544 def convert_flux(self, jyperk=None, eta=None, d=None, insitu=None):
1545 """\
1546 Return a scan where all spectra are converted to either
1547 Jansky or Kelvin depending upon the flux units of the scan table.
1548 By default the function tries to look the values up internally.
1549 If it can't find them (or if you want to over-ride), you must
1550 specify EITHER jyperk OR eta (and D which it will try to look up
1551 also if you don't set it). jyperk takes precedence if you set both.
1552
1553 Parameters:
1554
1555 jyperk: the Jy / K conversion factor
1556
1557 eta: the aperture efficiency
1558
1559 d: the geometric diameter (metres)
1560
1561 insitu: if False a new scantable is returned.
1562 Otherwise, the scaling is done in-situ
1563 The default is taken from .asaprc (False)
1564
1565 """
1566 if insitu is None: insitu = rcParams['insitu']
1567 self._math._setinsitu(insitu)
1568 varlist = vars()
1569 jyperk = jyperk or -1.0
1570 d = d or -1.0
1571 eta = eta or -1.0
1572 s = scantable(self._math._convertflux(self, d, eta, jyperk))
1573 s._add_history("convert_flux", varlist)
1574 if insitu: self._assign(s)
1575 else: return s
1576
1577 @asaplog_post_dec
1578 def gain_el(self, poly=None, filename="", method="linear", insitu=None):
1579 """\
1580 Return a scan after applying a gain-elevation correction.
1581 The correction can be made via either a polynomial or a
1582 table-based interpolation (and extrapolation if necessary).
1583 You specify polynomial coefficients, an ascii table or neither.
1584 If you specify neither, then a polynomial correction will be made
1585 with built in coefficients known for certain telescopes (an error
1586 will occur if the instrument is not known).
1587 The data and Tsys are *divided* by the scaling factors.
1588
1589 Parameters:
1590
1591 poly: Polynomial coefficients (default None) to compute a
1592 gain-elevation correction as a function of
1593 elevation (in degrees).
1594
1595 filename: The name of an ascii file holding correction factors.
1596 The first row of the ascii file must give the column
1597 names and these MUST include columns
1598 "ELEVATION" (degrees) and "FACTOR" (multiply data
1599 by this) somewhere.
1600 The second row must give the data type of the
1601 column. Use 'R' for Real and 'I' for Integer.
1602 An example file would be
1603 (actual factors are arbitrary) :
1604
1605 TIME ELEVATION FACTOR
1606 R R R
1607 0.1 0 0.8
1608 0.2 20 0.85
1609 0.3 40 0.9
1610 0.4 60 0.85
1611 0.5 80 0.8
1612 0.6 90 0.75
1613
1614 method: Interpolation method when correcting from a table.
1615 Values are "nearest", "linear" (default), "cubic"
1616 and "spline"
1617
1618 insitu: if False a new scantable is returned.
1619 Otherwise, the scaling is done in-situ
1620 The default is taken from .asaprc (False)
1621
1622 """
1623
1624 if insitu is None: insitu = rcParams['insitu']
1625 self._math._setinsitu(insitu)
1626 varlist = vars()
1627 poly = poly or ()
1628 from os.path import expandvars
1629 filename = expandvars(filename)
1630 s = scantable(self._math._gainel(self, poly, filename, method))
1631 s._add_history("gain_el", varlist)
1632 if insitu:
1633 self._assign(s)
1634 else:
1635 return s
1636
1637 @asaplog_post_dec
1638 def freq_align(self, reftime=None, method='cubic', insitu=None):
1639 """\
1640 Return a scan where all rows have been aligned in frequency/velocity.
1641 The alignment frequency frame (e.g. LSRK) is that set by function
1642 set_freqframe.
1643
1644 Parameters:
1645
1646 reftime: reference time to align at. By default, the time of
1647 the first row of data is used.
1648
1649 method: Interpolation method for regridding the spectra.
1650 Choose from "nearest", "linear", "cubic" (default)
1651 and "spline"
1652
1653 insitu: if False a new scantable is returned.
1654 Otherwise, the scaling is done in-situ
1655 The default is taken from .asaprc (False)
1656
1657 """
1658 if insitu is None: insitu = rcParams["insitu"]
1659 self._math._setinsitu(insitu)
1660 varlist = vars()
1661 reftime = reftime or ""
1662 s = scantable(self._math._freq_align(self, reftime, method))
1663 s._add_history("freq_align", varlist)
1664 if insitu: self._assign(s)
1665 else: return s
1666
1667 @asaplog_post_dec
1668 def opacity(self, tau=None, insitu=None):
1669 """\
1670 Apply an opacity correction. The data
1671 and Tsys are multiplied by the correction factor.
1672
1673 Parameters:
1674
1675 tau: (list of) opacity from which the correction factor is
1676 exp(tau*ZD)
1677 where ZD is the zenith-distance.
1678 If a list is provided, it has to be of length nIF,
1679 nIF*nPol or 1 and in order of IF/POL, e.g.
1680 [opif0pol0, opif0pol1, opif1pol0 ...]
1681 if tau is `None` the opacities are determined from a
1682 model.
1683
1684 insitu: if False a new scantable is returned.
1685 Otherwise, the scaling is done in-situ
1686 The default is taken from .asaprc (False)
1687
1688 """
1689 if insitu is None: insitu = rcParams['insitu']
1690 self._math._setinsitu(insitu)
1691 varlist = vars()
1692 if not hasattr(tau, "__len__"):
1693 tau = [tau]
1694 s = scantable(self._math._opacity(self, tau))
1695 s._add_history("opacity", varlist)
1696 if insitu: self._assign(s)
1697 else: return s
1698
1699 @asaplog_post_dec
1700 def bin(self, width=5, insitu=None):
1701 """\
1702 Return a scan where all spectra have been binned up.
1703
1704 Parameters:
1705
1706 width: The bin width (default=5) in pixels
1707
1708 insitu: if False a new scantable is returned.
1709 Otherwise, the scaling is done in-situ
1710 The default is taken from .asaprc (False)
1711
1712 """
1713 if insitu is None: insitu = rcParams['insitu']
1714 self._math._setinsitu(insitu)
1715 varlist = vars()
1716 s = scantable(self._math._bin(self, width))
1717 s._add_history("bin", varlist)
1718 if insitu:
1719 self._assign(s)
1720 else:
1721 return s
1722
1723 @asaplog_post_dec
1724 def resample(self, width=5, method='cubic', insitu=None):
1725 """\
1726 Return a scan where all spectra have been binned up.
1727
1728 Parameters:
1729
1730 width: The bin width (default=5) in pixels
1731
1732 method: Interpolation method when correcting from a table.
1733 Values are "nearest", "linear", "cubic" (default)
1734 and "spline"
1735
1736 insitu: if False a new scantable is returned.
1737 Otherwise, the scaling is done in-situ
1738 The default is taken from .asaprc (False)
1739
1740 """
1741 if insitu is None: insitu = rcParams['insitu']
1742 self._math._setinsitu(insitu)
1743 varlist = vars()
1744 s = scantable(self._math._resample(self, method, width))
1745 s._add_history("resample", varlist)
1746 if insitu: self._assign(s)
1747 else: return s
1748
1749 @asaplog_post_dec
1750 def average_pol(self, mask=None, weight='none'):
1751 """\
1752 Average the Polarisations together.
1753
1754 Parameters:
1755
1756 mask: An optional mask defining the region, where the
1757 averaging will be applied. The output will have all
1758 specified points masked.
1759
1760 weight: Weighting scheme. 'none' (default), 'var' (1/var(spec)
1761 weighted), or 'tsys' (1/Tsys**2 weighted)
1762
1763 """
1764 varlist = vars()
1765 mask = mask or ()
1766 s = scantable(self._math._averagepol(self, mask, weight.upper()))
1767 s._add_history("average_pol", varlist)
1768 return s
1769
1770 @asaplog_post_dec
1771 def average_beam(self, mask=None, weight='none'):
1772 """\
1773 Average the Beams together.
1774
1775 Parameters:
1776 mask: An optional mask defining the region, where the
1777 averaging will be applied. The output will have all
1778 specified points masked.
1779
1780 weight: Weighting scheme. 'none' (default), 'var' (1/var(spec)
1781 weighted), or 'tsys' (1/Tsys**2 weighted)
1782
1783 """
1784 varlist = vars()
1785 mask = mask or ()
1786 s = scantable(self._math._averagebeams(self, mask, weight.upper()))
1787 s._add_history("average_beam", varlist)
1788 return s
1789
1790 def parallactify(self, pflag):
1791 """\
1792 Set a flag to indicate whether this data should be treated as having
1793 been 'parallactified' (total phase == 0.0)
1794
1795 Parameters:
1796
1797 pflag: Bool indicating whether to turn this on (True) or
1798 off (False)
1799
1800 """
1801 varlist = vars()
1802 self._parallactify(pflag)
1803 self._add_history("parallactify", varlist)
1804
1805 @asaplog_post_dec
1806 def convert_pol(self, poltype=None):
1807 """\
1808 Convert the data to a different polarisation type.
1809 Note that you will need cross-polarisation terms for most conversions.
1810
1811 Parameters:
1812
1813 poltype: The new polarisation type. Valid types are:
1814 "linear", "circular", "stokes" and "linpol"
1815
1816 """
1817 varlist = vars()
1818 s = scantable(self._math._convertpol(self, poltype))
1819 s._add_history("convert_pol", varlist)
1820 return s
1821
1822 @asaplog_post_dec
1823 def smooth(self, kernel="hanning", width=5.0, order=2, plot=False, insitu=None):
1824 """\
1825 Smooth the spectrum by the specified kernel (conserving flux).
1826
1827 Parameters:
1828
1829 kernel: The type of smoothing kernel. Select from
1830 'hanning' (default), 'gaussian', 'boxcar', 'rmedian'
1831 or 'poly'
1832
1833 width: The width of the kernel in pixels. For hanning this is
1834 ignored otherwise it defauls to 5 pixels.
1835 For 'gaussian' it is the Full Width Half
1836 Maximum. For 'boxcar' it is the full width.
1837 For 'rmedian' and 'poly' it is the half width.
1838
1839 order: Optional parameter for 'poly' kernel (default is 2), to
1840 specify the order of the polnomial. Ignored by all other
1841 kernels.
1842
1843 plot: plot the original and the smoothed spectra.
1844 In this each indivual fit has to be approved, by
1845 typing 'y' or 'n'
1846
1847 insitu: if False a new scantable is returned.
1848 Otherwise, the scaling is done in-situ
1849 The default is taken from .asaprc (False)
1850
1851 """
1852 if insitu is None: insitu = rcParams['insitu']
1853 self._math._setinsitu(insitu)
1854 varlist = vars()
1855
1856 if plot: orgscan = self.copy()
1857
1858 s = scantable(self._math._smooth(self, kernel.lower(), width, order))
1859 s._add_history("smooth", varlist)
1860
1861 if plot:
1862 if rcParams['plotter.gui']:
1863 from asap.asaplotgui import asaplotgui as asaplot
1864 else:
1865 from asap.asaplot import asaplot
1866 self._p=asaplot()
1867 self._p.set_panels()
1868 ylab=s._get_ordinate_label()
1869 #self._p.palette(0,["#777777","red"])
1870 for r in xrange(s.nrow()):
1871 xsm=s._getabcissa(r)
1872 ysm=s._getspectrum(r)
1873 xorg=orgscan._getabcissa(r)
1874 yorg=orgscan._getspectrum(r)
1875 self._p.clear()
1876 self._p.hold()
1877 self._p.set_axes('ylabel',ylab)
1878 self._p.set_axes('xlabel',s._getabcissalabel(r))
1879 self._p.set_axes('title',s._getsourcename(r))
1880 self._p.set_line(label='Original',color="#777777")
1881 self._p.plot(xorg,yorg)
1882 self._p.set_line(label='Smoothed',color="red")
1883 self._p.plot(xsm,ysm)
1884 ### Ugly part for legend
1885 for i in [0,1]:
1886 self._p.subplots[0]['lines'].append([self._p.subplots[0]['axes'].lines[i]])
1887 self._p.release()
1888 ### Ugly part for legend
1889 self._p.subplots[0]['lines']=[]
1890 res = raw_input("Accept smoothing ([y]/n): ")
1891 if res.upper() == 'N':
1892 s._setspectrum(yorg, r)
1893 self._p.unmap()
1894 self._p = None
1895 del orgscan
1896
1897 if insitu: self._assign(s)
1898 else: return s
1899
1900 @asaplog_post_dec
1901 def old_poly_baseline(self, mask=None, order=0, plot=False, uselin=False, insitu=None, rows=None):
1902 """\
1903 Return a scan which has been baselined (all rows) by a polynomial.
1904
1905 Parameters:
1906
1907 mask: an optional mask
1908
1909 order: the order of the polynomial (default is 0)
1910
1911 plot: plot the fit and the residual. In this each
1912 indivual fit has to be approved, by typing 'y'
1913 or 'n'
1914
1915 uselin: use linear polynomial fit
1916
1917 insitu: if False a new scantable is returned.
1918 Otherwise, the scaling is done in-situ
1919 The default is taken from .asaprc (False)
1920
1921 rows: row numbers of spectra to be processed.
1922 (default is None: for all rows)
1923
1924 Example:
1925 # return a scan baselined by a third order polynomial,
1926 # not using a mask
1927 bscan = scan.poly_baseline(order=3)
1928
1929 """
1930 if insitu is None: insitu = rcParams['insitu']
1931 if not insitu:
1932 workscan = self.copy()
1933 else:
1934 workscan = self
1935 varlist = vars()
1936 if mask is None:
1937 mask = [True for i in xrange(self.nchan())]
1938
1939 try:
1940 f = fitter()
1941 if uselin:
1942 f.set_function(lpoly=order)
1943 else:
1944 f.set_function(poly=order)
1945
1946 if rows == None:
1947 rows = xrange(workscan.nrow())
1948 elif isinstance(rows, int):
1949 rows = [ rows ]
1950
1951 if len(rows) > 0:
1952 self.blpars = []
1953 self.masklists = []
1954 self.actualmask = []
1955
1956 for r in rows:
1957 f.x = workscan._getabcissa(r)
1958 f.y = workscan._getspectrum(r)
1959 f.mask = mask_and(mask, workscan._getmask(r)) # (CAS-1434)
1960 f.data = None
1961 f.fit()
1962 if plot:
1963 f.plot(residual=True)
1964 x = raw_input("Accept fit ( [y]/n ): ")
1965 if x.upper() == 'N':
1966 self.blpars.append(None)
1967 self.masklists.append(None)
1968 self.actualmask.append(None)
1969 continue
1970 workscan._setspectrum(f.fitter.getresidual(), r)
1971 self.blpars.append(f.get_parameters())
1972 self.masklists.append(workscan.get_masklist(f.mask, row=r, silent=True))
1973 self.actualmask.append(f.mask)
1974
1975 if plot:
1976 f._p.unmap()
1977 f._p = None
1978 workscan._add_history("poly_baseline", varlist)
1979 if insitu:
1980 self._assign(workscan)
1981 else:
1982 return workscan
1983 except RuntimeError:
1984 msg = "The fit failed, possibly because it didn't converge."
1985 raise RuntimeError(msg)
1986
1987 @asaplog_post_dec
1988 def poly_baseline(self, mask=None, order=0, plot=False, batch=False, insitu=None, rows=None):
1989 """\
1990 Return a scan which has been baselined (all rows) by a polynomial.
1991 Parameters:
1992 mask: an optional mask
1993 order: the order of the polynomial (default is 0)
1994 plot: plot the fit and the residual. In this each
1995 indivual fit has to be approved, by typing 'y'
1996 or 'n'. Ignored if batch = True.
1997 batch: if True a faster algorithm is used and logs
1998 including the fit results are not output
1999 (default is False)
2000 insitu: if False a new scantable is returned.
2001 Otherwise, the scaling is done in-situ
2002 The default is taken from .asaprc (False)
2003 rows: row numbers of spectra to be baselined.
2004 (default is None: for all rows)
2005 Example:
2006 # return a scan baselined by a third order polynomial,
2007 # not using a mask
2008 bscan = scan.poly_baseline(order=3)
2009 """
2010
2011 varlist = vars()
2012
2013 if insitu is None: insitu = rcParams["insitu"]
2014 if insitu:
2015 workscan = self
2016 else:
2017 workscan = self.copy()
2018
2019 nchan = workscan.nchan()
2020
2021 if mask is None:
2022 mask = [True for i in xrange(nchan)]
2023
2024 try:
2025 if rows == None:
2026 rows = xrange(workscan.nrow())
2027 elif isinstance(rows, int):
2028 rows = [ rows ]
2029
2030 if len(rows) > 0:
2031 workscan.blpars = []
2032 workscan.masklists = []
2033 workscan.actualmask = []
2034
2035 if batch:
2036 workscan._poly_baseline_batch(mask, order)
2037 elif plot:
2038 f = fitter()
2039 f.set_function(lpoly=order)
2040 for r in rows:
2041 f.x = workscan._getabcissa(r)
2042 f.y = workscan._getspectrum(r)
2043 f.mask = mask_and(mask, workscan._getmask(r)) # (CAS-1434)
2044 f.data = None
2045 f.fit()
2046
2047 f.plot(residual=True)
2048 accept_fit = raw_input("Accept fit ( [y]/n ): ")
2049 if accept_fit.upper() == "N":
2050 self.blpars.append(None)
2051 self.masklists.append(None)
2052 self.actualmask.append(None)
2053 continue
2054 workscan._setspectrum(f.fitter.getresidual(), r)
2055 workscan.blpars.append(f.get_parameters())
2056 workscan.masklists.append(workscan.get_masklist(f.mask, row=r))
2057 workscan.actualmask.append(f.mask)
2058
2059 f._p.unmap()
2060 f._p = None
2061 else:
2062 for r in rows:
2063 fitparams = workscan._poly_baseline(mask, order, r)
2064 params = fitparams.getparameters()
2065 fmtd = ", ".join(["p%d = %3.6f" % (i, v) for i, v in enumerate(params)])
2066 errors = fitparams.geterrors()
2067 fmask = mask_and(mask, workscan._getmask(r))
2068
2069 workscan.blpars.append({"params":params,
2070 "fixed": fitparams.getfixedparameters(),
2071 "formatted":fmtd, "errors":errors})
2072 workscan.masklists.append(workscan.get_masklist(fmask, r, silent=True))
2073 workscan.actualmask.append(fmask)
2074
2075 asaplog.push(fmtd)
2076
2077 workscan._add_history("poly_baseline", varlist)
2078
2079 if insitu:
2080 self._assign(workscan)
2081 else:
2082 return workscan
2083
2084 except RuntimeError, e:
2085 msg = "The fit failed, possibly because it didn't converge."
2086 if rcParams["verbose"]:
2087 asaplog.push(str(e))
2088 asaplog.push(str(msg))
2089 return
2090 else:
2091 raise RuntimeError(str(e)+'\n'+msg)
2092
2093
2094 def auto_poly_baseline(self, mask=None, edge=(0, 0), order=0,
2095 threshold=3, chan_avg_limit=1, plot=False,
2096 insitu=None, rows=None):
2097 """\
2098 Return a scan which has been baselined (all rows) by a polynomial.
2099 Spectral lines are detected first using linefinder and masked out
2100 to avoid them affecting the baseline solution.
2101
2102 Parameters:
2103
2104 mask: an optional mask retreived from scantable
2105
2106 edge: an optional number of channel to drop at the edge of
2107 spectrum. If only one value is
2108 specified, the same number will be dropped from
2109 both sides of the spectrum. Default is to keep
2110 all channels. Nested tuples represent individual
2111 edge selection for different IFs (a number of spectral
2112 channels can be different)
2113
2114 order: the order of the polynomial (default is 0)
2115
2116 threshold: the threshold used by line finder. It is better to
2117 keep it large as only strong lines affect the
2118 baseline solution.
2119
2120 chan_avg_limit:
2121 a maximum number of consequtive spectral channels to
2122 average during the search of weak and broad lines.
2123 The default is no averaging (and no search for weak
2124 lines). If such lines can affect the fitted baseline
2125 (e.g. a high order polynomial is fitted), increase this
2126 parameter (usually values up to 8 are reasonable). Most
2127 users of this method should find the default value
2128 sufficient.
2129
2130 plot: plot the fit and the residual. In this each
2131 indivual fit has to be approved, by typing 'y'
2132 or 'n'
2133
2134 insitu: if False a new scantable is returned.
2135 Otherwise, the scaling is done in-situ
2136 The default is taken from .asaprc (False)
2137 rows: row numbers of spectra to be processed.
2138 (default is None: for all rows)
2139
2140
2141 Example::
2142
2143 scan2 = scan.auto_poly_baseline(order=7, insitu=False)
2144
2145 """
2146 if insitu is None: insitu = rcParams['insitu']
2147 varlist = vars()
2148 from asap.asaplinefind import linefinder
2149 from asap import _is_sequence_or_number as _is_valid
2150
2151 # check whether edge is set up for each IF individually
2152 individualedge = False;
2153 if len(edge) > 1:
2154 if isinstance(edge[0], list) or isinstance(edge[0], tuple):
2155 individualedge = True;
2156
2157 if not _is_valid(edge, int) and not individualedge:
2158 raise ValueError, "Parameter 'edge' has to be an integer or a \
2159 pair of integers specified as a tuple. Nested tuples are allowed \
2160 to make individual selection for different IFs."
2161
2162 curedge = (0, 0)
2163 if individualedge:
2164 for edgepar in edge:
2165 if not _is_valid(edgepar, int):
2166 raise ValueError, "Each element of the 'edge' tuple has \
2167 to be a pair of integers or an integer."
2168 else:
2169 curedge = edge;
2170
2171 if not insitu:
2172 workscan = self.copy()
2173 else:
2174 workscan = self
2175
2176 # setup fitter
2177 f = fitter()
2178 f.set_function(lpoly=order)
2179
2180 # setup line finder
2181 fl = linefinder()
2182 fl.set_options(threshold=threshold,avg_limit=chan_avg_limit)
2183
2184 fl.set_scan(workscan)
2185
2186 if mask is None:
2187 mask = _n_bools(workscan.nchan(), True)
2188
2189 if rows is None:
2190 rows = xrange(workscan.nrow())
2191 elif isinstance(rows, int):
2192 rows = [ rows ]
2193
2194 # Save parameters of baseline fits & masklists as a class attribute.
2195 # NOTICE: It does not reflect changes in scantable!
2196 if len(rows) > 0:
2197 self.blpars=[]
2198 self.masklists=[]
2199 self.actualmask=[]
2200 asaplog.push("Processing:")
2201 for r in rows:
2202 msg = " Scan[%d] Beam[%d] IF[%d] Pol[%d] Cycle[%d]" % \
2203 (workscan.getscan(r), workscan.getbeam(r), workscan.getif(r), \
2204 workscan.getpol(r), workscan.getcycle(r))
2205 asaplog.push(msg, False)
2206
2207 # figure out edge parameter
2208 if individualedge:
2209 if len(edge) >= workscan.getif(r):
2210 raise RuntimeError, "Number of edge elements appear to " \
2211 "be less than the number of IFs"
2212 curedge = edge[workscan.getif(r)]
2213
2214 actualmask = mask_and(mask, workscan._getmask(r)) # (CAS-1434)
2215
2216 # setup line finder
2217 fl.find_lines(r, actualmask, curedge)
2218
2219 f.x = workscan._getabcissa(r)
2220 f.y = workscan._getspectrum(r)
2221 f.mask = fl.get_mask()
2222 f.data = None
2223 f.fit()
2224
2225 # Show mask list
2226 masklist=workscan.get_masklist(f.mask, row=r, silent=True)
2227 msg = "mask range: "+str(masklist)
2228 asaplog.push(msg, False)
2229
2230 if plot:
2231 f.plot(residual=True)
2232 x = raw_input("Accept fit ( [y]/n ): ")
2233 if x.upper() == 'N':
2234 self.blpars.append(None)
2235 self.masklists.append(None)
2236 self.actualmask.append(None)
2237 continue
2238
2239 workscan._setspectrum(f.fitter.getresidual(), r)
2240 self.blpars.append(f.get_parameters())
2241 self.masklists.append(masklist)
2242 self.actualmask.append(f.mask)
2243 if plot:
2244 f._p.unmap()
2245 f._p = None
2246 workscan._add_history("auto_poly_baseline", varlist)
2247 if insitu:
2248 self._assign(workscan)
2249 else:
2250 return workscan
2251
2252 @asaplog_post_dec
2253 def rotate_linpolphase(self, angle):
2254 """\
2255 Rotate the phase of the complex polarization O=Q+iU correlation.
2256 This is always done in situ in the raw data. So if you call this
2257 function more than once then each call rotates the phase further.
2258
2259 Parameters:
2260
2261 angle: The angle (degrees) to rotate (add) by.
2262
2263 Example::
2264
2265 scan.rotate_linpolphase(2.3)
2266
2267 """
2268 varlist = vars()
2269 self._math._rotate_linpolphase(self, angle)
2270 self._add_history("rotate_linpolphase", varlist)
2271 return
2272
2273 @asaplog_post_dec
2274 def rotate_xyphase(self, angle):
2275 """\
2276 Rotate the phase of the XY correlation. This is always done in situ
2277 in the data. So if you call this function more than once
2278 then each call rotates the phase further.
2279
2280 Parameters:
2281
2282 angle: The angle (degrees) to rotate (add) by.
2283
2284 Example::
2285
2286 scan.rotate_xyphase(2.3)
2287
2288 """
2289 varlist = vars()
2290 self._math._rotate_xyphase(self, angle)
2291 self._add_history("rotate_xyphase", varlist)
2292 return
2293
2294 @asaplog_post_dec
2295 def swap_linears(self):
2296 """\
2297 Swap the linear polarisations XX and YY, or better the first two
2298 polarisations as this also works for ciculars.
2299 """
2300 varlist = vars()
2301 self._math._swap_linears(self)
2302 self._add_history("swap_linears", varlist)
2303 return
2304
2305 @asaplog_post_dec
2306 def invert_phase(self):
2307 """\
2308 Invert the phase of the complex polarisation
2309 """
2310 varlist = vars()
2311 self._math._invert_phase(self)
2312 self._add_history("invert_phase", varlist)
2313 return
2314
2315 @asaplog_post_dec
2316 def add(self, offset, insitu=None):
2317 """\
2318 Return a scan where all spectra have the offset added
2319
2320 Parameters:
2321
2322 offset: the offset
2323
2324 insitu: if False a new scantable is returned.
2325 Otherwise, the scaling is done in-situ
2326 The default is taken from .asaprc (False)
2327
2328 """
2329 if insitu is None: insitu = rcParams['insitu']
2330 self._math._setinsitu(insitu)
2331 varlist = vars()
2332 s = scantable(self._math._unaryop(self, offset, "ADD", False))
2333 s._add_history("add", varlist)
2334 if insitu:
2335 self._assign(s)
2336 else:
2337 return s
2338
2339 @asaplog_post_dec
2340 def scale(self, factor, tsys=True, insitu=None):
2341 """\
2342
2343 Return a scan where all spectra are scaled by the given 'factor'
2344
2345 Parameters:
2346
2347 factor: the scaling factor (float or 1D float list)
2348
2349 insitu: if False a new scantable is returned.
2350 Otherwise, the scaling is done in-situ
2351 The default is taken from .asaprc (False)
2352
2353 tsys: if True (default) then apply the operation to Tsys
2354 as well as the data
2355
2356 """
2357 if insitu is None: insitu = rcParams['insitu']
2358 self._math._setinsitu(insitu)
2359 varlist = vars()
2360 s = None
2361 import numpy
2362 if isinstance(factor, list) or isinstance(factor, numpy.ndarray):
2363 if isinstance(factor[0], list) or isinstance(factor[0], numpy.ndarray):
2364 from asapmath import _array2dOp
2365 s = _array2dOp( self.copy(), factor, "MUL", tsys )
2366 else:
2367 s = scantable( self._math._arrayop( self.copy(), factor, "MUL", tsys ) )
2368 else:
2369 s = scantable(self._math._unaryop(self.copy(), factor, "MUL", tsys))
2370 s._add_history("scale", varlist)
2371 if insitu:
2372 self._assign(s)
2373 else:
2374 return s
2375
2376 def set_sourcetype(self, match, matchtype="pattern",
2377 sourcetype="reference"):
2378 """\
2379 Set the type of the source to be an source or reference scan
2380 using the provided pattern.
2381
2382 Parameters:
2383
2384 match: a Unix style pattern, regular expression or selector
2385
2386 matchtype: 'pattern' (default) UNIX style pattern or
2387 'regex' regular expression
2388
2389 sourcetype: the type of the source to use (source/reference)
2390
2391 """
2392 varlist = vars()
2393 basesel = self.get_selection()
2394 stype = -1
2395 if sourcetype.lower().startswith("r"):
2396 stype = 1
2397 elif sourcetype.lower().startswith("s"):
2398 stype = 0
2399 else:
2400 raise ValueError("Illegal sourcetype use s(ource) or r(eference)")
2401 if matchtype.lower().startswith("p"):
2402 matchtype = "pattern"
2403 elif matchtype.lower().startswith("r"):
2404 matchtype = "regex"
2405 else:
2406 raise ValueError("Illegal matchtype, use p(attern) or r(egex)")
2407 sel = selector()
2408 if isinstance(match, selector):
2409 sel = match
2410 else:
2411 sel.set_query("SRCNAME == %s('%s')" % (matchtype, match))
2412 self.set_selection(basesel+sel)
2413 self._setsourcetype(stype)
2414 self.set_selection(basesel)
2415 self._add_history("set_sourcetype", varlist)
2416
2417 @asaplog_post_dec
2418 @preserve_selection
2419 def auto_quotient(self, preserve=True, mode='paired', verify=False):
2420 """\
2421 This function allows to build quotients automatically.
2422 It assumes the observation to have the same number of
2423 "ons" and "offs"
2424
2425 Parameters:
2426
2427 preserve: you can preserve (default) the continuum or
2428 remove it. The equations used are
2429
2430 preserve: Output = Toff * (on/off) - Toff
2431
2432 remove: Output = Toff * (on/off) - Ton
2433
2434 mode: the on/off detection mode
2435 'paired' (default)
2436 identifies 'off' scans by the
2437 trailing '_R' (Mopra/Parkes) or
2438 '_e'/'_w' (Tid) and matches
2439 on/off pairs from the observing pattern
2440 'time'
2441 finds the closest off in time
2442
2443 .. todo:: verify argument is not implemented
2444
2445 """
2446 varlist = vars()
2447 modes = ["time", "paired"]
2448 if not mode in modes:
2449 msg = "please provide valid mode. Valid modes are %s" % (modes)
2450 raise ValueError(msg)
2451 s = None
2452 if mode.lower() == "paired":
2453 sel = self.get_selection()
2454 sel.set_query("SRCTYPE==psoff")
2455 self.set_selection(sel)
2456 offs = self.copy()
2457 sel.set_query("SRCTYPE==pson")
2458 self.set_selection(sel)
2459 ons = self.copy()
2460 s = scantable(self._math._quotient(ons, offs, preserve))
2461 elif mode.lower() == "time":
2462 s = scantable(self._math._auto_quotient(self, mode, preserve))
2463 s._add_history("auto_quotient", varlist)
2464 return s
2465
2466 @asaplog_post_dec
2467 def mx_quotient(self, mask = None, weight='median', preserve=True):
2468 """\
2469 Form a quotient using "off" beams when observing in "MX" mode.
2470
2471 Parameters:
2472
2473 mask: an optional mask to be used when weight == 'stddev'
2474
2475 weight: How to average the off beams. Default is 'median'.
2476
2477 preserve: you can preserve (default) the continuum or
2478 remove it. The equations used are:
2479
2480 preserve: Output = Toff * (on/off) - Toff
2481
2482 remove: Output = Toff * (on/off) - Ton
2483
2484 """
2485 mask = mask or ()
2486 varlist = vars()
2487 on = scantable(self._math._mx_extract(self, 'on'))
2488 preoff = scantable(self._math._mx_extract(self, 'off'))
2489 off = preoff.average_time(mask=mask, weight=weight, scanav=False)
2490 from asapmath import quotient
2491 q = quotient(on, off, preserve)
2492 q._add_history("mx_quotient", varlist)
2493 return q
2494
2495 @asaplog_post_dec
2496 def freq_switch(self, insitu=None):
2497 """\
2498 Apply frequency switching to the data.
2499
2500 Parameters:
2501
2502 insitu: if False a new scantable is returned.
2503 Otherwise, the swictching is done in-situ
2504 The default is taken from .asaprc (False)
2505
2506 """
2507 if insitu is None: insitu = rcParams['insitu']
2508 self._math._setinsitu(insitu)
2509 varlist = vars()
2510 s = scantable(self._math._freqswitch(self))
2511 s._add_history("freq_switch", varlist)
2512 if insitu:
2513 self._assign(s)
2514 else:
2515 return s
2516
2517 @asaplog_post_dec
2518 def recalc_azel(self):
2519 """Recalculate the azimuth and elevation for each position."""
2520 varlist = vars()
2521 self._recalcazel()
2522 self._add_history("recalc_azel", varlist)
2523 return
2524
2525 @asaplog_post_dec
2526 def __add__(self, other):
2527 varlist = vars()
2528 s = None
2529 if isinstance(other, scantable):
2530 s = scantable(self._math._binaryop(self, other, "ADD"))
2531 elif isinstance(other, float):
2532 s = scantable(self._math._unaryop(self, other, "ADD", False))
2533 else:
2534 raise TypeError("Other input is not a scantable or float value")
2535 s._add_history("operator +", varlist)
2536 return s
2537
2538 @asaplog_post_dec
2539 def __sub__(self, other):
2540 """
2541 implicit on all axes and on Tsys
2542 """
2543 varlist = vars()
2544 s = None
2545 if isinstance(other, scantable):
2546 s = scantable(self._math._binaryop(self, other, "SUB"))
2547 elif isinstance(other, float):
2548 s = scantable(self._math._unaryop(self, other, "SUB", False))
2549 else:
2550 raise TypeError("Other input is not a scantable or float value")
2551 s._add_history("operator -", varlist)
2552 return s
2553
2554 @asaplog_post_dec
2555 def __mul__(self, other):
2556 """
2557 implicit on all axes and on Tsys
2558 """
2559 varlist = vars()
2560 s = None
2561 if isinstance(other, scantable):
2562 s = scantable(self._math._binaryop(self, other, "MUL"))
2563 elif isinstance(other, float):
2564 s = scantable(self._math._unaryop(self, other, "MUL", False))
2565 else:
2566 raise TypeError("Other input is not a scantable or float value")
2567 s._add_history("operator *", varlist)
2568 return s
2569
2570
2571 @asaplog_post_dec
2572 def __div__(self, other):
2573 """
2574 implicit on all axes and on Tsys
2575 """
2576 varlist = vars()
2577 s = None
2578 if isinstance(other, scantable):
2579 s = scantable(self._math._binaryop(self, other, "DIV"))
2580 elif isinstance(other, float):
2581 if other == 0.0:
2582 raise ZeroDivisionError("Dividing by zero is not recommended")
2583 s = scantable(self._math._unaryop(self, other, "DIV", False))
2584 else:
2585 raise TypeError("Other input is not a scantable or float value")
2586 s._add_history("operator /", varlist)
2587 return s
2588
2589 @asaplog_post_dec
2590 def get_fit(self, row=0):
2591 """\
2592 Print or return the stored fits for a row in the scantable
2593
2594 Parameters:
2595
2596 row: the row which the fit has been applied to.
2597
2598 """
2599 if row > self.nrow():
2600 return
2601 from asap.asapfit import asapfit
2602 fit = asapfit(self._getfit(row))
2603 asaplog.push( '%s' %(fit) )
2604 return fit.as_dict()
2605
2606 def flag_nans(self):
2607 """\
2608 Utility function to flag NaN values in the scantable.
2609 """
2610 import numpy
2611 basesel = self.get_selection()
2612 for i in range(self.nrow()):
2613 sel = self.get_row_selector(i)
2614 self.set_selection(basesel+sel)
2615 nans = numpy.isnan(self._getspectrum(0))
2616 if numpy.any(nans):
2617 bnans = [ bool(v) for v in nans]
2618 self.flag(bnans)
2619 self.set_selection(basesel)
2620
2621 def get_row_selector(self, rowno):
2622 #return selector(beams=self.getbeam(rowno),
2623 # ifs=self.getif(rowno),
2624 # pols=self.getpol(rowno),
2625 # scans=self.getscan(rowno),
2626 # cycles=self.getcycle(rowno))
2627 return selector(rows=[rowno])
2628
2629 def _add_history(self, funcname, parameters):
2630 if not rcParams['scantable.history']:
2631 return
2632 # create date
2633 sep = "##"
2634 from datetime import datetime
2635 dstr = datetime.now().strftime('%Y/%m/%d %H:%M:%S')
2636 hist = dstr+sep
2637 hist += funcname+sep#cdate+sep
2638 if parameters.has_key('self'): del parameters['self']
2639 for k, v in parameters.iteritems():
2640 if type(v) is dict:
2641 for k2, v2 in v.iteritems():
2642 hist += k2
2643 hist += "="
2644 if isinstance(v2, scantable):
2645 hist += 'scantable'
2646 elif k2 == 'mask':
2647 if isinstance(v2, list) or isinstance(v2, tuple):
2648 hist += str(self._zip_mask(v2))
2649 else:
2650 hist += str(v2)
2651 else:
2652 hist += str(v2)
2653 else:
2654 hist += k
2655 hist += "="
2656 if isinstance(v, scantable):
2657 hist += 'scantable'
2658 elif k == 'mask':
2659 if isinstance(v, list) or isinstance(v, tuple):
2660 hist += str(self._zip_mask(v))
2661 else:
2662 hist += str(v)
2663 else:
2664 hist += str(v)
2665 hist += sep
2666 hist = hist[:-2] # remove trailing '##'
2667 self._addhistory(hist)
2668
2669
2670 def _zip_mask(self, mask):
2671 mask = list(mask)
2672 i = 0
2673 segments = []
2674 while mask[i:].count(1):
2675 i += mask[i:].index(1)
2676 if mask[i:].count(0):
2677 j = i + mask[i:].index(0)
2678 else:
2679 j = len(mask)
2680 segments.append([i, j])
2681 i = j
2682 return segments
2683
2684 def _get_ordinate_label(self):
2685 fu = "("+self.get_fluxunit()+")"
2686 import re
2687 lbl = "Intensity"
2688 if re.match(".K.", fu):
2689 lbl = "Brightness Temperature "+ fu
2690 elif re.match(".Jy.", fu):
2691 lbl = "Flux density "+ fu
2692 return lbl
2693
2694 def _check_ifs(self):
2695 #nchans = [self.nchan(i) for i in range(self.nif(-1))]
2696 nchans = [self.nchan(i) for i in self.getifnos()]
2697 nchans = filter(lambda t: t > 0, nchans)
2698 return (sum(nchans)/len(nchans) == nchans[0])
2699
2700 @asaplog_post_dec
2701 #def _fill(self, names, unit, average, getpt, antenna):
2702 def _fill(self, names, unit, average, opts={}):
2703 first = True
2704 fullnames = []
2705 for name in names:
2706 name = os.path.expandvars(name)
2707 name = os.path.expanduser(name)
2708 if not os.path.exists(name):
2709 msg = "File '%s' does not exists" % (name)
2710 raise IOError(msg)
2711 fullnames.append(name)
2712 if average:
2713 asaplog.push('Auto averaging integrations')
2714 stype = int(rcParams['scantable.storage'].lower() == 'disk')
2715 for name in fullnames:
2716 tbl = Scantable(stype)
2717 if is_ms( name ):
2718 r = msfiller( tbl )
2719 else:
2720 r = filler( tbl )
2721 rx = rcParams['scantable.reference']
2722 r.setreferenceexpr(rx)
2723 #r = filler(tbl)
2724 #rx = rcParams['scantable.reference']
2725 #r.setreferenceexpr(rx)
2726 msg = "Importing %s..." % (name)
2727 asaplog.push(msg, False)
2728 #opts = {'ms': {'antenna' : antenna, 'getpt': getpt} }
2729 r.open(name, opts)# antenna, -1, -1, getpt)
2730 r.fill()
2731 if average:
2732 tbl = self._math._average((tbl, ), (), 'NONE', 'SCAN')
2733 if not first:
2734 tbl = self._math._merge([self, tbl])
2735 Scantable.__init__(self, tbl)
2736 r.close()
2737 del r, tbl
2738 first = False
2739 #flush log
2740 asaplog.post()
2741 if unit is not None:
2742 self.set_fluxunit(unit)
2743 if not is_casapy():
2744 self.set_freqframe(rcParams['scantable.freqframe'])
2745
2746 def __getitem__(self, key):
2747 if key < 0:
2748 key += self.nrow()
2749 if key >= self.nrow():
2750 raise IndexError("Row index out of range.")
2751 return self._getspectrum(key)
2752
2753 def __setitem__(self, key, value):
2754 if key < 0:
2755 key += self.nrow()
2756 if key >= self.nrow():
2757 raise IndexError("Row index out of range.")
2758 if not hasattr(value, "__len__") or \
2759 len(value) > self.nchan(self.getif(key)):
2760 raise ValueError("Spectrum length doesn't match.")
2761 return self._setspectrum(value, key)
2762
2763 def __len__(self):
2764 return self.nrow()
2765
2766 def __iter__(self):
2767 for i in range(len(self)):
2768 yield self[i]
Note: See TracBrowser for help on using the repository browser.