source: branches/casa-prerelease/pre-asap/python/scantable.py@ 2066

Last change on this file since 2066 was 2047, checked in by WataruKawasaki, 14 years ago

New Development: Yes

JIRA Issue: Yes CAS-2847

Ready for Test: Yes

Interface Changes: No

What Interface Changed:

Test Programs:

Put in Release Notes: Yes

Module(s): scantable

Description: added {auto_}sinusoid_baseline() for sinusoidal baseline fitting. also minor bug fixes for asapfitter.


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