"Proper" mosaicing

[KshitijT]

@gigjozsa's formula:

[KshitijT]

Not sure if "Linmos" task from Miriad is upto doing this for MeerKat; for one, MeerKat is (more or less obviously) one of the antennas which Miriad "recognizes" (see pbplot documentation, while this can be circumvented by manually putting beam parameters using "puthd" or something similar, one can only specify a Gaussian beam (parameters limited to beam FWHM) (see Primary beams in Miriad).

[gigjozsa]

Above formula describes how we can find the proper weight w_i for a given pixel, where i is the index of the image it is from. \sigma_i is the rms of the pixel in the ith image and depends on the primary beam A in that image. sigma_i scales like 1/A(r_i), where A(r_i) is the primary beam attenuation at the position of the pixel (at radius r). xi is the value of that pixel after primary beam correction (also scaling like 1/A) \widebar{x} is the weighted mean and using Gaussian error propagation we can calculate the corresponding standard deviation \sigma\widebar{x}. If we want to choose the weights such that the \sigma_\widebar{x} is minimized, the gradient with respect to the weights has to vanish, which is written down explicitly above. The solution is w_i = \frac{1}{\sigma_i}.

[bennahugo]

Yea a gaussian beam probably won't cut the mustard. How are we going to apply beam corrections? First need to build spectral deconvolution support into DDFacet? Since the beams rotate and the primary lobe is not circularly symmetric we should correct for it while imaging and not on the observation-averaged map. I've already added capabilities to it to handle MeerKAT duelcorr data and beams (assuming they are properly centred) into DDFacet, but currently the imager only supports various MFS modes. So we need to roll in a fast spectral deconvolution mode - hardest part will be prediction. The imager requires full resolution (or some way to interpolate) to full resolution when predicting in the major cycle.

Alternatively we can make snapshots, apply the beam in image space and reproject the snapshots before mosaicking.

[bennahugo]

Noting the montage steps here for reference of a quick and dirty mosaicing implementation - this is not correct as it simply weights the intrinsic maps by the square of the beam individually before addition.

> $ more mosaic.sh                                                                                    
#!/bin/bash

rm -f imglist
rm -f images.tbl
rm -f base.hdr
outdir="output"
imlist=("G330.89-0.36-MFS-image.fits.int.fits" "G331.59-0.36-MFS-image.fits.int.fits" "G332.31-0.36-MFS
-image.fits.int.fits" "G332.65+0.35-MFS-image.fits.int.fits" "G331.95+0.35-MFS-image.fits.int.fits" "G3
31.24+0.35-MFS-image.fits.int.fits")

echo "|                                               fname|" >> imglist
echo "|                                                char|" >> imglist
for i in ${imlist[*]}
do
  appmap="$(sed 's/.int.fits//g' <<< $i)"
  echo " $appmap" >> imglist
done

mImgtbl -t imglist ${outdir} images.tbl
mMakeHdr -p 4.444e-4 images.tbl base.hdr EQUJ 2000
for i in ${imlist[*]}
do
  noisemap="$(sed 's/.int.fits//g' <<< $i).noisescale.fits"
  mProjectPP -t 0.01 -w "${outdir}/${noisemap}" "${outdir}/$i" "${outdir}/${i}.projected.fits" base.hdr
  #mProject "${outdir}/$i" "${outdir}/${i}.projected.fits" base.hdr
done
rm -f imglist
echo "|                                               fname|" >> imglist
echo "|                                                char|" >> imglist
for i in ${imlist[*]}
do
  echo " $i.projected.fits" >> imglist
done
mImgtbl -t imglist ${outdir} imagesproj.tbl
echo "Overlap"
mOverlaps imagesproj.tbl diffs.tbl
echo "Diff"
mDiffExec -p ${outdir} diffs.tbl base.hdr ${outdir}
echo "Fit"
mFitExec diffs.tbl fits.tbl ${outdir}
echo "Background model"
mBgModel imagesproj.tbl fits.tbl corrections.tbl
echo "Background correct"
mBgExec -p ${outdir} imagesproj.tbl corrections.tbl ${outdir}
echo "Nan correct"
#for i in ${imlist[*]}
#do
#       mFixNaN -v 0.0 "$i.projected.fits" "$i.projected.fixed.fits"
#done
echo "Add"
mAdd -p ${outdir} imagesproj.tbl base.hdr "${outdir}/mosiac.fits"

It does an okay job, but the noise is not corrected. CASA supposedly implements this formula (applause to @gigjozsa for deriving this from first principals!) with a wide range of approximate responses see: https://casa.nrao.edu/docs/CasaRef/imager.linearmosaic.html https://casa.nrao.edu/docs/CasaRef/vpmanager-Tool.html

Using an AIRY disk response and running the mosaicer does not work - a blank image is returned.

I've just discovered there is another mosaicing tool in casa (wtf...), so I'll try this tomorrow: https://casa.nrao.edu/docs/CasaRef/linearmosaic.makemosaic.html

[IanHeywood]

Can you elaborate as to what's up with Montage? It's pretty much all I use for mosaicking so if it's doing something bogus I'd like to know!

Cheers.

[bennahugo]

@IanHeywood - not complete bogus but not quite the radio version as I understand it - the normalization by the sum of contributing beam weights should to be done per pixel after the addition step.

I also see patches where there are edge effects visible in the mosaics it produces, but maybe the projected images can be smoothed by convolved mask like is done in ddfacet

[IanHeywood]

I thought the -w switch took care of this. Are you seeing sharp boundaries where there's obviously a cut PB edge?

[bennahugo]

it is more a sub image edges - maybe that is just a smoothing issue - the leftover correlated noise on one map doesn't quite line up with the noise on the other map. I'm going to play around a little more with this

[IanHeywood]

Ah ok. What you're mosaicking and how you're going about it is probably a secret but you can tell me all about it one day. :)

I don't think I've ever done anything than stitch together entire PB-corrected images.

[paoloserra]

I've tried to mosaic MeerKATHI HI cubes in montage following this tutorial. A couple of comments:

• Montage fails if the cubes have a 4th axis (Stokes), which is the case for WSClean products. For this reason I've implemented a "remove 4th axis" option in MeerKATHI, see https://github.com/ska-sa/meerkathi/commit/2f2d92e38fabf6ef33bcc0252a107859536d1531.
• Montage fails if the cubes are on different channel grids. Maybe there are ways to make Montage regrid also along the frequency/velocity axis but I wasn't able to find them. For this reason I've implemented a new option in MeerKATHI that allows one to regrid to their preferred channelisation. This is done when Doppler correcting the .MS files, see https://github.com/ska-sa/meerkathi/commit/3787ed9727418855419c68242c0cc6c97cb95643. With this option, MeerKATHI can be used to make cubes on a same channel grid, and Montage will mosaic them properly.
• It isn't clear to me whether Montage is able to make the PB correction based on some input settings, or whether one needs to PB-correct first, and then use Montage. Just in case, I've implemented an option in MeerKATHI to create the necessary PB cube each time you create an HI cube, see https://github.com/ska-sa/meerkathi/commit/2f2d92e38fabf6ef33bcc0252a107859536d1531. For now this is Gaussian PB. Happy to have your advice on how to make it more realistic/sophisticated.

Given these new MeerKATHI options, my current approach to HI mosaicking is to:

• Use MeerKATHI to make HI cubes and PB cubes that are already on the same channel grid.
• Regrid the HI cubes and input PB cubes on the sky to a common mosaic header using Montage.
• Mosaic the regridded cubes in Python using the regridded PB cubes and applying the equation:

taken from this page.

The whole thing is a fairly short Python script that I'm happy to share if somebody is interested (once I've commented the code a little better though). The script produces also a noise cube and a weights cube to be used in SoFiA.

