source: branches/Release12/python/scantable.py@ 2468

Last change on this file since 2468 was 800, checked in by mar637, 19 years ago

Request: added channel based flagging
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 63.9 KB
RevLine 
[102]1from asap._asap import sdtable
[226]2from asap import rcParams
[714]3from asap import print_log
[189]4from numarray import ones,zeros
[102]5import sys
6
7class scantable(sdtable):
8 """
9 The ASAP container for scans
10 """
[710]11
[484]12 def __init__(self, filename, average=None, unit=None):
[102]13 """
14 Create a scantable from a saved one or make a reference
15 Parameters:
[181]16 filename: the name of an asap table on disk
17 or
18 the name of a rpfits/sdfits/ms file
19 (integrations within scans are auto averaged
20 and the whole file is read)
21 or
22 [advanced] a reference to an existing
[102]23 scantable
[484]24 average: average all integrations withinb a scan on read.
25 The default (True) is taken from .asaprc.
26 unit: brightness unit; must be consistent with K or Jy.
[340]27 Over-rides the default selected by the reader
28 (input rpfits/sdfits/ms) or replaces the value
29 in existing scantables
[710]30 """
[489]31 if average is None or type(average) is not bool:
[710]32 average = rcParams['scantable.autoaverage']
[489]33
[484]34 varlist = vars()
[113]35 self._p = None
[710]36 from asap import asaplog
[181]37 if isinstance(filename,sdtable):
[710]38 sdtable.__init__(self, filename)
[340]39 if unit is not None:
[710]40 self.set_fluxunit(unit)
[181]41 else:
[411]42 import os.path
43 if not os.path.exists(filename):
[710]44 s = "File '%s' not found." % (filename)
45 if rcParams['verbose']:
46 asaplog.push(s)
47 print asaplog.pop().strip()
48 return
49 raise IOError(s)
[411]50 filename = os.path.expandvars(filename)
51 if os.path.isdir(filename):
[181]52 # crude check if asap table
[411]53 if os.path.exists(filename+'/table.info'):
[181]54 sdtable.__init__(self, filename)
[340]55 if unit is not None:
[710]56 self.set_fluxunit(unit)
[181]57 else:
[718]58 msg = "The given file '%s'is not a valid asap table." % (filename)
59 if rcParams['verbose']:
60 print msg
61 return
62 else:
63 raise IOError(msg)
[710]64 else:
[181]65 from asap._asap import sdreader
[340]66 ifSel = -1
67 beamSel = -1
[710]68 r = sdreader()
69 r._open(filename,ifSel,beamSel)
70 asaplog.push('Importing data...')
[760]71 print_log()
[411]72 r._read([-1])
73 tbl = r._getdata()
[346]74 if unit is not None:
75 tbl.set_fluxunit(unit)
[489]76 if average:
77 from asap._asap import average as _av
[710]78 asaplog.push('Auto averaging integrations...')
[760]79 print_log()
[513]80 tbl2 = _av((tbl,),(),True,'none')
[226]81 sdtable.__init__(self,tbl2)
[513]82 del tbl2
[226]83 else:
84 sdtable.__init__(self,tbl)
[513]85 del r,tbl
[484]86 self._add_history("scantable", varlist)
[714]87 print_log()
[102]88
[451]89 def save(self, name=None, format=None, stokes=False, overwrite=False):
[116]90 """
[710]91 Store the scantable on disk. This can be an asap (aips++) Table, SDFITS,
[280]92 Image FITS or MS2 format.
[116]93 Parameters:
[196]94 name: the name of the outputfile. For format="FITS" this
[489]95 is the directory file name into which all the files
[497]96 will be written (default is 'asap_FITS'). For format
97 "ASCII" this is the root file name (data in 'name'.txt
98 and header in 'name'_header.txt)
[116]99 format: an optional file format. Default is ASAP.
[280]100 Allowed are - 'ASAP' (save as ASAP [aips++] Table),
[194]101 'SDFITS' (save as SDFITS file)
[710]102 'FITS' (saves each row as a FITS Image)
[200]103 'ASCII' (saves as ascii text file)
[226]104 'MS2' (saves as an aips++
105 MeasurementSet V2)
[489]106 stokes: Convert to Stokes parameters (only available
107 currently with FITS and ASCII formats.
108 Default is False.
[411]109 overwrite: If the file should be overwritten if it exists.
[256]110 The default False is to return with warning
[411]111 without writing the output. USE WITH CARE.
[116]112 Example:
113 scan.save('myscan.asap')
114 scan.save('myscan.sdfits','SDFITS')
115 """
[411]116 from os import path
[226]117 if format is None: format = rcParams['scantable.save']
[256]118 suffix = '.'+format.lower()
119 if name is None or name =="":
120 name = 'scantable'+suffix
[718]121 from asap import asaplog
122 msg = "No filename given. Using default name %s..." % name
123 asaplog.push(msg)
[411]124 name = path.expandvars(name)
[256]125 if path.isfile(name) or path.isdir(name):
126 if not overwrite:
[718]127 msg = "File %s exists." % name
128 if rcParams['verbose']:
129 print msg
130 return
131 else:
132 raise IOError(msg)
[451]133 format2 = format.upper()
134 if format2 == 'ASAP':
[116]135 self._save(name)
136 else:
137 from asap._asap import sdwriter as _sw
[451]138 w = _sw(format2)
[445]139 w.write(self, name, stokes)
[718]140
141 print_log()
[116]142 return
143
[102]144 def copy(self):
145 """
146 Return a copy of this scantable.
147 Parameters:
[113]148 none
[102]149 Example:
150 copiedscan = scan.copy()
151 """
[113]152 sd = scantable(sdtable._copy(self))
153 return sd
154
[102]155 def get_scan(self, scanid=None):
156 """
157 Return a specific scan (by scanno) or collection of scans (by
158 source name) in a new scantable.
159 Parameters:
[513]160 scanid: a (list of) scanno or a source name, unix-style
161 patterns are accepted for source name matching, e.g.
162 '*_R' gets all 'ref scans
[102]163 Example:
[513]164 # get all scans containing the source '323p459'
165 newscan = scan.get_scan('323p459')
166 # get all 'off' scans
167 refscans = scan.get_scan('*_R')
168 # get a susbset of scans by scanno (as listed in scan.summary())
169 newscan = scan.get_scan([0,2,7,10])
[102]170 """
171 if scanid is None:
[718]172 if rcParams['verbose']:
173 print "Please specify a scan no or name to retrieve from the scantable"
174 return
175 else:
176 raise RuntimeError("No scan given")
177
[102]178 try:
179 if type(scanid) is str:
[113]180 s = sdtable._getsource(self,scanid)
[102]181 return scantable(s)
182 elif type(scanid) is int:
[381]183 s = sdtable._getscan(self,[scanid])
184 return scantable(s)
185 elif type(scanid) is list:
[113]186 s = sdtable._getscan(self,scanid)
[102]187 return scantable(s)
[381]188 else:
[718]189 msg = "Illegal scanid type, use 'int' or 'list' if ints."
190 if rcParams['verbose']:
191 print msg
192 else:
193 raise TypeError(msg)
[102]194 except RuntimeError:
[718]195 if rcParams['verbose']: print "Couldn't find any match."
196 else: raise
[102]197
198 def __str__(self):
[381]199 return sdtable._summary(self,True)
[102]200
[381]201 def summary(self,filename=None, verbose=None):
[102]202 """
203 Print a summary of the contents of this scantable.
204 Parameters:
205 filename: the name of a file to write the putput to
206 Default - no file output
[381]207 verbose: print extra info such as the frequency table
208 The default (False) is taken from .asaprc
[102]209 """
[381]210 info = sdtable._summary(self, verbose)
211 if verbose is None: verbose = rcParams['scantable.verbosesummary']
[102]212 if filename is not None:
[256]213 if filename is "":
214 filename = 'scantable_summary.txt'
[415]215 from os.path import expandvars, isdir
[411]216 filename = expandvars(filename)
[415]217 if not isdir(filename):
[413]218 data = open(filename, 'w')
219 data.write(info)
220 data.close()
221 else:
[718]222 msg = "Illegal file name '%s'." % (filename)
223 if rcParams['verbose']:
224 print msg
225 else:
226 raise IOError(msg)
227 if rcParams['verbose']:
[791]228 try:
229 from IPython.genutils import page as pager
230 except ImportError:
231 from pydoc import pager
232 pager(info)
[718]233 else:
234 return info
[710]235
[553]236 def set_cursor(self, beam=0, IF=0, pol=0):
[102]237 """
238 Set the spectrum for individual operations.
239 Parameters:
[553]240 beam, IF, pol: a number
[102]241 Example:
[256]242 scan.set_cursor(0,0,1)
[135]243 pol1sig = scan.stats(all=False) # returns std dev for beam=0
[181]244 # if=0, pol=1
[102]245 """
[484]246 varlist = vars()
[553]247 self.setbeam(beam)
248 self.setpol(pol)
249 self.setif(IF)
[484]250 self._add_history("set_cursor",varlist)
[102]251 return
252
[256]253 def get_cursor(self):
[113]254 """
255 Return/print a the current 'cursor' into the Beam/IF/Pol cube.
256 Parameters:
257 none
258 Returns:
259 a list of values (currentBeam,currentIF,currentPol)
260 Example:
261 none
262 """
[102]263 i = self.getbeam()
264 j = self.getif()
265 k = self.getpol()
[721]266 from asap import asaplog
267 out = "--------------------------------------------------\n"
268 out += " Cursor position\n"
269 out += "--------------------------------------------------\n"
270 out += 'Beam=%d IF=%d Pol=%d ' % (i,j,k)
271 asaplog.push(out)
272 print_log()
273 return i,j,k
[102]274
[436]275 def stats(self, stat='stddev', mask=None, allaxes=None):
[102]276 """
[135]277 Determine the specified statistic of the current beam/if/pol
[102]278 Takes a 'mask' as an optional parameter to specify which
279 channels should be excluded.
280 Parameters:
[135]281 stat: 'min', 'max', 'sumsq', 'sum', 'mean'
282 'var', 'stddev', 'avdev', 'rms', 'median'
283 mask: an optional mask specifying where the statistic
[102]284 should be determined.
[436]285 allaxes: if True apply to all spectra. Otherwise
286 apply only to the selected (beam/pol/if)spectra only.
287 The default is taken from .asaprc (True if none)
[102]288 Example:
[113]289 scan.set_unit('channel')
[102]290 msk = scan.create_mask([100,200],[500,600])
[135]291 scan.stats(stat='mean', mask=m)
[102]292 """
[436]293 if allaxes is None: allaxes = rcParams['scantable.allaxes']
[135]294 from asap._asap import stats as _stats
[256]295 from numarray import array,zeros,Float
[102]296 if mask == None:
297 mask = ones(self.nchan())
[256]298 axes = ['Beam','IF','Pol','Time']
299
[436]300 beamSel,IFSel,polSel = (self.getbeam(),self.getif(),self.getpol())
301 if allaxes:
[256]302 n = self.nbeam()*self.nif()*self.npol()*self.nrow()
303 shp = [self.nbeam(),self.nif(),self.npol(),self.nrow()]
304 arr = array(zeros(n),shape=shp,type=Float)
305
306 for i in range(self.nbeam()):
307 self.setbeam(i)
308 for j in range(self.nif()):
309 self.setif(j)
310 for k in range(self.npol()):
311 self.setpol(k)
312 arr[i,j,k,:] = _stats(self,mask,stat,-1)
313 retval = {'axes': axes, 'data': arr, 'cursor':None}
314 tm = [self._gettime(val) for val in range(self.nrow())]
[718]315 if rcParams['verbose']:
[256]316 self._print_values(retval,stat,tm)
[436]317 self.setbeam(beamSel)
318 self.setif(IFSel)
319 self.setpol(polSel)
[256]320 return retval
[102]321
322 else:
[241]323 statval = _stats(self,mask,stat,-1)
[181]324 out = ''
325 for l in range(self.nrow()):
326 tm = self._gettime(l)
327 out += 'Time[%s]:\n' % (tm)
[436]328 if self.nbeam() > 1: out += ' Beam[%d] ' % (beamSel)
329 if self.nif() > 1: out += ' IF[%d] ' % (IFSel)
330 if self.npol() > 1: out += ' Pol[%d] ' % (polSel)
[181]331 out += '= %3.3f\n' % (statval[l])
332 out += "--------------------------------------------------\n"
[226]333
[718]334 if rcParams['verbose']:
[181]335 print "--------------------------------------------------"
336 print " ",stat
337 print "--------------------------------------------------"
338 print out
[256]339 retval = {'axes': axes, 'data': array(statval), 'cursor':(i,j,k)}
340 return retval
[102]341
[436]342 def stddev(self,mask=None, allaxes=None):
[135]343 """
344 Determine the standard deviation of the current beam/if/pol
345 Takes a 'mask' as an optional parameter to specify which
346 channels should be excluded.
347 Parameters:
348 mask: an optional mask specifying where the standard
349 deviation should be determined.
[436]350 allaxes: optional flag to show all or a cursor selected
[226]351 spectrum of Beam/IF/Pol. Default is all or taken
352 from .asaprc
[135]353
354 Example:
355 scan.set_unit('channel')
356 msk = scan.create_mask([100,200],[500,600])
357 scan.stddev(mask=m)
358 """
[436]359 if allaxes is None: allaxes = rcParams['scantable.allaxes']
360 return self.stats(stat='stddev',mask=mask, allaxes=allaxes);
[135]361
[436]362 def get_tsys(self, allaxes=None):
[113]363 """
364 Return the System temperatures.
365 Parameters:
[436]366 allaxes: if True apply to all spectra. Otherwise
367 apply only to the selected (beam/pol/if)spectra only.
368 The default is taken from .asaprc (True if none)
[113]369 Returns:
370 a list of Tsys values.
371 """
[436]372 if allaxes is None: allaxes = rcParams['scantable.allaxes']
[256]373 from numarray import array,zeros,Float
374 axes = ['Beam','IF','Pol','Time']
375
[436]376 if allaxes:
[256]377 n = self.nbeam()*self.nif()*self.npol()*self.nrow()
378 shp = [self.nbeam(),self.nif(),self.npol(),self.nrow()]
379 arr = array(zeros(n),shape=shp,type=Float)
380
381 for i in range(self.nbeam()):
382 self.setbeam(i)
383 for j in range(self.nif()):
384 self.setif(j)
385 for k in range(self.npol()):
386 self.setpol(k)
387 arr[i,j,k,:] = self._gettsys()
388 retval = {'axes': axes, 'data': arr, 'cursor':None}
389 tm = [self._gettime(val) for val in range(self.nrow())]
[718]390 if rcParams['verbose']:
[256]391 self._print_values(retval,'Tsys',tm)
392 return retval
393
[102]394 else:
[256]395 i,j,k = (self.getbeam(),self.getif(),self.getpol())
396 statval = self._gettsys()
[181]397 out = ''
398 for l in range(self.nrow()):
399 tm = self._gettime(l)
400 out += 'Time[%s]:\n' % (tm)
401 if self.nbeam() > 1: out += ' Beam[%d] ' % (i)
402 if self.nif() > 1: out += ' IF[%d] ' % (j)
403 if self.npol() > 1: out += ' Pol[%d] ' % (k)
[256]404 out += '= %3.3f\n' % (statval[l])
405 out += "--------------------------------------------------\n"
406
[718]407 if rcParams['verbose']:
[181]408 print "--------------------------------------------------"
[256]409 print " TSys"
[181]410 print "--------------------------------------------------"
411 print out
[256]412 retval = {'axes': axes, 'data': array(statval), 'cursor':(i,j,k)}
413 return retval
[710]414
[407]415 def get_time(self, row=-1):
[113]416 """
417 Get a list of time stamps for the observations.
[181]418 Return a string for each integration in the scantable.
[113]419 Parameters:
[407]420 row: row no of integration. Default -1 return all rows
[113]421 Example:
422 none
423 """
424 out = []
[407]425 if row == -1:
426 for i in range(self.nrow()):
427 out.append(self._gettime(i))
428 return out
429 else:
430 if row < self.nrow():
431 return self._gettime(row)
[102]432
[714]433 def get_sourcename(self, row=-1):
434 """
[787]435 Get a list source names for the observations.
[714]436 Return a string for each integration in the scantable.
437 Parameters:
438 row: row no of integration. Default -1 return all rows
439 Example:
440 none
441 """
442 out = []
443 if row == -1:
444 return [self._getsourcename(i) for i in range(self.nrow())]
445 else:
446 if 0 <= row < self.nrow():
447 return self._getsourcename(row)
448
[787]449 def get_elevation(self, row=-1):
450 """
451 Get a list of elevations for the observations.
452 Return a float for each integration in the scantable.
453 Parameters:
454 row: row no of integration. Default -1 return all rows
455 Example:
456 none
457 """
458 out = []
459 if row == -1:
460 return [self._getelevation(i) for i in range(self.nrow())]
461 else:
462 if 0 <= row < self.nrow():
463 return self._getelevation(row)
464
[788]465 def get_azimuth(self, row=-1):
[787]466 """
467 Get a list of azimuths for the observations.
468 Return a float for each integration in the scantable.
469 Parameters:
470 row: row no of integration. Default -1 return all rows
471 Example:
472 none
473 """
474 out = []
475 if row == -1:
476 return [self._getazimuth(i) for i in range(self.nrow())]
477 else:
478 if 0 <= row < self.nrow():
479 return self._getazimuth(row)
480
481 def get_parangle(self, row=-1):
482 """
483 Get a list of parallactic angles for the observations.
484 Return a float for each integration in the scantable.
485 Parameters:
486 row: row no of integration. Default -1 return all rows
487 Example:
488 none
489 """
490 out = []
491 if row == -1:
492 return [self._getparangle(i) for i in range(self.nrow())]
493 else:
494 if 0 <= row < self.nrow():
495 return self._getparangle(row)
496
[102]497 def set_unit(self, unit='channel'):
498 """
499 Set the unit for all following operations on this scantable
500 Parameters:
501 unit: optional unit, default is 'channel'
[113]502 one of '*Hz','km/s','channel', ''
[102]503 """
[484]504 varlist = vars()
[113]505 if unit in ['','pixel', 'channel']:
506 unit = ''
507 inf = list(self._getcoordinfo())
508 inf[0] = unit
509 self._setcoordinfo(inf)
510 if self._p: self.plot()
[484]511 self._add_history("set_unit",varlist)
[113]512
[484]513 def set_instrument(self, instr):
[358]514 """
515 Set the instrument for subsequent processing
516 Parameters:
[710]517 instr: Select from 'ATPKSMB', 'ATPKSHOH', 'ATMOPRA',
[407]518 'DSS-43' (Tid), 'CEDUNA', and 'HOBART'
[358]519 """
520 self._setInstrument(instr)
[484]521 self._add_history("set_instument",vars())
[718]522 print_log()
[358]523
[276]524 def set_doppler(self, doppler='RADIO'):
525 """
526 Set the doppler for all following operations on this scantable.
527 Parameters:
528 doppler: One of 'RADIO', 'OPTICAL', 'Z', 'BETA', 'GAMMA'
529 """
[484]530 varlist = vars()
[276]531 inf = list(self._getcoordinfo())
532 inf[2] = doppler
533 self._setcoordinfo(inf)
534 if self._p: self.plot()
[484]535 self._add_history("set_doppler",vars())
[718]536 print_log()
[710]537
[226]538 def set_freqframe(self, frame=None):
[113]539 """
540 Set the frame type of the Spectral Axis.
541 Parameters:
[591]542 frame: an optional frame type, default 'LSRK'. Valid frames are:
543 'REST','TOPO','LSRD','LSRK','BARY',
544 'GEO','GALACTO','LGROUP','CMB'
[113]545 Examples:
546 scan.set_freqframe('BARY')
547 """
[484]548 if frame is None: frame = rcParams['scantable.freqframe']
549 varlist = vars()
[113]550 valid = ['REST','TOPO','LSRD','LSRK','BARY', \
551 'GEO','GALACTO','LGROUP','CMB']
[591]552
[113]553 if frame in valid:
554 inf = list(self._getcoordinfo())
555 inf[1] = frame
556 self._setcoordinfo(inf)
[484]557 self._add_history("set_freqframe",varlist)
[102]558 else:
[718]559 msg = "Please specify a valid freq type. Valid types are:\n",valid
560 if rcParams['verbose']:
561 print msg
562 else:
563 raise TypeError(msg)
564 print_log()
[710]565
[113]566 def get_unit(self):
567 """
568 Get the default unit set in this scantable
569 Parameters:
570 Returns:
571 A unit string
572 """
573 inf = self._getcoordinfo()
574 unit = inf[0]
575 if unit == '': unit = 'channel'
576 return unit
[102]577
[158]578 def get_abcissa(self, rowno=0):
[102]579 """
[158]580 Get the abcissa in the current coordinate setup for the currently
[113]581 selected Beam/IF/Pol
582 Parameters:
[226]583 rowno: an optional row number in the scantable. Default is the
584 first row, i.e. rowno=0
[113]585 Returns:
[407]586 The abcissa values and it's format string (as a dictionary)
[113]587 """
[256]588 abc = self._getabcissa(rowno)
[710]589 lbl = self._getabcissalabel(rowno)
[718]590 print_log()
[158]591 return abc, lbl
[113]592
593 def create_mask(self, *args, **kwargs):
594 """
[102]595 Compute and return a mask based on [min,max] windows.
[189]596 The specified windows are to be INCLUDED, when the mask is
[113]597 applied.
[102]598 Parameters:
599 [min,max],[min2,max2],...
600 Pairs of start/end points specifying the regions
601 to be masked
[189]602 invert: optional argument. If specified as True,
603 return an inverted mask, i.e. the regions
604 specified are EXCLUDED
[513]605 row: create the mask using the specified row for
606 unit conversions, default is row=0
607 only necessary if frequency varies over rows.
[102]608 Example:
[113]609 scan.set_unit('channel')
610
611 a)
[102]612 msk = scan.set_mask([400,500],[800,900])
[189]613 # masks everything outside 400 and 500
[113]614 # and 800 and 900 in the unit 'channel'
615
616 b)
617 msk = scan.set_mask([400,500],[800,900], invert=True)
[189]618 # masks the regions between 400 and 500
[113]619 # and 800 and 900 in the unit 'channel'
[710]620
[102]621 """
[513]622 row = 0
623 if kwargs.has_key("row"):
624 row = kwargs.get("row")
625 data = self._getabcissa(row)
[113]626 u = self._getcoordinfo()[0]
[718]627 if rcParams['verbose']:
[113]628 if u == "": u = "channel"
[718]629 from asap import asaplog
630 msg = "The current mask window unit is %s" % u
631 asaplog.push(msg)
[102]632 n = self.nchan()
[189]633 msk = zeros(n)
[710]634 # test if args is a 'list' or a 'normal *args - UGLY!!!
635
636 ws = (isinstance(args[-1][-1],int) or isinstance(args[-1][-1],float)) and args or args[0]
637 for window in ws:
[102]638 if (len(window) != 2 or window[0] > window[1] ):
[718]639 raise TypeError("A window needs to be defined as [min,max]")
[102]640 for i in range(n):
[113]641 if data[i] >= window[0] and data[i] < window[1]:
[189]642 msk[i] = 1
[113]643 if kwargs.has_key('invert'):
644 if kwargs.get('invert'):
645 from numarray import logical_not
[710]646 msk = logical_not(msk)
[718]647 print_log()
[102]648 return msk
[710]649
[256]650 def get_restfreqs(self):
651 """
652 Get the restfrequency(s) stored in this scantable.
653 The return value(s) are always of unit 'Hz'
654 Parameters:
655 none
656 Returns:
657 a list of doubles
658 """
659 return list(self._getrestfreqs())
[102]660
[402]661 def lines(self):
[391]662 """
[402]663 Print the list of known spectral lines
664 """
[714]665 l = sdtable._lines(self)
666 if rcParams['verbose']:
667 print l
668 else:
669 return l
[402]670
[484]671 def set_restfreqs(self, freqs=None, unit='Hz', lines=None, source=None,
672 theif=None):
[402]673 """
[393]674 Select the restfrequency for the specified source and IF OR
675 replace for all IFs. If the 'freqs' argument holds a scalar,
676 then that rest frequency will be applied to the selected
677 data (and added to the list of available rest frequencies).
678 In this way, you can set a rest frequency for each
679 source and IF combination. If the 'freqs' argument holds
680 a vector, then it MUST be of length the number of IFs
681 (and the available restfrequencies will be replaced by
[710]682 this vector). In this case, *all* data ('source' and
683 'theif' are ignored) have the restfrequency set per IF according
684 to the corresponding value you give in the 'freqs' vector.
685 E.g. 'freqs=[1e9,2e9]' would mean IF 0 gets restfreq 1e9 and
[393]686 IF 1 gets restfreq 2e9.
[402]687
688 You can also specify the frequencies via known line names
689 in the argument 'lines'. Use 'freqs' or 'lines'. 'freqs'
690 takes precedence. See the list of known names via function
691 scantable.lines()
[391]692 Parameters:
[402]693 freqs: list of rest frequencies
[393]694 unit: unit for rest frequency (default 'Hz')
[402]695 lines: list of known spectral lines names (alternative to freqs).
696 See possible list via scantable.lines()
[391]697 source: Source name (blank means all)
698 theif: IF (-1 means all)
699 Example:
[393]700 scan.set_restfreqs(freqs=1.4e9, source='NGC253', theif=2)
701 scan.set_restfreqs(freqs=[1.4e9,1.67e9])
[391]702 """
[484]703 varlist = vars()
[391]704 if source is None:
705 source = ""
706 if theif is None:
707 theif = -1
[393]708 t = type(freqs)
709 if t is int or t is float:
710 freqs = [freqs]
[402]711 if freqs is None:
712 freqs = []
713 t = type(lines)
714 if t is str:
715 lines = [lines]
716 if lines is None:
717 lines = []
718 sdtable._setrestfreqs(self, freqs, unit, lines, source, theif)
[484]719 self._add_history("set_restfreqs", varlist)
[391]720
721
[710]722
[800]723 def flag_spectrum(self, flags, row=-1, allaxes=None):
[102]724 """
725 This flags a selected spectrum in the scan 'for good'.
726 USE WITH CARE - not reversible.
727 Use masks for non-permanent exclusion of channels.
728 Parameters:
[800]729 flags: a mask created with scantable.create_mask
730 row: the number of the row (intergration) to flag,
731 the deafult is -1, which is ALL rows.
732 allaxes: if True apply to all spectra. Otherwise
733 apply only to the selected (beam/pol/if)spectra only
734 The default is taken from .asaprc (True if none)
[102]735 Example:
[800]736 scan.flag_spectrum() # flags everything, not very useful
737 msk = scan.create_mask([511,513])
738 scan,set_cursor(IF=1)
739 scan.flag_spectrum(msk, allaxes=False)
740 # flags the central three channels for all rows in the table,
741 # using the selected IF only (e.g. a birdie)
[102]742 """
[800]743 if allaxes is None: allaxes = rcParams['scantable.allaxes']
744 if len(flags) != self.nchan():
745 msg = "Length of mask invalid.",valid
746 if rcParams['verbose']:
747 print msg
748 return
749 else:
750 raise AssertionError(msg)
751 rows = range(self.nrow())
752 if row > -1: rows = [row]
753 for r in rows:
754 if allaxes:
755 beamSel,IFSel,polSel = (self.getbeam(),self.getif(),self.getpol())
756 for i in range(self.nbeam()):
757 self.setbeam(i)
758 for j in range(self.nif()):
759 self.setif(j)
760 for k in range(self.npol()):
761 self.setpol(k)
762 sdtable._flagspectrum(self, flags, r)
763 self.setbeam(beamSel)
764 self.setif(IFSel)
765 self.setpol(polSel)
766 else:
767 sdtable._flagspectrum(self, flags, r)
768 self._add_history("flag_spectrum", vars())
[102]769 return
[113]770
[256]771 def _print_values(self, dat, label='', timestamps=[]):
772 d = dat['data']
773 a = dat['axes']
774 shp = d.getshape()
775 out = ''
776 for i in range(shp[3]):
777 out += '%s [%s]:\n' % (a[3],timestamps[i])
778 t = d[:,:,:,i]
779 for j in range(shp[0]):
780 if shp[0] > 1: out += ' %s[%d] ' % (a[0],j)
781 for k in range(shp[1]):
782 if shp[1] > 1: out += ' %s[%d] ' % (a[1],k)
783 for l in range(shp[2]):
784 if shp[2] > 1: out += ' %s[%d] ' % (a[2],l)
785 out += '= %3.3f\n' % (t[j,k,l])
[484]786 out += "-"*80
787 out += "\n"
788 print "-"*80
[256]789 print " ", label
[484]790 print "-"*80
[710]791 print out
[484]792
793 def history(self):
794 hist = list(self._gethistory())
[791]795 out = "-"*80
[484]796 for h in hist:
[489]797 if h.startswith("---"):
[791]798 out += "\n"+h
[489]799 else:
800 items = h.split("##")
801 date = items[0]
802 func = items[1]
803 items = items[2:]
[791]804 out += "\n"+date+"\n"
805 out += "Function: %s\n Parameters:" % (func)
[489]806 for i in items:
807 s = i.split("=")
[791]808 out += "\n %s = %s" % (s[0],s[1])
809 out += "\n"+"-"*80
810 try:
811 from IPython.genutils import page as pager
812 except ImportError:
813 from pydoc import pager
814 pager(out)
[484]815 return
816
[513]817 #
818 # Maths business
819 #
820
[558]821 def average_time(self, mask=None, scanav=False, weight='tint'):
[513]822 """
823 Return the (time) average of a scan, or apply it 'insitu'.
824 Note:
825 in channels only
826 The cursor of the output scan is set to 0.
827 Parameters:
828 one scan or comma separated scans
829 mask: an optional mask (only used for 'var' and 'tsys'
830 weighting)
[558]831 scanav: True averages each scan separately
832 False (default) averages all scans together,
[524]833 weight: Weighting scheme. 'none', 'var' (1/var(spec)
[535]834 weighted), 'tsys' (1/Tsys**2 weighted), 'tint'
835 (integration time weighted) or 'tintsys' (Tint/Tsys**2).
836 The default is 'tint'
[513]837 Example:
838 # time average the scantable without using a mask
[710]839 newscan = scan.average_time()
[513]840 """
841 varlist = vars()
[524]842 if weight is None: weight = 'tint'
[513]843 if mask is None: mask = ()
[710]844 from asap._asap import average as _av
[513]845 s = scantable(_av((self,), mask, scanav, weight))
846 s._add_history("average_time",varlist)
[718]847 print_log()
[513]848 return s
[710]849
[513]850 def convert_flux(self, jyperk=None, eta=None, d=None, insitu=None,
851 allaxes=None):
852 """
853 Return a scan where all spectra are converted to either
854 Jansky or Kelvin depending upon the flux units of the scan table.
855 By default the function tries to look the values up internally.
856 If it can't find them (or if you want to over-ride), you must
857 specify EITHER jyperk OR eta (and D which it will try to look up
858 also if you don't set it). jyperk takes precedence if you set both.
859 Parameters:
860 jyperk: the Jy / K conversion factor
861 eta: the aperture efficiency
862 d: the geomtric diameter (metres)
863 insitu: if False a new scantable is returned.
864 Otherwise, the scaling is done in-situ
865 The default is taken from .asaprc (False)
[800]866 allaxes: if True apply to all spectra. Otherwise
[513]867 apply only to the selected (beam/pol/if)spectra only
868 The default is taken from .asaprc (True if none)
869 """
870 if allaxes is None: allaxes = rcParams['scantable.allaxes']
871 if insitu is None: insitu = rcParams['insitu']
872 varlist = vars()
873 if jyperk is None: jyperk = -1.0
874 if d is None: d = -1.0
875 if eta is None: eta = -1.0
876 if not insitu:
877 from asap._asap import convertflux as _convert
878 s = scantable(_convert(self, d, eta, jyperk, allaxes))
879 s._add_history("convert_flux", varlist)
[718]880 print_log()
[513]881 return s
882 else:
883 from asap._asap import convertflux_insitu as _convert
884 _convert(self, d, eta, jyperk, allaxes)
885 self._add_history("convert_flux", varlist)
[718]886 print_log()
[513]887 return
888
889 def gain_el(self, poly=None, filename="", method="linear",
890 insitu=None, allaxes=None):
891 """
892 Return a scan after applying a gain-elevation correction.
893 The correction can be made via either a polynomial or a
894 table-based interpolation (and extrapolation if necessary).
895 You specify polynomial coefficients, an ascii table or neither.
896 If you specify neither, then a polynomial correction will be made
897 with built in coefficients known for certain telescopes (an error
898 will occur if the instrument is not known).
899 The data and Tsys are *divided* by the scaling factors.
900 Parameters:
901 poly: Polynomial coefficients (default None) to compute a
902 gain-elevation correction as a function of
903 elevation (in degrees).
904 filename: The name of an ascii file holding correction factors.
905 The first row of the ascii file must give the column
906 names and these MUST include columns
907 "ELEVATION" (degrees) and "FACTOR" (multiply data
908 by this) somewhere.
909 The second row must give the data type of the
910 column. Use 'R' for Real and 'I' for Integer.
911 An example file would be
912 (actual factors are arbitrary) :
913
914 TIME ELEVATION FACTOR
915 R R R
916 0.1 0 0.8
917 0.2 20 0.85
918 0.3 40 0.9
919 0.4 60 0.85
920 0.5 80 0.8
921 0.6 90 0.75
922 method: Interpolation method when correcting from a table.
923 Values are "nearest", "linear" (default), "cubic"
924 and "spline"
925 insitu: if False a new scantable is returned.
926 Otherwise, the scaling is done in-situ
927 The default is taken from .asaprc (False)
928 allaxes: If True apply to all spectra. Otherwise
929 apply only to the selected (beam/pol/if) spectra only
930 The default is taken from .asaprc (True if none)
931 """
932
933 if allaxes is None: allaxes = rcParams['scantable.allaxes']
934 if insitu is None: insitu = rcParams['insitu']
935 varlist = vars()
936 if poly is None:
937 poly = ()
938 from os.path import expandvars
939 filename = expandvars(filename)
940 if not insitu:
941 from asap._asap import gainel as _gainEl
942 s = scantable(_gainEl(self, poly, filename, method, allaxes))
943 s._add_history("gain_el", varlist)
[718]944 print_log()
[513]945 return s
946 else:
947 from asap._asap import gainel_insitu as _gainEl
948 _gainEl(self, poly, filename, method, allaxes)
949 self._add_history("gain_el", varlist)
[718]950 print_log()
[513]951 return
[710]952
[513]953 def freq_align(self, reftime=None, method='cubic', perif=False,
954 insitu=None):
955 """
956 Return a scan where all rows have been aligned in frequency/velocity.
957 The alignment frequency frame (e.g. LSRK) is that set by function
958 set_freqframe.
959 Parameters:
960 reftime: reference time to align at. By default, the time of
961 the first row of data is used.
962 method: Interpolation method for regridding the spectra.
963 Choose from "nearest", "linear", "cubic" (default)
964 and "spline"
965 perif: Generate aligners per freqID (no doppler tracking) or
966 per IF (scan-based doppler tracking)
967 insitu: if False a new scantable is returned.
968 Otherwise, the scaling is done in-situ
969 The default is taken from .asaprc (False)
970 """
971 if insitu is None: insitu = rcParams['insitu']
972 varlist = vars()
973 if reftime is None: reftime = ''
974 perfreqid = not perif
975 if not insitu:
976 from asap._asap import freq_align as _align
977 s = scantable(_align(self, reftime, method, perfreqid))
978 s._add_history("freq_align", varlist)
[718]979 print_log()
[513]980 return s
981 else:
982 from asap._asap import freq_align_insitu as _align
983 _align(self, reftime, method, perfreqid)
984 self._add_history("freq_align", varlist)
[718]985 print_log()
[513]986 return
987
988 def opacity(self, tau, insitu=None, allaxes=None):
989 """
990 Apply an opacity correction. The data
991 and Tsys are multiplied by the correction factor.
992 Parameters:
993 tau: Opacity from which the correction factor is
994 exp(tau*ZD)
995 where ZD is the zenith-distance
996 insitu: if False a new scantable is returned.
997 Otherwise, the scaling is done in-situ
998 The default is taken from .asaprc (False)
999 allaxes: if True apply to all spectra. Otherwise
1000 apply only to the selected (beam/pol/if)spectra only
1001 The default is taken from .asaprc (True if none)
1002 """
1003 if allaxes is None: allaxes = rcParams['scantable.allaxes']
1004 if insitu is None: insitu = rcParams['insitu']
1005 varlist = vars()
1006 if not insitu:
1007 from asap._asap import opacity as _opacity
1008 s = scantable(_opacity(self, tau, allaxes))
1009 s._add_history("opacity", varlist)
[718]1010 print_log()
[513]1011 return s
1012 else:
1013 from asap._asap import opacity_insitu as _opacity
1014 _opacity(self, tau, allaxes)
1015 self._add_history("opacity", varlist)
[718]1016 print_log()
[513]1017 return
1018
1019 def bin(self, width=5, insitu=None):
1020 """
1021 Return a scan where all spectra have been binned up.
1022 width: The bin width (default=5) in pixels
1023 insitu: if False a new scantable is returned.
1024 Otherwise, the scaling is done in-situ
1025 The default is taken from .asaprc (False)
1026 """
1027 if insitu is None: insitu = rcParams['insitu']
1028 varlist = vars()
1029 if not insitu:
1030 from asap._asap import bin as _bin
1031 s = scantable(_bin(self, width))
1032 s._add_history("bin",varlist)
[718]1033 print_log()
[513]1034 return s
1035 else:
1036 from asap._asap import bin_insitu as _bin
1037 _bin(self, width)
1038 self._add_history("bin",varlist)
[718]1039 print_log()
[513]1040 return
1041
[710]1042
[513]1043 def resample(self, width=5, method='cubic', insitu=None):
1044 """
1045 Return a scan where all spectra have been binned up
1046 width: The bin width (default=5) in pixels
1047 method: Interpolation method when correcting from a table.
1048 Values are "nearest", "linear", "cubic" (default)
1049 and "spline"
1050 insitu: if False a new scantable is returned.
1051 Otherwise, the scaling is done in-situ
1052 The default is taken from .asaprc (False)
1053 """
1054 if insitu is None: insitu = rcParams['insitu']
1055 varlist = vars()
1056 if not insitu:
1057 from asap._asap import resample as _resample
1058 s = scantable(_resample(self, method, width))
1059 s._add_history("resample",varlist)
[718]1060 print_log()
[513]1061 return s
1062 else:
1063 from asap._asap import resample_insitu as _resample
1064 _resample(self, method, width)
1065 self._add_history("resample",varlist)
[718]1066 print_log()
[513]1067 return
1068
1069 def average_pol(self, mask=None, weight='none', insitu=None):
1070 """
1071 Average the Polarisations together.
1072 The polarisation cursor of the output scan is set to 0
1073 Parameters:
1074 mask: An optional mask defining the region, where the
1075 averaging will be applied. The output will have all
1076 specified points masked.
[532]1077 weight: Weighting scheme. 'none' (default), 'var' (1/var(spec)
1078 weighted), or 'tsys' (1/Tsys**2 weighted)
[513]1079 insitu: if False a new scantable is returned.
1080 Otherwise, the scaling is done in-situ
1081 The default is taken from .asaprc (False)
1082 """
1083 if insitu is None: insitu = rcParams['insitu']
1084 varlist = vars()
1085
1086 if mask is None:
1087 mask = ()
1088 if not insitu:
1089 from asap._asap import averagepol as _avpol
1090 s = scantable(_avpol(self, mask, weight))
1091 s._add_history("average_pol",varlist)
[718]1092 print_log()
[513]1093 return s
1094 else:
1095 from asap._asap import averagepol_insitu as _avpol
1096 _avpol(self, mask, weight)
1097 self._add_history("average_pol",varlist)
[718]1098 print_log()
[513]1099 return
1100
1101 def smooth(self, kernel="hanning", width=5.0, insitu=None, allaxes=None):
1102 """
1103 Smooth the spectrum by the specified kernel (conserving flux).
1104 Parameters:
1105 scan: The input scan
1106 kernel: The type of smoothing kernel. Select from
1107 'hanning' (default), 'gaussian' and 'boxcar'.
1108 The first three characters are sufficient.
1109 width: The width of the kernel in pixels. For hanning this is
1110 ignored otherwise it defauls to 5 pixels.
1111 For 'gaussian' it is the Full Width Half
1112 Maximum. For 'boxcar' it is the full width.
1113 insitu: if False a new scantable is returned.
1114 Otherwise, the scaling is done in-situ
1115 The default is taken from .asaprc (False)
1116 allaxes: If True (default) apply to all spectra. Otherwise
1117 apply only to the selected (beam/pol/if)spectra only
1118 The default is taken from .asaprc (True if none)
1119 Example:
1120 none
1121 """
1122 if allaxes is None: allaxes = rcParams['scantable.allaxes']
1123 if insitu is None: insitu = rcParams['insitu']
1124 varlist = vars()
1125 if not insitu:
1126 from asap._asap import smooth as _smooth
1127 s = scantable(_smooth(self,kernel,width,allaxes))
1128 s._add_history("smooth", varlist)
[718]1129 print_log()
[513]1130 return s
1131 else:
1132 from asap._asap import smooth_insitu as _smooth
1133 _smooth(self,kernel,width,allaxes)
1134 self._add_history("smooth", varlist)
[718]1135 print_log()
[513]1136 return
1137
[785]1138 def poly_baseline(self, mask=None, order=0, insitu=None, allaxes=None):
[513]1139 """
1140 Return a scan which has been baselined (all rows) by a polynomial.
1141 Parameters:
[785]1142 scan: a scantable
1143 mask: an optional mask
1144 order: the order of the polynomial (default is 0)
1145 insitu: if False a new scantable is returned.
1146 Otherwise, the scaling is done in-situ
1147 The default is taken from .asaprc (False)
1148 allaxes: If True (default) apply to all spectra. Otherwise
1149 apply only to the selected (beam/pol/if)spectra only
1150 The default is taken from .asaprc (True if none)
[513]1151 Example:
1152 # return a scan baselined by a third order polynomial,
1153 # not using a mask
1154 bscan = scan.poly_baseline(order=3)
[579]1155 """
[785]1156 if allaxes is None: allaxes = rcParams['scantable.allaxes']
[513]1157 if insitu is None: insitu = rcParams['insitu']
1158 varlist = vars()
1159 if mask is None:
1160 from numarray import ones
[626]1161 mask = list(ones(self.nchan()))
[513]1162 from asap.asapfitter import fitter
1163 f = fitter()
1164 f.set_scan(self, mask)
1165 f.set_function(poly=order)
[786]1166 sf = f.auto_fit(insitu,allaxes)
[513]1167 if insitu:
1168 self._add_history("poly_baseline", varlist)
[718]1169 print_log()
[513]1170 return
1171 else:
1172 sf._add_history("poly_baseline", varlist)
[718]1173 print_log()
[513]1174 return sf
1175
1176 def auto_poly_baseline(self, mask=None, edge=(0,0), order=0,
1177 threshold=3,insitu=None):
1178 """
1179 Return a scan which has been baselined (all rows) by a polynomial.
1180 Spectral lines are detected first using linefinder and masked out
1181 to avoid them affecting the baseline solution.
1182
1183 Parameters:
1184 mask: an optional mask retreived from scantable
1185 edge: an optional number of channel to drop at
1186 the edge of spectrum. If only one value is
1187 specified, the same number will be dropped from
1188 both sides of the spectrum. Default is to keep
1189 all channels
1190 order: the order of the polynomial (default is 0)
1191 threshold: the threshold used by line finder. It is better to
1192 keep it large as only strong lines affect the
1193 baseline solution.
1194 insitu: if False a new scantable is returned.
1195 Otherwise, the scaling is done in-situ
1196 The default is taken from .asaprc (False)
1197
1198 Example:
[516]1199 scan2=scan.auto_poly_baseline(order=7)
[513]1200 """
1201 if insitu is None: insitu = rcParams['insitu']
1202 varlist = vars()
1203 from asap.asapfitter import fitter
1204 from asap.asaplinefind import linefinder
1205 from asap import _is_sequence_or_number as _is_valid
[710]1206
[513]1207 if not _is_valid(edge, int):
[718]1208
[516]1209 raise RuntimeError, "Parameter 'edge' has to be an integer or a \
1210 pair of integers specified as a tuple"
[710]1211
[513]1212 # setup fitter
1213 f = fitter()
1214 f.set_function(poly=order)
1215
1216 # setup line finder
1217 fl=linefinder()
1218 fl.set_options(threshold=threshold)
1219
1220 if not insitu:
1221 workscan=self.copy()
1222 else:
1223 workscan=self
1224
1225 sel=workscan.get_cursor()
1226 rows=range(workscan.nrow())
[718]1227 from asap import asaplog
[513]1228 for i in range(workscan.nbeam()):
1229 workscan.setbeam(i)
1230 for j in range(workscan.nif()):
1231 workscan.setif(j)
1232 for k in range(workscan.npol()):
1233 workscan.setpol(k)
[718]1234 asaplog.push("Processing:")
1235 msg = 'Beam[%d], IF[%d], Pol[%d]' % (i,j,k)
1236 asaplog.push(msg)
[513]1237 for iRow in rows:
1238 fl.set_scan(workscan,mask,edge)
1239 fl.find_lines(iRow)
1240 f.set_scan(workscan, fl.get_mask())
1241 f.x=workscan._getabcissa(iRow)
1242 f.y=workscan._getspectrum(iRow)
1243 f.data=None
1244 f.fit()
1245 x=f.get_parameters()
1246 workscan._setspectrum(f.fitter.getresidual(),iRow)
1247 workscan.set_cursor(sel[0],sel[1],sel[2])
1248 if not insitu:
[718]1249 return workscan
[513]1250
1251 def rotate_linpolphase(self, angle, allaxes=None):
1252 """
[710]1253 Rotate the phase of the complex polarization O=Q+iU correlation.
1254 This is always done in situ in the raw data. So if you call this
[513]1255 function more than once then each call rotates the phase further.
1256 Parameters:
1257 angle: The angle (degrees) to rotate (add) by.
1258 allaxes: If True apply to all spectra. Otherwise
1259 apply only to the selected (beam/pol/if)spectra only.
1260 The default is taken from .asaprc (True if none)
1261 Examples:
1262 scan.rotate_linpolphase(2.3)
1263 """
1264 if allaxes is None: allaxes = rcParams['scantable.allaxes']
1265 varlist = vars()
1266 from asap._asap import _rotate_linpolphase as _rotate
1267 _rotate(self, angle, allaxes)
1268 self._add_history("rotate_linpolphase", varlist)
[718]1269 print_log()
[513]1270 return
1271
[710]1272
[513]1273 def rotate_xyphase(self, angle, allaxes=None):
1274 """
1275 Rotate the phase of the XY correlation. This is always done in situ
1276 in the data. So if you call this function more than once
1277 then each call rotates the phase further.
1278 Parameters:
1279 angle: The angle (degrees) to rotate (add) by.
1280 allaxes: If True apply to all spectra. Otherwise
1281 apply only to the selected (beam/pol/if)spectra only.
1282 The default is taken from .asaprc (True if none)
1283 Examples:
1284 scan.rotate_xyphase(2.3)
1285 """
1286 if allaxes is None: allaxes = rcParams['scantable.allaxes']
1287 varlist = vars()
[670]1288 from asap._asap import _rotate_xyphase
[513]1289 _rotate_xyphase(self, angle, allaxes)
1290 self._add_history("rotate_xyphase", varlist)
[718]1291 print_log()
[513]1292 return
1293
[798]1294 def invert_phase(self, allaxes=None):
1295 """
1296 Take the complex conjugate of the complex polarisation cross
1297 correlation. This is always done in situ in the data. If
1298 you call this function twice, you end up with the original phase
1299 Parameters:
1300 allaxes: If True apply to all spectra. Otherwise
1301 apply only to the selected (beam/pol/if)spectra only.
1302 The default is taken from .asaprc (True if none)
1303 Examples:
1304 scan.invert_phase()
1305 """
1306 if allaxes is None: allaxes = rcParams['scantable.allaxes']
1307 varlist = vars()
1308 from asap._asap import _invert_phase
1309 _invert_phase(self, allaxes)
1310 self._add_history("invert_phase", varlist)
1311 print_log()
1312 return
[513]1313
[798]1314 def swap_ppol(self, allaxes=None):
1315 """
1316 Swaps the order of parallel polarisation products (ie AA, BB
1317 becomes BB, AA) . This is always done in situ in the data.
1318 If you call this function twice, you end up with the original
1319 order
1320 Parameters:
1321 allaxes: If True apply to all spectra. Otherwise
1322 apply only to the selected (beam/pol/if)spectra only.
1323 The default is taken from .asaprc (True if none)
1324 Examples:
1325 scan.swap_pol()
1326 """
1327 if allaxes is None: allaxes = rcParams['scantable.allaxes']
1328 varlist = vars()
1329 from asap._asap import _swap_pol
1330 _swap_pol(self, allaxes)
1331 self._add_history("swap_phase", varlist)
1332 print_log()
1333 return
1334
1335
[513]1336 def add(self, offset, insitu=None, allaxes=None):
1337 """
1338 Return a scan where all spectra have the offset added
1339 Parameters:
1340 offset: the offset
1341 insitu: if False a new scantable is returned.
1342 Otherwise, the scaling is done in-situ
1343 The default is taken from .asaprc (False)
1344 allaxes: if True apply to all spectra. Otherwise
1345 apply only to the selected (beam/pol/if)spectra only
1346 The default is taken from .asaprc (True if none)
1347 """
1348 if allaxes is None: allaxes = rcParams['scantable.allaxes']
1349 if insitu is None: insitu = rcParams['insitu']
1350 varlist = vars()
1351 if not insitu:
1352 from asap._asap import add as _add
1353 s = scantable(_add(self, offset, allaxes))
1354 s._add_history("add",varlist)
[718]1355 print_log()
[513]1356 return s
1357 else:
1358 from asap._asap import add_insitu as _add
1359 _add(self, offset, allaxes)
1360 self._add_history("add",varlist)
[718]1361 print_log()
[513]1362 return
1363
1364 def scale(self, factor, insitu=None, allaxes=None, tsys=True):
1365 """
1366 Return a scan where all spectra are scaled by the give 'factor'
1367 Parameters:
1368 factor: the scaling factor
1369 insitu: if False a new scantable is returned.
1370 Otherwise, the scaling is done in-situ
1371 The default is taken from .asaprc (False)
1372 allaxes: if True apply to all spectra. Otherwise
1373 apply only to the selected (beam/pol/if)spectra only.
1374 The default is taken from .asaprc (True if none)
1375 tsys: if True (default) then apply the operation to Tsys
1376 as well as the data
1377 """
1378 if allaxes is None: allaxes = rcParams['scantable.allaxes']
1379 if insitu is None: insitu = rcParams['insitu']
1380 varlist = vars()
1381 if not insitu:
1382 from asap._asap import scale as _scale
1383 s = scantable(_scale(self, factor, allaxes, tsys))
1384 s._add_history("scale",varlist)
[718]1385 print_log()
[513]1386 return s
1387 else:
1388 from asap._asap import scale_insitu as _scale
1389 _scale(self, factor, allaxes, tsys)
1390 self._add_history("scale",varlist)
[718]1391 print_log()
[513]1392 return
1393
[670]1394 def auto_quotient(self, mode='suffix', preserve=True):
1395 """
1396 This function allows to build quotients automatically.
1397 It assumes the observation to have the same numer of
1398 "ons" and "offs"
1399 It will support "closest off in time" in the future
1400 Parameters:
1401 mode: the on/off detection mode; 'suffix' (default)
1402 'suffix' identifies 'off' scans by the
1403 trailing '_R' (Mopra/Parkes) or
1404 '_e'/'_w' (Tid)
[710]1405 preserve: you can preserve (default) the continuum or
1406 remove it. The equations used are
[670]1407 preserve: Output = Toff * (on/off) - Toff
1408 remove: Output = Tref * (on/off) - Ton
1409 """
1410 modes = ["suffix","time"]
1411 if not mode in modes:
1412 print "please provide valid mode. Valid modes are %s" % (modes)
1413 return None
1414 from asap._asap import quotient as _quot
1415 if mode == "suffix":
1416 srcs = self.get_scan("*[^_ewR]")
1417 refs = self.get_scan("*[_ewR]")
[710]1418 if isinstance(srcs,scantable) and isinstance(refs,scantable):
[784]1419 from asap import asaplog
[710]1420 ns,nr = srcs.nrow(),refs.nrow()
[784]1421 msg = "Found %i Off and %i On scans" % (ns,nr)
1422 asaplog.push(msg)
[710]1423 if nr > ns:
[784]1424 asaplog("Found more Off integrations than On scans - dropping excess Offs.")
[710]1425 refs = refs.get_scan(range(ns))
[718]1426 print_log()
[710]1427 return scantable(_quot(srcs,refs, preserve))
1428 else:
[718]1429 msg = "Couldn't find any on/off pairs"
1430 if rcParams['verbose']:
1431 print msg
1432 return
1433 else:
1434 raise RuntimeError()
[670]1435 else:
[718]1436 if rcParams['verbose']: print "not yet implemented"
[670]1437 return None
[710]1438
[513]1439 def quotient(self, other, isreference=True, preserve=True):
1440 """
[582]1441 Return the quotient of a 'source' (on) scan and a 'reference' (off)
1442 scan.
[513]1443 The reference can have just one row, even if the signal has many.
1444 Otherwise they must have the same number of rows.
1445 The cursor of the output scan is set to 0
1446 Parameters:
[582]1447 other: the 'other' scan
[513]1448 isreference: if the 'other' scan is the reference (default)
1449 or source
[710]1450 preserve: you can preserve (default) the continuum or
1451 remove it. The equations used are
[582]1452 preserve: Output = Toff * (on/off) - Toff
1453 remove: Output = Tref * (on/off) - Ton
[513]1454 Example:
1455 # src is a scantable for an 'on' scan, ref for an 'off' scantable
1456 q1 = src.quotient(ref)
1457 q2 = ref.quotient(src, isreference=False)
1458 # gives the same result
1459 """
1460 from asap._asap import quotient as _quot
[718]1461 try:
1462 s = None
1463 if isreference:
1464 s = scantable(_quot(self, other, preserve))
1465 else:
1466 s = scantable(_quot(other, self, preserve))
1467 print_log()
1468 return s
1469 except RuntimeError,e:
1470 if rcParams['verbose']:
1471 print e
1472 return
1473 else: raise
[513]1474
[718]1475 def freq_switch(self, insitu=None):
1476 """
1477 Apply frequency switching to the data.
1478 Parameters:
1479 insitu: if False a new scantable is returned.
1480 Otherwise, the swictching is done in-situ
1481 The default is taken from .asaprc (False)
1482 Example:
1483 none
1484 """
1485 if insitu is None: insitu = rcParams['insitu']
1486 varlist = vars()
1487 try:
1488 if insitu:
1489 from asap._asap import _frequency_switch_insitu
1490 _frequency_switch_insitu(self)
1491 self._add_history("freq_switch", varlist)
1492 print_log()
1493 return
1494 else:
1495 from asap._asap import _frequency_switch
1496 sf = scantable(_frequency_switch(self))
1497 sf._add_history("freq_switch", varlist)
1498 print_log()
1499 return sf
1500 except RuntimeError,e:
1501 if rcParams['verbose']: print e
1502 else: raise
1503
[773]1504 def recalc_azel(self):
1505 """
1506 Recalculate the azimuth and elevation for each position.
1507 Parameters:
1508 none
1509 Example:
1510 """
[778]1511 varlist = vars()
[773]1512 self._recalc_azel()
1513 self._add_history("recalc_azel", varlist)
[778]1514 print_log()
[773]1515 return
1516
[513]1517 def __add__(self, other):
1518 varlist = vars()
1519 s = None
1520 if isinstance(other, scantable):
[710]1521 from asap._asap import b_operate as _bop
[513]1522 s = scantable(_bop(self, other, 'add', True))
1523 elif isinstance(other, float):
1524 from asap._asap import add as _add
1525 s = scantable(_add(self, other, True))
1526 else:
[718]1527 raise TypeError("Other input is not a scantable or float value")
[513]1528 s._add_history("operator +", varlist)
[718]1529 print_log()
[513]1530 return s
1531
1532 def __sub__(self, other):
1533 """
1534 implicit on all axes and on Tsys
1535 """
1536 varlist = vars()
1537 s = None
1538 if isinstance(other, scantable):
[710]1539 from asap._asap import b_operate as _bop
[513]1540 s = scantable(_bop(self, other, 'sub', True))
1541 elif isinstance(other, float):
1542 from asap._asap import add as _add
1543 s = scantable(_add(self, -other, True))
1544 else:
[718]1545 raise TypeError("Other input is not a scantable or float value")
[513]1546 s._add_history("operator -", varlist)
[718]1547 print_log()
[513]1548 return s
[710]1549
[513]1550 def __mul__(self, other):
1551 """
1552 implicit on all axes and on Tsys
1553 """
1554 varlist = vars()
1555 s = None
1556 if isinstance(other, scantable):
[710]1557 from asap._asap import b_operate as _bop
[513]1558 s = scantable(_bop(self, other, 'mul', True))
1559 elif isinstance(other, float):
1560 if other == 0.0:
[718]1561 raise ZeroDivisionError("Dividing by zero is not recommended")
[513]1562 from asap._asap import scale as _sca
1563 s = scantable(_sca(self, other, True))
1564 else:
[718]1565 raise TypeError("Other input is not a scantable or float value")
[513]1566 s._add_history("operator *", varlist)
[718]1567 print_log()
[513]1568 return s
1569
[710]1570
[513]1571 def __div__(self, other):
1572 """
1573 implicit on all axes and on Tsys
1574 """
1575 varlist = vars()
1576 s = None
1577 if isinstance(other, scantable):
[710]1578 from asap._asap import b_operate as _bop
[513]1579 s = scantable(_bop(self, other, 'div', True))
1580 elif isinstance(other, float):
1581 if other == 0.0:
[718]1582 raise ZeroDivisionError("Dividing by zero is not recommended")
[513]1583 from asap._asap import scale as _sca
1584 s = scantable(_sca(self, 1.0/other, True))
1585 else:
[718]1586 raise TypeError("Other input is not a scantable or float value")
[513]1587 s._add_history("operator /", varlist)
[718]1588 print_log()
[513]1589 return s
1590
[530]1591 def get_fit(self, row=0):
1592 """
1593 Print or return the stored fits for a row in the scantable
1594 Parameters:
1595 row: the row which the fit has been applied to.
1596 """
1597 if row > self.nrow():
1598 return
1599 from asap import asapfit
1600 fit = asapfit(self._getfit(row))
[718]1601 if rcParams['verbose']:
[530]1602 print fit
1603 return
1604 else:
1605 return fit.as_dict()
1606
[484]1607 def _add_history(self, funcname, parameters):
1608 # create date
1609 sep = "##"
1610 from datetime import datetime
1611 dstr = datetime.now().strftime('%Y/%m/%d %H:%M:%S')
1612 hist = dstr+sep
1613 hist += funcname+sep#cdate+sep
1614 if parameters.has_key('self'): del parameters['self']
1615 for k,v in parameters.iteritems():
1616 if type(v) is dict:
1617 for k2,v2 in v.iteritems():
1618 hist += k2
1619 hist += "="
1620 if isinstance(v2,scantable):
1621 hist += 'scantable'
1622 elif k2 == 'mask':
[513]1623 if isinstance(v2,list) or isinstance(v2,tuple):
1624 hist += str(self._zip_mask(v2))
1625 else:
1626 hist += str(v2)
[484]1627 else:
[513]1628 hist += str(v2)
[484]1629 else:
1630 hist += k
1631 hist += "="
1632 if isinstance(v,scantable):
1633 hist += 'scantable'
1634 elif k == 'mask':
[513]1635 if isinstance(v,list) or isinstance(v,tuple):
1636 hist += str(self._zip_mask(v))
1637 else:
1638 hist += str(v)
[484]1639 else:
1640 hist += str(v)
1641 hist += sep
1642 hist = hist[:-2] # remove trailing '##'
1643 self._addhistory(hist)
1644
[710]1645
[484]1646 def _zip_mask(self, mask):
1647 mask = list(mask)
1648 i = 0
1649 segments = []
1650 while mask[i:].count(1):
1651 i += mask[i:].index(1)
1652 if mask[i:].count(0):
1653 j = i + mask[i:].index(0)
1654 else:
[710]1655 j = len(mask)
[484]1656 segments.append([i,j])
[710]1657 i = j
[484]1658 return segments
[714]1659
[626]1660 def _get_ordinate_label(self):
1661 fu = "("+self.get_fluxunit()+")"
1662 import re
1663 lbl = "Intensity"
1664 if re.match(".K.",fu):
1665 lbl = "Brightness Temperature "+ fu
1666 elif re.match(".Jy.",fu):
1667 lbl = "Flux density "+ fu
1668 return lbl
[710]1669
Note: See TracBrowser for help on using the repository browser.