Changeset 1826 for trunk/python
- Timestamp:
- 08/03/10 11:41:13 (14 years ago)
- Location:
- trunk/python
- Files:
-
- 15 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/python/asapfit.py
r1232 r1826 1 1 from _asap import fitentry 2 from asap import rcParams2 from asap.parameters import rcParams 3 3 4 4 class asapfit(fitentry): -
trunk/python/asapfitter.py
r1819 r1826 1 1 import _asap 2 from asap import rcParams 3 from asap import print_log, print_log_dec 4 from asap import _n_bools 5 from asap import mask_and 6 from asap import asaplog 2 from asap.parameters import rcParams 3 from asap.logging import asaplog, print_log, print_log_dec 4 from asap.utils import _n_bools, mask_and 5 7 6 8 7 class fitter: … … 697 696 scan = self.data 698 697 rows = xrange(scan.nrow()) 699 # Save parameters of baseline fits as a class attribute. 698 # Save parameters of baseline fits as a class attribute. 700 699 # NOTICE: This does not reflect changes in scantable! 701 700 if len(rows) > 0: self.blpars=[] -
trunk/python/asaplinefind.py
r1644 r1826 71 71 in the noise box (default) 72 72 'median' median of deviations in the noise box 73 73 74 74 Note: For bad baselines threshold should be increased, 75 75 and avg_limit decreased (or even switched off completely by -
trunk/python/asaplot.py
r1819 r1826 13 13 """ 14 14 def __init__(self, rows=1, cols=0, title='', size=None, buffering=False): 15 16 15 """ 16 Create a new instance of the ASAPlot plotting class. 17 17 18 19 20 18 If rows < 1 then a separate call to set_panels() is required to define 19 the panel layout; refer to the doctext for set_panels(). 20 """ 21 21 v = vars() 22 22 del v['self'] -
trunk/python/asaplotbase.py
r1819 r1826 5 5 import sys 6 6 from re import match 7 8 7 import matplotlib 9 8 … … 12 11 from numpy import sqrt 13 12 from matplotlib import rc, rcParams 14 from asap import rcParams as asaprcParams15 13 from matplotlib.ticker import OldScalarFormatter 14 15 from asap.parameters import rcParams as asaprcParams 16 from asap.logging import asaplog 16 17 17 18 # API change in mpl >= 0.98 … … 21 22 from matplotlib.transforms import blend_xy_sep_transform as blended_transform_factory 22 23 23 from asap import asaplog 24 25 if int(matplotlib.__version__.split(".")[1]) < 87: 24 if int(matplotlib.__version__.split(".")[1]) < 99: 26 25 #print "Warning: matplotlib version < 0.87. This might cause errors. Please upgrade." 27 asaplog.push( "matplotlib version < 0. 87. This might cause errors. Please upgrade." )26 asaplog.push( "matplotlib version < 0.99. This might cause errors. Please upgrade." ) 28 27 print_log( 'WARN' ) 29 28 … … 316 315 self.figmgr.toolbar.draw_rubberband(event, event.x, event.y, 317 316 self.rect['x'], self.rect['y']) 318 317 319 318 def region_disable(event): 320 319 self.register('motion_notify', None) … … 655 654 cols, i+1) 656 655 if asaprcParams['plotter.axesformatting'] != 'mpl': 657 656 658 657 self.subplots[i]['axes'].xaxis.set_major_formatter(OldScalarFormatter()) 659 658 else: -
trunk/python/asaplotgui.py
r1819 r1826 17 17 18 18 def __init__(self, rows=1, cols=0, title='', size=None, buffering=False): 19 20 19 """ 20 Create a new instance of the ASAPlot plotting class. 21 21 22 23 24 22 If rows < 1 then a separate call to set_panels() is required to define 23 the panel layout; refer to the doctext for set_panels(). 24 """ 25 25 v = vars() 26 26 del v['self'] … … 49 49 50 50 def map(self): 51 52 53 54 55 56 51 """ 52 Reveal the ASAPlot graphics window and bring it to the top of the 53 window stack. 54 """ 55 self.window.wm_deiconify() 56 self.window.lift() 57 57 58 58 def quit(self): 59 60 61 62 59 """ 60 Destroy the ASAPlot graphics window. 61 """ 62 self.window.destroy() 63 63 64 64 def show(self, hardrefresh=True): 65 66 67 68 65 """ 66 Show graphics dependent on the current buffering state. 67 """ 68 if not self.buffering: 69 69 if hardrefresh: 70 70 asaplotbase.show(self) 71 72 71 self.window.wm_deiconify() 72 self.canvas.show() 73 73 74 74 def terminate(self): 75 76 77 78 75 """ 76 Clear the figure. 77 """ 78 self.window.destroy() 79 79 80 80 def unmap(self): 81 82 83 84 81 """ 82 Hide the ASAPlot graphics window. 83 """ 84 self.window.wm_withdraw() -
trunk/python/asapmath.py
r1819 r1826 1 1 from asap.scantable import scantable 2 from asap import rcParams 3 from asap import print_log, print_log_dec 4 from asap import selector 5 from asap import asaplog 2 from asap.paramaters import rcParams 3 from asap.logging import asaplog, print_log, print_log_dec 4 from asap.selector import selector 6 5 from asap import asaplotgui 7 6 … … 937 936 """ 938 937 Calibrate data. 939 938 940 939 Parameters: 941 940 scantab: scantable 942 941 scannos: list of scan number 943 942 calmode: calibration mode 944 verify: verify calibration 943 verify: verify calibration 945 944 """ 946 945 antname = scantab.get_antennaname() … … 979 978 scal = scantab.copy() 980 979 981 return scal 980 return scal 982 981 983 982 def apexcal( scantab, scannos=[], calmode='none', verify=False ): … … 990 989 calmode: calibration mode 991 990 992 verify: verify calibration 991 verify: verify calibration 993 992 """ 994 993 from asap._asap import stmath … … 1008 1007 calmode: calibration mode 1009 1008 1010 verify: verify calibration 1009 verify: verify calibration 1011 1010 """ 1012 1011 from asap._asap import stmath … … 1020 1019 Split Measurement set by antenna name, save data as a scantables, 1021 1020 and return a list of filename. 1022 Notice this method can only be available from CASA. 1021 Notice this method can only be available from CASA. 1023 1022 Prameter 1024 filename: the name of Measurement set to be read. 1023 filename: the name of Measurement set to be read. 1025 1024 outprefix: the prefix of output scantable name. 1026 1025 the names of output scantable will be … … 1030 1029 The default False is to return with warning 1031 1030 without writing the output. USE WITH CARE. 1032 1031 1033 1032 """ 1034 1033 # Import the table toolkit from CASA … … 1071 1070 if not os.path.isdir(filename) \ 1072 1071 or not os.path.exists(filename+'/ANTENNA') \ 1073 or not os.path.exists(filename+'/table.f1'): 1072 or not os.path.exists(filename+'/table.f1'): 1074 1073 s = "File '%s' is not a Measurement set." % (filename) 1075 1074 if rcParams['verbose']: … … 1119 1118 value: float list operand 1120 1119 mode: operation mode (ADD, SUB, MUL, DIV) 1121 tsys: if True, operate tsys as well 1120 tsys: if True, operate tsys as well 1122 1121 """ 1123 1122 nrow = scan.nrow() … … 1134 1133 from asap._asap import stmath 1135 1134 stm = stmath() 1136 # insitu must be True 1135 # insitu must be True 1137 1136 stm._setinsitu( True ) 1138 1137 s = scan.copy() … … 1150 1149 del stm 1151 1150 return s 1152 1153 1154 -
trunk/python/asapreader.py
r1819 r1826 1 1 from asap._asap import stfiller 2 from asap import print_log, print_log_dec2 from asap.logging import print_log, print_log_dec 3 3 4 4 class reader(stfiller): -
trunk/python/casatoolbar.py
r1819 r1826 10 10 11 11 ### select the nearest spectrum in pick radius 12 ### and display spectral value on the toolbar. 12 ### and display spectral value on the toolbar. 13 13 def _select_spectrum(self,event): 14 14 # Do not fire event when in zooming/panning mode 15 15 mode = self.figmgr.toolbar.mode 16 if not mode == '':17 18 # When selected point is out of panels16 if not mode == '': 17 return 18 # When selected point is out of panels 19 19 if event.inaxes == None: 20 20 return 21 21 # If not left button 22 22 if event.button != 1: 23 23 return 24 24 25 25 xclick=event.xdata … … 27 27 dist2=1000. 28 28 pickline=None 29 30 29 # If the pannel has picable objects 30 pflag=False 31 31 for lin in event.inaxes.lines: 32 if not lin.pickable(): continue 33 pflag=True 34 flag,pind = lin.contains(event) 35 if not flag: continue 36 # Get nearest point 37 inds = pind['ind'] 38 xlin = lin.get_xdata() 39 ylin = lin.get_ydata() 40 for i in inds: 41 d2=(xlin[i]-xclick)**2+(ylin[i]-yclick)**2 42 if dist2 >= d2: 43 dist2 = d2 44 pickline = lin 45 # No pickcable line in the pannel 46 if not pflag: return 47 # Pickable but too far from mouse position 48 elif pickline is None: 49 picked='No line selected.' 50 self.figmgr.toolbar.set_message(picked) 51 return 32 if not lin.pickable(): 33 continue 34 pflag=True 35 flag,pind = lin.contains(event) 36 if not flag: 37 continue 38 # Get nearest point 39 inds = pind['ind'] 40 xlin = lin.get_xdata() 41 ylin = lin.get_ydata() 42 for i in inds: 43 d2=(xlin[i]-xclick)**2+(ylin[i]-yclick)**2 44 if dist2 >= d2: 45 dist2 = d2 46 pickline = lin 47 # No pickcable line in the pannel 48 if not pflag: 49 return 50 # Pickable but too far from mouse position 51 elif pickline is None: 52 picked='No line selected.' 53 self.figmgr.toolbar.set_message(picked) 54 return 52 55 del pind, inds, xlin, ylin 53 56 # Spectra are Picked 54 57 theplot = self.plotter._plotter 55 56 57 58 59 58 thetoolbar = self.figmgr.toolbar 59 thecanvas = self.figmgr.canvas 60 # Disconnect the default motion notify event 61 # Notice! the other buttons are also diabled!!! 62 thecanvas.mpl_disconnect(thetoolbar._idDrag) 60 63 # Get picked spectrum 61 64 xdata = pickline.get_xdata() … … 64 67 titp=event.inaxes.title.get_text() 65 68 panel0=event.inaxes 66 67 68 69 picked="Selected: '"+titl+"' in panel '"+titp+"'." 70 thetoolbar.set_message(picked) 71 # Generate a navigation window 69 72 #naviwin=Navigationwindow(titp,titl) 70 73 #------------------------------------------------------# 71 74 # Show spectrum data at mouse position 72 75 def spec_data(event): 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 76 # Getting spectrum data of neiboring point 77 xclick=event.xdata 78 if event.inaxes != panel0: 79 return 80 ipoint=len(xdata)-1 81 for i in range(len(xdata)-1): 82 xl=xclick-xdata[i] 83 xr=xclick-xdata[i+1] 84 if xl*xr <= 0.: 85 ipoint = i 86 break 87 # Output spectral value on the navigation window 88 posi='[ %s, %s ]: x = %.2f value = %.2f'\ 89 %(titl,titp,xdata[ipoint],ydata[ipoint]) 90 #naviwin.posi.set(posi) 91 thetoolbar.set_message(posi) 89 92 #------------------------------------------------------# 90 93 # Disconnect from mouse events 91 94 def discon(event): 92 #naviwin.window.destroy() 93 theplot.register('motion_notify',None) 94 # Re-activate the default motion_notify_event 95 thetoolbar._idDrag=thecanvas.mpl_connect('motion_notify_event', thetoolbar.mouse_move) 96 theplot.register('button_release',None) 97 return 95 #naviwin.window.destroy() 96 theplot.register('motion_notify',None) 97 # Re-activate the default motion_notify_event 98 thetoolbar._idDrag=thecanvas.mpl_connect('motion_notify_event', 99 thetoolbar.mouse_move) 100 theplot.register('button_release',None) 101 return 98 102 #------------------------------------------------------# 99 103 # Show data value along with mouse movement 100 104 theplot.register('motion_notify',spec_data) 101 105 # Finish events when mouse button is released 102 106 theplot.register('button_release',discon) 103 107 104 108 105 ### Calculate statistics of the selected area. 109 ### Calculate statistics of the selected area. 106 110 def _single_mask(self,event): 107 111 # Do not fire event when in zooming/panning mode 108 if not self.figmgr.toolbar.mode == '': return 112 if not self.figmgr.toolbar.mode == '': 113 return 109 114 # When selected point is out of panels 110 115 if event.inaxes == None: 111 return 112 if event.button ==1: baseinv=True 113 elif event.button == 3: baseinv=False 114 else: return 115 116 def _calc_stats(): 117 msk=mymask.get_mask() 118 mymask.scan.stats(stat='max',mask=msk) 119 mymask.scan.stats(stat='min',mask=msk) 120 mymask.scan.stats(stat='sum',mask=msk) 121 mymask.scan.stats(stat='mean',mask=msk) 122 mymask.scan.stats(stat='median',mask=msk) 123 mymask.scan.stats(stat='rms',mask=msk) 124 mymask.scan.stats(stat='stddev',mask=msk) 125 126 # Interactive mask definition 116 return 117 if event.button ==1: 118 baseinv=True 119 elif event.button == 3: 120 baseinv=False 121 else: 122 return 123 124 def _calc_stats(): 125 msk=mymask.get_mask() 126 mymask.scan.stats(stat='max',mask=msk) 127 mymask.scan.stats(stat='min',mask=msk) 128 mymask.scan.stats(stat='sum',mask=msk) 129 mymask.scan.stats(stat='mean',mask=msk) 130 mymask.scan.stats(stat='median',mask=msk) 131 mymask.scan.stats(stat='rms',mask=msk) 132 mymask.scan.stats(stat='stddev',mask=msk) 133 134 # Interactive mask definition 127 135 from asap.interactivemask import interactivemask 128 129 130 131 132 133 134 135 136 136 mymask=interactivemask(plotter=self.plotter,scan=self.plotter._data) 137 # Create initial mask 138 mymask.set_basemask(invert=baseinv) 139 # Inherit event 140 mymask.set_startevent(event) 141 # Set callback func 142 mymask.set_callback(_calc_stats) 143 # Selected mask 144 mymask.select_mask(once=True,showmask=False) 137 145 138 146 ##################################### … … 140 148 ##################################### 141 149 ### TkAgg 142 if matplotlib.get_backend() == 'TkAgg': import Tkinter as Tk 150 if matplotlib.get_backend() == 'TkAgg': 151 import Tkinter as Tk 152 143 153 class CustomToolbarTkAgg(CustomToolbarCommon, Tk.Frame): 144 154 def __init__(self,parent): 145 155 from asap.asapplotter import asapplotter 146 if not isinstance(parent,asapplotter): return False 147 if not parent._plotter: return False 156 if not isinstance(parent,asapplotter): 157 return False 158 if not parent._plotter: 159 return False 148 160 self._p=parent._plotter 149 161 self.figmgr=self._p.figmgr … … 172 184 173 185 def _NewButton(self, master, text, command, side=Tk.LEFT): 174 if (os.uname()[0] == 'Darwin'):186 if os.uname()[0] == 'Darwin': 175 187 b = Tk.Button(master=master, text=text, command=command) 176 188 else: 177 b = Tk.Button(master=master, text=text, padx=2, pady=2, command=command) 189 b = Tk.Button(master=master, text=text, padx=2, pady=2, 190 command=command) 178 191 b.pack(side=side) 179 192 return b 180 193 181 194 def spec_show(self): 182 195 if not self.figmgr.toolbar.mode == '' or not self.button: return … … 212 225 self.button=True 213 226 self.spec_show() 214 227 215 228 def disable_button(self): 216 229 if not self.button: return 217 self.bStat.config(relief='raised', state=Tk.DISABLED)218 self.bSpec.config(relief='raised', state=Tk.DISABLED)230 self.bStat.config(relief='raised', state=Tk.DISABLED) 231 self.bSpec.config(relief='raised', state=Tk.DISABLED) 219 232 self.button=False 220 233 self.mode='' -
trunk/python/interactivemask.py
r1819 r1826 1 from asap import rcParams2 from asap import _n_bools, mask_and, mask_or1 from asap.parameters import rcParams 2 from asap.utils import _n_bools, mask_and, mask_or 3 3 from asap.scantable import scantable 4 4 … … 15 15 my_mask.finish_selection(callback=func) 16 16 mask=my_mask.get_mask() 17 18 Modify mask region by selecting a region on a plot with mouse. 17 18 Modify mask region by selecting a region on a plot with mouse. 19 19 """ 20 20 … … 22 22 """ 23 23 Create a interactive masking object. 24 Either or both 'plotter' or/and 'scan' should be defined. 24 Either or both 'plotter' or/and 'scan' should be defined. 25 25 26 26 Parameters: … … 63 63 self.ydataold=None 64 64 self._polygons=[] 65 65 66 66 67 67 def set_basemask(self,masklist=[],invert=False): 68 68 """ 69 69 Set initial channel mask. 70 70 71 71 Parameters: 72 72 masklist: [[min, max], [min2, max2], ...] 73 A list of pairs of start/end points (inclusive) 73 A list of pairs of start/end points (inclusive) 74 74 specifying the regions to be masked 75 75 invert: optional argument. If specified as True, … … 77 77 specified are excluded 78 78 You can reset the mask selection by running this method with 79 the default parameters. 79 the default parameters. 80 80 """ 81 81 # Verify input parameters … … 97 97 """ 98 98 Inherit an event from the parent function. 99 99 100 100 Parameters: 101 101 event: 'button_press_event' object to be inherited to 102 start interactive region selection . 102 start interactive region selection . 103 103 """ 104 104 from matplotlib.backend_bases import MouseEvent 105 if isinstance(event,MouseEvent) and event.name=='button_press_event': 105 if isinstance(event,MouseEvent) and event.name=='button_press_event': 106 106 self.event=event 107 107 else: 108 108 msg="Invalid event." 109 raise TypeError(msg) 109 raise TypeError(msg) 110 110 111 111 def set_callback(self,callback): 112 112 """ 113 Set callback function to run when finish_selection() is executed. 113 Set callback function to run when finish_selection() is executed. 114 114 callback: The post processing function to run after 115 115 the mask selections are completed. … … 123 123 Do interactive mask selection. 124 124 Modify masks interactively by adding/deleting regions with 125 mouse drawing.(left-button: mask; right-button: UNmask) 126 Note that the interactive region selection is available only 127 when GUI plotter is active. 125 mouse drawing.(left-button: mask; right-button: UNmask) 126 Note that the interactive region selection is available only 127 when GUI plotter is active. 128 128 129 129 Parameters: 130 130 once: If specified as True, you can modify masks only 131 once. Else if False, you can modify them repeatedly. 131 once. Else if False, you can modify them repeatedly. 132 132 showmask: If specified as True, the masked regions are plotted 133 133 on the plotter. 134 134 Note this parameter is valid only when once=True. 135 Otherwise, maskes are forced to be plotted for reference. 135 Otherwise, maskes are forced to be plotted for reference. 136 136 """ 137 137 # Return if GUI is not active … … 181 181 if self.event != None: 182 182 self._region_start(self.event) 183 else: 183 else: 184 184 self.p._plotter.register('button_press',None) 185 185 self.p._plotter.register('button_press',self._region_start) … … 205 205 def _region_draw(self,event): 206 206 sameaxes=(event.inaxes == self.rect['axes']) 207 if sameaxes: 207 if sameaxes: 208 208 xnow=event.x 209 209 ynow=event.y … … 224 224 # Delete the rubber band 225 225 self.p._plotter.figmgr.toolbar.release(event) 226 227 if event.inaxes == self.rect['axes']: 226 227 if event.inaxes == self.rect['axes']: 228 228 xend=event.x 229 229 yend=event.y … … 235 235 xdataend=self.xdataold 236 236 ydataend=self.ydataold 237 237 238 238 self.rect['world'][2:4] = [xdataend, ydataend] 239 239 self.rect['pixel'][2:4] = [xend, yend] … … 300 300 if len(ifs) > 1: 301 301 for k in xrange(len(ifs)-1): 302 self._polygons.append(self.p._plotter.subplots[j]['axes'].axvspan(projs[k][i][0],projs[k][i][1],facecolor='#aaaaaa')) 302 self._polygons.append(self.p._plotter.subplots[j]['axes'].axvspan(projs[k][i][0],projs[k][i][1],facecolor='#aaaaaa')) 303 303 self._polygons.append(self.p._plotter.subplots[j]['axes'].axvspan(msks[i][0],msks[i][1],facecolor='yellow')) 304 304 self.p._plotter.canvas.draw() … … 314 314 Specifying the callback function here will overwrite 315 315 the one set by set_callback(func) 316 316 317 317 Note this function is automatically called at the end of 318 select_mask() if once=True. 318 select_mask() if once=True. 319 319 """ 320 320 if callback: self.callback=callback … … 324 324 if not self.newplot: 325 325 self.clear_polygon() 326 else: 326 else: 327 327 self.p._plotter.unmap() 328 328 self.p._plotter = None … … 334 334 def clear_polygon(self): 335 335 """ 336 Erase masks plots from the plotter. 336 Erase masks plots from the plotter. 337 337 """ 338 338 if len(self._polygons)>0: … … 347 347 Get the interactively selected channel mask. 348 348 Returns: 349 A list of channel mask. 349 A list of channel mask. 350 350 """ 351 351 return self.mask 352 353 -
trunk/python/lagflagger.py
r1823 r1826 1 from asap import scantable 1 __all_ = ['lagplotter'] 2 3 from asap.scantable import scantable 2 4 from matplotlib.pylab import * 3 5 from numpy import array, ma, logical_not -
trunk/python/linecatalog.py
r1819 r1826 7 7 __revision__ = "$Revision$" 8 8 from asap._asap import linecatalog as lcbase 9 from asap import rcParams10 from asap import asaplog9 from asap.parameters import rcParams 10 from asap.logging import asaplog 11 11 import os 12 12 … … 130 130 if k < 0: k = self.nrow()-k 131 131 return self.get_row(k) 132 133 134 135 136 -
trunk/python/opacity.py
r1823 r1826 1 __all__ = ["model", "skydip"] 1 2 import os 2 3 import math 3 from asap import scantable4 from asap import merge5 from asap import fitter6 from asap import selector7 from asap import rcParams4 from asap.scantable import scantable 5 from asap.asapmath import merge 6 from asap.asapfitter import fitter 7 from asap.selector import selector 8 from asap.parameters import rcParams 8 9 from asap._asap import atmosphere 9 10 … … 45 46 elevation: observatory elevation about sea level (in meters) 46 47 """ 47 self._atm = atmosphere(temperature, self._to_pascals(pressure), 48 self._atm = atmosphere(temperature, self._to_pascals(pressure), 48 49 humidity) 49 50 self.set_observatory_elevation(elevation) -
trunk/python/selector.py
r1819 r1826 1 1 from asap._asap import selector as _selector 2 from asap import unique, _to_list2 from asap.utils import unique, _to_list 3 3 4 4 class selector(_selector): … … 165 165 Set a sequence of row numbers (0-based). Power users Only! 166 166 NOTICE row numbers can be changed easily by sorting, 167 prior selection, etc. 167 prior selection, etc. 168 168 Parameters: 169 169 rows: a list of integers. Default [] is to unset the selection. 170 170 """ 171 vec = _to_list(rows, int) 171 vec = _to_list(rows, int) 172 172 if isinstance(vec,list): 173 173 self._setrows(vec) … … 177 177 def set_types(self, types=[]): 178 178 """ 179 Set a sequence of source types. 179 Set a sequence of source types. 180 180 Parameters: 181 181 types: a list of integers. Default [] is to unset the selection. 182 182 """ 183 vec = _to_list(types, int) 183 vec = _to_list(types, int) 184 184 if isinstance(vec,list): 185 185 self._settypes(vec) -
trunk/python/simplelinefinder.py
r1823 r1826 1 1 import math 2 from asap import asaplog2 from asap.logging import asaplog 3 3 4 4 class StatCalculator: … … 7 7 self.s2=0. 8 8 self.cnt=0 9 9 10 10 def mean(self): 11 11 if self.cnt<=0: 12 12 raise RuntimeError, "At least one data point has to be defined" 13 13 return self.s/float(self.cnt) 14 14 15 15 def variance(self): 16 16 if self.cnt<=1: 17 17 raise RuntimeError, "At least two data points has to be defined" 18 18 return math.sqrt((self.s2/self.cnt)-(self.s/self.cnt)**2+1e-12) 19 19 20 20 def rms(self): 21 21 """ … … 25 25 raise RuntimeError, "At least one data point has to be defined" 26 26 return math.sqrt(self.s2/self.cnt) 27 27 28 28 def add(self, pt): 29 29 self.s = self.s + pt … … 36 36 A simplified class to search for spectral features. The algorithm assumes that the bandpass 37 37 is taken out perfectly and no spectral channels are flagged (except some edge channels). 38 It works with a list or tuple rather than a scantable and returns the channel pairs. 38 It works with a list or tuple rather than a scantable and returns the channel pairs. 39 39 There is an optional feature to attempt to split the detected lines into components, although 40 it should be used with caution. This class is largely intended to be used with scripts. 41 40 it should be used with caution. This class is largely intended to be used with scripts. 41 42 42 The fully featured version of the algorithm working with scantables is called linefinder. 43 43 ''' 44 44 45 45 def __init__(self): 46 46 ''' 47 Initialize the class. 47 Initialize the class. 48 48 ''' 49 49 self._median = None 50 50 self._rms = None 51 51 52 52 def writeLog(self, str): 53 53 """ 54 54 Write user defined string into log file 55 """ 55 """ 56 56 asaplog.push(str) 57 57 58 58 def invertChannelSelection(self, nchan, chans, edge = (0,0)): 59 59 """ 60 This method converts a tuple with channel ranges to a tuple which covers all channels 60 This method converts a tuple with channel ranges to a tuple which covers all channels 61 61 not selected by the original tuple (optionally edge channels can be discarded) 62 63 nchan - number of channels in the spectrum. 62 63 nchan - number of channels in the spectrum. 64 64 chans - tuple (with even number of elements) containing start and stop channel for all selected ranges 65 65 edge - one number or two element tuple (separate values for two ends) defining how many channels to reject … … 69 69 """ 70 70 if nchan<=1: 71 raise RuntimeError, "Number of channels is supposed to be at least 2, you have %i"% nchan 71 raise RuntimeError, "Number of channels is supposed to be at least 2, you have %i"% nchan 72 72 if len(chans)%2!=0: 73 73 raise RuntimeError, "chans is supposed to be a tuple with even number of elements" … … 106 106 spc - tuple with the spectrum (first is the vector of abcissa values, second is the spectrum itself) 107 107 vel_range - a 2-element tuple with start and stop velocity of the range 108 108 109 109 return: a 2-element tuple with channels 110 110 Note, if supplied range is completely outside the spectrum, an empty tuple will be returned … … 145 145 """ 146 146 A helper method to compare x and y by absolute value (to do sorting) 147 147 148 148 x - first value 149 149 y - second value 150 150 151 151 return -1,0 or 1 depending on the result of comparison 152 152 """ … … 157 157 else: 158 158 return 0 159 159 160 160 def rms(self): 161 161 """ … … 167 167 raise RuntimeError, "call find_lines before using the rms method" 168 168 return self._rms 169 169 170 170 def median(self): 171 171 """ … … 176 176 raise RuntimeError, "call find_lines before using the median method" 177 177 return self._median 178 178 179 179 def _mergeIntervals(self, lines, spc): 180 180 """ 181 181 A helper method to merge intervals. 182 182 183 183 lines - list of tuples with first and last channels of all intervals 184 184 spc - spectrum (to be able to test whether adjacent intervals have the … … 196 196 raise RuntimeError, "this shouldn't have happened!" 197 197 lines.pop(toberemoved[i]) 198 198 199 199 def _splitIntervals(self,lines,spc,threshold=3,minchan=3): 200 200 """ … … 203 203 Noise is dealt with by taking into account only those extrema, where a difference with 204 204 respect to surrounding spectral points exceeds threshold times rms (stored inside this 205 class, so the main body of the line detection should be executed first) and there are 206 at least minchan such points. 205 class, so the main body of the line detection should be executed first) and there are 206 at least minchan such points. 207 207 """ 208 208 if minchan<1: … … 219 219 if wasIncreasing != None: 220 220 if isIncreasing != wasIncreasing: 221 derivSignReversals.append((ch,isIncreasing)) 221 derivSignReversals.append((ch,isIncreasing)) 222 222 wasIncreasing = isIncreasing 223 223 if len(derivSignReversals)==0: … … 237 237 newlines.append((startchan,line[1])) 238 238 return newlines 239 239 240 240 def find_lines(self,spc,threshold=3,edge=0,minchan=3, tailsearch = True, splitFeatures = False): 241 241 """ … … 243 243 taken out perfectly and no channels are flagged within the spectrum. A detection 244 244 is reported if consequtive minchan number of channels is consistently above or 245 below the median value. The threshold is given in terms of the rms calculated 246 using 80% of the lowest data points by the absolute value (with respect to median) 247 245 below the median value. The threshold is given in terms of the rms calculated 246 using 80% of the lowest data points by the absolute value (with respect to median) 247 248 248 spc - a list or tuple with the spectrum, no default 249 249 threshold - detection threshold, default is 3 sigma, see above for the definition … … 254 254 tailsearch - if True (default), the algorithm attempts to widen each line until 255 255 its flux crosses the median. It merges lines if necessary. Set this 256 option off if you need to split the lines according to some criterion 256 option off if you need to split the lines according to some criterion 257 257 splitFeatures - if True, the algorithm attempts to split each detected spectral feature into 258 258 a number of spectral lines (just one local extremum). The default action is 259 259 not to do it (may give an adverse results if the noise is high) 260 260 261 261 This method returns a list of tuples each containing start and stop 0-based channel 262 262 number of every detected line. Empty list if nothing has been detected. 263 263 264 264 Note. The median and rms about this median is stored inside this class and can 265 265 be obtained with rms and median methods. … … 270 270 raise RuntimeError, "edge is too high (%i), you rejected all channels (%i)" % (edge, len(spc)) 271 271 if threshold<=0: 272 raise RuntimeError, "threshold parameter of find_lines should be positive, you have %s" % threshold 272 raise RuntimeError, "threshold parameter of find_lines should be positive, you have %s" % threshold 273 273 if minchan<=0: 274 274 raise RuntimeError, "minchan parameter of find_lines should be positive, you have %s" % minchan 275 276 # temporary storage to get statistics, apply edge rejection here 275 276 # temporary storage to get statistics, apply edge rejection here 277 277 tmpspc = spc[edge:len(spc)-edge+1] 278 278 if len(tmpspc)<2: 279 raise RuntimeError, "Too many channels are rejected. Decrease edge parameter or provide a longer spectrum." 279 raise RuntimeError, "Too many channels are rejected. Decrease edge parameter or provide a longer spectrum." 280 280 tmpspc.sort() 281 281 self._median=tmpspc[len(tmpspc)/2] … … 286 286 sc = StatCalculator() 287 287 for i in tmpspc[:-int(0.2*len(tmpspc))]: 288 sc.add(i) 288 sc.add(i) 289 289 self._rms=sc.rms() 290 291 self.writeLog("Spectral line detection with edge=%i, threshold=%f, minchan=%i and tailsearch=%s" % (edge,threshold, minchan, tailsearch)) 292 self.writeLog("statistics: median=%f, rms=%f" % (self._median, self._rms)) 293 290 291 self.writeLog("Spectral line detection with edge=%i, threshold=%f, minchan=%i and tailsearch=%s" % (edge,threshold, minchan, tailsearch)) 292 self.writeLog("statistics: median=%f, rms=%f" % (self._median, self._rms)) 293 294 294 #actual line detection 295 295 lines=[] … … 319 319 if nchan>=minchan: 320 320 lines.append((startchan,len(spc)-edge-1)) 321 321 322 322 if tailsearch: 323 323 for i in range(len(lines)): … … 325 325 curRange = list(lines[i]) 326 326 for x in range(curRange[0],edge,-1): 327 isAbove=(spc[x] > self._median) 327 isAbove=(spc[x] > self._median) 328 328 if wasAbove == None: 329 329 wasAbove = isAbove … … 332 332 break 333 333 for x in range(curRange[1],len(spc)-edge): 334 isAbove=(spc[x] > self._median) 334 isAbove=(spc[x] > self._median) 335 335 if isAbove!=wasAbove: 336 336 curRange[1]=x-1 337 337 break 338 338 lines[i]=tuple(curRange) 339 self._mergeIntervals(lines,spc) 339 self._mergeIntervals(lines,spc) 340 340 if splitFeatures: 341 341 return self._splitIntervals(lines,spc,threshold,minchan) 342 return lines 342 return lines
Note:
See TracChangeset
for help on using the changeset viewer.