source: trunk/python/scantable.py@ 915

Last change on this file since 915 was 914, checked in by mar637, 19 years ago

Stokes polarimetry handling and conversion
bug fix: forgot self on _check_ifs

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