source: trunk/python/scantable.py@ 1672

Last change on this file since 1672 was 1617, checked in by Malte Marquarding, 15 years ago

added help for parallactify

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 70.4 KB
RevLine 
[876]1from asap._asap import Scantable
[226]2from asap import rcParams
[1589]3from asap import print_log, print_log_dec
[1118]4from asap import asaplog
[946]5from asap import selector
[1153]6from asap import linecatalog
[1600]7from asap.coordinate import coordinate
[1295]8from asap import _n_bools, mask_not, mask_and, mask_or
[102]9
[876]10class scantable(Scantable):
[102]11 """
12 The ASAP container for scans
13 """
[1589]14 @print_log_dec
[1586]15 def __init__(self, filename, average=None, unit=None, parallactify=None):
[102]16 """
17 Create a scantable from a saved one or make a reference
18 Parameters:
[181]19 filename: the name of an asap table on disk
20 or
21 the name of a rpfits/sdfits/ms file
22 (integrations within scans are auto averaged
23 and the whole file is read)
24 or
25 [advanced] a reference to an existing
[102]26 scantable
[484]27 average: average all integrations withinb a scan on read.
28 The default (True) is taken from .asaprc.
29 unit: brightness unit; must be consistent with K or Jy.
[340]30 Over-rides the default selected by the reader
31 (input rpfits/sdfits/ms) or replaces the value
32 in existing scantables
[1586]33 parallactify: Indcicate that the data had been parallatified.
34 Default is taken form rc file.
[710]35 """
[976]36 if average is None:
[710]37 average = rcParams['scantable.autoaverage']
[1593]38 parallactify = parallactify or rcParams['scantable.parallactify']
[1259]39 varlist = vars()
[876]40 from asap._asap import stmath
41 self._math = stmath()
42 if isinstance(filename, Scantable):
43 Scantable.__init__(self, filename)
[181]44 else:
[1435]45 if isinstance(filename, str):# or \
46# (isinstance(filename, list) or isinstance(filename, tuple)) \
47# and isinstance(filename[-1], str):
[976]48 import os.path
49 filename = os.path.expandvars(filename)
50 filename = os.path.expanduser(filename)
51 if not os.path.exists(filename):
52 s = "File '%s' not found." % (filename)
[718]53 if rcParams['verbose']:
[976]54 asaplog.push(s)
[1589]55 #print asaplog.pop().strip()
[718]56 return
[976]57 raise IOError(s)
[1115]58 if os.path.isdir(filename) \
[1118]59 and not os.path.exists(filename+'/table.f1'):
[976]60 # crude check if asap table
61 if os.path.exists(filename+'/table.info'):
[1118]62 ondisk = rcParams['scantable.storage'] == 'disk'
63 Scantable.__init__(self, filename, ondisk)
[976]64 if unit is not None:
65 self.set_fluxunit(unit)
66 self.set_freqframe(rcParams['scantable.freqframe'])
[718]67 else:
[1118]68 msg = "The given file '%s'is not a valid " \
69 "asap table." % (filename)
[976]70 if rcParams['verbose']:
71 print msg
72 return
73 else:
74 raise IOError(msg)
[226]75 else:
[1118]76 self._fill([filename], unit, average)
77 elif (isinstance(filename, list) or isinstance(filename, tuple)) \
[976]78 and isinstance(filename[-1], str):
79 self._fill(filename, unit, average)
[1586]80 self.parallactify(parallactify)
[1259]81 self._add_history("scantable", varlist)
[102]82
[1589]83 @print_log_dec
[876]84 def save(self, name=None, format=None, overwrite=False):
[116]85 """
[1280]86 Store the scantable on disk. This can be an asap (aips++) Table,
87 SDFITS or MS2 format.
[116]88 Parameters:
[1093]89 name: the name of the outputfile. For format "ASCII"
90 this is the root file name (data in 'name'.txt
[497]91 and header in 'name'_header.txt)
[116]92 format: an optional file format. Default is ASAP.
[280]93 Allowed are - 'ASAP' (save as ASAP [aips++] Table),
[194]94 'SDFITS' (save as SDFITS file)
[200]95 'ASCII' (saves as ascii text file)
[226]96 'MS2' (saves as an aips++
97 MeasurementSet V2)
[1573]98 'FITS' (save as image FITS - not
[1443]99 readable by class)
100 'CLASS' (save as FITS readable by CLASS)
[411]101 overwrite: If the file should be overwritten if it exists.
[256]102 The default False is to return with warning
[411]103 without writing the output. USE WITH CARE.
[116]104 Example:
105 scan.save('myscan.asap')
[1118]106 scan.save('myscan.sdfits', 'SDFITS')
[116]107 """
[411]108 from os import path
[1593]109 format = format or rcParams['scantable.save']
[256]110 suffix = '.'+format.lower()
[1118]111 if name is None or name == "":
[256]112 name = 'scantable'+suffix
[718]113 msg = "No filename given. Using default name %s..." % name
114 asaplog.push(msg)
[411]115 name = path.expandvars(name)
[256]116 if path.isfile(name) or path.isdir(name):
117 if not overwrite:
[718]118 msg = "File %s exists." % name
119 if rcParams['verbose']:
120 print msg
121 return
122 else:
123 raise IOError(msg)
[451]124 format2 = format.upper()
125 if format2 == 'ASAP':
[116]126 self._save(name)
127 else:
[989]128 from asap._asap import stwriter as stw
[1118]129 writer = stw(format2)
130 writer.write(self, name)
[116]131 return
132
[102]133 def copy(self):
134 """
135 Return a copy of this scantable.
[1348]136 Note:
137 This makes a full (deep) copy. scan2 = scan1 makes a reference.
[102]138 Parameters:
[113]139 none
[102]140 Example:
141 copiedscan = scan.copy()
142 """
[876]143 sd = scantable(Scantable._copy(self))
[113]144 return sd
145
[1093]146 def drop_scan(self, scanid=None):
147 """
148 Return a new scantable where the specified scan number(s) has(have)
149 been dropped.
150 Parameters:
151 scanid: a (list of) scan number(s)
152 """
153 from asap import _is_sequence_or_number as _is_valid
154 from asap import _to_list
155 from asap import unique
156 if not _is_valid(scanid):
157 if rcParams['verbose']:
158 print "Please specify a scanno to drop from the scantable"
159 return
160 else:
161 raise RuntimeError("No scan given")
162 try:
163 scanid = _to_list(scanid)
164 allscans = unique([ self.getscan(i) for i in range(self.nrow())])
165 for sid in scanid: allscans.remove(sid)
[1118]166 if len(allscans) == 0:
167 raise ValueError("Can't remove all scans")
[1093]168 except ValueError:
169 if rcParams['verbose']:
170 print "Couldn't find any match."
171 return
172 else: raise
173 try:
[1593]174 sel = selector(scans=allscans)
[1594]175 return self._select_copy(sel)
[1093]176 except RuntimeError:
[1118]177 if rcParams['verbose']:
178 print "Couldn't find any match."
179 else:
180 raise
[1093]181
[1594]182 def _select_copy(self, selection):
183 orig = self.get_selection()
184 self.set_selection(orig+selection)
185 cp = self.copy()
186 self.set_selection(orig)
187 return cp
188
[102]189 def get_scan(self, scanid=None):
190 """
191 Return a specific scan (by scanno) or collection of scans (by
192 source name) in a new scantable.
[1348]193 Note:
194 See scantable.drop_scan() for the inverse operation.
[102]195 Parameters:
[513]196 scanid: a (list of) scanno or a source name, unix-style
197 patterns are accepted for source name matching, e.g.
198 '*_R' gets all 'ref scans
[102]199 Example:
[513]200 # get all scans containing the source '323p459'
201 newscan = scan.get_scan('323p459')
202 # get all 'off' scans
203 refscans = scan.get_scan('*_R')
204 # get a susbset of scans by scanno (as listed in scan.summary())
[1118]205 newscan = scan.get_scan([0, 2, 7, 10])
[102]206 """
207 if scanid is None:
[718]208 if rcParams['verbose']:
[1118]209 print "Please specify a scan no or name to " \
210 "retrieve from the scantable"
[718]211 return
212 else:
213 raise RuntimeError("No scan given")
214
[102]215 try:
[946]216 bsel = self.get_selection()
217 sel = selector()
[102]218 if type(scanid) is str:
[946]219 sel.set_name(scanid)
[1594]220 return self._select_copy(sel)
[102]221 elif type(scanid) is int:
[946]222 sel.set_scans([scanid])
[1594]223 return self._select_copy(sel)
[381]224 elif type(scanid) is list:
[946]225 sel.set_scans(scanid)
[1594]226 return self._select_copy(sel)
[381]227 else:
[718]228 msg = "Illegal scanid type, use 'int' or 'list' if ints."
229 if rcParams['verbose']:
230 print msg
231 else:
232 raise TypeError(msg)
[102]233 except RuntimeError:
[718]234 if rcParams['verbose']: print "Couldn't find any match."
235 else: raise
[102]236
237 def __str__(self):
[1118]238 return Scantable._summary(self, True)
[102]239
[976]240 def summary(self, filename=None):
[102]241 """
242 Print a summary of the contents of this scantable.
243 Parameters:
244 filename: the name of a file to write the putput to
245 Default - no file output
[381]246 verbose: print extra info such as the frequency table
247 The default (False) is taken from .asaprc
[102]248 """
[976]249 info = Scantable._summary(self, True)
[102]250 if filename is not None:
[256]251 if filename is "":
252 filename = 'scantable_summary.txt'
[415]253 from os.path import expandvars, isdir
[411]254 filename = expandvars(filename)
[415]255 if not isdir(filename):
[413]256 data = open(filename, 'w')
257 data.write(info)
258 data.close()
259 else:
[718]260 msg = "Illegal file name '%s'." % (filename)
261 if rcParams['verbose']:
262 print msg
263 else:
264 raise IOError(msg)
265 if rcParams['verbose']:
[794]266 try:
267 from IPython.genutils import page as pager
268 except ImportError:
269 from pydoc import pager
270 pager(info)
[718]271 else:
272 return info
[710]273
[1512]274 def get_spectrum(self, rowno):
[1471]275 """Return the spectrum for the current row in the scantable as a list.
276 Parameters:
[1573]277 rowno: the row number to retrieve the spectrum from
[1471]278 """
279 return self._getspectrum(rowno)
[946]280
[1471]281 def get_mask(self, rowno):
282 """Return the mask for the current row in the scantable as a list.
283 Parameters:
[1573]284 rowno: the row number to retrieve the mask from
[1471]285 """
286 return self._getmask(rowno)
287
288 def set_spectrum(self, spec, rowno):
289 """Return the spectrum for the current row in the scantable as a list.
290 Parameters:
291 spec: the spectrum
[1573]292 rowno: the row number to set the spectrum for
[1471]293 """
294 assert(len(spec) == self.nchan())
295 return self._setspectrum(spec, rowno)
296
[1600]297 def get_coordinate(self, rowno):
298 """Return the (spectral) coordinate for a a given 'rowno'.
299 NOTE:
300 * This coordinate is only valid until a scantable method modifies
301 the frequency axis.
302 * This coordinate does contain the original frequency set-up
303 NOT the new frame. The conversions however are done using the user
304 specified frame (e.g. LSRK/TOPO). To get the 'real' coordinate,
305 use scantable.freq_align first. Without it there is no closure,
306 i.e.
307 c = myscan.get_coordinate(0)
308 c.to_frequency(c.get_reference_pixel()) != c.get_reference_value()
309
310 Parameters:
311 rowno: the row number for the spectral coordinate
312
313 """
314 return coordinate(Scantable.get_coordinate(self, rowno))
315
[946]316 def get_selection(self):
317 """
[1005]318 Get the selection object currently set on this scantable.
319 Parameters:
320 none
321 Example:
322 sel = scan.get_selection()
323 sel.set_ifs(0) # select IF 0
324 scan.set_selection(sel) # apply modified selection
[946]325 """
326 return selector(self._getselection())
327
[1576]328 def set_selection(self, selection=None, **kw):
[946]329 """
[1005]330 Select a subset of the data. All following operations on this scantable
331 are only applied to thi selection.
332 Parameters:
333 selection: a selector object (default unset the selection)
334 Examples:
335 sel = selector() # create a selection object
[1118]336 self.set_scans([0, 3]) # select SCANNO 0 and 3
[1005]337 scan.set_selection(sel) # set the selection
338 scan.summary() # will only print summary of scanno 0 an 3
339 scan.set_selection() # unset the selection
[946]340 """
[1576]341 if selection is None:
342 # reset
343 if len(kw) == 0:
344 selection = selector()
345 else:
346 # try keywords
347 for k in kw:
348 if k not in selector.fields:
349 raise KeyError("Invalid selection key '%s', valid keys are %s" % (k, selector.fields))
350 selection = selector(**kw)
[946]351 self._setselection(selection)
352
[876]353 def stats(self, stat='stddev', mask=None):
[102]354 """
[135]355 Determine the specified statistic of the current beam/if/pol
[102]356 Takes a 'mask' as an optional parameter to specify which
357 channels should be excluded.
358 Parameters:
[135]359 stat: 'min', 'max', 'sumsq', 'sum', 'mean'
360 'var', 'stddev', 'avdev', 'rms', 'median'
361 mask: an optional mask specifying where the statistic
[102]362 should be determined.
363 Example:
[113]364 scan.set_unit('channel')
[1118]365 msk = scan.create_mask([100, 200], [500, 600])
[135]366 scan.stats(stat='mean', mask=m)
[102]367 """
[1593]368 mask = mask or []
[876]369 if not self._check_ifs():
[1118]370 raise ValueError("Cannot apply mask as the IFs have different "
371 "number of channels. Please use setselection() "
372 "to select individual IFs")
[256]373
[876]374 statvals = self._math._stats(self, mask, stat)
[1590]375 def cb(i):
376 return statvals[i]
[256]377
[1590]378 return self._row_callback(cb, stat)
[102]379
[1118]380 def stddev(self, mask=None):
[135]381 """
382 Determine the standard deviation of the current beam/if/pol
383 Takes a 'mask' as an optional parameter to specify which
384 channels should be excluded.
385 Parameters:
386 mask: an optional mask specifying where the standard
387 deviation should be determined.
388
389 Example:
390 scan.set_unit('channel')
[1118]391 msk = scan.create_mask([100, 200], [500, 600])
[135]392 scan.stddev(mask=m)
393 """
[1118]394 return self.stats(stat='stddev', mask=mask);
[135]395
[1003]396
[1259]397 def get_column_names(self):
[1003]398 """
399 Return a list of column names, which can be used for selection.
400 """
[1259]401 return list(Scantable.get_column_names(self))
[1003]402
[876]403 def get_tsys(self):
[113]404 """
405 Return the System temperatures.
406 Returns:
[876]407 a list of Tsys values for the current selection
[113]408 """
[256]409
[876]410 return self._row_callback(self._gettsys, "Tsys")
[256]411
[876]412 def _row_callback(self, callback, label):
413 out = ""
[1118]414 outvec = []
[1590]415 sep = '-'*50
[876]416 for i in range(self.nrow()):
417 tm = self._gettime(i)
418 src = self._getsourcename(i)
[1590]419 out += 'Scan[%d] (%s) ' % (self.getscan(i), src)
[876]420 out += 'Time[%s]:\n' % (tm)
[1590]421 if self.nbeam(-1) > 1:
422 out += ' Beam[%d] ' % (self.getbeam(i))
423 if self.nif(-1) > 1: out += ' IF[%d] ' % (self.getif(i))
424 if self.npol(-1) > 1: out += ' Pol[%d] ' % (self.getpol(i))
[876]425 outvec.append(callback(i))
426 out += '= %3.3f\n' % (outvec[i])
[1590]427 out += sep+'\n'
[876]428 if rcParams['verbose']:
[1590]429 print sep
[876]430 print " %s" % (label)
[1590]431 print sep
[876]432 print out
[1175]433 return outvec
[256]434
[1070]435 def _get_column(self, callback, row=-1):
436 """
437 """
438 if row == -1:
439 return [callback(i) for i in range(self.nrow())]
440 else:
441 if 0 <= row < self.nrow():
442 return callback(row)
[256]443
[1070]444
[1348]445 def get_time(self, row=-1, asdatetime=False):
[113]446 """
447 Get a list of time stamps for the observations.
[1348]448 Return a datetime object for each integration time stamp in the scantable.
[113]449 Parameters:
[1348]450 row: row no of integration. Default -1 return all rows
451 asdatetime: return values as datetime objects rather than strings
[113]452 Example:
453 none
454 """
[1175]455 from time import strptime
456 from datetime import datetime
[1392]457 times = self._get_column(self._gettime, row)
[1348]458 if not asdatetime:
[1392]459 return times
[1175]460 format = "%Y/%m/%d/%H:%M:%S"
461 if isinstance(times, list):
462 return [datetime(*strptime(i, format)[:6]) for i in times]
463 else:
464 return datetime(*strptime(times, format)[:6])
[102]465
[1348]466
467 def get_inttime(self, row=-1):
468 """
469 Get a list of integration times for the observations.
470 Return a time in seconds for each integration in the scantable.
471 Parameters:
472 row: row no of integration. Default -1 return all rows.
473 Example:
474 none
475 """
[1573]476 return self._get_column(self._getinttime, row)
[1348]477
[1573]478
[714]479 def get_sourcename(self, row=-1):
480 """
[794]481 Get a list source names for the observations.
[714]482 Return a string for each integration in the scantable.
483 Parameters:
[1348]484 row: row no of integration. Default -1 return all rows.
[714]485 Example:
486 none
487 """
[1070]488 return self._get_column(self._getsourcename, row)
[714]489
[794]490 def get_elevation(self, row=-1):
491 """
492 Get a list of elevations for the observations.
493 Return a float for each integration in the scantable.
494 Parameters:
[1348]495 row: row no of integration. Default -1 return all rows.
[794]496 Example:
497 none
498 """
[1070]499 return self._get_column(self._getelevation, row)
[794]500
501 def get_azimuth(self, row=-1):
502 """
503 Get a list of azimuths for the observations.
504 Return a float for each integration in the scantable.
505 Parameters:
[1348]506 row: row no of integration. Default -1 return all rows.
[794]507 Example:
508 none
509 """
[1070]510 return self._get_column(self._getazimuth, row)
[794]511
512 def get_parangle(self, row=-1):
513 """
514 Get a list of parallactic angles for the observations.
515 Return a float for each integration in the scantable.
516 Parameters:
[1348]517 row: row no of integration. Default -1 return all rows.
[794]518 Example:
519 none
520 """
[1070]521 return self._get_column(self._getparangle, row)
[794]522
[1070]523 def get_direction(self, row=-1):
524 """
525 Get a list of Positions on the sky (direction) for the observations.
[1594]526 Return a string for each integration in the scantable.
[1070]527 Parameters:
528 row: row no of integration. Default -1 return all rows
529 Example:
530 none
531 """
532 return self._get_column(self._getdirection, row)
533
[1391]534 def get_directionval(self, row=-1):
535 """
536 Get a list of Positions on the sky (direction) for the observations.
537 Return a float for each integration in the scantable.
538 Parameters:
539 row: row no of integration. Default -1 return all rows
540 Example:
541 none
542 """
543 return self._get_column(self._getdirectionvec, row)
544
545
[102]546 def set_unit(self, unit='channel'):
547 """
548 Set the unit for all following operations on this scantable
549 Parameters:
550 unit: optional unit, default is 'channel'
[1118]551 one of '*Hz', 'km/s', 'channel', ''
[102]552 """
[484]553 varlist = vars()
[1118]554 if unit in ['', 'pixel', 'channel']:
[113]555 unit = ''
556 inf = list(self._getcoordinfo())
557 inf[0] = unit
558 self._setcoordinfo(inf)
[1118]559 self._add_history("set_unit", varlist)
[113]560
[1589]561 @print_log_dec
[484]562 def set_instrument(self, instr):
[358]563 """
[1348]564 Set the instrument for subsequent processing.
[358]565 Parameters:
[710]566 instr: Select from 'ATPKSMB', 'ATPKSHOH', 'ATMOPRA',
[407]567 'DSS-43' (Tid), 'CEDUNA', and 'HOBART'
[358]568 """
569 self._setInstrument(instr)
[1118]570 self._add_history("set_instument", vars())
[358]571
[1589]572 @print_log_dec
[1190]573 def set_feedtype(self, feedtype):
574 """
575 Overwrite the feed type, which might not be set correctly.
576 Parameters:
577 feedtype: 'linear' or 'circular'
578 """
579 self._setfeedtype(feedtype)
580 self._add_history("set_feedtype", vars())
581
[1589]582 @print_log_dec
[276]583 def set_doppler(self, doppler='RADIO'):
584 """
585 Set the doppler for all following operations on this scantable.
586 Parameters:
587 doppler: One of 'RADIO', 'OPTICAL', 'Z', 'BETA', 'GAMMA'
588 """
[484]589 varlist = vars()
[276]590 inf = list(self._getcoordinfo())
591 inf[2] = doppler
592 self._setcoordinfo(inf)
[1118]593 self._add_history("set_doppler", vars())
[710]594
[1589]595 @print_log_dec
[226]596 def set_freqframe(self, frame=None):
[113]597 """
598 Set the frame type of the Spectral Axis.
599 Parameters:
[591]600 frame: an optional frame type, default 'LSRK'. Valid frames are:
[1118]601 'REST', 'TOPO', 'LSRD', 'LSRK', 'BARY',
602 'GEO', 'GALACTO', 'LGROUP', 'CMB'
[113]603 Examples:
604 scan.set_freqframe('BARY')
605 """
[1593]606 frame = frame or rcParams['scantable.freqframe']
[484]607 varlist = vars()
[1118]608 valid = ['REST', 'TOPO', 'LSRD', 'LSRK', 'BARY', \
609 'GEO', 'GALACTO', 'LGROUP', 'CMB']
[591]610
[989]611 if frame in valid:
[113]612 inf = list(self._getcoordinfo())
613 inf[1] = frame
614 self._setcoordinfo(inf)
[1118]615 self._add_history("set_freqframe", varlist)
[102]616 else:
[1118]617 msg = "Please specify a valid freq type. Valid types are:\n", valid
[718]618 if rcParams['verbose']:
619 print msg
620 else:
621 raise TypeError(msg)
[710]622
[989]623 def set_dirframe(self, frame=""):
624 """
625 Set the frame type of the Direction on the sky.
626 Parameters:
627 frame: an optional frame type, default ''. Valid frames are:
628 'J2000', 'B1950', 'GALACTIC'
629 Examples:
630 scan.set_dirframe('GALACTIC')
631 """
632 varlist = vars()
633 try:
634 Scantable.set_dirframe(self, frame)
[1118]635 except RuntimeError, msg:
[989]636 if rcParams['verbose']:
637 print msg
638 else:
639 raise
[1118]640 self._add_history("set_dirframe", varlist)
[989]641
[113]642 def get_unit(self):
643 """
644 Get the default unit set in this scantable
645 Returns:
646 A unit string
647 """
648 inf = self._getcoordinfo()
649 unit = inf[0]
650 if unit == '': unit = 'channel'
651 return unit
[102]652
[158]653 def get_abcissa(self, rowno=0):
[102]654 """
[158]655 Get the abcissa in the current coordinate setup for the currently
[113]656 selected Beam/IF/Pol
657 Parameters:
[226]658 rowno: an optional row number in the scantable. Default is the
659 first row, i.e. rowno=0
[113]660 Returns:
[1348]661 The abcissa values and the format string (as a dictionary)
[113]662 """
[256]663 abc = self._getabcissa(rowno)
[710]664 lbl = self._getabcissalabel(rowno)
[158]665 return abc, lbl
[113]666
[1118]667 def flag(self, mask=None):
[1001]668 """
669 Flag the selected data using an optional channel mask.
670 Parameters:
671 mask: an optional channel mask, created with create_mask. Default
672 (no mask) is all channels.
673 """
674 varlist = vars()
[1593]675 mask = mask or []
[1001]676 try:
677 self._flag(mask)
[1118]678 except RuntimeError, msg:
[1001]679 if rcParams['verbose']:
680 print msg
681 return
682 else: raise
683 self._add_history("flag", varlist)
684
[1589]685 @print_log_dec
[1584]686 def lag_flag(self, start, end, unit="MHz", insitu=None):
[1192]687 """
688 Flag the data in 'lag' space by providing a frequency to remove.
[1584]689 Flagged data in the scantable gets interpolated over the region.
[1192]690 No taper is applied.
691 Parameters:
[1579]692 start: the start frequency (really a period within the
693 bandwidth) or period to remove
694 end: the end frequency or period to remove
[1584]695 unit: the frequency unit (default "MHz") or "" for
[1579]696 explicit lag channels
[1203]697 Notes:
[1579]698 It is recommended to flag edges of the band or strong
[1348]699 signals beforehand.
[1192]700 """
701 if insitu is None: insitu = rcParams['insitu']
702 self._math._setinsitu(insitu)
703 varlist = vars()
[1579]704 base = { "GHz": 1000000000., "MHz": 1000000., "kHz": 1000., "Hz": 1.}
705 if not (unit == "" or base.has_key(unit)):
[1192]706 raise ValueError("%s is not a valid unit." % unit)
707 try:
[1579]708 if unit == "":
709 s = scantable(self._math._lag_flag(self, start, end, "lags"))
710 else:
711 s = scantable(self._math._lag_flag(self, start*base[unit],
712 end*base[unit], "frequency"))
[1192]713 except RuntimeError, msg:
714 if rcParams['verbose']:
715 print msg
716 return
717 else: raise
718 s._add_history("lag_flag", varlist)
719 if insitu:
720 self._assign(s)
721 else:
722 return s
[1001]723
[1589]724 @print_log_dec
[113]725 def create_mask(self, *args, **kwargs):
726 """
[1118]727 Compute and return a mask based on [min, max] windows.
[189]728 The specified windows are to be INCLUDED, when the mask is
[113]729 applied.
[102]730 Parameters:
[1118]731 [min, max], [min2, max2], ...
[1024]732 Pairs of start/end points (inclusive)specifying the regions
[102]733 to be masked
[189]734 invert: optional argument. If specified as True,
735 return an inverted mask, i.e. the regions
736 specified are EXCLUDED
[513]737 row: create the mask using the specified row for
738 unit conversions, default is row=0
739 only necessary if frequency varies over rows.
[102]740 Example:
[113]741 scan.set_unit('channel')
742 a)
[1118]743 msk = scan.create_mask([400, 500], [800, 900])
[189]744 # masks everything outside 400 and 500
[113]745 # and 800 and 900 in the unit 'channel'
746
747 b)
[1118]748 msk = scan.create_mask([400, 500], [800, 900], invert=True)
[189]749 # masks the regions between 400 and 500
[113]750 # and 800 and 900 in the unit 'channel'
[1024]751 c)
752 mask only channel 400
[1554]753 msk = scan.create_mask([400])
[102]754 """
[1554]755 row = kwargs.get("row", 0)
[513]756 data = self._getabcissa(row)
[113]757 u = self._getcoordinfo()[0]
[718]758 if rcParams['verbose']:
[113]759 if u == "": u = "channel"
[718]760 msg = "The current mask window unit is %s" % u
[1118]761 i = self._check_ifs()
762 if not i:
[876]763 msg += "\nThis mask is only valid for IF=%d" % (self.getif(i))
[718]764 asaplog.push(msg)
[102]765 n = self.nchan()
[1295]766 msk = _n_bools(n, False)
[710]767 # test if args is a 'list' or a 'normal *args - UGLY!!!
768
[1118]769 ws = (isinstance(args[-1][-1], int) or isinstance(args[-1][-1], float)) \
770 and args or args[0]
[710]771 for window in ws:
[1554]772 if len(window) == 1:
773 window = [window[0], window[0]]
774 if len(window) == 0 or len(window) > 2:
775 raise ValueError("A window needs to be defined as [start(, end)]")
[1545]776 if window[0] > window[1]:
777 tmp = window[0]
778 window[0] = window[1]
779 window[1] = tmp
[102]780 for i in range(n):
[1024]781 if data[i] >= window[0] and data[i] <= window[1]:
[1295]782 msk[i] = True
[113]783 if kwargs.has_key('invert'):
784 if kwargs.get('invert'):
[1295]785 msk = mask_not(msk)
[102]786 return msk
[710]787
[256]788 def get_restfreqs(self):
789 """
790 Get the restfrequency(s) stored in this scantable.
791 The return value(s) are always of unit 'Hz'
792 Parameters:
793 none
794 Returns:
795 a list of doubles
796 """
797 return list(self._getrestfreqs())
[102]798
[402]799
[931]800 def set_restfreqs(self, freqs=None, unit='Hz'):
801 """
802 Set or replace the restfrequency specified and
803 If the 'freqs' argument holds a scalar,
804 then that rest frequency will be applied to all the selected
805 data. If the 'freqs' argument holds
806 a vector, then it MUST be of equal or smaller length than
807 the number of IFs (and the available restfrequencies will be
808 replaced by this vector). In this case, *all* data have
809 the restfrequency set per IF according
810 to the corresponding value you give in the 'freqs' vector.
[1118]811 E.g. 'freqs=[1e9, 2e9]' would mean IF 0 gets restfreq 1e9 and
[931]812 IF 1 gets restfreq 2e9.
[1395]813 You can also specify the frequencies via a linecatalog.
[1153]814
[931]815 Parameters:
816 freqs: list of rest frequency values or string idenitfiers
817 unit: unit for rest frequency (default 'Hz')
[402]818
[931]819 Example:
820 # set the given restfrequency for the whole table
821 scan.set_restfreqs(freqs=1.4e9)
[1348]822 # If thee number of IFs in the data is >= 2 IF0 gets the first
[931]823 # value IF1 the second...
[1118]824 scan.set_restfreqs(freqs=[1.4e9, 1.67e9])
[931]825 #set the given restfrequency for the whole table (by name)
826 scan.set_restfreqs(freqs="OH1667")
[391]827
[931]828 Note:
829 To do more sophisticate Restfrequency setting, e.g. on a
830 source and IF basis, use scantable.set_selection() before using
831 this function.
[1593]832 # provide your scantable is called scan
[931]833 selection = selector()
834 selection.set_name("ORION*")
835 selection.set_ifs([1])
836 scan.set_selection(selection)
837 scan.set_restfreqs(freqs=86.6e9)
838
839 """
840 varlist = vars()
[1157]841 from asap import linecatalog
842 # simple value
[1118]843 if isinstance(freqs, int) or isinstance(freqs, float):
[1153]844 self._setrestfreqs(freqs, "",unit)
[1157]845 # list of values
[1118]846 elif isinstance(freqs, list) or isinstance(freqs, tuple):
[1157]847 # list values are scalars
[1118]848 if isinstance(freqs[-1], int) or isinstance(freqs[-1], float):
[931]849 sel = selector()
850 savesel = self._getselection()
[1322]851 iflist = self.getifnos()
[931]852 for i in xrange(len(freqs)):
[1322]853 sel.set_ifs(iflist[i])
[931]854 self._setselection(sel)
[1153]855 self._setrestfreqs(freqs[i], "",unit)
[931]856 self._setselection(savesel)
[1157]857 # list values are tuples, (value, name)
858 elif isinstance(freqs[-1], dict):
859 sel = selector()
860 savesel = self._getselection()
[1322]861 iflist = self.getifnos()
[1157]862 for i in xrange(len(freqs)):
[1322]863 sel.set_ifs(iflist[i])
[1259]864 self._setselection(sel)
[1157]865 self._setrestfreqs(freqs[i]["value"],
866 freqs[i]["name"], "MHz")
867 self._setselection(savesel)
868 # freqs are to be taken from a linecatalog
[1153]869 elif isinstance(freqs, linecatalog):
870 sel = selector()
871 savesel = self._getselection()
[1322]872 iflist = self.getifnos()
[1153]873 for i in xrange(freqs.nrow()):
[1322]874 sel.set_ifs(iflist[i])
[1153]875 self._setselection(sel)
876 self._setrestfreqs(freqs.get_frequency(i),
877 freqs.get_name(i), "MHz")
878 # ensure that we are not iterating past nIF
879 if i == self.nif()-1: break
880 self._setselection(savesel)
[931]881 else:
882 return
883 self._add_history("set_restfreqs", varlist)
884
[1360]885 def shift_refpix(self, delta):
[1589]886 """
887 Shift the reference pixel of the Spectra Coordinate by an
888 integer amount.
889 Parameters:
890 delta: the amount to shift by
[1360]891 Note:
[1589]892 Be careful using this with broadband data.
[1360]893 """
[1589]894 Scantable.shift(self, delta)
[931]895
[1259]896 def history(self, filename=None):
897 """
898 Print the history. Optionally to a file.
[1348]899 Parameters:
900 filename: The name of the file to save the history to.
[1259]901 """
[484]902 hist = list(self._gethistory())
[794]903 out = "-"*80
[484]904 for h in hist:
[489]905 if h.startswith("---"):
[794]906 out += "\n"+h
[489]907 else:
908 items = h.split("##")
909 date = items[0]
910 func = items[1]
911 items = items[2:]
[794]912 out += "\n"+date+"\n"
913 out += "Function: %s\n Parameters:" % (func)
[489]914 for i in items:
915 s = i.split("=")
[1118]916 out += "\n %s = %s" % (s[0], s[1])
[794]917 out += "\n"+"-"*80
[1259]918 if filename is not None:
919 if filename is "":
920 filename = 'scantable_history.txt'
921 import os
922 filename = os.path.expandvars(os.path.expanduser(filename))
923 if not os.path.isdir(filename):
924 data = open(filename, 'w')
925 data.write(out)
926 data.close()
927 else:
928 msg = "Illegal file name '%s'." % (filename)
929 if rcParams['verbose']:
930 print msg
931 else:
932 raise IOError(msg)
933 if rcParams['verbose']:
934 try:
935 from IPython.genutils import page as pager
936 except ImportError:
937 from pydoc import pager
938 pager(out)
939 else:
940 return out
[484]941 return
[513]942 #
943 # Maths business
944 #
[1589]945 @print_log_dec
[931]946 def average_time(self, mask=None, scanav=False, weight='tint', align=False):
[513]947 """
[1070]948 Return the (time) weighted average of a scan.
[513]949 Note:
[1070]950 in channels only - align if necessary
[513]951 Parameters:
952 mask: an optional mask (only used for 'var' and 'tsys'
953 weighting)
[558]954 scanav: True averages each scan separately
955 False (default) averages all scans together,
[1099]956 weight: Weighting scheme.
957 'none' (mean no weight)
958 'var' (1/var(spec) weighted)
959 'tsys' (1/Tsys**2 weighted)
960 'tint' (integration time weighted)
961 'tintsys' (Tint/Tsys**2)
962 'median' ( median averaging)
[535]963 The default is 'tint'
[931]964 align: align the spectra in velocity before averaging. It takes
965 the time of the first spectrum as reference time.
[513]966 Example:
967 # time average the scantable without using a mask
[710]968 newscan = scan.average_time()
[513]969 """
970 varlist = vars()
[1593]971 weight = weight or 'TINT'
972 mask = mask or ()
973 scanav = (scanav and 'SCAN') or 'NONE'
[1118]974 scan = (self, )
[989]975 try:
[1118]976 if align:
977 scan = (self.freq_align(insitu=False), )
978 s = None
979 if weight.upper() == 'MEDIAN':
980 s = scantable(self._math._averagechannel(scan[0], 'MEDIAN',
981 scanav))
982 else:
983 s = scantable(self._math._average(scan, mask, weight.upper(),
984 scanav))
985 except RuntimeError, msg:
[989]986 if rcParams['verbose']:
987 print msg
988 return
989 else: raise
[1099]990 s._add_history("average_time", varlist)
[513]991 return s
[710]992
[1589]993 @print_log_dec
[876]994 def convert_flux(self, jyperk=None, eta=None, d=None, insitu=None):
[513]995 """
996 Return a scan where all spectra are converted to either
997 Jansky or Kelvin depending upon the flux units of the scan table.
998 By default the function tries to look the values up internally.
999 If it can't find them (or if you want to over-ride), you must
1000 specify EITHER jyperk OR eta (and D which it will try to look up
1001 also if you don't set it). jyperk takes precedence if you set both.
1002 Parameters:
1003 jyperk: the Jy / K conversion factor
1004 eta: the aperture efficiency
1005 d: the geomtric diameter (metres)
1006 insitu: if False a new scantable is returned.
1007 Otherwise, the scaling is done in-situ
1008 The default is taken from .asaprc (False)
1009 """
1010 if insitu is None: insitu = rcParams['insitu']
[876]1011 self._math._setinsitu(insitu)
[513]1012 varlist = vars()
[1593]1013 jyperk = jyperk or -1.0
1014 d = d or -1.0
1015 eta = eta or -1.0
[876]1016 s = scantable(self._math._convertflux(self, d, eta, jyperk))
1017 s._add_history("convert_flux", varlist)
1018 if insitu: self._assign(s)
1019 else: return s
[513]1020
[1589]1021 @print_log_dec
[876]1022 def gain_el(self, poly=None, filename="", method="linear", insitu=None):
[513]1023 """
1024 Return a scan after applying a gain-elevation correction.
1025 The correction can be made via either a polynomial or a
1026 table-based interpolation (and extrapolation if necessary).
1027 You specify polynomial coefficients, an ascii table or neither.
1028 If you specify neither, then a polynomial correction will be made
1029 with built in coefficients known for certain telescopes (an error
1030 will occur if the instrument is not known).
1031 The data and Tsys are *divided* by the scaling factors.
1032 Parameters:
1033 poly: Polynomial coefficients (default None) to compute a
1034 gain-elevation correction as a function of
1035 elevation (in degrees).
1036 filename: The name of an ascii file holding correction factors.
1037 The first row of the ascii file must give the column
1038 names and these MUST include columns
1039 "ELEVATION" (degrees) and "FACTOR" (multiply data
1040 by this) somewhere.
1041 The second row must give the data type of the
1042 column. Use 'R' for Real and 'I' for Integer.
1043 An example file would be
1044 (actual factors are arbitrary) :
1045
1046 TIME ELEVATION FACTOR
1047 R R R
1048 0.1 0 0.8
1049 0.2 20 0.85
1050 0.3 40 0.9
1051 0.4 60 0.85
1052 0.5 80 0.8
1053 0.6 90 0.75
1054 method: Interpolation method when correcting from a table.
1055 Values are "nearest", "linear" (default), "cubic"
1056 and "spline"
1057 insitu: if False a new scantable is returned.
1058 Otherwise, the scaling is done in-situ
1059 The default is taken from .asaprc (False)
1060 """
1061
1062 if insitu is None: insitu = rcParams['insitu']
[876]1063 self._math._setinsitu(insitu)
[513]1064 varlist = vars()
[1593]1065 poly = poly or ()
[513]1066 from os.path import expandvars
1067 filename = expandvars(filename)
[876]1068 s = scantable(self._math._gainel(self, poly, filename, method))
1069 s._add_history("gain_el", varlist)
[1593]1070 if insitu:
1071 self._assign(s)
1072 else:
1073 return s
[710]1074
[1589]1075 @print_log_dec
[931]1076 def freq_align(self, reftime=None, method='cubic', insitu=None):
[513]1077 """
1078 Return a scan where all rows have been aligned in frequency/velocity.
1079 The alignment frequency frame (e.g. LSRK) is that set by function
1080 set_freqframe.
1081 Parameters:
1082 reftime: reference time to align at. By default, the time of
1083 the first row of data is used.
1084 method: Interpolation method for regridding the spectra.
1085 Choose from "nearest", "linear", "cubic" (default)
1086 and "spline"
1087 insitu: if False a new scantable is returned.
1088 Otherwise, the scaling is done in-situ
1089 The default is taken from .asaprc (False)
1090 """
[931]1091 if insitu is None: insitu = rcParams["insitu"]
[876]1092 self._math._setinsitu(insitu)
[513]1093 varlist = vars()
[1593]1094 reftime = reftime or ""
[931]1095 s = scantable(self._math._freq_align(self, reftime, method))
[876]1096 s._add_history("freq_align", varlist)
1097 if insitu: self._assign(s)
1098 else: return s
[513]1099
[1589]1100 @print_log_dec
[876]1101 def opacity(self, tau, insitu=None):
[513]1102 """
1103 Apply an opacity correction. The data
1104 and Tsys are multiplied by the correction factor.
1105 Parameters:
1106 tau: Opacity from which the correction factor is
1107 exp(tau*ZD)
1108 where ZD is the zenith-distance
1109 insitu: if False a new scantable is returned.
1110 Otherwise, the scaling is done in-situ
1111 The default is taken from .asaprc (False)
1112 """
1113 if insitu is None: insitu = rcParams['insitu']
[876]1114 self._math._setinsitu(insitu)
[513]1115 varlist = vars()
[876]1116 s = scantable(self._math._opacity(self, tau))
1117 s._add_history("opacity", varlist)
1118 if insitu: self._assign(s)
1119 else: return s
[513]1120
[1589]1121 @print_log_dec
[513]1122 def bin(self, width=5, insitu=None):
1123 """
1124 Return a scan where all spectra have been binned up.
[1348]1125 Parameters:
[513]1126 width: The bin width (default=5) in pixels
1127 insitu: if False a new scantable is returned.
1128 Otherwise, the scaling is done in-situ
1129 The default is taken from .asaprc (False)
1130 """
1131 if insitu is None: insitu = rcParams['insitu']
[876]1132 self._math._setinsitu(insitu)
[513]1133 varlist = vars()
[876]1134 s = scantable(self._math._bin(self, width))
[1118]1135 s._add_history("bin", varlist)
[1589]1136 if insitu:
1137 self._assign(s)
1138 else:
1139 return s
[513]1140
[1589]1141 @print_log_dec
[513]1142 def resample(self, width=5, method='cubic', insitu=None):
1143 """
[1348]1144 Return a scan where all spectra have been binned up.
[1573]1145
[1348]1146 Parameters:
[513]1147 width: The bin width (default=5) in pixels
1148 method: Interpolation method when correcting from a table.
1149 Values are "nearest", "linear", "cubic" (default)
1150 and "spline"
1151 insitu: if False a new scantable is returned.
1152 Otherwise, the scaling is done in-situ
1153 The default is taken from .asaprc (False)
1154 """
1155 if insitu is None: insitu = rcParams['insitu']
[876]1156 self._math._setinsitu(insitu)
[513]1157 varlist = vars()
[876]1158 s = scantable(self._math._resample(self, method, width))
[1118]1159 s._add_history("resample", varlist)
[876]1160 if insitu: self._assign(s)
1161 else: return s
[513]1162
[1589]1163 @print_log_dec
[946]1164 def average_pol(self, mask=None, weight='none'):
1165 """
1166 Average the Polarisations together.
1167 Parameters:
1168 mask: An optional mask defining the region, where the
1169 averaging will be applied. The output will have all
1170 specified points masked.
1171 weight: Weighting scheme. 'none' (default), 'var' (1/var(spec)
1172 weighted), or 'tsys' (1/Tsys**2 weighted)
1173 """
1174 varlist = vars()
[1593]1175 mask = mask or ()
[1010]1176 s = scantable(self._math._averagepol(self, mask, weight.upper()))
[1118]1177 s._add_history("average_pol", varlist)
[992]1178 return s
[513]1179
[1589]1180 @print_log_dec
[1145]1181 def average_beam(self, mask=None, weight='none'):
1182 """
1183 Average the Beams together.
1184 Parameters:
1185 mask: An optional mask defining the region, where the
1186 averaging will be applied. The output will have all
1187 specified points masked.
1188 weight: Weighting scheme. 'none' (default), 'var' (1/var(spec)
1189 weighted), or 'tsys' (1/Tsys**2 weighted)
1190 """
1191 varlist = vars()
[1593]1192 mask = mask or ()
[1145]1193 s = scantable(self._math._averagebeams(self, mask, weight.upper()))
1194 s._add_history("average_beam", varlist)
1195 return s
1196
[1586]1197 def parallactify(self, pflag):
[1617]1198 """
1199 Set a flag to inidcate whether this data should be treated as having
1200 been 'parallactified' (total phase == 0.0)
1201 Parameters:
1202 pflag: Bool inidcating whether to turn this on (True) or
1203 off (False)
1204 """
[1586]1205 varlist = vars()
1206 self._parallactify(pflag)
1207 self._add_history("parallactify", varlist)
1208
[1589]1209 @print_log_dec
[992]1210 def convert_pol(self, poltype=None):
1211 """
1212 Convert the data to a different polarisation type.
[1565]1213 Note that you will need cross-polarisation terms for most conversions.
[992]1214 Parameters:
1215 poltype: The new polarisation type. Valid types are:
[1565]1216 "linear", "circular", "stokes" and "linpol"
[992]1217 """
1218 varlist = vars()
1219 try:
1220 s = scantable(self._math._convertpol(self, poltype))
[1118]1221 except RuntimeError, msg:
[992]1222 if rcParams['verbose']:
[1118]1223 print msg
1224 return
[992]1225 else:
1226 raise
[1118]1227 s._add_history("convert_pol", varlist)
[992]1228 return s
1229
[1589]1230 @print_log_dec
[1574]1231 def smooth(self, kernel="hanning", width=5.0, order=2, insitu=None):
[513]1232 """
1233 Smooth the spectrum by the specified kernel (conserving flux).
1234 Parameters:
1235 kernel: The type of smoothing kernel. Select from
[1574]1236 'hanning' (default), 'gaussian', 'boxcar', 'rmedian'
1237 or 'poly'
[513]1238 width: The width of the kernel in pixels. For hanning this is
1239 ignored otherwise it defauls to 5 pixels.
1240 For 'gaussian' it is the Full Width Half
1241 Maximum. For 'boxcar' it is the full width.
[1574]1242 For 'rmedian' and 'poly' it is the half width.
1243 order: Optional parameter for 'poly' kernel (default is 2), to
1244 specify the order of the polnomial. Ignored by all other
1245 kernels.
[513]1246 insitu: if False a new scantable is returned.
1247 Otherwise, the scaling is done in-situ
1248 The default is taken from .asaprc (False)
1249 Example:
1250 none
1251 """
1252 if insitu is None: insitu = rcParams['insitu']
[876]1253 self._math._setinsitu(insitu)
[513]1254 varlist = vars()
[1574]1255 s = scantable(self._math._smooth(self, kernel.lower(), width, order))
[876]1256 s._add_history("smooth", varlist)
1257 if insitu: self._assign(s)
1258 else: return s
[513]1259
[1589]1260 @print_log_dec
[1617]1261 def poly_baseline(self, mask=None, order=0, plot=False, uselin=False,
1262 insitu=None):
[513]1263 """
1264 Return a scan which has been baselined (all rows) by a polynomial.
1265 Parameters:
[794]1266 mask: an optional mask
1267 order: the order of the polynomial (default is 0)
[1061]1268 plot: plot the fit and the residual. In this each
1269 indivual fit has to be approved, by typing 'y'
1270 or 'n'
[1391]1271 uselin: use linear polynomial fit
[794]1272 insitu: if False a new scantable is returned.
1273 Otherwise, the scaling is done in-situ
1274 The default is taken from .asaprc (False)
[513]1275 Example:
1276 # return a scan baselined by a third order polynomial,
1277 # not using a mask
1278 bscan = scan.poly_baseline(order=3)
[579]1279 """
[513]1280 if insitu is None: insitu = rcParams['insitu']
1281 varlist = vars()
1282 if mask is None:
[1295]1283 mask = [True for i in xrange(self.nchan(-1))]
[513]1284 from asap.asapfitter import fitter
[1217]1285 try:
1286 f = fitter()
1287 f.set_scan(self, mask)
[1391]1288 #f.set_function(poly=order)
1289 if uselin:
1290 f.set_function(lpoly=order)
1291 else:
1292 f.set_function(poly=order)
[1217]1293 s = f.auto_fit(insitu, plot=plot)
1294 s._add_history("poly_baseline", varlist)
1295 if insitu: self._assign(s)
1296 else: return s
1297 except RuntimeError:
1298 msg = "The fit failed, possibly because it didn't converge."
1299 if rcParams['verbose']:
1300 print msg
1301 return
1302 else:
1303 raise RuntimeError(msg)
[513]1304
[1118]1305 def auto_poly_baseline(self, mask=[], edge=(0, 0), order=0,
[1280]1306 threshold=3, chan_avg_limit=1, plot=False,
1307 insitu=None):
[880]1308 """
1309 Return a scan which has been baselined (all rows) by a polynomial.
1310 Spectral lines are detected first using linefinder and masked out
1311 to avoid them affecting the baseline solution.
1312
1313 Parameters:
1314 mask: an optional mask retreived from scantable
1315 edge: an optional number of channel to drop at
1316 the edge of spectrum. If only one value is
1317 specified, the same number will be dropped from
1318 both sides of the spectrum. Default is to keep
[907]1319 all channels. Nested tuples represent individual
[976]1320 edge selection for different IFs (a number of spectral
1321 channels can be different)
[880]1322 order: the order of the polynomial (default is 0)
1323 threshold: the threshold used by line finder. It is better to
1324 keep it large as only strong lines affect the
1325 baseline solution.
[1280]1326 chan_avg_limit:
1327 a maximum number of consequtive spectral channels to
1328 average during the search of weak and broad lines.
1329 The default is no averaging (and no search for weak
1330 lines). If such lines can affect the fitted baseline
1331 (e.g. a high order polynomial is fitted), increase this
1332 parameter (usually values up to 8 are reasonable). Most
1333 users of this method should find the default value
1334 sufficient.
[1061]1335 plot: plot the fit and the residual. In this each
1336 indivual fit has to be approved, by typing 'y'
1337 or 'n'
[880]1338 insitu: if False a new scantable is returned.
1339 Otherwise, the scaling is done in-situ
1340 The default is taken from .asaprc (False)
1341
1342 Example:
1343 scan2=scan.auto_poly_baseline(order=7)
1344 """
1345 if insitu is None: insitu = rcParams['insitu']
1346 varlist = vars()
1347 from asap.asapfitter import fitter
1348 from asap.asaplinefind import linefinder
1349 from asap import _is_sequence_or_number as _is_valid
1350
[976]1351 # check whether edge is set up for each IF individually
[1118]1352 individualedge = False;
1353 if len(edge) > 1:
1354 if isinstance(edge[0], list) or isinstance(edge[0], tuple):
1355 individualedge = True;
[907]1356
[1118]1357 if not _is_valid(edge, int) and not individualedge:
[909]1358 raise ValueError, "Parameter 'edge' has to be an integer or a \
[907]1359 pair of integers specified as a tuple. Nested tuples are allowed \
1360 to make individual selection for different IFs."
[919]1361
[1118]1362 curedge = (0, 0)
1363 if individualedge:
1364 for edgepar in edge:
1365 if not _is_valid(edgepar, int):
1366 raise ValueError, "Each element of the 'edge' tuple has \
1367 to be a pair of integers or an integer."
[907]1368 else:
[1118]1369 curedge = edge;
[880]1370
1371 # setup fitter
1372 f = fitter()
1373 f.set_function(poly=order)
1374
1375 # setup line finder
[1118]1376 fl = linefinder()
[1268]1377 fl.set_options(threshold=threshold,avg_limit=chan_avg_limit)
[880]1378
1379 if not insitu:
[1118]1380 workscan = self.copy()
[880]1381 else:
[1118]1382 workscan = self
[880]1383
[907]1384 fl.set_scan(workscan)
1385
[1118]1386 rows = range(workscan.nrow())
[880]1387 asaplog.push("Processing:")
1388 for r in rows:
[1118]1389 msg = " Scan[%d] Beam[%d] IF[%d] Pol[%d] Cycle[%d]" % \
1390 (workscan.getscan(r), workscan.getbeam(r), workscan.getif(r), \
1391 workscan.getpol(r), workscan.getcycle(r))
[880]1392 asaplog.push(msg, False)
[907]1393
[976]1394 # figure out edge parameter
[1118]1395 if individualedge:
1396 if len(edge) >= workscan.getif(r):
1397 raise RuntimeError, "Number of edge elements appear to " \
1398 "be less than the number of IFs"
1399 curedge = edge[workscan.getif(r)]
[919]1400
[976]1401 # setup line finder
[1118]1402 fl.find_lines(r, mask, curedge)
[1536]1403 f.set_data(workscan._getabcissa(r), workscan._getspectrum(r),
1404 mask_and(workscan._getmask(r), fl.get_mask()))
[880]1405 f.fit()
1406 x = f.get_parameters()
[1061]1407 if plot:
1408 f.plot(residual=True)
1409 x = raw_input("Accept fit ( [y]/n ): ")
1410 if x.upper() == 'N':
1411 continue
[880]1412 workscan._setspectrum(f.fitter.getresidual(), r)
[1061]1413 if plot:
1414 f._p.unmap()
1415 f._p = None
1416 workscan._add_history("auto_poly_baseline", varlist)
[880]1417 if insitu:
1418 self._assign(workscan)
1419 else:
1420 return workscan
1421
[1589]1422 @print_log_dec
[914]1423 def rotate_linpolphase(self, angle):
1424 """
1425 Rotate the phase of the complex polarization O=Q+iU correlation.
1426 This is always done in situ in the raw data. So if you call this
1427 function more than once then each call rotates the phase further.
1428 Parameters:
1429 angle: The angle (degrees) to rotate (add) by.
1430 Examples:
1431 scan.rotate_linpolphase(2.3)
1432 """
1433 varlist = vars()
[936]1434 self._math._rotate_linpolphase(self, angle)
[914]1435 self._add_history("rotate_linpolphase", varlist)
1436 return
[710]1437
[1589]1438 @print_log_dec
[914]1439 def rotate_xyphase(self, angle):
1440 """
1441 Rotate the phase of the XY correlation. This is always done in situ
1442 in the data. So if you call this function more than once
1443 then each call rotates the phase further.
1444 Parameters:
1445 angle: The angle (degrees) to rotate (add) by.
1446 Examples:
1447 scan.rotate_xyphase(2.3)
1448 """
1449 varlist = vars()
[936]1450 self._math._rotate_xyphase(self, angle)
[914]1451 self._add_history("rotate_xyphase", varlist)
1452 return
1453
[1589]1454 @print_log_dec
[914]1455 def swap_linears(self):
1456 """
[1573]1457 Swap the linear polarisations XX and YY, or better the first two
[1348]1458 polarisations as this also works for ciculars.
[914]1459 """
1460 varlist = vars()
[936]1461 self._math._swap_linears(self)
[914]1462 self._add_history("swap_linears", varlist)
1463 return
1464
[1589]1465 @print_log_dec
[914]1466 def invert_phase(self):
1467 """
1468 Invert the phase of the complex polarisation
1469 """
1470 varlist = vars()
[936]1471 self._math._invert_phase(self)
[914]1472 self._add_history("invert_phase", varlist)
1473 return
1474
[1589]1475 @print_log_dec
[876]1476 def add(self, offset, insitu=None):
[513]1477 """
1478 Return a scan where all spectra have the offset added
1479 Parameters:
1480 offset: the offset
1481 insitu: if False a new scantable is returned.
1482 Otherwise, the scaling is done in-situ
1483 The default is taken from .asaprc (False)
1484 """
1485 if insitu is None: insitu = rcParams['insitu']
[876]1486 self._math._setinsitu(insitu)
[513]1487 varlist = vars()
[876]1488 s = scantable(self._math._unaryop(self, offset, "ADD", False))
[1118]1489 s._add_history("add", varlist)
[876]1490 if insitu:
1491 self._assign(s)
1492 else:
[513]1493 return s
1494
[1589]1495 @print_log_dec
[1308]1496 def scale(self, factor, tsys=True, insitu=None):
[513]1497 """
1498 Return a scan where all spectra are scaled by the give 'factor'
1499 Parameters:
1500 factor: the scaling factor
1501 insitu: if False a new scantable is returned.
1502 Otherwise, the scaling is done in-situ
1503 The default is taken from .asaprc (False)
1504 tsys: if True (default) then apply the operation to Tsys
1505 as well as the data
1506 """
1507 if insitu is None: insitu = rcParams['insitu']
[876]1508 self._math._setinsitu(insitu)
[513]1509 varlist = vars()
[876]1510 s = scantable(self._math._unaryop(self, factor, "MUL", tsys))
[1118]1511 s._add_history("scale", varlist)
[876]1512 if insitu:
1513 self._assign(s)
1514 else:
[513]1515 return s
1516
[1504]1517 def set_sourcetype(self, match, matchtype="pattern",
1518 sourcetype="reference"):
[1502]1519 """
1520 Set the type of the source to be an source or reference scan
1521 using the provided pattern:
1522 Parameters:
[1504]1523 match: a Unix style pattern, regular expression or selector
1524 matchtype: 'pattern' (default) UNIX style pattern or
1525 'regex' regular expression
[1502]1526 sourcetype: the type of the source to use (source/reference)
1527 """
1528 varlist = vars()
1529 basesel = self.get_selection()
1530 stype = -1
1531 if sourcetype.lower().startswith("r"):
1532 stype = 1
1533 elif sourcetype.lower().startswith("s"):
1534 stype = 0
[1504]1535 else:
[1502]1536 raise ValueError("Illegal sourcetype use s(ource) or r(eference)")
[1504]1537 if matchtype.lower().startswith("p"):
1538 matchtype = "pattern"
1539 elif matchtype.lower().startswith("r"):
1540 matchtype = "regex"
1541 else:
1542 raise ValueError("Illegal matchtype, use p(attern) or r(egex)")
[1502]1543 sel = selector()
1544 if isinstance(match, selector):
1545 sel = match
1546 else:
[1504]1547 sel.set_query("SRCNAME == %s('%s')" % (matchtype, match))
[1502]1548 self.set_selection(basesel+sel)
1549 self._setsourcetype(stype)
1550 self.set_selection(basesel)
[1573]1551 self._add_history("set_sourcetype", varlist)
[1502]1552
[1589]1553 @print_log_dec
[1348]1554 def auto_quotient(self, preserve=True, mode='paired'):
[670]1555 """
1556 This function allows to build quotients automatically.
1557 It assumes the observation to have the same numer of
1558 "ons" and "offs"
1559 Parameters:
[710]1560 preserve: you can preserve (default) the continuum or
1561 remove it. The equations used are
[670]1562 preserve: Output = Toff * (on/off) - Toff
[1070]1563 remove: Output = Toff * (on/off) - Ton
[1573]1564 mode: the on/off detection mode
[1348]1565 'paired' (default)
1566 identifies 'off' scans by the
1567 trailing '_R' (Mopra/Parkes) or
1568 '_e'/'_w' (Tid) and matches
1569 on/off pairs from the observing pattern
[1502]1570 'time'
1571 finds the closest off in time
[1348]1572
[670]1573 """
[1348]1574 modes = ["time", "paired"]
[670]1575 if not mode in modes:
[876]1576 msg = "please provide valid mode. Valid modes are %s" % (modes)
1577 raise ValueError(msg)
1578 varlist = vars()
[1348]1579 s = None
1580 if mode.lower() == "paired":
1581 basesel = self.get_selection()
[1356]1582 sel = selector()+basesel
1583 sel.set_query("SRCTYPE==1")
1584 self.set_selection(sel)
[1348]1585 offs = self.copy()
1586 sel.set_query("SRCTYPE==0")
[1356]1587 self.set_selection(sel)
[1348]1588 ons = self.copy()
1589 s = scantable(self._math._quotient(ons, offs, preserve))
1590 self.set_selection(basesel)
1591 elif mode.lower() == "time":
1592 s = scantable(self._math._auto_quotient(self, mode, preserve))
[1118]1593 s._add_history("auto_quotient", varlist)
[876]1594 return s
[710]1595
[1589]1596 @print_log_dec
[1145]1597 def mx_quotient(self, mask = None, weight='median', preserve=True):
[1141]1598 """
[1143]1599 Form a quotient using "off" beams when observing in "MX" mode.
1600 Parameters:
[1145]1601 mask: an optional mask to be used when weight == 'stddev'
[1143]1602 weight: How to average the off beams. Default is 'median'.
[1145]1603 preserve: you can preserve (default) the continuum or
1604 remove it. The equations used are
1605 preserve: Output = Toff * (on/off) - Toff
1606 remove: Output = Toff * (on/off) - Ton
[1217]1607 """
[1593]1608 mask = mask or ()
[1141]1609 varlist = vars()
1610 on = scantable(self._math._mx_extract(self, 'on'))
[1143]1611 preoff = scantable(self._math._mx_extract(self, 'off'))
1612 off = preoff.average_time(mask=mask, weight=weight, scanav=False)
[1217]1613 from asapmath import quotient
[1145]1614 q = quotient(on, off, preserve)
[1143]1615 q._add_history("mx_quotient", varlist)
[1217]1616 return q
[513]1617
[1589]1618 @print_log_dec
[718]1619 def freq_switch(self, insitu=None):
1620 """
1621 Apply frequency switching to the data.
1622 Parameters:
1623 insitu: if False a new scantable is returned.
1624 Otherwise, the swictching is done in-situ
1625 The default is taken from .asaprc (False)
1626 Example:
1627 none
1628 """
1629 if insitu is None: insitu = rcParams['insitu']
[876]1630 self._math._setinsitu(insitu)
[718]1631 varlist = vars()
[876]1632 s = scantable(self._math._freqswitch(self))
[1118]1633 s._add_history("freq_switch", varlist)
[876]1634 if insitu: self._assign(s)
1635 else: return s
[718]1636
[1589]1637 @print_log_dec
[780]1638 def recalc_azel(self):
1639 """
1640 Recalculate the azimuth and elevation for each position.
1641 Parameters:
1642 none
1643 Example:
1644 """
1645 varlist = vars()
[876]1646 self._recalcazel()
[780]1647 self._add_history("recalc_azel", varlist)
1648 return
1649
[1589]1650 @print_log_dec
[513]1651 def __add__(self, other):
1652 varlist = vars()
1653 s = None
1654 if isinstance(other, scantable):
[1573]1655 s = scantable(self._math._binaryop(self, other, "ADD"))
[513]1656 elif isinstance(other, float):
[876]1657 s = scantable(self._math._unaryop(self, other, "ADD", False))
[513]1658 else:
[718]1659 raise TypeError("Other input is not a scantable or float value")
[513]1660 s._add_history("operator +", varlist)
1661 return s
1662
[1589]1663 @print_log_dec
[513]1664 def __sub__(self, other):
1665 """
1666 implicit on all axes and on Tsys
1667 """
1668 varlist = vars()
1669 s = None
1670 if isinstance(other, scantable):
[1588]1671 s = scantable(self._math._binaryop(self, other, "SUB"))
[513]1672 elif isinstance(other, float):
[876]1673 s = scantable(self._math._unaryop(self, other, "SUB", False))
[513]1674 else:
[718]1675 raise TypeError("Other input is not a scantable or float value")
[513]1676 s._add_history("operator -", varlist)
1677 return s
[710]1678
[1589]1679 @print_log_dec
[513]1680 def __mul__(self, other):
1681 """
1682 implicit on all axes and on Tsys
1683 """
1684 varlist = vars()
1685 s = None
1686 if isinstance(other, scantable):
[1588]1687 s = scantable(self._math._binaryop(self, other, "MUL"))
[513]1688 elif isinstance(other, float):
[876]1689 s = scantable(self._math._unaryop(self, other, "MUL", False))
[513]1690 else:
[718]1691 raise TypeError("Other input is not a scantable or float value")
[513]1692 s._add_history("operator *", varlist)
1693 return s
1694
[710]1695
[1589]1696 @print_log_dec
[513]1697 def __div__(self, other):
1698 """
1699 implicit on all axes and on Tsys
1700 """
1701 varlist = vars()
1702 s = None
1703 if isinstance(other, scantable):
[1589]1704 s = scantable(self._math._binaryop(self, other, "DIV"))
[513]1705 elif isinstance(other, float):
1706 if other == 0.0:
[718]1707 raise ZeroDivisionError("Dividing by zero is not recommended")
[876]1708 s = scantable(self._math._unaryop(self, other, "DIV", False))
[513]1709 else:
[718]1710 raise TypeError("Other input is not a scantable or float value")
[513]1711 s._add_history("operator /", varlist)
1712 return s
1713
[530]1714 def get_fit(self, row=0):
1715 """
1716 Print or return the stored fits for a row in the scantable
1717 Parameters:
1718 row: the row which the fit has been applied to.
1719 """
1720 if row > self.nrow():
1721 return
[976]1722 from asap.asapfit import asapfit
[530]1723 fit = asapfit(self._getfit(row))
[718]1724 if rcParams['verbose']:
[530]1725 print fit
1726 return
1727 else:
1728 return fit.as_dict()
1729
[1483]1730 def flag_nans(self):
1731 """
1732 Utility function to flag NaN values in the scantable.
1733 """
1734 import numpy
1735 basesel = self.get_selection()
1736 for i in range(self.nrow()):
[1589]1737 sel = self.get_row_selector(i)
1738 self.set_selection(basesel+sel)
[1483]1739 nans = numpy.isnan(self._getspectrum(0))
1740 if numpy.any(nans):
1741 bnans = [ bool(v) for v in nans]
1742 self.flag(bnans)
1743 self.set_selection(basesel)
1744
[1588]1745 def get_row_selector(self, rowno):
1746 return selector(beams=self.getbeam(rowno),
1747 ifs=self.getif(rowno),
1748 pols=self.getpol(rowno),
1749 scans=self.getscan(rowno),
1750 cycles=self.getcycle(rowno))
[1573]1751
[484]1752 def _add_history(self, funcname, parameters):
[1435]1753 if not rcParams['scantable.history']:
1754 return
[484]1755 # create date
1756 sep = "##"
1757 from datetime import datetime
1758 dstr = datetime.now().strftime('%Y/%m/%d %H:%M:%S')
1759 hist = dstr+sep
1760 hist += funcname+sep#cdate+sep
1761 if parameters.has_key('self'): del parameters['self']
[1118]1762 for k, v in parameters.iteritems():
[484]1763 if type(v) is dict:
[1118]1764 for k2, v2 in v.iteritems():
[484]1765 hist += k2
1766 hist += "="
[1118]1767 if isinstance(v2, scantable):
[484]1768 hist += 'scantable'
1769 elif k2 == 'mask':
[1118]1770 if isinstance(v2, list) or isinstance(v2, tuple):
[513]1771 hist += str(self._zip_mask(v2))
1772 else:
1773 hist += str(v2)
[484]1774 else:
[513]1775 hist += str(v2)
[484]1776 else:
1777 hist += k
1778 hist += "="
[1118]1779 if isinstance(v, scantable):
[484]1780 hist += 'scantable'
1781 elif k == 'mask':
[1118]1782 if isinstance(v, list) or isinstance(v, tuple):
[513]1783 hist += str(self._zip_mask(v))
1784 else:
1785 hist += str(v)
[484]1786 else:
1787 hist += str(v)
1788 hist += sep
1789 hist = hist[:-2] # remove trailing '##'
1790 self._addhistory(hist)
1791
[710]1792
[484]1793 def _zip_mask(self, mask):
1794 mask = list(mask)
1795 i = 0
1796 segments = []
1797 while mask[i:].count(1):
1798 i += mask[i:].index(1)
1799 if mask[i:].count(0):
1800 j = i + mask[i:].index(0)
1801 else:
[710]1802 j = len(mask)
[1118]1803 segments.append([i, j])
[710]1804 i = j
[484]1805 return segments
[714]1806
[626]1807 def _get_ordinate_label(self):
1808 fu = "("+self.get_fluxunit()+")"
1809 import re
1810 lbl = "Intensity"
[1118]1811 if re.match(".K.", fu):
[626]1812 lbl = "Brightness Temperature "+ fu
[1118]1813 elif re.match(".Jy.", fu):
[626]1814 lbl = "Flux density "+ fu
1815 return lbl
[710]1816
[876]1817 def _check_ifs(self):
1818 nchans = [self.nchan(i) for i in range(self.nif(-1))]
[889]1819 nchans = filter(lambda t: t > 0, nchans)
[876]1820 return (sum(nchans)/len(nchans) == nchans[0])
[976]1821
1822 def _fill(self, names, unit, average):
1823 import os
1824 from asap._asap import stfiller
1825 first = True
1826 fullnames = []
1827 for name in names:
1828 name = os.path.expandvars(name)
1829 name = os.path.expanduser(name)
1830 if not os.path.exists(name):
1831 msg = "File '%s' does not exists" % (name)
1832 if rcParams['verbose']:
1833 asaplog.push(msg)
1834 print asaplog.pop().strip()
1835 return
1836 raise IOError(msg)
1837 fullnames.append(name)
1838 if average:
1839 asaplog.push('Auto averaging integrations')
[1079]1840 stype = int(rcParams['scantable.storage'].lower() == 'disk')
[976]1841 for name in fullnames:
[1073]1842 tbl = Scantable(stype)
1843 r = stfiller(tbl)
[1504]1844 rx = rcParams['scantable.reference']
1845 r._setreferenceexpr(rx)
[976]1846 msg = "Importing %s..." % (name)
[1118]1847 asaplog.push(msg, False)
[976]1848 print_log()
[1118]1849 r._open(name, -1, -1)
[976]1850 r._read()
1851 if average:
[1118]1852 tbl = self._math._average((tbl, ), (), 'NONE', 'SCAN')
[976]1853 if not first:
1854 tbl = self._math._merge([self, tbl])
1855 Scantable.__init__(self, tbl)
1856 r._close()
[1118]1857 del r, tbl
[976]1858 first = False
1859 if unit is not None:
1860 self.set_fluxunit(unit)
1861 self.set_freqframe(rcParams['scantable.freqframe'])
1862
[1402]1863 def __getitem__(self, key):
1864 if key < 0:
1865 key += self.nrow()
1866 if key >= self.nrow():
1867 raise IndexError("Row index out of range.")
1868 return self._getspectrum(key)
1869
1870 def __setitem__(self, key, value):
1871 if key < 0:
1872 key += self.nrow()
1873 if key >= self.nrow():
1874 raise IndexError("Row index out of range.")
1875 if not hasattr(value, "__len__") or \
1876 len(value) > self.nchan(self.getif(key)):
1877 raise ValueError("Spectrum length doesn't match.")
1878 return self._setspectrum(value, key)
1879
1880 def __len__(self):
1881 return self.nrow()
1882
1883 def __iter__(self):
1884 for i in range(len(self)):
1885 yield self[i]
Note: See TracBrowser for help on using the repository browser.