source: trunk/python/scantable.py@ 2599

Last change on this file since 2599 was 2599, checked in by Kana Sugimoto, 12 years ago

New Development: No

JIRA Issue: No (a bug fix)

Ready for Test: Yes

Interface Changes: No

What Interface Changed:

Test Programs:

Put in Release Notes: No

Module(s): scantable, and sdtasks with restfreq parameter

Description:

Fixed a bug in scantable.set_restfreq(), i.e., the method ignored user selection
of scantable and caused unexpected modification of restvfrequencies.
It now modifies rest frequencies of selected data only.


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