source: trunk/python/scantable.py@ 908

Last change on this file since 908 was 907, checked in by vor010, 19 years ago

LineFinder & auto_poly_baseline: a support of
new scantable format has been added. Now findLines accept a mask and edge
parameters, which can therefore be different for different rows. Constructor
need now just a scan table. Boost types removed from STLineFinder. auto_poly_baseline can accept a nested tuple of edges, which would be interpreted as
different edge parameters for different IFs

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 52.0 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 _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 RuntimeError, "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 RuntimeError, "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# from asap._asap import _rotate_linpolphase as _rotate
1130# _rotate(self, angle, allaxes)
1131# self._add_history("rotate_linpolphase", varlist)
1132# print_log()
1133# return
1134#
1135#
1136# def rotate_xyphase(self, angle):
1137# """
1138# Rotate the phase of the XY correlation. This is always done in situ
1139# in the data. So if you call this function more than once
1140# then each call rotates the phase further.
1141# Parameters:
1142# angle: The angle (degrees) to rotate (add) by.
1143# Examples:
1144# scan.rotate_xyphase(2.3)
1145# """
1146# varlist = vars()
1147# from asap._asap import _rotate_xyphase
1148# _rotate_xyphase(self, angle, allaxes)
1149# self._add_history("rotate_xyphase", varlist)
1150# print_log()
1151# return
1152
1153
1154 def add(self, offset, insitu=None):
1155 """
1156 Return a scan where all spectra have the offset added
1157 Parameters:
1158 offset: the offset
1159 insitu: if False a new scantable is returned.
1160 Otherwise, the scaling is done in-situ
1161 The default is taken from .asaprc (False)
1162 """
1163 if insitu is None: insitu = rcParams['insitu']
1164 self._math._setinsitu(insitu)
1165 varlist = vars()
1166 s = scantable(self._math._unaryop(self, offset, "ADD", False))
1167 s._add_history("add",varlist)
1168 print_log()
1169 if insitu:
1170 self._assign(s)
1171 else:
1172 return s
1173
1174 def scale(self, factor, tsys=True, insitu=None,):
1175 """
1176 Return a scan where all spectra are scaled by the give 'factor'
1177 Parameters:
1178 factor: the scaling factor
1179 insitu: if False a new scantable is returned.
1180 Otherwise, the scaling is done in-situ
1181 The default is taken from .asaprc (False)
1182 tsys: if True (default) then apply the operation to Tsys
1183 as well as the data
1184 """
1185 if insitu is None: insitu = rcParams['insitu']
1186 self._math._setinsitu(insitu)
1187 varlist = vars()
1188 s = scantable(self._math._unaryop(self, factor, "MUL", tsys))
1189 s._add_history("scale",varlist)
1190 print_log()
1191 if insitu:
1192 self._assign(s)
1193 else:
1194 return s
1195
1196 def auto_quotient(self, mode='time', preserve=True):
1197 """
1198 This function allows to build quotients automatically.
1199 It assumes the observation to have the same numer of
1200 "ons" and "offs"
1201 It will support "closest off in time" in the future
1202 Parameters:
1203 mode: the on/off detection mode; 'suffix' (default)
1204 'suffix' identifies 'off' scans by the
1205 trailing '_R' (Mopra/Parkes) or
1206 '_e'/'_w' (Tid)
1207 preserve: you can preserve (default) the continuum or
1208 remove it. The equations used are
1209 preserve: Output = Toff * (on/off) - Toff
1210 remove: Output = Tref * (on/off) - Ton
1211 """
1212 modes = ["time"]
1213 if not mode in modes:
1214 msg = "please provide valid mode. Valid modes are %s" % (modes)
1215 raise ValueError(msg)
1216 varlist = vars()
1217 s = scantable(self._math._quotient(self, mode, preserve))
1218 s._add_history("auto_quotient",varlist)
1219 print_log()
1220 return s
1221
1222
1223
1224
1225 def freq_switch(self, insitu=None):
1226 """
1227 Apply frequency switching to the data.
1228 Parameters:
1229 insitu: if False a new scantable is returned.
1230 Otherwise, the swictching is done in-situ
1231 The default is taken from .asaprc (False)
1232 Example:
1233 none
1234 """
1235 if insitu is None: insitu = rcParams['insitu']
1236 self._math._setinsitu(insitu)
1237 varlist = vars()
1238 s = scantable(self._math._freqswitch(self))
1239 s._add_history("freq_switch",varlist)
1240 print_log()
1241 if insitu: self._assign(s)
1242 else: return s
1243
1244 def recalc_azel(self):
1245 """
1246 Recalculate the azimuth and elevation for each position.
1247 Parameters:
1248 none
1249 Example:
1250 """
1251 varlist = vars()
1252 self._recalcazel()
1253 self._add_history("recalc_azel", varlist)
1254 print_log()
1255 return
1256
1257 def __add__(self, other):
1258 varlist = vars()
1259 s = None
1260 if isinstance(other, scantable):
1261 print "scantable + scantable NYI"
1262 return
1263 elif isinstance(other, float):
1264 s = scantable(self._math._unaryop(self, other, "ADD", False))
1265 else:
1266 raise TypeError("Other input is not a scantable or float value")
1267 s._add_history("operator +", varlist)
1268 print_log()
1269 return s
1270
1271 def __sub__(self, other):
1272 """
1273 implicit on all axes and on Tsys
1274 """
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, "SUB", 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 __mul__(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, "MUL", 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
1306 def __div__(self, other):
1307 """
1308 implicit on all axes and on Tsys
1309 """
1310 varlist = vars()
1311 s = None
1312 if isinstance(other, scantable):
1313 print "scantable / scantable NYI"
1314 return
1315 elif isinstance(other, float):
1316 if other == 0.0:
1317 raise ZeroDivisionError("Dividing by zero is not recommended")
1318 s = scantable(self._math._unaryop(self, other, "DIV", False))
1319 else:
1320 raise TypeError("Other input is not a scantable or float value")
1321 s._add_history("operator /", varlist)
1322 print_log()
1323 return s
1324
1325 def get_fit(self, row=0):
1326 """
1327 Print or return the stored fits for a row in the scantable
1328 Parameters:
1329 row: the row which the fit has been applied to.
1330 """
1331 if row > self.nrow():
1332 return
1333 from asap import asapfit
1334 fit = asapfit(self._getfit(row))
1335 if rcParams['verbose']:
1336 print fit
1337 return
1338 else:
1339 return fit.as_dict()
1340
1341 def _add_history(self, funcname, parameters):
1342 # create date
1343 sep = "##"
1344 from datetime import datetime
1345 dstr = datetime.now().strftime('%Y/%m/%d %H:%M:%S')
1346 hist = dstr+sep
1347 hist += funcname+sep#cdate+sep
1348 if parameters.has_key('self'): del parameters['self']
1349 for k,v in parameters.iteritems():
1350 if type(v) is dict:
1351 for k2,v2 in v.iteritems():
1352 hist += k2
1353 hist += "="
1354 if isinstance(v2,scantable):
1355 hist += 'scantable'
1356 elif k2 == 'mask':
1357 if isinstance(v2,list) or isinstance(v2,tuple):
1358 hist += str(self._zip_mask(v2))
1359 else:
1360 hist += str(v2)
1361 else:
1362 hist += str(v2)
1363 else:
1364 hist += k
1365 hist += "="
1366 if isinstance(v,scantable):
1367 hist += 'scantable'
1368 elif k == 'mask':
1369 if isinstance(v,list) or isinstance(v,tuple):
1370 hist += str(self._zip_mask(v))
1371 else:
1372 hist += str(v)
1373 else:
1374 hist += str(v)
1375 hist += sep
1376 hist = hist[:-2] # remove trailing '##'
1377 self._addhistory(hist)
1378
1379
1380 def _zip_mask(self, mask):
1381 mask = list(mask)
1382 i = 0
1383 segments = []
1384 while mask[i:].count(1):
1385 i += mask[i:].index(1)
1386 if mask[i:].count(0):
1387 j = i + mask[i:].index(0)
1388 else:
1389 j = len(mask)
1390 segments.append([i,j])
1391 i = j
1392 return segments
1393
1394 def _get_ordinate_label(self):
1395 fu = "("+self.get_fluxunit()+")"
1396 import re
1397 lbl = "Intensity"
1398 if re.match(".K.",fu):
1399 lbl = "Brightness Temperature "+ fu
1400 elif re.match(".Jy.",fu):
1401 lbl = "Flux density "+ fu
1402 return lbl
1403
1404 def _check_ifs(self):
1405 nchans = [self.nchan(i) for i in range(self.nif(-1))]
1406 nchans = filter(lambda t: t > 0, nchans)
1407 return (sum(nchans)/len(nchans) == nchans[0])
Note: See TracBrowser for help on using the repository browser.