imstat.py

Go to the documentation of this file.

#
# This file was generated using xslt from its XML file
#
# Copyright 2009, Associated Universities Inc., Washington DC
#
import sys
import os
from  casac import *
import string
from taskinit import casalog
#from taskmanager import tm
import task_imstat
def imstat(imagename='', axes=-1, region='', box='', chans='', stokes='', listit=True, verbose=True, mask='', stretch=False, logfile='', append=True):

        """Displays statistical information from an image or image region

     Many parameters are determined from the specified region of an image.
     For this version, the region can be specified by a set of rectangular
     pixel coordinates, the channel ranges and the Stokes.

     For directed output, run as 
                    myoutput = imstat()
   

Keyword arguments:
imagename    Name of input image
                     Default: none; Example: imagename='ngc5921_task.im'
axes         axes to compute statistics over. -1 => all axes.
region       Region of interest. See help par.region.
box          A box region specified in pixels on the directional plane
                     Default: none (whole 2-D plane); 
             Example: box='10,10,50,50'
             box = '10,10,30,30,35,35,50,50' (two boxes)
chans        Zero based channel numbers
                 Range of channel numbers to include in statistics
             All spectral windows are included
                     Default:''= all;  Example: chans='3~20'    
stokes       Stokes parameters to analyze.
             Default: all; Example: stokes='IQUV';  
             Example:stokes='I,Q'
             Options: 'I','Q','U','V','RR','RL','LR','LL','XX','YX','XY','YY', ...
listit       Print stats and bounding box to logger? 
verbose      Print additional messages to logger?
mask         Mask to use. See help par.mask. Default is none.
stretch      Stretch the mask if necessary and possible? See help par.stretch 
logfile      Name of file to write fit results.
append       If logfile exists, append to it (True) or overwrite it (False).


      General procedure:

         1.  Specify inputs, then

         2.  myoutput = imstat()
               or specify inputs directly in calling sequence to task
             myoutput = imstat(imagename='image.im', etc)

         3.  myoutput['KEYS'] will contain the result associated with any
               of the keys given below
        
        KEYS CURRENTLY AVAILABLE
        blc          - absolute PIXEL coordinate of the bottom left corner of 
                       the bounding box surrounding the selected region
        blcf         - Same as blc, but uses WORLD coordinates instead of pixels
        trc          - the absolute PIXEL coordinate of the top right corner 
                       of the bounding box surrounding the selected region
        trcf         - Same as trc, but uses WORLD coordinates instead of pixels
        flux         - the integrated flux density if the beam is defined and 
                       the if brightness units are $Jy/beam$
        npts         - the number of unmasked points used
        max          - the maximum pixel value
        min          - minimum pixel value
        maxpos       - absolute PIXEL coordinate of maximum pixel value
        maxposf      - Same as maxpos, but uses WORLD coordinates instead of pixels
        minpos       - absolute pixel coordinate of minimum pixel value
        minposf      - Same as minpos, but uses WORLD coordinates instead of pixels
        sum          - the sum of the pixel values: $\sum I_i$
        sumsq        - the sum of the squares of the pixel values: $\sum I_i^2$
        mean         - the mean of pixel values: 
                       $\bar{I} = \sum I_i / n$
        sigma        - the standard deviation about the mean: 
                       $\sigma^2 = (\sum I_i - \bar{I})^2 / (n-1)$
        rms          - the root mean square: 
                       $\sqrt {\sum I_i^2 / n}$
        median       - the median pixel value
        medabsdevmed - the median of the absolute deviations from the 
                       median
        quartile     - the inter-quartile range. Find the points 
                       which are 25% largest and 75% largest (the median is 
                       50% largest), find their difference and divide that 
                       difference by 2.

        Additional Examples
        # Selected two box region
        # box 1, bottom-left coord is 2,3 and top-right coord is 14,15
        # box 2, bottom-left coord is 30,31 and top-right coord is 42,43
        imstat( 'myImage', box='2,3,14,15;30,31,42,43' )

        # Select the same two box regions but only channels 4 and 5
        imstat( 'myImage', box='2,3,14,15;30,31,42,43', chan='4~5' )

        # Select all channels greater the 20 as well as channel 0.
        # Then the mean and standard deviation are printed
        results = imstat( 'myImage', chans='>20;0' )
        print "Mean is: ", results['mean'], "  s.d. ", results['sigma']

        # Find statistical information for the Q stokes value only
        # then the I stokes values only, and printing out the statistical
        # values that we are interested in.
        s1 = imstat( 'myimage', stokes='Q' )
        s2 = imstat( 'myimage', stokes='I' )
        print "       |  MIN  |   MAX  | MEAN"
        print "  Q    | ",s1['min'][0],"  |  ",s1['max'][0],"  |  ",,"  |  ",s1['mean'][0]
        print "  I    | ",s2['min'][0],"  |  ",s2['max'][0],"  |  ",,"  |  ",s2['mean'][0]

# evaluate statistics for each spectral plane in an ra x dec x frequency image
ia.fromshape("", [20,30,40])
# give pixels non-zero values
ia.addnoise()
# These are the display axes, the calculation of statistics occurs
# for each (hyper)plane along axes not listed in the axes parameter,
# in this case axis 2 (the frequency axis)
# display the rms for each frequency plane (your mileage will vary with
# the values).
stats = ia.statistics(axes=[0,1])
 stats["rms"]
  Out[10]: 
array([ 0.99576014,  1.03813124,  0.97749186,  0.97587883,  1.04189885,
        1.03784776,  1.03371549,  1.03153074,  1.00841606,  0.947155  ,
        0.97335404,  0.94389403,  1.0010221 ,  0.97151822,  1.03942156,
        1.01158476,  0.96957082,  1.04212773,  1.00589049,  0.98696715,
        1.00451481,  1.02307892,  1.03102005,  0.97334671,  0.95209879,
        1.02088714,  0.96999902,  0.98661619,  1.01039267,  0.96842754,
        0.99464947,  1.01536798,  1.02466023,  0.96956468,  0.98090756,
        0.9835844 ,  0.95698935,  1.05487967,  0.99846411,  0.99634868])



        """

#
#    The following is work around to avoid a bug with current python translation
#
        mytmp = {}

        mytmp['imagename'] = imagename
        mytmp['axes'] = axes
        mytmp['region'] = region
        mytmp['box'] = box
        mytmp['chans'] = chans
        mytmp['stokes'] = stokes
        mytmp['listit'] = listit
        mytmp['verbose'] = verbose
        mytmp['mask'] = mask
        mytmp['stretch'] = stretch
        mytmp['logfile'] = logfile
        mytmp['append'] = append
        pathname='file:///'+os.environ.get('CASAPATH').split()[0]+'/share/xml/'
        trec = casac.utils().torecord(pathname+'imstat.xml')

        casalog.origin('imstat')
        if trec.has_key('imstat') and casac.utils().verify(mytmp, trec['imstat']) :
            result = task_imstat.imstat(imagename, axes, region, box, chans, stokes, listit, verbose, mask, stretch, logfile, append)

        else :
          result = False
        return result