source: trunk/python/__init__.py@ 1662

Last change on this file since 1662 was 1645, checked in by Max Voronkov, 15 years ago

added a simplified version of the linefinder (largely to be used with some scripts than on real single-dish data)

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 24.6 KB
RevLine 
[100]1"""
2This is the ATNF Single Dish Analysis package.
3
4"""
[1080]5import os,sys,shutil, platform
[1589]6import functools
[226]7
[1080]8# Set up AIPSPATH and first time use of asap i.e. ~/.asap/*
9plf = None
10if sys.platform == "linux2":
11 if platform.architecture()[0] == '64bit':
12 plf = 'linux_64b'
13 else:
14 plf = 'linux_gnu'
15elif sys.platform == 'darwin':
16 plf = 'darwin'
17else:
18 # Shouldn't happen - default to linux
19 plf = 'linux'
20asapdata = __path__[-1]
[1171]21# Allow user defined data location
22if os.environ.has_key("ASAPDATA"):
23 if os.path.exists(os.environ["ASAPDATA"]):
24 asapdata = os.environ["ASAPDATA"]
[1371]25# use AIPSPATH if defined and "data" dir present
26if not os.environ.has_key("AIPSPATH") or \
27 not os.path.exists(os.environ["AIPSPATH"].split()[0]+"/data"):
[1316]28 os.environ["AIPSPATH"] = "%s %s somwhere" % ( asapdata, plf)
[1171]29# set up user space
[1080]30userdir = os.environ["HOME"]+"/.asap"
31if not os.path.exists(userdir):
32 print 'First time ASAP use. Setting up ~/.asap'
33 os.mkdir(userdir)
34 shutil.copyfile(asapdata+"/data/ipythonrc-asap", userdir+"/ipythonrc-asap")
[1588]35 shutil.copyfile(asapdata+"/data/ipy_user_conf.py",
[1434]36 userdir+"/ipy_user_conf.py")
[1080]37 f = file(userdir+"/asapuserfuncs.py", "w")
38 f.close()
39 f = file(userdir+"/ipythonrc", "w")
40 f.close()
[1472]41else:
42 # upgrade to support later ipython versions
43 if not os.path.exists(userdir+"/ipy_user_conf.py"):
[1588]44 shutil.copyfile(asapdata+"/data/ipy_user_conf.py",
[1472]45 userdir+"/ipy_user_conf.py")
46
[1171]47# remove from namespace
[1080]48del asapdata, userdir, shutil, platform
49
[513]50def _validate_bool(b):
[226]51 'Convert b to a boolean or raise'
52 bl = b.lower()
53 if bl in ('f', 'no', 'false', '0', 0): return False
54 elif bl in ('t', 'yes', 'true', '1', 1): return True
55 else:
56 raise ValueError('Could not convert "%s" to boolean' % b)
57
[513]58def _validate_int(s):
[226]59 'convert s to int or raise'
60 try: return int(s)
61 except ValueError:
62 raise ValueError('Could not convert "%s" to int' % s)
63
[513]64def _asap_fname():
[226]65 """
66 Return the path to the rc file
67
68 Search order:
69
70 * current working dir
71 * environ var ASAPRC
[274]72 * HOME/.asaprc
[706]73
[226]74 """
75
76 fname = os.path.join( os.getcwd(), '.asaprc')
77 if os.path.exists(fname): return fname
78
79 if os.environ.has_key('ASAPRC'):
80 path = os.environ['ASAPRC']
81 if os.path.exists(path):
82 fname = os.path.join(path, '.asaprc')
83 if os.path.exists(fname):
84 return fname
85
86 if os.environ.has_key('HOME'):
87 home = os.environ['HOME']
88 fname = os.path.join(home, '.asaprc')
89 if os.path.exists(fname):
90 return fname
91 return None
92
[706]93
[226]94defaultParams = {
95 # general
[513]96 'verbose' : [True, _validate_bool],
97 'useplotter' : [True, _validate_bool],
[542]98 'insitu' : [True, _validate_bool],
[706]99
[226]100 # plotting
[706]101 'plotter.gui' : [True, _validate_bool],
[226]102 'plotter.stacking' : ['p', str],
103 'plotter.panelling' : ['s', str],
[700]104 'plotter.colours' : ['', str],
105 'plotter.linestyles' : ['', str],
[710]106 'plotter.decimate' : [False, _validate_bool],
107 'plotter.ganged' : [True, _validate_bool],
[1022]108 'plotter.histogram' : [False, _validate_bool],
[1097]109 'plotter.papertype' : ['A4', str],
[1560]110 'plotter.axesformatting' : ['mpl', str],
[710]111
[226]112 # scantable
113 'scantable.save' : ['ASAP', str],
[513]114 'scantable.autoaverage' : [True, _validate_bool],
[226]115 'scantable.freqframe' : ['LSRK', str], #default frequency frame
[1076]116 'scantable.verbosesummary' : [False, _validate_bool],
[1422]117 'scantable.storage' : ['memory', str],
[1504]118 'scantable.history' : [True, _validate_bool],
[1586]119 'scantable.reference' : ['.*(e|w|_R)$', str],
120 'scantable.parallactify' : [False, _validate_bool]
[226]121 # fitter
122 }
123
[255]124def list_rcparameters():
[706]125
[255]126 print """
[737]127# general
128# print verbose output
129verbose : True
[255]130
[737]131# preload a default plotter
132useplotter : True
[255]133
[737]134# apply operations on the input scantable or return new one
135insitu : True
[706]136
[737]137# plotting
[710]138
[737]139# do we want a GUI or plot to a file
140plotter.gui : True
[710]141
[737]142# default mode for colour stacking
143plotter.stacking : Pol
[255]144
[737]145# default mode for panelling
146plotter.panelling : scan
[255]147
[1560]148# push panels together, to share axis labels
[737]149plotter.ganged : True
[710]150
[1513]151# decimate the number of points plotted by a factor of
[737]152# nchan/1024
153plotter.decimate : False
[733]154
[737]155# default colours/linestyles
156plotter.colours :
157plotter.linestyles :
[700]158
[1022]159# enable/disable histogram plotting
160plotter.histogram : False
161
[1097]162# ps paper type
163plotter.papertype : A4
164
[1513]165# The formatting style of the xaxis
[1560]166plotter.axesformatting : 'mpl' (default) or 'asap' (for old versions of matplotlib)
[1513]167
[737]168# scantable
[1076]169
[1259]170# default storage of scantable ('memory'/'disk')
[1076]171scantable.storage : memory
[1422]172
173# write history of each call to scantable
174scantable.history : True
175
[737]176# default ouput format when saving
177scantable.save : ASAP
[1422]178
[737]179# auto averaging on read
180scantable.autoaverage : True
[255]181
[737]182# default frequency frame to set when function
[1097]183# scantable.set_freqframe is called
[737]184scantable.freqframe : LSRK
[255]185
[737]186# Control the level of information printed by summary
187scantable.verbosesummary : False
[706]188
[1504]189# Control the identification of reference (off) scans
190# This is has to be a regular expression
[1595]191scantable.reference : .*(e|w|_R)$
[1560]192
[1595]193# Indicate whether the data was parallactified (total phase offest == 0.0)
194scantable.parallactify : False
195
[737]196# Fitter
197"""
[706]198
[226]199def rc_params():
200 'Return the default params updated from the values in the rc file'
[706]201
[513]202 fname = _asap_fname()
[706]203
[226]204 if fname is None or not os.path.exists(fname):
205 message = 'could not find rc file; returning defaults'
206 ret = dict([ (key, tup[0]) for key, tup in defaultParams.items()])
207 #print message
208 return ret
[706]209
[226]210 cnt = 0
211 for line in file(fname):
212 cnt +=1
213 line = line.strip()
214 if not len(line): continue
215 if line.startswith('#'): continue
216 tup = line.split(':',1)
217 if len(tup) !=2:
218 print ('Illegal line #%d\n\t%s\n\tin file "%s"' % (cnt, line, fname))
219 continue
[706]220
[226]221 key, val = tup
222 key = key.strip()
223 if not defaultParams.has_key(key):
224 print ('Bad key "%s" on line %d in %s' % (key, cnt, fname))
225 continue
[706]226
[226]227 default, converter = defaultParams[key]
228
229 ind = val.find('#')
230 if ind>=0: val = val[:ind] # ignore trailing comments
231 val = val.strip()
232 try: cval = converter(val) # try to convert to proper type or raise
[1080]233 except ValueError, msg:
[226]234 print ('Bad val "%s" on line #%d\n\t"%s"\n\tin file "%s"\n\t%s' % (val, cnt, line, fname, msg))
235 continue
236 else:
237 # Alles Klar, update dict
238 defaultParams[key][0] = cval
239
240 # strip the conveter funcs and return
241 ret = dict([ (key, tup[0]) for key, tup in defaultParams.items()])
[466]242 print ('loaded rc file %s'%fname)
[226]243
244 return ret
245
246
247# this is the instance used by the asap classes
[706]248rcParams = rc_params()
[226]249
250rcParamsDefault = dict(rcParams.items()) # a copy
251
252def rc(group, **kwargs):
253 """
254 Set the current rc params. Group is the grouping for the rc, eg
[379]255 for scantable.save the group is 'scantable', for plotter.stacking, the
256 group is 'plotter', and so on. kwargs is a list of attribute
[226]257 name/value pairs, eg
258
[379]259 rc('scantable', save='SDFITS')
[226]260
261 sets the current rc params and is equivalent to
[706]262
[379]263 rcParams['scantable.save'] = 'SDFITS'
[226]264
265 Use rcdefaults to restore the default rc params after changes.
266 """
267
[379]268 aliases = {}
[706]269
[226]270 for k,v in kwargs.items():
271 name = aliases.get(k) or k
[1422]272 if len(group):
273 key = '%s.%s' % (group, name)
274 else:
275 key = name
[226]276 if not rcParams.has_key(key):
277 raise KeyError('Unrecognized key "%s" for group "%s" and name "%s"' % (key, group, name))
[706]278
[226]279 rcParams[key] = v
280
281
282def rcdefaults():
283 """
284 Restore the default rc params - the ones that were created at
285 asap load time
286 """
287 rcParams.update(rcParamsDefault)
288
[1295]289def _n_bools(n, val):
[1588]290 return [ val for i in xrange(n) ]
[513]291
292def _is_sequence_or_number(param, ptype=int):
293 if isinstance(param,tuple) or isinstance(param,list):
[928]294 if len(param) == 0: return True # empty list
[513]295 out = True
296 for p in param:
297 out &= isinstance(p,ptype)
298 return out
299 elif isinstance(param, ptype):
300 return True
301 return False
302
[928]303def _to_list(param, ptype=int):
304 if isinstance(param, ptype):
305 if ptype is str: return param.split()
306 else: return [param]
307 if _is_sequence_or_number(param, ptype):
308 return param
309 return None
[715]310
[944]311def unique(x):
[992]312 """
313 Return the unique values in a list
314 Parameters:
315 x: the list to reduce
316 Examples:
317 x = [1,2,3,3,4]
318 print unique(x)
319 [1,2,3,4]
320 """
[944]321 return dict([ (val, 1) for val in x]).keys()
322
[992]323def list_files(path=".",suffix="rpf"):
324 """
325 Return a list files readable by asap, such as rpf, sdfits, mbf, asap
326 Parameters:
327 path: The directory to list (default '.')
328 suffix: The file extension (default rpf)
329 Example:
330 files = list_files("data/","sdfits")
331 print files
332 ['data/2001-09-01_0332_P363.sdfits',
333 'data/2003-04-04_131152_t0002.sdfits',
334 'data/Sgr_86p262_best_SPC.sdfits']
335 """
336 if not os.path.isdir(path):
337 return None
[1295]338 valid = "rpf rpf.1 rpf.2 sdf sdfits mbf asap".split()
[992]339 if not suffix in valid:
340 return None
341 files = [os.path.expanduser(os.path.expandvars(path+"/"+f)) for f in os.listdir(path)]
342 return filter(lambda x: x.endswith(suffix),files)
343
[715]344# workaround for ipython, which redirects this if banner=0 in ipythonrc
345sys.stdout = sys.__stdout__
346sys.stderr = sys.__stderr__
347
348# Logging
349from asap._asap import Log as _asaplog
350global asaplog
[710]351asaplog=_asaplog()
[715]352if rcParams['verbose']:
353 asaplog.enable()
354else:
355 asaplog.disable()
356
[1589]357
358def print_log_dec(f):
359 @functools.wraps(f)
360 def wrap_it(*args, **kw):
361 val = f(*args, **kw)
362 print_log()
363 return val
364 return wrap_it
365
[715]366def print_log():
[1589]367 log = asaplog.pop().strip()
[715]368 if len(log) and rcParams['verbose']: print log
369 return
370
[1295]371def mask_and(a, b):
372 assert(len(a)==len(b))
373 return [ a[i] & b[i] for i in xrange(len(a)) ]
[1134]374
[1295]375def mask_or(a, b):
376 assert(len(a)==len(b))
377 return [ a[i] | b[i] for i in xrange(len(a)) ]
378
379def mask_not(a):
380 return [ not i for i in a ]
381
[1117]382from asapfitter import fitter
[895]383from asapreader import reader
[944]384from selector import selector
[710]385
[100]386from asapmath import *
[1080]387from scantable import scantable
[1097]388from asaplinefind import linefinder
[1645]389from simplelinefinder import simplelinefinder
[1134]390from linecatalog import linecatalog
[285]391
[928]392if rcParams['useplotter']:
[1295]393 try:
[1423]394 from asapplotter import asapplotter
395 gui = os.environ.has_key('DISPLAY') and rcParams['plotter.gui']
396 if gui:
[1422]397 import matplotlib
398 matplotlib.use("TkAgg")
[1423]399 import pylab
[1422]400 xyplotter = pylab
[1423]401 plotter = asapplotter(gui)
402 del gui
[1295]403 except ImportError:
[1423]404 print "Matplotlib not installed. No plotting available"
[285]405
[574]406__date__ = '$Date: 2009-10-04 12:24:37 +0000 (Sun, 04 Oct 2009) $'.split()[1]
[1560]407__version__ = '$Revision: 1645 $'
[100]408
[1193]409def is_ipython():
[1259]410 return '__IP' in dir(sys.modules["__main__"])
[1589]411
[1193]412if is_ipython():
[1014]413 def version(): print "ASAP %s(%s)"% (__version__, __date__)
[1588]414
415 def list_scans(t = scantable):
[1560]416 import inspect
417 print "The user created scantables are: ",
418 globs=inspect.currentframe().f_back.f_locals.copy()
419 out = [ k for k,v in globs.iteritems() \
420 if isinstance(v, scantable) and not k.startswith("_") ]
421 print out
422 return out
[100]423
[715]424 def commands():
425 x = """
[113]426 [The scan container]
427 scantable - a container for integrations/scans
[182]428 (can open asap/rpfits/sdfits and ms files)
[113]429 copy - returns a copy of a scan
430 get_scan - gets a specific scan out of a scantable
[984]431 (by name or number)
[1093]432 drop_scan - drops a specific scan out of a scantable
433 (by number)
[984]434 set_selection - set a new subselection of the data
435 get_selection - get the current selection object
[113]436 summary - print info about the scantable contents
[182]437 stats - get specified statistic of the spectra in
438 the scantable
439 stddev - get the standard deviation of the spectra
440 in the scantable
[113]441 get_tsys - get the TSys
442 get_time - get the timestamps of the integrations
[1351]443 get_inttime - get the integration time
[733]444 get_sourcename - get the source names of the scans
[794]445 get_azimuth - get the azimuth of the scans
446 get_elevation - get the elevation of the scans
447 get_parangle - get the parallactic angle of the scans
[876]448 get_unit - get the current unit
[513]449 set_unit - set the abcissa unit to be used from this
450 point on
[255]451 get_abcissa - get the abcissa values and name for a given
452 row (time)
[1259]453 get_column_names - get the names of the columns in the scantable
454 for use with selector.set_query
[113]455 set_freqframe - set the frame info for the Spectral Axis
456 (e.g. 'LSRK')
[276]457 set_doppler - set the doppler to be used from this point on
[984]458 set_dirframe - set the frame for the direction on the sky
[240]459 set_instrument - set the instrument name
[1190]460 set_feedtype - set the feed type
[255]461 get_fluxunit - get the brightness flux unit
[240]462 set_fluxunit - set the brightness flux unit
[1508]463 set_sourcetype - set the type of the source - source or reference
[188]464 create_mask - return an mask in the current unit
465 for the given region. The specified regions
466 are NOT masked
[255]467 get_restfreqs - get the current list of rest frequencies
468 set_restfreqs - set a list of rest frequencies
[1471]469 shift_refpix - shift the reference pixel of the IFs
470 set_spectrum - overwrite the spectrum for a given row
471 get_spectrum - retrieve the spectrum for a given
472 get_mask - retrieve the mask for a given
[1012]473 flag - flag selected channels in the data
[1192]474 lag_flag - flag specified frequency in the data
[1151]475 save - save the scantable to disk as either 'ASAP',
476 'SDFITS' or 'ASCII'
[486]477 nbeam,nif,nchan,npol - the number of beams/IFs/Pols/Chans
[733]478 nscan - the number of scans in the scantable
[1422]479 nrow - the number of spectra in the scantable
[486]480 history - print the history of the scantable
[530]481 get_fit - get a fit which has been stored witnh the data
[706]482 average_time - return the (weighted) time average of a scan
[513]483 or a list of scans
484 average_pol - average the polarisations together.
[1144]485 average_beam - average the beams together.
[992]486 convert_pol - convert to a different polarisation type
[690]487 auto_quotient - return the on/off quotient with
[1069]488 automatic detection of the on/off scans (closest
489 in time off is selected)
[1144]490 mx_quotient - Form a quotient using MX data (off beams)
[1014]491 scale, *, / - return a scan scaled by a given factor
[1560]492 add, + - return a scan with given value added
493 sub, - - return a scan with given value subtracted
[513]494 bin - return a scan with binned channels
495 resample - return a scan with resampled channels
496 smooth - return the spectrally smoothed scan
497 poly_baseline - fit a polynomial baseline to all Beams/IFs/Pols
[706]498 auto_poly_baseline - automatically fit a polynomial baseline
[780]499 recalc_azel - recalculate azimuth and elevation based on
500 the pointing
[513]501 gain_el - apply gain-elevation correction
502 opacity - apply opacity correction
503 convert_flux - convert to and from Jy and Kelvin brightness
[255]504 units
[513]505 freq_align - align spectra in frequency frame
[1014]506 invert_phase - Invert the phase of the cross-correlation
[1351]507 swap_linears - Swap XX and YY (or RR LL)
[513]508 rotate_xyphase - rotate XY phase of cross correlation
509 rotate_linpolphase - rotate the phase of the complex
510 polarization O=Q+iU correlation
[733]511 freq_switch - perform frequency switching on the data
512 stats - Determine the specified statistic, e.g. 'min'
513 'max', 'rms' etc.
514 stddev - Determine the standard deviation of the current
515 beam/if/pol
[1588]516 get_row_selector - get the selection object for a specified row
517 number
[1014]518 [Selection]
519 selector - a selection object to set a subset of a scantable
520 set_scans - set (a list of) scans by index
521 set_cycles - set (a list of) cycles by index
522 set_beams - set (a list of) beamss by index
523 set_ifs - set (a list of) ifs by index
524 set_polarisations - set (a list of) polarisations by name
525 or by index
526 set_names - set a selection by name (wildcards allowed)
527 set_tsys - set a selection by tsys thresholds
[1259]528 set_query - set a selection by SQL-like query, e.g. BEAMNO==1
[1351]529 ( also get_ functions for all these )
[1014]530 reset - unset all selections
[1351]531 + - merge two selections
[733]532
[513]533 [Math] Mainly functions which operate on more than one scantable
[100]534
[706]535 average_time - return the (weighted) time average
[513]536 of a list of scans
537 quotient - return the on/off quotient
538 simple_math - simple mathematical operations on two scantables,
539 'add', 'sub', 'mul', 'div'
[1069]540 quotient - build quotient of the given on and off scans
[1351]541 (matched pairs and 1 off - n on are valid)
[1117]542 merge - merge a list of scantables
[1069]543
[1151]544 [Line Catalog]
545 linecatalog - a linecatalog wrapper, taking an ASCII or
546 internal format table
547 summary - print a summary of the current selection
548 set_name - select a subset by name pattern, e.g. '*OH*'
549 set_strength_limits - select a subset by line strength limits
550 set_frequency_limits - select a subset by frequency limits
551 reset - unset all selections
552 save - save the current subset to a table (internal
553 format)
554 get_row - get the name and frequency from a specific
555 row in the table
[513]556 [Fitting]
[113]557 fitter
558 auto_fit - return a scan where the function is
559 applied to all Beams/IFs/Pols.
560 commit - return a new scan where the fits have been
561 commited.
562 fit - execute the actual fitting process
[984]563 store_fit - store the fit parameters in the data (scantable)
[113]564 get_chi2 - get the Chi^2
565 set_scan - set the scantable to be fit
566 set_function - set the fitting function
567 set_parameters - set the parameters for the function(s), and
568 set if they should be held fixed during fitting
[513]569 set_gauss_parameters - same as above but specialised for individual
570 gaussian components
[113]571 get_parameters - get the fitted parameters
[513]572 plot - plot the resulting fit and/or components and
573 residual
[210]574 [Plotter]
575 asapplotter - a plotter for asap, default plotter is
576 called 'plotter'
[984]577 plot - plot a scantable
[1151]578 plot_lines - plot a linecatalog overlay
[378]579 save - save the plot to a file ('png' ,'ps' or 'eps')
[210]580 set_mode - set the state of the plotter, i.e.
581 what is to be plotted 'colour stacked'
582 and what 'panelled'
[984]583 set_selection - only plot a selected part of the data
[733]584 set_range - set a 'zoom' window [xmin,xmax,ymin,ymax]
[255]585 set_legend - specify user labels for the legend indeces
586 set_title - specify user labels for the panel indeces
[733]587 set_abcissa - specify a user label for the abcissa
[255]588 set_ordinate - specify a user label for the ordinate
[378]589 set_layout - specify the multi-panel layout (rows,cols)
[733]590 set_colors - specify a set of colours to use
591 set_linestyles - specify a set of linestyles to use if only
592 using one color
[1151]593 set_font - set general font properties, e.g. 'family'
[1054]594 set_histogram - plot in historam style
[733]595 set_mask - set a plotting mask for a specific polarization
[1171]596 text - draw text annotations either in data or relative
597 coordinates
598 arrow - draw arrow annotations either in data or relative
599 coordinates
600 axhline,axvline - draw horizontal/vertical lines
601 axhspan,axvspan - draw horizontal/vertical regions
[1560]602 annotate - draw an arrow with label
603 create_mask - create a scnatble mask interactively
[1175]604
605 xyplotter - matplotlib/pylab plotting functions
606
[182]607 [Reading files]
608 reader - access rpfits/sdfits files
[984]609 open - attach reader to a file
610 close - detach reader from file
[182]611 read - read in integrations
612 summary - list info about all integrations
613
[113]614 [General]
615 commands - this command
616 print - print details about a variable
[1560]617 list_scans - list all scantables created by the user
[992]618 list_files - list all files readable by asap (default rpf)
[113]619 del - delete the given variable from memory
620 range - create a list of values, e.g.
621 range(3) = [0,1,2], range(2,5) = [2,3,4]
622 help - print help for one of the listed functions
623 execfile - execute an asap script, e.g. execfile('myscript')
[255]624 list_rcparameters - print out a list of possible values to be
[274]625 put into $HOME/.asaprc
[1171]626 rc - set rc parameters from within asap
[466]627 mask_and,mask_or,
628 mask_not - boolean operations on masks created with
629 scantable.create_mask
[706]630
[210]631 Note:
632 How to use this with help:
633 # function 'summary'
634 [xxx] is just a category
635 Every 'sub-level' in this list should be replaces by a '.' Period when
[706]636 using help
[210]637 Example:
638 ASAP> help scantable # to get info on ths scantable
639 ASAP> help scantable.summary # to get help on the scantable's
640 ASAP> help average_time
641
[715]642 """
[1151]643 if rcParams['verbose']:
644 try:
645 from IPython.genutils import page as pager
646 except ImportError:
647 from pydoc import pager
648 pager(x)
649 else:
650 print x
[715]651 return
[113]652
[706]653def welcome():
654 return """Welcome to ASAP v%s (%s) - the ATNF Spectral Analysis Package
[100]655
[1035]656Please report any bugs via:
[1378]657http://svn.atnf.csiro.au/trac/asap/simpleticket
[100]658
[378]659[IMPORTANT: ASAP is 0-based]
[706]660Type commands() to get a list of all available ASAP commands.""" % (__version__, __date__)
Note: See TracBrowser for help on using the repository browser.