widefield.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_widefield
def widefield(vis=[''], imagename='', outlierfile='', field='', spw='', selectdata=False, timerange='', uvrange='', antenna='', scan='', mode='mfs', niter=500, gain=0.1, threshold='0.0Jy', psfmode='clark', ftmachine='', facets=3, wprojplanes=64, multiscale=[], negcomponent=0, interactive=False, mask=[], nchan=1, start=0, width=1, imsize=[256, 256], cell=['1.0arcsec', '1.0arcsec'], phasecenter='', restfreq='', stokes='I', weighting='natural', robust=0.0, npixels=0, noise='1.0Jy', cyclefactor=1.5, cyclespeedup=-1, npercycle=100, uvtaper=False, outertaper=[''], innertaper=['1.0'], restoringbeam=[''], calready=False):

        """Wide-field imaging and deconvolution with selected algorithm

        Wide-field imaging and deconvolution with selected algorithm:

        This is the main wide-field imaging/deconvolution task.  It
        uses the wprojection method for a large field of view, can
        make many facets, and can include outlier fields.  Several
        deconvolution algorithms are supported.  Interactive cleaning
        is also supported.

        For making large images (>2000 on a size), see hints at the
        end of the descriptions.  For making images larger than about
        5000x5000, the available memory must be larger than 2 Gbytes. For such 
        images therefore  a computer with a 64-bit operating system may be
        needed.


        Keyword arguments:
        vis -- Name of all input visibility files
                default: none; example: vis='ngc5921.ms'
                example: vis=['data01.ms', 'data02.ms']
        imagename -- Pre-name of output images:
                default: none; example: imagename='n5921'
                if outlier fields are included, then
                   imagename=['n5921', 'outlier1', outlier2']
                   and the first imagename is the wide-field image
                output images names are: n5921.clean, n5921.residual,
                n5921.model, n5921.interactive.mask
        mode -- Type of selection
                default: 'mfs'; example: mode='channel';
                Options: 'mfs', channel, velocity, frequency'
        alg -- Algorithm to use
                default: 'clark';
                Options: 'clark', 'hogbom','multiscale','entropy'
                    Strongly advise 'clark'.  multiscale and entropy
                    well-tested.
        imsize -- Image pixel size (x,y)
                default = [256,256]; example: imsize=[500,500], or imsize=500
                example for multiple fields: imsize=[(1000, 1000), (100, 100)]
        cell -- Cell size (x,y)
                default=['1arcsec,'1arcsec']
                example: cell=['0.5arcsec,'0.5arcsec'], or cell='0.5arcsec'
        phasecenter -- direction position or the field for the image center
                A list of the above is needed for multiple-fields
                default: '' -->field='0' as center; example: phasecenter='6'
                   phasecenter='J2000 19h30m00 -40d00m00'
                   phasecenter=['J2000 19h30m00 -40d00m00', 'J2000 19h57m00 40d00m00']
                      for wide-field, plus one outlier field.
        stokes -- Stokes parameters to image
                default='I'; example: stokes='IQUV';
                Options: 'I','IV','IQU','IQUV'
        niter -- Number iterations, set to zero for no CLEANing
                default: 500; example: niter=500
        gain -- Loop gain for CLEANing
                default: 0.1; example: gain=0.1
        threshold -- Flux level at which to stop CLEANing (units=mJy)
                default: 0.0; example: threshold=0.0
        mask -- Name(s) of mask image(s) used for CLEANing
                default: ''  example: mask='orion.mask'
                Number of mask fields must equal number of imaged fields
        cleanbox -- List of [blc-x,blc-y,trc-x,trc-y] values
                default: []; example: cleanbox=[110,110,150,145]
                Note: This can also be a filename with clean values:
                fieldindex blc-x blc-y trc-x trc-y
                cleanbox = 'interactive' is very useful.
        --- Data Selection
        nchan -- Number of channels to select
                default: 1; example: nchan=45
        start -- Start channel, 0-relative
                default=0; example: start=5
                if mode='frequency' then a frequency value e.g start='1.4GHz'
        width -- Channel width (value > 1 indicates channel averaging)
                default=1; example: width=5
                if mode='frequency' then a frequency value e.g  width='10kHz'
        step -- Step in channel number
                default=1; example: step=2
        field -- Select field using field id(s) or field name(s).
                  [run listobs to obtain the list id's or names]
               default: ''=all fields
               If field string is a non-negative integer, it is assumed a field index
                 otherwise, it is assumed a field name
               field='0~2'; field ids 0,1,2
               field='0,4,5~7'; field ids 0,4,5,6,7
               field='3C286,3C295'; field named 3C286 adn 3C295
               field = '3,4C*'; field id 3, all names starting with 4C
               example for multiple ms in vis parameter:
               field=['0~2', '1,2']
        spw -- Select spectral window/channels
               default: ''=all spectral windows and channels
               spw='0~2,4'; spectral windows 0,1,2,4 (all channels)
               spw='<2';  spectral windows less than 2 (i.e. 0,1)
               spw='0:5~61'; spw 0, channels 5 to 61
               spw='0,10,3:3~45'; spw 0,10 all channels, spw 3, channels 3 to 45.
               spw='0~2:2~6'; spw 0,1,2 with channels 2 through 6 in each.
               spw='0:0~10;15~60'; spectral window 0 with channels 0-10,15-60
               spw='0:0~10,1:20~30,2:1;2;3'; spw 0, channels 0-10,
                        spw 1, channels 20-30, and spw 2, channels, 1,2 and 3
               For multiple ms in vis parameter:
               spw=['0,10,3:3~45', '<2']
        timerange -- Select time range subset of data (not implemented yet)
            default='' meaning no time selection
            example: timerange='YYYY/MM/DD/HH:MM:SS.sss'
            timerange='< YYYY/MM/DD/HH:MM:SS.sss'
            timerange='> YYYY/MM/DD/HH:MM:SS.sss'
            timerange='ddd/HH:MM:SS.sss'
            timerange='< ddd/HH:MM:SS.sss'
            timerange='> ddd/HH:MM:SS.sss'
        restfreq -- Specify rest frequency to use for image
            default='' (i.e., try to use the restfreq specified in the visibility data)

        --- Weighting
        weighting -- Weighting to apply to visibilities
                default='natural'; example: weighting='uniform';
                Options: 'natural','uniform','briggs','briggsabs','radial', 'superuniform'
        robust -- 'briggs' and 'brigssabs' robustness parameter
                default=0.0; example: robust=0.5;
                Options: -2.0 to 2.0; -2 (uniform)/+2 (natural)
        npixels -- number of pixels to determine uv-cell size for weight calculation
                 -- Used with superuniform or briggs weighting schemes
                  example: npixels=3

        --- widefield controls
        ftmachine -- Gridding method for the image;
                ft (standard interferometric gridding).
                wproject (wprojection algorithm for gridding)
                default: wproject
        wprojplanes -- Number w-projection planes to use for gridding
                default: 256
                example: wprojplanes=64
                   Good value = BMAX(klambda) * Map width(arcmin)^2 / 600
        facets   -- Number of facets along one axis on central image
                image is divided in facets x facets rectangles.
                default: 1
                example: facets=3 makes 3x3 images to cover the field
                if ftmachine = 'ft', only faceting is used
                if ftmachine = 'wproject', both wplanes and faceting  
                         can be used  (see below).

        cyclefactor -- Change the threshold at which the deconvolution cycle will
                stop and degrid and subtract from the visibilities. For bad PSFs,
                reconcile often (cyclefactor=4 or 5); For good PSFs, use
                cyclefactor 1.5 to 2.0.
                default: 2.5; example: cyclefactor=4, but decreases speed considerably.
                <cycle threshold = cyclefactor * max sidelobe * max residual>
        cyclespeedup -- Cycle threshold doubles in this number of iterations
                default: -1; example: cyclespeedup=500

        --- MEM parameters (Experimental, not well-tested)
        sigma -- Target image sigma
                default: '0.001Jy'; example: sigma='0.1Jy'
        targetflux -- Target flux for final image
                default: '1.0Jy'; example: targetflux='200Jy'
        constrainflux -- Constrain image to match target flux;
                otherwise, targetflux is used to initialize model only.
                default: False; example: constrainflux=True
        prior -- Name of MEM prior images
                default: ['']; example: prior='source_mem.image'

        --- Multi-scale parameters (Experimental, not well-tested)
        negcomponent -- Stop component search when the largest scale has found this
                number of negative components; -1 means continue component search
                even if the largest component is negative.
                default: 2; example: negcomponent=-1
        scales -- Used for alg='multiscale'; set a number of scales or a vector
                default: [0,3,10]; example: scales=[0.0,3.0,10.0, 30]
        --  interactive masking
        npercycle -- when cleanbox is set to 'interactive',
           this is the number of iterations between each clean to update mask
           interactively. However, this number can be adjusted during execution.

        uvtaper -- Apply additional uv tapering of the visibilities.
               default: uvtaper=False; example: uvtaper=True
                  uvtaper=True expandable parameters
                     outertaper -- uv-taper on outer baselines in uv-plane
                           [bmaj, bmin, bpa] taper Gaussian scale in uv or
                            angular units. NOTE: uv taper in (klambda) is  
                            roughly on-sky FWHM(arcsec/200)
                         default: outertaper=[]; no outer taper applied
                            example: outertaper=['5klambda']  circular taper
                                FWHM=5 kilo-lambda
                                outertaper=['5klambda','3klambda','45.0deg']
                                outertaper=['10arcsec'] on-sky FWHM 10"
                                outertaper=['300.0'] default units are meters
                                   in aperture plane
                     innertaper -- uv-taper in center of uv-plane
                             NOT YET IMPLEMENTED

        restoringbeam -- Output Gaussian restoring beam for CLEAN image
                [bmaj, bmin, bpa] elliptical Gaussian restoring beam
                default units are in arc-seconds for bmaj,bmin, degrees
                for bpa default: restoringbeam=[]; Use PSF calculated
                from dirty beam.
                example: restoringbeam=['10arcsec'] circular Gaussian
                       FWHM 10" example:
                       restoringbeam=['10.0','5.0','45.0deg'] 10"x5"
                       at 45 degrees

        calready -- if True will create scratch columns if they are 
                not there. And after clean completes the predicted model 
                visibility is from the clean components are
                written to the ms.        

        async --  Run asynchronously
                default = False; do not run asychronously

 ======================================================================

                      HINTS ON RUNNING WIDEFIELD

      1.  Decide if the images will be specified directly in the  
          inputs or with an outlier file.  For more than a few fields,
          an outlier file more convenient.

         Direct Method:

            cell = ['1.0arcsec', '1.0arcsec']
            imagename = ['M1_0','M1_1','M1_2]
            imsize = [[1024,1024],[128,128],[128,128]]
            phasecenter = ['J2000 13h27m20.98 43d26m28.0', 
                      'J2000 13h30m52.159 43d23m08.02', 'J2000 13h24m08.16 43d09m48.0']

          Text file method  (in outlier.txt)

            imagename = 'M1'
            outlierfile = 'outlier.txt'
               [phasecenter, imsize ignored]

            Contents of outlier.txt
            C   0   1024 1024   13 27 20.98     43 26 28.0
            C   1    128  128   13 30 52.158    43 23 08.00
            C   2    128  128   13 24 08.163    43 09 48.00

         In both cases the following images will be made:
             M1_0.image, M1_1.image, M1_2.image     cleaned images
             M1.0.model, M1_1.model, M1_2.model     model images
             M1.0.residual, M1_1.residual, M1_2.residual     residual images

       2.  Wprojection:  It is fastest to use wprojection without faceting.
             ftmachine = 'wproject'
             wprojplane = NN

           The value of NN should be chosen as small as possible to reduce
           execution time.  The algorithm
               NN = BMAX(klambda) * imagewidth (arcmin)^2 / 600, with a minimum
                    of 16, should be adequate.

       3.  Depending on the memory of the computer, a limit of about
       5000x5000 may occur for example if a computer has 2Gbyte of
       RAM. Also a 32-bit computer has a maximum limit of 2Gbyte
       memory usable per process, irrespective of how much physical
       RAM is present. Hence it is recommended to move to a 64-bit
       computer with more than 2 GByte of RAM for >5000x5000 images
           

       4. For data with extremely large 'w' values, i.e low frequency,
       long baseline and very widefield image, the wprojection
       convolution can be very large and either not fit in memory or
       slow for processing.  In these cases you should consider using
       both ftmachine='wproject' and facets=xx where is 3.

  
        """
        if type(vis)==str: vis=[vis]
        if type(multiscale)==int: multiscale=[multiscale]
        if type(imsize)==int: imsize=[imsize]
        if type(cell)==float: cell=[cell]
        if type(outertaper)==str: outertaper=[outertaper]
        if type(innertaper)==str: innertaper=[innertaper]
        if type(restoringbeam)==str: restoringbeam=[restoringbeam]

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

        mytmp['vis'] = vis
        mytmp['imagename'] = imagename
        mytmp['outlierfile'] = outlierfile
        mytmp['field'] = field
        mytmp['spw'] = spw
        mytmp['selectdata'] = selectdata
        mytmp['timerange'] = timerange
        mytmp['uvrange'] = uvrange
        mytmp['antenna'] = antenna
        mytmp['scan'] = scan
        mytmp['mode'] = mode
        mytmp['niter'] = niter
        mytmp['gain'] = gain
        mytmp['threshold'] = threshold
        mytmp['psfmode'] = psfmode
        mytmp['ftmachine'] = ftmachine
        mytmp['facets'] = facets
        mytmp['wprojplanes'] = wprojplanes
        mytmp['multiscale'] = multiscale
        mytmp['negcomponent'] = negcomponent
        mytmp['interactive'] = interactive
        mytmp['mask'] = mask
        mytmp['nchan'] = nchan
        mytmp['start'] = start
        mytmp['width'] = width
        mytmp['imsize'] = imsize
        if type(cell) == str :
           mytmp['cell'] = casac.quanta().quantity(cell)
        else :
           mytmp['cell'] = cell
        mytmp['phasecenter'] = phasecenter
        mytmp['restfreq'] = restfreq
        mytmp['stokes'] = stokes
        mytmp['weighting'] = weighting
        mytmp['robust'] = robust
        mytmp['npixels'] = npixels
        mytmp['noise'] = noise
        mytmp['cyclefactor'] = cyclefactor
        mytmp['cyclespeedup'] = cyclespeedup
        mytmp['npercycle'] = npercycle
        mytmp['uvtaper'] = uvtaper
        mytmp['outertaper'] = outertaper
        mytmp['innertaper'] = innertaper
        mytmp['restoringbeam'] = restoringbeam
        mytmp['calready'] = calready
        pathname='file:///'+os.environ.get('CASAPATH').split()[0]+'/share/xml/'
        trec = casac.utils().torecord(pathname+'widefield.xml')

        casalog.origin('widefield')
        if trec.has_key('widefield') and casac.utils().verify(mytmp, trec['widefield']) :
            result = task_widefield.widefield(vis, imagename, outlierfile, field, spw, selectdata, timerange, uvrange, antenna, scan, mode, niter, gain, threshold, psfmode, ftmachine, facets, wprojplanes, multiscale, negcomponent, interactive, mask, nchan, start, width, imsize, cell, phasecenter, restfreq, stokes, weighting, robust, npixels, noise, cyclefactor, cyclespeedup, npercycle, uvtaper, outertaper, innertaper, restoringbeam, calready)

        else :
          result = False
        return result