Changeset 113


Ignore:
Timestamp:
12/01/04 11:29:39 (20 years ago)
Author:
mar637
Message:

version 0.1a

Location:
trunk/python
Files:
1 added
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/python/__init__.py

    r100 r113  
    55#import _asap
    66#from asaplot import ASAPlot
    7 #from asapfitter import AsapFitter
     7from asapfitter import *
    88from asapreader import reader
    99from asapmath import *
    … …  
    2222    #    t = meta_t
    2323    globs = sys.modules['__main__'].__dict__.iteritems()
    24     return map(lambda x: x[0], filter(lambda x: isinstance(x[1], t), globs))
     24    print "The user created scantables are:"
     25    x = map(lambda x: x[0], filter(lambda x: isinstance(x[1], t), globs))
     26    print x
    2527
    26 def list_all():
    27      return ['reader','scantable','ASAPlot','AsapFitter']
     28def commands():
     29    x = """
     30    [Reading files]
     31        reader              - access rpfits/sdfits files
     32            read            - read in integrations
     33            summary         - list info about all integrations
     34    [The scan container]
     35        scantable           - a container for integrations/scans
     36            copy            - returns a copy of a scan
     37            get_scan        - gets a specific scan out of a scantable
     38            summary         - print info about the scantable contents
     39            set_selection   - set a specific Beam/IF/Pol for furthrt use
     40            get_selection   - print out the current selection
     41            rms             - get the rms of the spectra in the scantable
     42            get_tsys        - get the TSys
     43            get_time        - get the timestamps of the integrations
     44            set_unit        - set the units to be used from this point on
     45            set_freqframe   - set the frame info for the Spectral Axis
     46                              (e.g. 'LSRK')
     47            create_mask     - return a mask in thecurrent unit
     48            set_restfreqs   - give a list of rest frequencies
     49            flag_spectrum   - flag a whole Beam/IF/Pol
     50            save            - save the scantable to disk
     51            nbeam,nif,nchan,npol - the number of beams/IFs/Pols/Chans
     52    [Math]
     53        average_scans       - return the rms-weighted (time) average of
     54                              a scan or a list of scans
     55        average_pol         - averge the polarisations together.
     56                              The dimension won't be reduced and
     57                              all polarisations will contain the
     58                              averaged spectrum.
     59        quotient            - return the on/off quotient
     60        scale               - returns a scan scaled by a given factor
     61        bin                 - return a scan with binned channels
     62        hanning             - return the hanning smoothed scan
     63        poly_baseline       - fit a polynomial baseline to all Beams/IFs/Pols
    2864
    29 print """
    30 Welcome to ASAP - the ATNF Single Dish Analysis Package
     65        fitter
     66            auto_fit        - return a scan where the function is
     67                              applied to all Beams/IFs/Pols.
     68            commit          - return a new scan where the fits have been
     69                              commited.
     70            fit             - execute the actual fitting process
     71            get_chi2        - get the Chi^2
     72            set_scan        - set the scantable to be fit
     73            set_function    - set the fitting function
     74            set_parameters  - set the parameters for the function(s), and
     75                              set if they should be held fixed during fitting
     76            get_parameters  - get the fitted parameters
     77           
     78    [General]
     79        commands            - this command
     80        print               - print details about a variable
     81        list_scans          - list all scantables created bt the user
     82        del                 - delete the given variable from memory
     83        range               - create a list of values, e.g.
     84                              range(3) = [0,1,2], range(2,5) = [2,3,4]
     85        help                - print help for one of the listed functions
     86        execfile            - execute an asap script, e.g. execfile('myscript')
     87    """
     88    print x
     89    return
     90
     91print """Welcome to ASAP - the ATNF Single Dish Analysis Package
    3192This is a testing pre-release v0.1a
    3293
    3394Please report any bugs to:
    34                  Malte.Marquarding@atnf.csiro.au
     95Malte.Marquarding@atnf.csiro.au
    3596
    36 NOTE: ASAP is 0-based
     97[NOTE: ASAP is 0-based]
     98Type commands() to get a list of all available ASAP commands.
    3799"""

  • trunk/python/asapmath.py

    r101 r113  
    11from scantable import scantable
    2 def average_scan(scan):
     2
     3def average_time(*args, **kwargs):
    34    """
    4         Return a (time) averaged a scan, i.e. all correlator cycles
    5         are averaged into one "scan".
     5    Return the (time) average of a scan or list of scans. [in channels only]
     6    Parameters:
     7        one scan or comma separated  scans
     8        mask:     an optional mask
     9    Example:
     10        # return a time averaged scan from scana and scanb
     11        # without using a mask
     12        scanav = average_scans(scana,scanb)
     13        # return the (time) averaged scan, i.e. the average of
     14        # all correlator cycles
     15        scanav = average_time(scan)
     16       
    617    """
    7     from asap._asap import average as _av
    8     return scantable(_av(scan))
    9 
    10 
    11 def average_scans(*args, **kwargs):
    12     """
    13         Return the (time) average of a list of scans. [in channels only]
    14         Parameters:
    15             a comma separated list of scans
    16             mask:     an optional mask
    17         Example:
    18             scanav = average_scans(scana,scanb)
    19             # return a time averaged scan from scan and scanb
    20             # without using a mask
    21     """
    22 
     18    lst = args
    2319    if len(args) < 2:
    24         print "Please give at least two scantables"
    25         return
    26 
    27     from asap._asap import averages as _av
    28     d = [args[0].nbeam(),args[0].nif(),args[0].npol(),args[0].nchan()]
    29     for s in args:
     20        if type(args[0]) is list:
     21            if  len(args[0]) < 2:
     22                print "Please give at least two scantables"
     23                return
     24        else:
     25            s = args[0]
     26            if s.nrow() > 1:
     27                from asap._asap import average as _av
     28                return scantable(_av(s))
     29            else:
     30                print "Given scantable is already time averaged"
     31                return
     32        lst = tuple(args[0])
     33    else:
     34        lst = tuple(args)
     35    from asap._asap import averages as _avs
     36    d = [lst[0].nbeam(),lst[0].nif(),lst[0].npol(),lst[0].nchan()]
     37    for s in lst:
    3038        if not isinstance(s,scantable):
    3139            print "Please give a list of scantables"
    … …  
    3644            return
    3745    if kwargs.has_key('mask'):
    38         return scantable(_av(args, kwargs.get('mask')))
     46        return scantable(_avs(lst, kwargs.get('mask')))
    3947    else:
    4048        from numarray import ones
    41         mask = list(ones(d[3]))
    42         return scantable(_av((args), mask))
     49        mask = tuple(ones(d[3]))
     50        return scantable(_avs(lst, mask))
    4351
    4452def quotient(source, reference):
    … …  
    5765    Parameters:
    5866        scan:        a scantable
    59         factor:      the sclaing factor
     67        factor:      the scaling factor
    6068    Note:
    6169        This currently applies the all beams/IFs/pols
    … …  
    6371    from asap._asap import scale as _scale
    6472    return scantable(_scale(scan, factor))
     73
     74def add(scan, offset):
     75    """
     76    Return a scan where the offset is added.
     77    Parameters:
     78        scan:        a scantable
     79        offset:      the value to add
     80    Note:
     81        This currently applies the all beams/IFs/pols
     82    """
     83    from asap._asap import add as _add
     84    return scantable(_add(scan, offset))
    6585
    6686
    … …  
    7090    from asap._asap import bin as _bin
    7191    return scantable(_bin(scan, binwidth))
     92
     93def average_pol(scan, mask=None):
     94    """
     95    Average the Polarisations together.
     96    Parameters:
     97        scan   - a scantable
     98        mask   - an optional mask defining the region, where
     99                 the averaging will be applied. The output
     100                 will have all specified points masked.
     101                 The dimension won't be reduced and
     102                 all polarisations will contain the
     103                 averaged spectrum.
     104    Example:
     105        polav = average_pols(myscan)
     106    """
     107    from asap._asap import averagepol as _avpol
     108    from numarray import ones
     109    if mask is None:
     110        mask = tuple(ones(scan.nchan()))
     111    return scantable(_avpol(scan, mask))
     112   
     113def hanning(scan):
     114    """
     115    Hanning smooth the channels.
     116    Parameters:
     117         scan    - the input scan
     118    Example:
     119         none
     120    """
     121    from asap._asap import hanning as _han
     122    return scantable(_han(scan))
     123
     124   
     125def poly_baseline(scan, mask=None, order=0):
     126    """
     127    Return a scan which has been baselined by a polynomial.
     128    Parameters:
     129        scan:    a scantable
     130        mask:    an optional mask
     131        order:   the order of the polynomial (default is 0)
     132    Example:
     133        # return a scan baselined by a third order polynomial,
     134        # not using a mask
     135        bscan = poly_baseline(scan, order=3)
     136    """
     137    from asap.asapfitter import fitter
     138    if mask is None:
     139        from numarray import ones
     140        mask = tuple(ones(scan.nchan()))
     141    f = fitter()
     142    f._verbose(True)
     143    f.set_scan(scan, mask)
     144    f.set_function(poly=order)   
     145    sf = f.auto_fit()
     146    return sf

  • trunk/python/asapreader.py

    r103 r113  
    55    This class allows the user to import single dish files
    66    (rpfits,sdfits,ms).
     7    The reader reads in integrations from the file and reamins at
     8    the fileposition afterwards.
    79    Available functions are:
    810
    … …  
    3335            integrations:    a 'range' of integration numbers, e.g.
    3436                             range(100) or [0,1,2,3,4,10,11,100]
     37                             If not given (default) all integrations
     38                             are read in
    3539        Example:
    3640            r.read([0,1,2,3,4])    # reads in the first 5 integatrions
    … …  
    4145        if integrations is None:
    4246            integrations = [-1]
     47        print "Reading intergrations from disk..."
    4348        sdreader.read(self,integrations)
    4449        tbl = sdreader.getdata(self)
    … …  
    4651
    4752    def summary(self):
    48         print "Disabled"
     53        """
     54        Print a summary of all scans/integrations. This reads through the
     55        whole file once.
     56        Parameters:
     57             None
     58        Example:
     59             r.summary()
     60        """
     61        sdreader.reset(self)
     62        sdreader.read(self,[-1])
     63        tbl = sdreader.getdata(self)
     64        sdreader.reset(self)
     65        print tbl.summary()
    4966        return
    50         sdreader.reset(self)
    51         sdreader.read([-1])
    52         sdreader.reset(self)
    53         tbl = sdreader.getdata(self)
    54         print tbl.summary()

  • trunk/python/scantable.py

    r102 r113  
    77        The ASAP container for scans
    88    """
    9     def functions(self):
    10         return ['copy','get_scan','summary','rms','set_unit','create_mask']
    11 
     9   
    1210    def __init__(self, filename):
    1311        """
    … …  
    1816                         scantable
    1917        """
     18        self._vb = True
     19        self._p = None
    2020        sdtable.__init__(self, filename)
    21         self._unit = 'channel'
     21
     22    def _verbose(self, *args):
     23        """
     24        Set the verbose level True or False, to indicate if output
     25        should be printed as well as returned.
     26        """
     27        if type(args[0]) is bool:
     28            self._vb = args[0]
     29            return
     30        elif len(args) == 0:
     31            return self._vb
     32
    2233
    2334    def copy(self):
    … …  
    2536        Return a copy of this scantable.
    2637        Parameters:
    27 
     38            none
    2839        Example:
    2940            copiedscan = scan.copy()
    3041        """
    31         sd = sdtable._copy(self)
    32         return scantable(sd)
     42        sd = scantable(sdtable._copy(self))
     43        return sd
     44
    3345    def get_scan(self, scanid=None):
    3446        """
    … …  
    3850            scanid:    a scanno or a source name
    3951        Example:
    40             scan.getscan('323p459')
     52            scan.get_scan('323p459')
    4153            # gets all scans containing the source '323p459'
    4254        """
    … …  
    4557        try:
    4658            if type(scanid) is str:
    47                 s = sdtable.getsource(self,scanid)
     59                s = sdtable._getsource(self,scanid)
    4860                return scantable(s)
    4961            elif type(scanid) is int:
    50                 s = sdtable.getscan(self,scanid)
     62                s = sdtable._getscan(self,scanid)
    5163                return scantable(s)
    5264        except RuntimeError:
    … …  
    8698
    8799    def get_selection(self):
     100        """
     101        Return/print a the current 'cursor' into the Beam/IF/Pol cube.
     102        Parameters:
     103            none
     104        Returns:
     105            a list of values (currentBeam,currentIF,currentPol)
     106        Example:
     107            none
     108        """
    88109        i = self.getbeam()
    89110        j = self.getif()
    90111        k = self.getpol()
    91         out = 'Beam=%d IF=%d Pol=%d '% (i,j,k)
    92         print out
     112        if self._vb:
     113            out = 'Beam=%d IF=%d Pol=%d '% (i,j,k)
     114            print out
    93115        return i,j,k
    94116
    … …  
    105127
    106128        Example:
    107             scan.setuint('channel')
     129            scan.set_unit('channel')
    108130            msk = scan.create_mask([100,200],[500,600])
    109131            scan.rms(mask=m)
    … …  
    124146                        tmp.append(rmsval)
    125147                        out += 'Beam[%d], IF[%d], Pol[%d] = %3.3f\n' % (i,j,k,rmsval)
    126             print out
     148            if self._vb:
     149                print out
    127150            return tmp
    128151
    … …  
    133156            rmsval = _rms(self,mask)
    134157            out = 'Beam[%d], IF[%d], Pol[%d] = %3.3f' % (i,j,k,rmsval)
    135             print out
     158            if self._vb:
     159                print out
    136160            return rmsval
    137161
    138162    def get_tsys(self, all=True):
     163        """
     164        Return the System temperatures.
     165        Parameters:
     166            all:    optional parameter to get the Tsys values for all
     167                    Beams/IFs/Pols (default) or just the one selected
     168                    with scantable.set_selection()
     169                    [True or False]
     170        Returns:
     171            a list of Tsys values.
     172        """
    139173        if all:
    140174            tmp = []
    … …  
    146180                    for k in range(self.npol()):
    147181                        self.setpol(k)
    148                         ts = self.gettsys()
     182                        ts = self._gettsys()
    149183                        tmp.append(ts)
    150184                        out += 'TSys: Beam[%d], IF[%d], Pol[%d] = %3.3f\n' % (i,j,k,ts)
    151             print out
     185            if self._vb:
     186                print out
    152187            return tmp
    153188        else:
    … …  
    155190            j = self.getif()
    156191            k = self.getpol()
    157             ts = self.gettsys()
     192            ts = self._gettsys()
    158193            out = 'TSys: Beam[%d], IF[%d], Pol[%d] = %3.3f' % (i,j,k,ts)
    159             print out
     194            if self._vb:
     195                print out
    160196            return ts
     197       
     198    def get_time(self):
     199        """
     200        Get a list of time stamps for the observations.
     201        Return a string for each intergration in the scantable.
     202        Parameters:
     203            none
     204        Example:
     205            none
     206        """
     207        out = []
     208        for i in range(self.nrow()):
     209            out.append(self._gettime(i))
     210        return out
    161211
    162212    def set_unit(self, unit='channel'):
    … …  
    165215        Parameters:
    166216            unit:    optional unit, default is 'channel'
    167                      one of 'GHz','MHz','km/s','channel'
    168         """
    169         units = ['GHz','MHz','km/s','channel','pixel','']
    170         if unit in units:
    171             if unit in ['','pixel']:
    172                 unit = 'channel'
    173             self._unit = unit
     217                     one of '*Hz','km/s','channel', ''
     218        """
     219
     220        if unit in ['','pixel', 'channel']:
     221            unit = ''
     222        inf = list(self._getcoordinfo())
     223        inf[0] = unit
     224        self._setcoordinfo(inf)
     225        if self._p: self.plot()
     226
     227    def set_freqframe(self, frame='LSRK'):
     228        """
     229        Set the frame type of the Spectral Axis.
     230        Parameters:
     231            frame:   an optional frame type, default 'LSRK'.
     232        Examples:
     233            scan.set_freqframe('BARY')
     234        """
     235        valid = ['REST','TOPO','LSRD','LSRK','BARY', \
     236                   'GEO','GALACTO','LGROUP','CMB']
     237        if frame in valid:
     238            inf = list(self._getcoordinfo())
     239            inf[1] = frame
     240            self._setcoordinfo(inf)
    174241        else:
    175             print "Invalid unit given, please use one of:"
    176             print units
    177 
    178     def create_mask(self, *args):
     242            print "Please specify a valid freq type. Valid types are:\n",valid
     243           
     244    def get_unit(self):
     245        """
     246        Get the default unit set in this scantable
     247        Parameters:
     248        Returns:
     249            A unit string
     250        """
     251        inf = self._getcoordinfo()
     252        unit = inf[0]
     253        if unit == '': unit = 'channel'
     254        return unit
     255
     256    def get_abscissa(self, rowno=0):
     257        """
     258        Get the abscissa in the current coordinate setup for the currently
     259        selected Beam/IF/Pol
     260        Parameters:
     261            none
     262        Returns:
     263            The abscissa values and it's format string.
     264        """
     265        absc = self.getabscissa(rowno)
     266        lbl = self.getabscissalabel(rowno)
     267        return absc, lbl
     268
     269    def create_mask(self, *args, **kwargs):
    179270        """
    180271        Compute and return a mask based on [min,max] windows.
     272        The specified windows are to be EXCLUDED, when the mask is
     273        applied.
    181274        Parameters:
    182275            [min,max],[min2,max2],...
    183276                Pairs of start/end points specifying the regions
    184277                to be masked
    185         Example:
    186             scan.setunit('channel')
     278            invert:     return an inverted mask, i.e. the regions
     279                        specified are not masked (INCLUDED)
     280        Example:
     281            scan.set_unit('channel')
     282
     283            a)
    187284            msk = scan.set_mask([400,500],[800,900])
    188             masks the regions between 400 and 500
    189             and 800 and 900 in the unit 'channel'
    190         """
    191         print "The current mask window unit is", self._unit
     285            # masks the regions between 400 and 500
     286            # and 800 and 900 in the unit 'channel'
     287
     288            b)
     289            msk = scan.set_mask([400,500],[800,900], invert=True)
     290            # masks the regions outside 400 and 500
     291            # and 800 and 900 in the unit 'channel'
     292           
     293        """
     294        u = self._getcoordinfo()[0]
     295        if self._vb:
     296            if u == "": u = "channel"
     297            print "The current mask window unit is", u
    192298        n = self.nchan()
    193         if self._unit == 'channel':
    194             data = range(n)
    195         else:
    196             data = self.getabscissa(unit=self._unit)
     299        data = self.getabscissa()
    197300        msk = ones(n)
    198301        for  window in args:
    … …  
    201304                return
    202305            for i in range(n):
    203                 if data[i] > window[0] and data[i] < window[1]:
     306                if data[i] >= window[0] and data[i] < window[1]:
    204307                    msk[i] = 0
     308        if kwargs.has_key('invert'):
     309            if kwargs.get('invert'):
     310                from numarray import logical_not
     311                msk = logical_not(msk)
    205312        return msk
    206     def set_restfrequencies(self, freqs, unit='Hz'):
     313   
     314    def set_restfreqs(self, freqs, unit='Hz'):
    207315        """
    208316        Set the restfrequency(s) for this scantable.
    … …  
    211319            unit:     optional 'unit', default 'Hz'
    212320        Example:
    213             scan.set_restfrequencies([1000000000.0])
    214         """
    215         self.setrestfreqs(freqs, unit)
     321            scan.set_restfreqs([1000000000.0])
     322        """
     323        if type(freqs) is float or int:
     324            freqs = (freqs)
     325        sdtable._setrestfreqs(self,freqs, unit)
    216326        return
    217327
    … …  
    236346            print "Please specify a valid (Beam/IF/Pol)"
    237347        return
     348
     349    def plot(self, what='spectrum',col='Pol', panel=None):
     350        """
     351        Plot the spectra contained in the scan. Alternatively you can also
     352        Plot Tsys vs Time
     353        Parameters:
     354            what:     a choice of 'spectrum' (default) or 'tsys'
     355            col:      which out of Beams/IFs/Pols should be colour stacked
     356            panel:    set up multiple panels, currently not working.
     357        """
     358        validcol = {'Beam':self.nbeam(),'IF':self.nif(),'Pol':self.npol()}
     359        validyax = ['spectrum','tsys']
     360        if not self._p:
     361            from asap.asaplot import ASAPlot
     362            self._p = ASAPlot()
     363#            print "Plotting not enabled"
     364#            return
     365        npan = 1
     366        x = None
     367        if what == 'tsys':
     368            n = self.nrow()
     369            if n < 2:
     370                print "Only one integration. Can't plot."
     371                return
     372        if panel == 'Time':
     373            npan = self.nrow()
     374        self._p.hold()
     375        self._p.clear()
     376        xlab,ylab,tlab = None,None,None
     377        vb = self._verbose
     378        self._verbose(False)
     379        sel = self.get_selection()
     380        for i in range(npan):
     381            for j in range(validcol[col]):
     382                x = None
     383                y = None
     384                m = None
     385                tlab = self._getsourcename(i)
     386                if col == 'Beam':
     387                    self.setbeam(j)
     388                elif col == 'IF':
     389                    self.setif(j)
     390                elif col == 'Pol':
     391                    self.setpol(j)
     392                if what == 'tsys':
     393                    x = range(self.nrow())
     394                    xlab = 'Time [pixel]'
     395                    m = list(ones(len(x)))
     396                    y = []
     397                    ylab = r'$T_{sys}$'
     398                    for k in range(len(x)):
     399                        y.append(self._gettsys(k))
     400                else:
     401                    x,xlab = self.get_abscissa(i)
     402                    y = self.getspectrum(i)
     403                    ylab = r'Flux'
     404                    m = self.getmask(i)
     405                llab = col+' '+str(j)
     406                self._p.set_line(label=llab)
     407                self._p.plot(x,y,m)
     408        self._p.set_axes('xlabel',xlab)
     409        self._p.set_axes('ylabel',ylab)
     410        self._p.set_axes('title',tlab)
     411        self._p.release()
     412        self.set_selection(sel[0],sel[1],sel[2])
     413        self._verbose(vb)
     414        return
Note: See TracChangeset for help on using the changeset viewer.