source: trunk/python/__init__.py @ 1725

Last change on this file since 1725 was 1725, checked in by Malte Marquarding, 14 years ago

Finishing touches to opacity calculations, docs, plotting and model

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