Open in Colab: https://colab.research.google.com/github/casangi/casadocs/blob/master/docs/notebooks/UVTaper_Imaging_Weights.ipynb

---

UV-Taper Imaging Weights

UV-taper imaging weights are applied by calculating a 2D Gaussian in the UV-domain and multiplying it with the weightdensity grid (from previous stages of weighting schemes).

General formulae

• The relation between \(\sigma\) and the full-width-half-maximum \(\theta\) of a Gaussian is given by \(\theta = \sigma\sqrt{8\ln2}\) or \(\sigma = \frac{\theta}{\sqrt{8\ln2}}\). (Ref : https://brainder.org/2011/08/20/gaussian-kernels-convert-fwhm-to-sigma/ )

• For a Gaussian, the Fourier dual relationship (generalized to any number of dimensions) is given by \(\sigma_{uv} = \frac{1}{2 \pi \sigma_{lm}}\) (Ref : Bracewell)

These two expressions may be used to convert uvtaper inputs in either the image or uv-domain, to a \(\sigma_{uv}\) that may be used to evaluate the uv-domain Gaussian.

Evaluating the UV-Taper Gaussian

(1) Calculate the \(\sigma_{uv}\) of the UV-domain Gaussian

Convert the uvtaper input information into \(\sigma\) for the uv-domain Gaussian, separately for the major and minor axes. The input uvtaper information can be given either in the image domain or uv domain.

Input : FWHM in the image domain : beam_radian

The interpretation is that of convolving the existing beam with a ‘uvtaper beam’. The FWHM is specified in units of angle (deg,arcmin,arcsec), converted to radians. - $:nbsphinx-math:theta_{lm} = $ beam_radian - \(\sigma_{lm} = \frac{\theta_{lm}}{ 2 \sqrt{2 \ln2}}\) - \(\sigma_{uv} = \frac{1}{2 \pi \sigma_{lm}}\)

Input : HWHM in the uv domain : beam_lambda

The interpretation is ‘uv distance’ or ‘baseline length’ and is the half width of a Gaussian. The HWHM is specified in units of wavelength (lambda) - $:nbsphinx-math:theta_{uv} = 2 $ beam_lambda - \(\sigma_{uv}= \frac{\theta_{uv}}{ 2 \sqrt{2 \ln2}}\)

For each uv cell in the weightdensity grid, rotate the u and v vectors by -PA to get the effective distances along the major and minor axis of the uvtaper elliptical Gaussian. In this calculation, \(u,v\) are in units of meters, and \(u_{rot}, v_{rot}\) are in units of wavelength. This is a rotation by -PA. - \(u_{rot} = \frac{1}{\lambda_{obs}} \left( u~cos(pa) + v~sin(pa) \right)\) - \(v_{rot} = \frac{1}{\lambda_{obs}} \left( - u~sin(pa) + v~cos(pa) \right)\)

With this calculation, if the uvtaper Gaussian has PA=0, then \(u_{rot}=u, v_{rot}=v\). If the uvtaper Gaussian has PA=90, then \(u_{rot}=v, v_{rot}= -u\) .

Evaluate the Gaussian as \(e^{-\frac{1}{2} \left( \frac{u_{rot}^2}{\sigma_{uv,maj}^2} + \frac{v_{rot}^2}{\sigma_{uv,min}^2} \right)}\)

Each uv cell in the weightdensity grid is multiplied with this evaluated Gaussian.

Implementation in CASA : VisImagingWeight.cc

A quantity representing $:nbsphinx-math:frac{1}{ 2 sigma_{uv}^2} $ is computed for both ways of specifying uvtaper. This is stored in the variable ‘fact’. Following through the above math, the two input options result in the following. - For input as \(\theta_{lm}\), \(xx = \left( \frac{\pi^2}{4 \ln2} \right) beam_{radian}^2\) - For input as \(\theta_{uv}/2\), \(xx = \frac{\ln 2}{beam_{lambda}^2}\)

The Gaussian is then evaluated as \(e^{\left(- u_{rot}^2 xx_{maj} - r_{rot}^2 xx_{min} \right)}\)

Testing the calculations

Let us consider a dataset with a naturally roundish PSF of beam width approximately 45 arcsec (e.g. VLA D-config at L-Band). A uniformly weighted PSF has a beam width of about 28 arcsec. The imaging cell size is therefore chosen to be 5.0arcsec.

UV-taper values to test should be comparable to or larger than this natural PSF and so we choose 100 arcsec and 200 arcsec. We also test round and elliptical tapers, and start with naturally-weighted as well as uniformly-weighted PSFs. In all cases, inputs are supplied both as the FWHM in the image domain, as well as the HWHM in the UV-domain, and the output restoring beams must be equivalent with both these input methods (to well-within the pixel size of 5.0arcsec).

from casatasks import tclean,imhead
import os, numpy

def convert_lm_to_uv(theta_lm):
    """
    theta_lm : Full width at half maximum, in the image domain, in units of arcsec
    """
    theta_uv = ( 4*numpy.log(2)/numpy.pi ) / ( (theta_lm / 3600) * numpy.pi/180.0)
    print("FWHM of %3.3f arcsec maps to a FWHM of %3.3f lambda"%(theta_lm,theta_uv))
    return theta_uv

def convert_uv_to_lm(theta_uv):
    """
    theta_uv : Full width at half maximum, in the uv domain, in units of lambda
    """
    theta_lm  = 3600 * ( 4*numpy.log(2)/numpy.pi ) / (theta_uv * numpy.pi/180.0)


    print("FWHM of %3.3f arcsec maps to a FWHM of %3.3f lambda"%(theta_lm,theta_uv))

    return theta_lm

convert_uv_to_lm(convert_lm_to_uv(100.0));

FWHM of 100.000 arcsec maps to a FWHM of 1820.374 lambda
FWHM of 100.000 arcsec maps to a FWHM of 1820.374 lambda

def dispbeam(beam):
    """
    Print restoring beam info...
    """
    print("Restoring Beam : %3.4f %s  X  %3.4f %s ,  %3.4f %s"%(beam['major']['value'],
                                               beam['major']['unit'],
                                               beam['minor']['value'],
                                               beam['minor']['unit'],
                                               beam['positionangle']['value'],
                                               beam['positionangle']['unit'] ))

def run_im(imnames,uvtapers,weighting='natural'):
    os.system('rm -rf uvt*')
    vis = '/export/home/riya/rurvashi/CASADATA/casatestdata/unittest/tclean/refim_point.ms'

    for (imname,uvtaper) in zip(imnames,uvtapers):
        print("\n%s : uvtaper = %s"%(imname,uvtaper))
        tclean(vis=vis, spw='0:10', imagename=imname,
               uvtaper=uvtaper,
               weighting=weighting,
               imsize=200,
               cell='5.0arcsec',niter=0,
               calcpsf=True,calcres=False,
               restoration=False)
        beam =imhead(imname+'.psf')['restoringbeam']
        dispbeam(beam)

Theoretical Calculation :

When a UV-taper equivalent to \(\sigma_{taper,lm}\) in the image domain is applied to pre-existing beam of width \(\sigma_{orig}\), the two beams effectively convolve with each other, resulting in a new beam whose sigma is given by $ \sigma_{new} = \sqrt{\sigma_{orig}^2 + \sigma_{taper,lm}^2} $, where \(\sigma = \frac{\theta}{\sqrt{8\ln2}}\). This theoretical calculation is only an approximation (and usually a bad one) because the original beam is almost never an exact Gaussian. Therefore, it should be used with care.

def calc_convolve(theta_orig, theta_taper):
    """
    Calculate the width of a Gaussian resulting from the convolution of two Gaussians.
    This calculation is only for a circular Gaussian.
    Units of inputs : arcsec.
    """
    arcsec_to_radians = (1/3600.0)*numpy.pi/180.0
    sigma_orig = arcsec_to_radians * theta_orig/numpy.sqrt(8*numpy.log(2.0))
    sigma_taper = arcsec_to_radians * theta_taper/numpy.sqrt(8*numpy.log(2.0))

    sigma_new = numpy.sqrt(sigma_orig**2 + sigma_taper**2)
    theta_new = sigma_new * numpy.sqrt(8*numpy.log(2.0)) / arcsec_to_radians

    print("Convolution of FWHMs of %3.4f arcsec and %3.4f arcsec \
          yields %3.4f arcsec"%(theta_orig, theta_taper, theta_new))

Specify the uvtaper as a single quantity, either in the image domain or in the uv domain. Run tclean with no taper, taper specified in the image domain, and taper specified in the uv-domain. Extract and print out the restoring beam results.

imnames = ['uvt_orig' , 'uvt_taper_im' , 'uvt_taper_uv']

## Round tapers
imtaper = 100 #arcsec
uvtapers=['' , '%3.2farcsec'%(imtaper) , '%3.2flambda'%(convert_lm_to_uv(imtaper)/2.0)]

print("\nSettings for uvtaper in tclean : \n\
[ None,  FWHM in the image domain, HWHM in the uv-domain] ")
print(uvtapers)

FWHM of 100.000 arcsec maps to a FWHM of 1820.374 lambda

Settings for uvtaper in tclean :
[ None,  FWHM in the image domain, HWHM in the uv-domain]
['', '100.00arcsec', '910.19lambda']

run_im(imnames,uvtapers,weighting='natural')

uvt_orig : uvtaper =
Restoring Beam : 49.8538 arcsec  X  46.8010 arcsec ,  -87.6734 deg

uvt_taper_im : uvtaper = 100.00arcsec
Restoring Beam : 143.2574 arcsec  X  136.1750 arcsec ,  88.5291 deg

uvt_taper_uv : uvtaper = 910.19lambda
Restoring Beam : 143.2571 arcsec  X  136.1747 arcsec ,  88.5291 deg

calc_convolve(theta_orig=50.0, theta_taper=100.0)

Convolution of FWHMs of 50.0000 arcsec and 100.0000 arcsec           yields 111.8034 arcsec

run_im(imnames,uvtapers,weighting='uniform')

uvt_orig : uvtaper =
Restoring Beam : 28.2215 arcsec  X  27.3716 arcsec ,  -11.8527 deg

uvt_taper_im : uvtaper = 100.00arcsec
Restoring Beam : 100.5366 arcsec  X  100.1634 arcsec ,  55.6722 deg

uvt_taper_uv : uvtaper = 910.19lambda
Restoring Beam : 100.5363 arcsec  X  100.1631 arcsec ,  55.6722 deg

calc_convolve(theta_orig=28.0, theta_taper=100.0)

Convolution of FWHMs of 28.0000 arcsec and 100.0000 arcsec           yields 103.8460 arcsec

Specify the uvtaper as vector of quantities [bmaj.bmin,positionangle], either in the image domain or in the uv domain. Run tclean with no taper, taper specified in the image domain, and taper specified in the uv-domain. Extract and print out the restoring beam results.

For an elliptical Gaussian taper, we need to flip the inputs for the major and minor axis after conversion to the uv-domain, and rotate by 90deg.

imnames = ['uvt_orig' , 'uvt_taper_im' , 'uvt_taper_uv']

## Elliptical tapers
imtaper_maj=200 #arcsec
imtaper_min=100 #arcsec
#imtaper_pa=0 #deg  (A horizontally stretched taper in the image domain)
imtaper_pa=30 #deg  (A rotated Gaussian with PA=30deg)
uvtapers=['' , ['%3.2farcsec'%(imtaper_maj),
                '%3.2farcsec'%(imtaper_min) ,
                '%3.2fdeg'%(imtaper_pa) ] ,
              ['%3.2flambda'%(convert_lm_to_uv(imtaper_min)/2.0),
                '%3.2flambda'%(convert_lm_to_uv(imtaper_maj)/2.0) ,
                '%3.2fdeg'%(imtaper_pa+90.0) ] ]
## To get the uv-domain parameters, need to flip the major and minor axis,
## and rotate by 90deg.


print("\nSettings for uvtaper in tclean : \n\
[ None,  FWHM in the image domain, HWHM in the uv-domain] ")
print(uvtapers)

FWHM of 100.000 arcsec maps to a FWHM of 1820.374 lambda
FWHM of 200.000 arcsec maps to a FWHM of 910.187 lambda

Settings for uvtaper in tclean :
[ None,  FWHM in the image domain, HWHM in the uv-domain]
['', ['200.00arcsec', '100.00arcsec', '30.00deg'], ['910.19lambda', '455.09lambda', '120.00deg']]

run_im(imnames,uvtapers,weighting='natural')

uvt_orig : uvtaper =
Restoring Beam : 49.8538 arcsec  X  46.8010 arcsec ,  -87.6734 deg

uvt_taper_im : uvtaper = ['200.00arcsec', '100.00arcsec', '30.00deg']
Restoring Beam : 226.6692 arcsec  X  152.2435 arcsec ,  31.6937 deg

uvt_taper_uv : uvtaper = ['910.19lambda', '455.09lambda', '120.00deg']
Restoring Beam : 226.6705 arcsec  X  152.2433 arcsec ,  31.6936 deg

calc_convolve(theta_orig=50.0, theta_taper=100.0)
calc_convolve(theta_orig=50.0, theta_taper=200.0)

Convolution of FWHMs of 50.0000 arcsec and 100.0000 arcsec           yields 111.8034 arcsec
Convolution of FWHMs of 50.0000 arcsec and 200.0000 arcsec           yields 206.1553 arcsec

run_im(imnames,uvtapers,weighting='uniform')

uvt_orig : uvtaper =
Restoring Beam : 28.2215 arcsec  X  27.3716 arcsec ,  -11.8527 deg

uvt_taper_im : uvtaper = ['200.00arcsec', '100.00arcsec', '30.00deg']
Restoring Beam : 201.1375 arcsec  X  100.6616 arcsec ,  29.9237 deg

uvt_taper_uv : uvtaper = ['910.19lambda', '455.09lambda', '120.00deg']
Restoring Beam : 201.1391 arcsec  X  100.6613 arcsec ,  29.9237 deg

calc_convolve(theta_orig=28.0, theta_taper=100.0)
calc_convolve(theta_orig=28.0, theta_taper=200.0)

Convolution of FWHMs of 28.0000 arcsec and 100.0000 arcsec           yields 103.8460 arcsec
Convolution of FWHMs of 28.0000 arcsec and 200.0000 arcsec           yields 201.9505 arcsec

Conclusions

• In all cases, the resulting restoring beam is the same when specified either via the FWHM in the image domain, or as the HWHM in the uv-domain.

  • The difference in fitted restoring beam is \(< 0.01~arcsec\) for a pixel size of \(5.0~arcsec\) and beamwidths from 20 to 200 arcsec.

• Expected results are obtained for round versus elliptical beams, showing that the rotations are implemented accurately.

• Expected results are obtained when starting with a naturally weighted versus uniformly weighted PSF. This shows that the uvtaper is indeed being applied at the correct stage of weighting calculations.

• Note that the value of the resulting beam size is close to the theoretical calculation (with the understanding that the theoretical calculation is always only an approximation because the pre-existing uv-coverage can never be a true Gaussian).