source: trunk/python/scantable.py@ 987

Last change on this file since 987 was 981, checked in by mar637, 19 years ago

Removed align option from average in c++ as it is buggy.
python is now applying it.

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