Radio interferometry is amazing. Nobel prize winning amazing since it turns the weak floppy photons of the radio spectrum into powerhouses of resolution. However, there remains a sad fact of life that interferometers by themselves do not measure all the information in an image. Specifically, interferometers measure Fourier components of the image. But they are poorly suited to actually measuring Fourier components with small wave numbers ($\mathbf{k} \sim 0$). These spatial frequencies are typically measured by single dish telescopes and then… somehow… stitched together into single images. In theory, this is a well understood process, but there is a complication. In addition to only measuring large wave numbers, interferometers also only measure certain parts of the Fourier plane and so we need to… somehow… interpolate a complex function into places where it is not measured.

People gleefully point out fun facts like interferometers measure zero net flux by construction and that there are an infinite number of possible solutions for the deconvolution process. However, radio astronomers can stare at poorly deconvolved images and think “that isn’t right.” We have a great prior (more like a prejudice) on what the sky looks like and so we accept families of solutions consistent with our prejudice. To manage that, we tend to like minimum information deconvolution.

I’m trying to deconvolve a radio image from the VLA. The VLA is an utterly amazing instrument and it produces huge amounts of data. The supported software for the deconvolution is called CASA and the deconvolution and imaging of the data is handled through the clean or the affectionately named tclean where the t means test. The data I’m used are large and I have to produce large images (100 field mosaics, $2500\times 2500$ pixel images). This means that running clean once, for a single frequency channel is a multi-hour affair where I conduct a slow ballet around the various pitfalls of tclean to discover how to produces beautiful images. I’m writing now since I see some progress.

Here’s the full clean call that seems to be non-terrible:

tclean(vis=flist,specmode='cube',imagename='NH3_11',
       start="23.693764GHz", nchan=10, restfreq="23.69450GHz",
       width='10.416kHz', imsize=[2500,2500],
       gridder='mosaic', restoringbeam='common',
       cell=['0.5arcsec'], pbcor=True, outframe='LSRK',
       niter=100000, chanchunks=-1,
       weighting='briggs', robust=0.5,
       deconvolver='hogbom', 
       threshold='0.8mJy', mask='nh3_11_gas.mask', calcres=False,
       calcpsf=False, restart=True, 
       phasecenter='J2000 3h28m59.474s 31d16m11.85s')

First, I was seeing tons of streaking in the brightest channels of the image. At first, I worried these were RFI, but they showed up in every observation of the field. It was something wrong with the imager. I tried all kinds of stabilization tricks for tclean by restricting deconvolution to specific regions via a mask and feeding in model images which are an initial guess for what the true sky looks like before the effects of the interferometer are applied. What ended up working was simple: I was using the default natural weighting, which turns out to be non-ideal for imaging wide areas with lots of structure since it makes the beam large and blends together emission. Switching this to briggs weighting with robust=0.5 calmed down the deconvolution significantly.

Second, multiscale just isn’t nice and stable, even with heavy masking or even modeling. I have the excellent fortune of having single dish data to pair with the VLA so it’s not critical to get plausible information on large scales. Going with Clean Classic (now with 50% fewer calories) seems to be appropriate and thus the deconvolver=hogbom call.

I experimented with model images. I’m not sure what exactly is expected. However, it seems to want an image on exactly the same coordinate grid and it looks like it must be in Jy/pixel. The “exactly the same coordinate grid” seems to be hard because I lost a day on trying to make a model work that had completely indistinguishable coordinate information. I even went so far as to write over all the coordinate metadata in the model image with the template image. No luck. Still not on the same grid.

The general strategy of using a niter=0 call to initialize a template is great for messing around with masking, etc.

The MPI interface does not seem to work on my system yet. It merely spreads the crashes over multiple nodes. Forcing parallelization by splitting over frequency planes seems like the way to go.

Looking over it, the call is rather simple after a great deal of experimentation. It just seems that none of the advanced features produced stable imaging results.

I imagine a lot of things would go better with ALMA data rather than the VLA since there is less azimuthal structure in the PSF given it is a random array.

Note we also found out that the restoringbeam='common' doesn’t work properly either, so that’s just in there in the “wishful thinking category.” Expect an exciting blog post about ellipses soon!