So far I've used this script to mosaic a small numbers of ~2 GB HI cubes with good results.

MeerKATHI is currently set up to process one pointing at a time. I don't think there is any urgency to change this. In my opinion it is OK to mosaic outside MeerKATHI for the time being.

[bennahugo]

To point out:

DDFacet can give you a joint deconvolution mosaic (similar to CASA tclean, but with correct intrinsic corrections). You just specify the fields as F* .

• image plane beam corrections need to be done per subband before stitching, otherwise the flux scale is biased over wider bandwidths, and a new MFS mosaic created, which takes into account the sigma of the MFS slice. This needs an extra set of scripts to implement correctly.
• Montage does a linear fit - it is built for optical data. This is not correct - if you look at the mosaics carefully there are clear stitching edges visible because the residual noise per pointing differs between images. The mosaic recipe I've given above is there for eyecandy purposes and not for actual science purposes.

[paoloserra]

I think we'll need to implement both approaches in the pipeline. For example, I don't expect DDFacet is able to deconvolve large HI cubes with 3D clean masks (not for the moment at least). So for those we'll probably need to image each pointing independently and then mosaic. The mosaic could be used to define a better 3D clean mask to be used to re-image each pointing. Plenty to do/discuss at the busyweek!

[bennahugo]

@paoloserra unless we provide you with an HI deconvolution mode.... plenty to discuss at the DDF workweek :)

You get a nicely stitched image like the example above with DDF. The exact parset options you need to set are: https://github.com/saopicc/DDFacet/blob/master/DDFacet/Parset/DefaultParset.cfg#L5 and https://github.com/saopicc/DDFacet/blob/master/DDFacet/Parset/DefaultParset.cfg#L76

Use version 0.4.0 or 0.4.1 on my dockerhub: bhugo/ddfacet

[paoloserra]

nice! yes, if you guys are interested in HI imaging we can definitely discuss what we need.

[KshitijT]

Tagging @svw26 , since mosaicking.

[gigjozsa]

With #720 we do now have a working mosaic mode. Proper mosaic would be DDFacet ing with joint deconvolution, but that is not an issue but a whole project, which we should define at the next busy week. Therefore I am closing this as an issue. No one will forget it.