source: trunk/python/scantable.py@ 2090

Last change on this file since 2090 was 2081, checked in by WataruKawasaki, 14 years ago

New Development: No

JIRA Issue: Yes CAS-2847

Ready for Test: Yes

Interface Changes: No

What Interface Changed:

Test Programs:

Put in Release Notes: No

Module(s): Scantable

Description: Scantable::sinusoidBaseline(), Scantable::autoSinusoidBaseline()


  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 123.3 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 sinusoid_baseline(self, insitu=None, mask=None, nwave=None, maxwavelength=None,
2073 clipthresh=None, clipniter=None, plot=None, getresidual=None, outlog=None, blfile=None):
2074 """\
2075 Return a scan which has been baselined (all rows) by sinusoidal functions.
2076 Parameters:
2077 insitu: If False a new scantable is returned.
2078 Otherwise, the scaling is done in-situ
2079 The default is taken from .asaprc (False)
2080 mask: An optional mask
2081 nwave: the maximum wave number of sinusoids within
2082 maxwavelength * (spectral range).
2083 The default is 3 (i.e., sinusoids with wave
2084 number of 0(=constant), 1, 2, and 3 are
2085 used for fitting). Also it is possible to
2086 explicitly specify all the wave numbers to
2087 be used, by giving a list including them
2088 (e.g. [0,1,2,15,16]).
2089 maxwavelength: the longest sinusoidal wavelength. The
2090 default is 1.0 (unit: spectral range).
2091 clipthresh: Clipping threshold. (default is 3.0, unit: sigma)
2092 clipniter: maximum number of iteration of 'clipthresh'-sigma clipping (default is 1)
2093 plot: *** CURRENTLY UNAVAILABLE, ALWAYS FALSE ***
2094 plot the fit and the residual. In this each
2095 indivual fit has to be approved, by typing 'y'
2096 or 'n'
2097 getresidual: if False, returns best-fit values instead of
2098 residual. (default is True)
2099 outlog: Output the coefficients of the best-fit
2100 function to logger (default is False)
2101 blfile: Name of a text file in which the best-fit
2102 parameter values to be written
2103 (default is "": no file/logger output)
2104
2105 Example:
2106 # return a scan baselined by a combination of sinusoidal curves having
2107 # wave numbers in spectral window up to 10,
2108 # also with 3-sigma clipping, iteration up to 4 times
2109 bscan = scan.sinusoid_baseline(nwave=10,clipthresh=3.0,clipniter=4)
2110
2111 Note:
2112 The best-fit parameter values output in logger and/or blfile are now
2113 based on specunit of 'channel'.
2114 """
2115
2116 varlist = vars()
2117
2118 if insitu is None: insitu = rcParams["insitu"]
2119 if insitu:
2120 workscan = self
2121 else:
2122 workscan = self.copy()
2123
2124 nchan = workscan.nchan()
2125
2126 if mask is None: mask = [True for i in xrange(nchan)]
2127 if nwave is None: nwave = 3
2128 if maxwavelength is None: maxwavelength = 1.0
2129 if clipthresh is None: clipthresh = 3.0
2130 if clipniter is None: clipniter = 1
2131 if plot is None: plot = False
2132 if getresidual is None: getresidual = True
2133 if outlog is None: outlog = False
2134 if blfile is None: blfile = ""
2135
2136 if isinstance(nwave, int):
2137 in_nwave = nwave
2138 nwave = []
2139 for i in xrange(in_nwave+1): nwave.append(i)
2140
2141 outblfile = (blfile != "") and os.path.exists(os.path.expanduser(os.path.expandvars(blfile)))
2142
2143 try:
2144 #CURRENTLY, PLOT=true is UNAVAILABLE UNTIL sinusoidal fitting is implemented as a fitter method.
2145 workscan._sinusoid_baseline(mask, nwave, maxwavelength, clipthresh, clipniter, getresidual, outlog, blfile)
2146
2147 workscan._add_history("sinusoid_baseline", varlist)
2148
2149 if insitu:
2150 self._assign(workscan)
2151 else:
2152 return workscan
2153
2154 except RuntimeError, e:
2155 msg = "The fit failed, possibly because it didn't converge."
2156 if rcParams["verbose"]:
2157 asaplog.push(str(e))
2158 asaplog.push(str(msg))
2159 return
2160 else:
2161 raise RuntimeError(str(e)+'\n'+msg)
2162
2163
2164 def auto_sinusoid_baseline(self, insitu=None, mask=None, nwave=None, maxwavelength=None,
2165 clipthresh=None, clipniter=None, edge=None, threshold=None,
2166 chan_avg_limit=None, plot=None, getresidual=None, outlog=None, blfile=None):
2167 """\
2168 Return a scan which has been baselined (all rows) by cubic spline
2169 function (piecewise cubic polynomial).
2170 Spectral lines are detected first using linefinder and masked out
2171 to avoid them affecting the baseline solution.
2172
2173 Parameters:
2174 insitu: if False a new scantable is returned.
2175 Otherwise, the scaling is done in-situ
2176 The default is taken from .asaprc (False)
2177 mask: an optional mask retreived from scantable
2178 nwave: the maximum wave number of sinusoids within
2179 maxwavelength * (spectral range).
2180 The default is 3 (i.e., sinusoids with wave
2181 number of 0(=constant), 1, 2, and 3 are
2182 used for fitting). Also it is possible to
2183 explicitly specify all the wave numbers to
2184 be used, by giving a list including them
2185 (e.g. [0,1,2,15,16]).
2186 maxwavelength: the longest sinusoidal wavelength. The
2187 default is 1.0 (unit: spectral range).
2188 clipthresh: Clipping threshold. (default is 3.0, unit: sigma)
2189 clipniter: maximum number of iteration of 'clipthresh'-sigma clipping (default is 1)
2190 edge: an optional number of channel to drop at
2191 the edge of spectrum. If only one value is
2192 specified, the same number will be dropped
2193 from both sides of the spectrum. Default
2194 is to keep all channels. Nested tuples
2195 represent individual edge selection for
2196 different IFs (a number of spectral channels
2197 can be different)
2198 threshold: the threshold used by line finder. It is
2199 better to keep it large as only strong lines
2200 affect the baseline solution.
2201 chan_avg_limit:a maximum number of consequtive spectral
2202 channels to average during the search of
2203 weak and broad lines. The default is no
2204 averaging (and no search for weak lines).
2205 If such lines can affect the fitted baseline
2206 (e.g. a high order polynomial is fitted),
2207 increase this parameter (usually values up
2208 to 8 are reasonable). Most users of this
2209 method should find the default value sufficient.
2210 plot: *** CURRENTLY UNAVAILABLE, ALWAYS FALSE ***
2211 plot the fit and the residual. In this each
2212 indivual fit has to be approved, by typing 'y'
2213 or 'n'
2214 getresidual: if False, returns best-fit values instead of
2215 residual. (default is True)
2216 outlog: Output the coefficients of the best-fit
2217 function to logger (default is False)
2218 blfile: Name of a text file in which the best-fit
2219 parameter values to be written
2220 (default is "": no file/logger output)
2221
2222 Example:
2223 bscan = scan.auto_sinusoid_baseline(nwave=10, insitu=False)
2224
2225 Note:
2226 The best-fit parameter values output in logger and/or blfile are now
2227 based on specunit of 'channel'.
2228 """
2229
2230 varlist = vars()
2231
2232 if insitu is None: insitu = rcParams['insitu']
2233 if insitu:
2234 workscan = self
2235 else:
2236 workscan = self.copy()
2237
2238 nchan = workscan.nchan()
2239
2240 if mask is None: mask = [True for i in xrange(nchan)]
2241 if nwave is None: nwave = 3
2242 if maxwavelength is None: maxwavelength = 1.0
2243 if clipthresh is None: clipthresh = 3.0
2244 if clipniter is None: clipniter = 1
2245 if edge is None: edge = (0,0)
2246 if threshold is None: threshold = 3
2247 if chan_avg_limit is None: chan_avg_limit = 1
2248 if plot is None: plot = False
2249 if getresidual is None: getresidual = True
2250 if outlog is None: outlog = False
2251 if blfile is None: blfile = ""
2252
2253 outblfile = (blfile != "") and os.path.exists(os.path.expanduser(os.path.expandvars(blfile)))
2254
2255 from asap.asaplinefind import linefinder
2256 from asap import _is_sequence_or_number as _is_valid
2257
2258 if isinstance(nwave, int):
2259 in_nwave = nwave
2260 nwave = []
2261 for i in xrange(in_nwave+1): nwave.append(i)
2262
2263 if not (isinstance(edge, list) or isinstance(edge, tuple)): edge = [ edge ]
2264 individualedge = False;
2265 if len(edge) > 1: individualedge = isinstance(edge[0], list) or isinstance(edge[0], tuple)
2266
2267 if individualedge:
2268 for edgepar in edge:
2269 if not _is_valid(edgepar, int):
2270 raise ValueError, "Each element of the 'edge' tuple has \
2271 to be a pair of integers or an integer."
2272 else:
2273 if not _is_valid(edge, int):
2274 raise ValueError, "Parameter 'edge' has to be an integer or a \
2275 pair of integers specified as a tuple. \
2276 Nested tuples are allowed \
2277 to make individual selection for different IFs."
2278
2279 if len(edge) > 1:
2280 curedge = edge
2281 else:
2282 curedge = edge + edge
2283
2284 try:
2285 #CURRENTLY, PLOT=true is UNAVAILABLE UNTIL sinusoidal fitting is implemented as a fitter method.
2286 if individualedge:
2287 curedge = []
2288 for i in xrange(len(edge)):
2289 curedge += edge[i]
2290
2291 workscan._auto_sinusoid_baseline(mask, nwave, maxwavelength, clipthresh, clipniter, curedge, threshold, chan_avg_limit, getresidual, outlog, blfile)
2292
2293 workscan._add_history("auto_sinusoid_baseline", varlist)
2294
2295 if insitu:
2296 self._assign(workscan)
2297 else:
2298 return workscan
2299
2300 except RuntimeError, e:
2301 msg = "The fit failed, possibly because it didn't converge."
2302 if rcParams["verbose"]:
2303 asaplog.push(str(e))
2304 asaplog.push(str(msg))
2305 return
2306 else:
2307 raise RuntimeError(str(e)+'\n'+msg)
2308
2309
2310 @asaplog_post_dec
2311 def cspline_baseline(self, insitu=None, mask=None, npiece=None, clipthresh=None, clipniter=None, plot=None, outlog=None, blfile=None):
2312 """\
2313 Return a scan which has been baselined (all rows) by cubic spline function (piecewise cubic polynomial).
2314 Parameters:
2315 insitu: If False a new scantable is returned.
2316 Otherwise, the scaling is done in-situ
2317 The default is taken from .asaprc (False)
2318 mask: An optional mask
2319 npiece: Number of pieces. (default is 2)
2320 clipthresh: Clipping threshold. (default is 3.0, unit: sigma)
2321 clipniter: maximum number of iteration of 'clipthresh'-sigma clipping (default is 1)
2322 plot: *** CURRENTLY UNAVAILABLE, ALWAYS FALSE ***
2323 plot the fit and the residual. In this each
2324 indivual fit has to be approved, by typing 'y'
2325 or 'n'
2326 outlog: Output the coefficients of the best-fit
2327 function to logger (default is False)
2328 blfile: Name of a text file in which the best-fit
2329 parameter values to be written
2330 (default is "": no file/logger output)
2331
2332 Example:
2333 # return a scan baselined by a cubic spline consisting of 2 pieces (i.e., 1 internal knot),
2334 # also with 3-sigma clipping, iteration up to 4 times
2335 bscan = scan.cspline_baseline(npiece=2,clipthresh=3.0,clipniter=4)
2336
2337 Note:
2338 The best-fit parameter values output in logger and/or blfile are now
2339 based on specunit of 'channel'.
2340 """
2341
2342 varlist = vars()
2343
2344 if insitu is None: insitu = rcParams["insitu"]
2345 if insitu:
2346 workscan = self
2347 else:
2348 workscan = self.copy()
2349
2350 nchan = workscan.nchan()
2351
2352 if mask is None: mask = [True for i in xrange(nchan)]
2353 if npiece is None: npiece = 2
2354 if clipthresh is None: clipthresh = 3.0
2355 if clipniter is None: clipniter = 1
2356 if plot is None: plot = False
2357 if outlog is None: outlog = False
2358 if blfile is None: blfile = ""
2359
2360 outblfile = (blfile != "") and os.path.exists(os.path.expanduser(os.path.expandvars(blfile)))
2361
2362 try:
2363 #CURRENTLY, PLOT=true UNAVAILABLE UNTIL cubic spline fitting is implemented as a fitter method.
2364 workscan._cspline_baseline(mask, npiece, clipthresh, clipniter, outlog, blfile)
2365
2366 workscan._add_history("cspline_baseline", varlist)
2367
2368 if insitu:
2369 self._assign(workscan)
2370 else:
2371 return workscan
2372
2373 except RuntimeError, e:
2374 msg = "The fit failed, possibly because it didn't converge."
2375 if rcParams["verbose"]:
2376 asaplog.push(str(e))
2377 asaplog.push(str(msg))
2378 return
2379 else:
2380 raise RuntimeError(str(e)+'\n'+msg)
2381
2382
2383 def auto_cspline_baseline(self, insitu=None, mask=None, npiece=None, clipthresh=None,
2384 clipniter=None, edge=None, threshold=None,
2385 chan_avg_limit=None, plot=None, outlog=None, blfile=None):
2386 """\
2387 Return a scan which has been baselined (all rows) by cubic spline
2388 function (piecewise cubic polynomial).
2389 Spectral lines are detected first using linefinder and masked out
2390 to avoid them affecting the baseline solution.
2391
2392 Parameters:
2393 insitu: if False a new scantable is returned.
2394 Otherwise, the scaling is done in-situ
2395 The default is taken from .asaprc (False)
2396 mask: an optional mask retreived from scantable
2397 npiece: Number of pieces. (default is 2)
2398 clipthresh: Clipping threshold. (default is 3.0, unit: sigma)
2399 clipniter: maximum number of iteration of 'clipthresh'-sigma clipping (default is 1)
2400 edge: an optional number of channel to drop at
2401 the edge of spectrum. If only one value is
2402 specified, the same number will be dropped
2403 from both sides of the spectrum. Default
2404 is to keep all channels. Nested tuples
2405 represent individual edge selection for
2406 different IFs (a number of spectral channels
2407 can be different)
2408 threshold: the threshold used by line finder. It is
2409 better to keep it large as only strong lines
2410 affect the baseline solution.
2411 chan_avg_limit:
2412 a maximum number of consequtive spectral
2413 channels to average during the search of
2414 weak and broad lines. The default is no
2415 averaging (and no search for weak lines).
2416 If such lines can affect the fitted baseline
2417 (e.g. a high order polynomial is fitted),
2418 increase this parameter (usually values up
2419 to 8 are reasonable). Most users of this
2420 method should find the default value sufficient.
2421 plot: *** CURRENTLY UNAVAILABLE, ALWAYS FALSE ***
2422 plot the fit and the residual. In this each
2423 indivual fit has to be approved, by typing 'y'
2424 or 'n'
2425 outlog: Output the coefficients of the best-fit
2426 function to logger (default is False)
2427 blfile: Name of a text file in which the best-fit
2428 parameter values to be written
2429 (default is "": no file/logger output)
2430
2431 Example:
2432 bscan = scan.auto_cspline_baseline(npiece=3, insitu=False)
2433
2434 Note:
2435 The best-fit parameter values output in logger and/or blfile are now
2436 based on specunit of 'channel'.
2437 """
2438
2439 varlist = vars()
2440
2441 if insitu is None: insitu = rcParams['insitu']
2442 if insitu:
2443 workscan = self
2444 else:
2445 workscan = self.copy()
2446
2447 nchan = workscan.nchan()
2448
2449 if mask is None: mask = [True for i in xrange(nchan)]
2450 if npiece is None: npiece = 2
2451 if clipthresh is None: clipthresh = 3.0
2452 if clipniter is None: clipniter = 1
2453 if edge is None: edge = (0, 0)
2454 if threshold is None: threshold = 3
2455 if chan_avg_limit is None: chan_avg_limit = 1
2456 if plot is None: plot = False
2457 if outlog is None: outlog = False
2458 if blfile is None: blfile = ""
2459
2460 outblfile = (blfile != "") and os.path.exists(os.path.expanduser(os.path.expandvars(blfile)))
2461
2462 from asap.asaplinefind import linefinder
2463 from asap import _is_sequence_or_number as _is_valid
2464
2465 if not (isinstance(edge, list) or isinstance(edge, tuple)): edge = [ edge ]
2466 individualedge = False;
2467 if len(edge) > 1: individualedge = isinstance(edge[0], list) or isinstance(edge[0], tuple)
2468
2469 if individualedge:
2470 for edgepar in edge:
2471 if not _is_valid(edgepar, int):
2472 raise ValueError, "Each element of the 'edge' tuple has \
2473 to be a pair of integers or an integer."
2474 else:
2475 if not _is_valid(edge, int):
2476 raise ValueError, "Parameter 'edge' has to be an integer or a \
2477 pair of integers specified as a tuple. \
2478 Nested tuples are allowed \
2479 to make individual selection for different IFs."
2480
2481 if len(edge) > 1:
2482 curedge = edge
2483 else:
2484 curedge = edge + edge
2485
2486 try:
2487 #CURRENTLY, PLOT=true UNAVAILABLE UNTIL cubic spline fitting is implemented as a fitter method.
2488 if individualedge:
2489 curedge = []
2490 for i in xrange(len(edge)):
2491 curedge += edge[i]
2492
2493 workscan._auto_cspline_baseline(mask, npiece, clipthresh, clipniter, curedge, threshold, chan_avg_limit, outlog, blfile)
2494
2495 workscan._add_history("auto_cspline_baseline", varlist)
2496
2497 if insitu:
2498 self._assign(workscan)
2499 else:
2500 return workscan
2501
2502 except RuntimeError, e:
2503 msg = "The fit failed, possibly because it didn't converge."
2504 if rcParams["verbose"]:
2505 asaplog.push(str(e))
2506 asaplog.push(str(msg))
2507 return
2508 else:
2509 raise RuntimeError(str(e)+'\n'+msg)
2510
2511
2512 @asaplog_post_dec
2513 def poly_baseline(self, insitu=None, mask=None, order=None, plot=None, outlog=None, blfile=None):
2514 """\
2515 Return a scan which has been baselined (all rows) by a polynomial.
2516 Parameters:
2517 insitu: if False a new scantable is returned.
2518 Otherwise, the scaling is done in-situ
2519 The default is taken from .asaprc (False)
2520 mask: an optional mask
2521 order: the order of the polynomial (default is 0)
2522 plot: plot the fit and the residual. In this each
2523 indivual fit has to be approved, by typing 'y'
2524 or 'n'
2525 outlog: Output the coefficients of the best-fit
2526 function to logger (default is False)
2527 blfile: Name of a text file in which the best-fit
2528 parameter values to be written
2529 (default is "": no file/logger output)
2530
2531 Example:
2532 # return a scan baselined by a third order polynomial,
2533 # not using a mask
2534 bscan = scan.poly_baseline(order=3)
2535 """
2536
2537 varlist = vars()
2538
2539 if insitu is None: insitu = rcParams["insitu"]
2540 if insitu:
2541 workscan = self
2542 else:
2543 workscan = self.copy()
2544
2545 nchan = workscan.nchan()
2546
2547 if mask is None: mask = [True for i in xrange(nchan)]
2548 if order is None: order = 0
2549 if plot is None: plot = False
2550 if outlog is None: outlog = False
2551 if blfile is None: blfile = ""
2552
2553 outblfile = (blfile != "") and os.path.exists(os.path.expanduser(os.path.expandvars(blfile)))
2554
2555 try:
2556 rows = xrange(workscan.nrow())
2557
2558 #if len(rows) > 0: workscan._init_blinfo()
2559
2560 if plot:
2561 if outblfile: blf = open(blfile, "a")
2562
2563 f = fitter()
2564 f.set_function(lpoly=order)
2565 for r in rows:
2566 f.x = workscan._getabcissa(r)
2567 f.y = workscan._getspectrum(r)
2568 f.mask = mask_and(mask, workscan._getmask(r)) # (CAS-1434)
2569 f.data = None
2570 f.fit()
2571
2572 f.plot(residual=True)
2573 accept_fit = raw_input("Accept fit ( [y]/n ): ")
2574 if accept_fit.upper() == "N":
2575 #workscan._append_blinfo(None, None, None)
2576 continue
2577
2578 blpars = f.get_parameters()
2579 masklist = workscan.get_masklist(f.mask, row=r, silent=True)
2580 #workscan._append_blinfo(blpars, masklist, f.mask)
2581 workscan._setspectrum(f.fitter.getresidual(), r)
2582
2583 if outblfile:
2584 rms = workscan.get_rms(f.mask, r)
2585 dataout = workscan.format_blparams_row(blpars["params"], blpars["fixed"], rms, str(masklist), r, True)
2586 blf.write(dataout)
2587
2588 f._p.unmap()
2589 f._p = None
2590
2591 if outblfile: blf.close()
2592 else:
2593 workscan._poly_baseline(mask, order, outlog, blfile)
2594
2595 workscan._add_history("poly_baseline", varlist)
2596
2597 if insitu:
2598 self._assign(workscan)
2599 else:
2600 return workscan
2601
2602 except RuntimeError, e:
2603 msg = "The fit failed, possibly because it didn't converge."
2604 if rcParams["verbose"]:
2605 asaplog.push(str(e))
2606 asaplog.push(str(msg))
2607 return
2608 else:
2609 raise RuntimeError(str(e)+'\n'+msg)
2610
2611
2612 def auto_poly_baseline(self, insitu=None, mask=None, order=None, edge=None, threshold=None,
2613 chan_avg_limit=None, plot=None, outlog=None, blfile=None):
2614 """\
2615 Return a scan which has been baselined (all rows) by a polynomial.
2616 Spectral lines are detected first using linefinder and masked out
2617 to avoid them affecting the baseline solution.
2618
2619 Parameters:
2620 insitu: if False a new scantable is returned.
2621 Otherwise, the scaling is done in-situ
2622 The default is taken from .asaprc (False)
2623 mask: an optional mask retreived from scantable
2624 order: the order of the polynomial (default is 0)
2625 edge: an optional number of channel to drop at
2626 the edge of spectrum. If only one value is
2627 specified, the same number will be dropped
2628 from both sides of the spectrum. Default
2629 is to keep all channels. Nested tuples
2630 represent individual edge selection for
2631 different IFs (a number of spectral channels
2632 can be different)
2633 threshold: the threshold used by line finder. It is
2634 better to keep it large as only strong lines
2635 affect the baseline solution.
2636 chan_avg_limit:
2637 a maximum number of consequtive spectral
2638 channels to average during the search of
2639 weak and broad lines. The default is no
2640 averaging (and no search for weak lines).
2641 If such lines can affect the fitted baseline
2642 (e.g. a high order polynomial is fitted),
2643 increase this parameter (usually values up
2644 to 8 are reasonable). Most users of this
2645 method should find the default value sufficient.
2646 plot: plot the fit and the residual. In this each
2647 indivual fit has to be approved, by typing 'y'
2648 or 'n'
2649 outlog: Output the coefficients of the best-fit
2650 function to logger (default is False)
2651 blfile: Name of a text file in which the best-fit
2652 parameter values to be written
2653 (default is "": no file/logger output)
2654
2655 Example:
2656 bscan = scan.auto_poly_baseline(order=7, insitu=False)
2657 """
2658
2659 varlist = vars()
2660
2661 if insitu is None: insitu = rcParams['insitu']
2662 if insitu:
2663 workscan = self
2664 else:
2665 workscan = self.copy()
2666
2667 nchan = workscan.nchan()
2668
2669 if mask is None: mask = [True for i in xrange(nchan)]
2670 if order is None: order = 0
2671 if edge is None: edge = (0, 0)
2672 if threshold is None: threshold = 3
2673 if chan_avg_limit is None: chan_avg_limit = 1
2674 if plot is None: plot = False
2675 if outlog is None: outlog = False
2676 if blfile is None: blfile = ""
2677
2678 outblfile = (blfile != "") and os.path.exists(os.path.expanduser(os.path.expandvars(blfile)))
2679
2680 from asap.asaplinefind import linefinder
2681 from asap import _is_sequence_or_number as _is_valid
2682
2683 if not (isinstance(edge, list) or isinstance(edge, tuple)): edge = [ edge ]
2684 individualedge = False;
2685 if len(edge) > 1: individualedge = isinstance(edge[0], list) or isinstance(edge[0], tuple)
2686
2687 if individualedge:
2688 for edgepar in edge:
2689 if not _is_valid(edgepar, int):
2690 raise ValueError, "Each element of the 'edge' tuple has \
2691 to be a pair of integers or an integer."
2692 else:
2693 if not _is_valid(edge, int):
2694 raise ValueError, "Parameter 'edge' has to be an integer or a \
2695 pair of integers specified as a tuple. \
2696 Nested tuples are allowed \
2697 to make individual selection for different IFs."
2698
2699 if len(edge) > 1:
2700 curedge = edge
2701 else:
2702 curedge = edge + edge
2703
2704 try:
2705 rows = xrange(workscan.nrow())
2706
2707 #if len(rows) > 0: workscan._init_blinfo()
2708
2709 if plot:
2710 if outblfile: blf = open(blfile, "a")
2711
2712 fl = linefinder()
2713 fl.set_options(threshold=threshold,avg_limit=chan_avg_limit)
2714 fl.set_scan(workscan)
2715 f = fitter()
2716 f.set_function(lpoly=order)
2717
2718 for r in rows:
2719 if individualedge:
2720 if len(edge) <= workscan.getif(r):
2721 raise RuntimeError, "Number of edge elements appear to " \
2722 "be less than the number of IFs"
2723 else:
2724 curedge = edge[workscan.getif(r)]
2725
2726 fl.find_lines(r, mask_and(mask, workscan._getmask(r)), curedge) # (CAS-1434)
2727
2728 f.x = workscan._getabcissa(r)
2729 f.y = workscan._getspectrum(r)
2730 f.mask = fl.get_mask()
2731 f.data = None
2732 f.fit()
2733
2734 f.plot(residual=True)
2735 accept_fit = raw_input("Accept fit ( [y]/n ): ")
2736 if accept_fit.upper() == "N":
2737 #workscan._append_blinfo(None, None, None)
2738 continue
2739
2740 blpars = f.get_parameters()
2741 masklist = workscan.get_masklist(f.mask, row=r, silent=True)
2742 #workscan._append_blinfo(blpars, masklist, f.mask)
2743 workscan._setspectrum(f.fitter.getresidual(), r)
2744
2745 if outblfile:
2746 rms = workscan.get_rms(f.mask, r)
2747 dataout = workscan.format_blparams_row(blpars["params"], blpars["fixed"], rms, str(masklist), r, True)
2748 blf.write(dataout)
2749
2750 f._p.unmap()
2751 f._p = None
2752
2753 if outblfile: blf.close()
2754
2755 else:
2756 if individualedge:
2757 curedge = []
2758 for i in xrange(len(edge)):
2759 curedge += edge[i]
2760
2761 workscan._auto_poly_baseline(mask, order, curedge, threshold, chan_avg_limit, outlog, blfile)
2762
2763 workscan._add_history("auto_poly_baseline", varlist)
2764
2765 if insitu:
2766 self._assign(workscan)
2767 else:
2768 return workscan
2769
2770 except RuntimeError, e:
2771 msg = "The fit failed, possibly because it didn't converge."
2772 if rcParams["verbose"]:
2773 asaplog.push(str(e))
2774 asaplog.push(str(msg))
2775 return
2776 else:
2777 raise RuntimeError(str(e)+'\n'+msg)
2778
2779
2780 ### OBSOLETE ##################################################################
2781 @asaplog_post_dec
2782 def old_poly_baseline(self, mask=None, order=0, plot=False, uselin=False, insitu=None, rows=None):
2783 """
2784 Return a scan which has been baselined (all rows) by a polynomial.
2785
2786 Parameters:
2787
2788 mask: an optional mask
2789
2790 order: the order of the polynomial (default is 0)
2791
2792 plot: plot the fit and the residual. In this each
2793 indivual fit has to be approved, by typing 'y'
2794 or 'n'
2795
2796 uselin: use linear polynomial fit
2797
2798 insitu: if False a new scantable is returned.
2799 Otherwise, the scaling is done in-situ
2800 The default is taken from .asaprc (False)
2801
2802 rows: row numbers of spectra to be processed.
2803 (default is None: for all rows)
2804
2805 Example:
2806 # return a scan baselined by a third order polynomial,
2807 # not using a mask
2808 bscan = scan.poly_baseline(order=3)
2809
2810 """
2811 if insitu is None: insitu = rcParams['insitu']
2812 if not insitu:
2813 workscan = self.copy()
2814 else:
2815 workscan = self
2816 varlist = vars()
2817 if mask is None:
2818 mask = [True for i in xrange(self.nchan())]
2819
2820 try:
2821 f = fitter()
2822 if uselin:
2823 f.set_function(lpoly=order)
2824 else:
2825 f.set_function(poly=order)
2826
2827 if rows == None:
2828 rows = xrange(workscan.nrow())
2829 elif isinstance(rows, int):
2830 rows = [ rows ]
2831
2832 if len(rows) > 0:
2833 self.blpars = []
2834 self.masklists = []
2835 self.actualmask = []
2836
2837 for r in rows:
2838 f.x = workscan._getabcissa(r)
2839 f.y = workscan._getspectrum(r)
2840 f.mask = mask_and(mask, workscan._getmask(r)) # (CAS-1434)
2841 f.data = None
2842 f.fit()
2843 if plot:
2844 f.plot(residual=True)
2845 x = raw_input("Accept fit ( [y]/n ): ")
2846 if x.upper() == 'N':
2847 self.blpars.append(None)
2848 self.masklists.append(None)
2849 self.actualmask.append(None)
2850 continue
2851 workscan._setspectrum(f.fitter.getresidual(), r)
2852 self.blpars.append(f.get_parameters())
2853 self.masklists.append(workscan.get_masklist(f.mask, row=r, silent=True))
2854 self.actualmask.append(f.mask)
2855
2856 if plot:
2857 f._p.unmap()
2858 f._p = None
2859 workscan._add_history("poly_baseline", varlist)
2860 if insitu:
2861 self._assign(workscan)
2862 else:
2863 return workscan
2864 except RuntimeError:
2865 msg = "The fit failed, possibly because it didn't converge."
2866 raise RuntimeError(msg)
2867
2868 def _init_blinfo(self):
2869 """\
2870 Initialise the following three auxiliary members:
2871 blpars : parameters of the best-fit baseline,
2872 masklists : mask data (edge positions of masked channels) and
2873 actualmask : mask data (in boolean list),
2874 to keep for use later (including output to logger/text files).
2875 Used by poly_baseline() and auto_poly_baseline() in case of
2876 'plot=True'.
2877 """
2878 self.blpars = []
2879 self.masklists = []
2880 self.actualmask = []
2881 return
2882
2883 def _append_blinfo(self, data_blpars, data_masklists, data_actualmask):
2884 """\
2885 Append baseline-fitting related info to blpars, masklist and
2886 actualmask.
2887 """
2888 self.blpars.append(data_blpars)
2889 self.masklists.append(data_masklists)
2890 self.actualmask.append(data_actualmask)
2891 return
2892
2893 @asaplog_post_dec
2894 def rotate_linpolphase(self, angle):
2895 """\
2896 Rotate the phase of the complex polarization O=Q+iU correlation.
2897 This is always done in situ in the raw data. So if you call this
2898 function more than once then each call rotates the phase further.
2899
2900 Parameters:
2901
2902 angle: The angle (degrees) to rotate (add) by.
2903
2904 Example::
2905
2906 scan.rotate_linpolphase(2.3)
2907
2908 """
2909 varlist = vars()
2910 self._math._rotate_linpolphase(self, angle)
2911 self._add_history("rotate_linpolphase", varlist)
2912 return
2913
2914 @asaplog_post_dec
2915 def rotate_xyphase(self, angle):
2916 """\
2917 Rotate the phase of the XY correlation. This is always done in situ
2918 in the data. So if you call this function more than once
2919 then each call rotates the phase further.
2920
2921 Parameters:
2922
2923 angle: The angle (degrees) to rotate (add) by.
2924
2925 Example::
2926
2927 scan.rotate_xyphase(2.3)
2928
2929 """
2930 varlist = vars()
2931 self._math._rotate_xyphase(self, angle)
2932 self._add_history("rotate_xyphase", varlist)
2933 return
2934
2935 @asaplog_post_dec
2936 def swap_linears(self):
2937 """\
2938 Swap the linear polarisations XX and YY, or better the first two
2939 polarisations as this also works for ciculars.
2940 """
2941 varlist = vars()
2942 self._math._swap_linears(self)
2943 self._add_history("swap_linears", varlist)
2944 return
2945
2946 @asaplog_post_dec
2947 def invert_phase(self):
2948 """\
2949 Invert the phase of the complex polarisation
2950 """
2951 varlist = vars()
2952 self._math._invert_phase(self)
2953 self._add_history("invert_phase", varlist)
2954 return
2955
2956 @asaplog_post_dec
2957 def add(self, offset, insitu=None):
2958 """\
2959 Return a scan where all spectra have the offset added
2960
2961 Parameters:
2962
2963 offset: the offset
2964
2965 insitu: if False a new scantable is returned.
2966 Otherwise, the scaling is done in-situ
2967 The default is taken from .asaprc (False)
2968
2969 """
2970 if insitu is None: insitu = rcParams['insitu']
2971 self._math._setinsitu(insitu)
2972 varlist = vars()
2973 s = scantable(self._math._unaryop(self, offset, "ADD", False))
2974 s._add_history("add", varlist)
2975 if insitu:
2976 self._assign(s)
2977 else:
2978 return s
2979
2980 @asaplog_post_dec
2981 def scale(self, factor, tsys=True, insitu=None):
2982 """\
2983
2984 Return a scan where all spectra are scaled by the given 'factor'
2985
2986 Parameters:
2987
2988 factor: the scaling factor (float or 1D float list)
2989
2990 insitu: if False a new scantable is returned.
2991 Otherwise, the scaling is done in-situ
2992 The default is taken from .asaprc (False)
2993
2994 tsys: if True (default) then apply the operation to Tsys
2995 as well as the data
2996
2997 """
2998 if insitu is None: insitu = rcParams['insitu']
2999 self._math._setinsitu(insitu)
3000 varlist = vars()
3001 s = None
3002 import numpy
3003 if isinstance(factor, list) or isinstance(factor, numpy.ndarray):
3004 if isinstance(factor[0], list) or isinstance(factor[0], numpy.ndarray):
3005 from asapmath import _array2dOp
3006 s = _array2dOp( self.copy(), factor, "MUL", tsys )
3007 else:
3008 s = scantable( self._math._arrayop( self.copy(), factor, "MUL", tsys ) )
3009 else:
3010 s = scantable(self._math._unaryop(self.copy(), factor, "MUL", tsys))
3011 s._add_history("scale", varlist)
3012 if insitu:
3013 self._assign(s)
3014 else:
3015 return s
3016
3017 def set_sourcetype(self, match, matchtype="pattern",
3018 sourcetype="reference"):
3019 """\
3020 Set the type of the source to be an source or reference scan
3021 using the provided pattern.
3022
3023 Parameters:
3024
3025 match: a Unix style pattern, regular expression or selector
3026
3027 matchtype: 'pattern' (default) UNIX style pattern or
3028 'regex' regular expression
3029
3030 sourcetype: the type of the source to use (source/reference)
3031
3032 """
3033 varlist = vars()
3034 basesel = self.get_selection()
3035 stype = -1
3036 if sourcetype.lower().startswith("r"):
3037 stype = 1
3038 elif sourcetype.lower().startswith("s"):
3039 stype = 0
3040 else:
3041 raise ValueError("Illegal sourcetype use s(ource) or r(eference)")
3042 if matchtype.lower().startswith("p"):
3043 matchtype = "pattern"
3044 elif matchtype.lower().startswith("r"):
3045 matchtype = "regex"
3046 else:
3047 raise ValueError("Illegal matchtype, use p(attern) or r(egex)")
3048 sel = selector()
3049 if isinstance(match, selector):
3050 sel = match
3051 else:
3052 sel.set_query("SRCNAME == %s('%s')" % (matchtype, match))
3053 self.set_selection(basesel+sel)
3054 self._setsourcetype(stype)
3055 self.set_selection(basesel)
3056 self._add_history("set_sourcetype", varlist)
3057
3058 @asaplog_post_dec
3059 @preserve_selection
3060 def auto_quotient(self, preserve=True, mode='paired', verify=False):
3061 """\
3062 This function allows to build quotients automatically.
3063 It assumes the observation to have the same number of
3064 "ons" and "offs"
3065
3066 Parameters:
3067
3068 preserve: you can preserve (default) the continuum or
3069 remove it. The equations used are
3070
3071 preserve: Output = Toff * (on/off) - Toff
3072
3073 remove: Output = Toff * (on/off) - Ton
3074
3075 mode: the on/off detection mode
3076 'paired' (default)
3077 identifies 'off' scans by the
3078 trailing '_R' (Mopra/Parkes) or
3079 '_e'/'_w' (Tid) and matches
3080 on/off pairs from the observing pattern
3081 'time'
3082 finds the closest off in time
3083
3084 .. todo:: verify argument is not implemented
3085
3086 """
3087 varlist = vars()
3088 modes = ["time", "paired"]
3089 if not mode in modes:
3090 msg = "please provide valid mode. Valid modes are %s" % (modes)
3091 raise ValueError(msg)
3092 s = None
3093 if mode.lower() == "paired":
3094 sel = self.get_selection()
3095 sel.set_query("SRCTYPE==psoff")
3096 self.set_selection(sel)
3097 offs = self.copy()
3098 sel.set_query("SRCTYPE==pson")
3099 self.set_selection(sel)
3100 ons = self.copy()
3101 s = scantable(self._math._quotient(ons, offs, preserve))
3102 elif mode.lower() == "time":
3103 s = scantable(self._math._auto_quotient(self, mode, preserve))
3104 s._add_history("auto_quotient", varlist)
3105 return s
3106
3107 @asaplog_post_dec
3108 def mx_quotient(self, mask = None, weight='median', preserve=True):
3109 """\
3110 Form a quotient using "off" beams when observing in "MX" mode.
3111
3112 Parameters:
3113
3114 mask: an optional mask to be used when weight == 'stddev'
3115
3116 weight: How to average the off beams. Default is 'median'.
3117
3118 preserve: you can preserve (default) the continuum or
3119 remove it. The equations used are:
3120
3121 preserve: Output = Toff * (on/off) - Toff
3122
3123 remove: Output = Toff * (on/off) - Ton
3124
3125 """
3126 mask = mask or ()
3127 varlist = vars()
3128 on = scantable(self._math._mx_extract(self, 'on'))
3129 preoff = scantable(self._math._mx_extract(self, 'off'))
3130 off = preoff.average_time(mask=mask, weight=weight, scanav=False)
3131 from asapmath import quotient
3132 q = quotient(on, off, preserve)
3133 q._add_history("mx_quotient", varlist)
3134 return q
3135
3136 @asaplog_post_dec
3137 def freq_switch(self, insitu=None):
3138 """\
3139 Apply frequency switching to the data.
3140
3141 Parameters:
3142
3143 insitu: if False a new scantable is returned.
3144 Otherwise, the swictching is done in-situ
3145 The default is taken from .asaprc (False)
3146
3147 """
3148 if insitu is None: insitu = rcParams['insitu']
3149 self._math._setinsitu(insitu)
3150 varlist = vars()
3151 s = scantable(self._math._freqswitch(self))
3152 s._add_history("freq_switch", varlist)
3153 if insitu:
3154 self._assign(s)
3155 else:
3156 return s
3157
3158 @asaplog_post_dec
3159 def recalc_azel(self):
3160 """Recalculate the azimuth and elevation for each position."""
3161 varlist = vars()
3162 self._recalcazel()
3163 self._add_history("recalc_azel", varlist)
3164 return
3165
3166 @asaplog_post_dec
3167 def __add__(self, other):
3168 varlist = vars()
3169 s = None
3170 if isinstance(other, scantable):
3171 s = scantable(self._math._binaryop(self, other, "ADD"))
3172 elif isinstance(other, float):
3173 s = scantable(self._math._unaryop(self, other, "ADD", False))
3174 else:
3175 raise TypeError("Other input is not a scantable or float value")
3176 s._add_history("operator +", varlist)
3177 return s
3178
3179 @asaplog_post_dec
3180 def __sub__(self, other):
3181 """
3182 implicit on all axes and on Tsys
3183 """
3184 varlist = vars()
3185 s = None
3186 if isinstance(other, scantable):
3187 s = scantable(self._math._binaryop(self, other, "SUB"))
3188 elif isinstance(other, float):
3189 s = scantable(self._math._unaryop(self, other, "SUB", False))
3190 else:
3191 raise TypeError("Other input is not a scantable or float value")
3192 s._add_history("operator -", varlist)
3193 return s
3194
3195 @asaplog_post_dec
3196 def __mul__(self, other):
3197 """
3198 implicit on all axes and on Tsys
3199 """
3200 varlist = vars()
3201 s = None
3202 if isinstance(other, scantable):
3203 s = scantable(self._math._binaryop(self, other, "MUL"))
3204 elif isinstance(other, float):
3205 s = scantable(self._math._unaryop(self, other, "MUL", False))
3206 else:
3207 raise TypeError("Other input is not a scantable or float value")
3208 s._add_history("operator *", varlist)
3209 return s
3210
3211
3212 @asaplog_post_dec
3213 def __div__(self, other):
3214 """
3215 implicit on all axes and on Tsys
3216 """
3217 varlist = vars()
3218 s = None
3219 if isinstance(other, scantable):
3220 s = scantable(self._math._binaryop(self, other, "DIV"))
3221 elif isinstance(other, float):
3222 if other == 0.0:
3223 raise ZeroDivisionError("Dividing by zero is not recommended")
3224 s = scantable(self._math._unaryop(self, other, "DIV", False))
3225 else:
3226 raise TypeError("Other input is not a scantable or float value")
3227 s._add_history("operator /", varlist)
3228 return s
3229
3230 @asaplog_post_dec
3231 def get_fit(self, row=0):
3232 """\
3233 Print or return the stored fits for a row in the scantable
3234
3235 Parameters:
3236
3237 row: the row which the fit has been applied to.
3238
3239 """
3240 if row > self.nrow():
3241 return
3242 from asap.asapfit import asapfit
3243 fit = asapfit(self._getfit(row))
3244 asaplog.push( '%s' %(fit) )
3245 return fit.as_dict()
3246
3247 def flag_nans(self):
3248 """\
3249 Utility function to flag NaN values in the scantable.
3250 """
3251 import numpy
3252 basesel = self.get_selection()
3253 for i in range(self.nrow()):
3254 sel = self.get_row_selector(i)
3255 self.set_selection(basesel+sel)
3256 nans = numpy.isnan(self._getspectrum(0))
3257 if numpy.any(nans):
3258 bnans = [ bool(v) for v in nans]
3259 self.flag(bnans)
3260 self.set_selection(basesel)
3261
3262 def get_row_selector(self, rowno):
3263 #return selector(beams=self.getbeam(rowno),
3264 # ifs=self.getif(rowno),
3265 # pols=self.getpol(rowno),
3266 # scans=self.getscan(rowno),
3267 # cycles=self.getcycle(rowno))
3268 return selector(rows=[rowno])
3269
3270 def _add_history(self, funcname, parameters):
3271 if not rcParams['scantable.history']:
3272 return
3273 # create date
3274 sep = "##"
3275 from datetime import datetime
3276 dstr = datetime.now().strftime('%Y/%m/%d %H:%M:%S')
3277 hist = dstr+sep
3278 hist += funcname+sep#cdate+sep
3279 if parameters.has_key('self'): del parameters['self']
3280 for k, v in parameters.iteritems():
3281 if type(v) is dict:
3282 for k2, v2 in v.iteritems():
3283 hist += k2
3284 hist += "="
3285 if isinstance(v2, scantable):
3286 hist += 'scantable'
3287 elif k2 == 'mask':
3288 if isinstance(v2, list) or isinstance(v2, tuple):
3289 hist += str(self._zip_mask(v2))
3290 else:
3291 hist += str(v2)
3292 else:
3293 hist += str(v2)
3294 else:
3295 hist += k
3296 hist += "="
3297 if isinstance(v, scantable):
3298 hist += 'scantable'
3299 elif k == 'mask':
3300 if isinstance(v, list) or isinstance(v, tuple):
3301 hist += str(self._zip_mask(v))
3302 else:
3303 hist += str(v)
3304 else:
3305 hist += str(v)
3306 hist += sep
3307 hist = hist[:-2] # remove trailing '##'
3308 self._addhistory(hist)
3309
3310
3311 def _zip_mask(self, mask):
3312 mask = list(mask)
3313 i = 0
3314 segments = []
3315 while mask[i:].count(1):
3316 i += mask[i:].index(1)
3317 if mask[i:].count(0):
3318 j = i + mask[i:].index(0)
3319 else:
3320 j = len(mask)
3321 segments.append([i, j])
3322 i = j
3323 return segments
3324
3325 def _get_ordinate_label(self):
3326 fu = "("+self.get_fluxunit()+")"
3327 import re
3328 lbl = "Intensity"
3329 if re.match(".K.", fu):
3330 lbl = "Brightness Temperature "+ fu
3331 elif re.match(".Jy.", fu):
3332 lbl = "Flux density "+ fu
3333 return lbl
3334
3335 def _check_ifs(self):
3336 #nchans = [self.nchan(i) for i in range(self.nif(-1))]
3337 nchans = [self.nchan(i) for i in self.getifnos()]
3338 nchans = filter(lambda t: t > 0, nchans)
3339 return (sum(nchans)/len(nchans) == nchans[0])
3340
3341 @asaplog_post_dec
3342 #def _fill(self, names, unit, average, getpt, antenna):
3343 def _fill(self, names, unit, average, opts={}):
3344 first = True
3345 fullnames = []
3346 for name in names:
3347 name = os.path.expandvars(name)
3348 name = os.path.expanduser(name)
3349 if not os.path.exists(name):
3350 msg = "File '%s' does not exists" % (name)
3351 raise IOError(msg)
3352 fullnames.append(name)
3353 if average:
3354 asaplog.push('Auto averaging integrations')
3355 stype = int(rcParams['scantable.storage'].lower() == 'disk')
3356 for name in fullnames:
3357 tbl = Scantable(stype)
3358 if is_ms( name ):
3359 r = msfiller( tbl )
3360 else:
3361 r = filler( tbl )
3362 rx = rcParams['scantable.reference']
3363 r.setreferenceexpr(rx)
3364 #r = filler(tbl)
3365 #rx = rcParams['scantable.reference']
3366 #r.setreferenceexpr(rx)
3367 msg = "Importing %s..." % (name)
3368 asaplog.push(msg, False)
3369 #opts = {'ms': {'antenna' : antenna, 'getpt': getpt} }
3370 r.open(name, opts)# antenna, -1, -1, getpt)
3371 r.fill()
3372 if average:
3373 tbl = self._math._average((tbl, ), (), 'NONE', 'SCAN')
3374 if not first:
3375 tbl = self._math._merge([self, tbl])
3376 Scantable.__init__(self, tbl)
3377 r.close()
3378 del r, tbl
3379 first = False
3380 #flush log
3381 asaplog.post()
3382 if unit is not None:
3383 self.set_fluxunit(unit)
3384 if not is_casapy():
3385 self.set_freqframe(rcParams['scantable.freqframe'])
3386
3387
3388 def __getitem__(self, key):
3389 if key < 0:
3390 key += self.nrow()
3391 if key >= self.nrow():
3392 raise IndexError("Row index out of range.")
3393 return self._getspectrum(key)
3394
3395 def __setitem__(self, key, value):
3396 if key < 0:
3397 key += self.nrow()
3398 if key >= self.nrow():
3399 raise IndexError("Row index out of range.")
3400 if not hasattr(value, "__len__") or \
3401 len(value) > self.nchan(self.getif(key)):
3402 raise ValueError("Spectrum length doesn't match.")
3403 return self._setspectrum(value, key)
3404
3405 def __len__(self):
3406 return self.nrow()
3407
3408 def __iter__(self):
3409 for i in range(len(self)):
3410 yield self[i]
Note: See TracBrowser for help on using the repository browser.