source: trunk/python/scantable.py@ 1002

Last change on this file since 1002 was 1001, checked in by mar637, 18 years ago

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