Revised workflow using new wsclean options (and some musings on residuals)

[IanHeywood]

A recent version (3.1.0 (2022-04-01)) of wsclean has some options that may simplify the solar imaging workflow. To recap, the plan so far has been:

1. Reference calibration / flagging / selfcal (likely using an automated pipeline of some kind)

2. Deconvolution of the self-calibrated visibilities in CORRECTED_DATA to construct a sky model (wsclean)

3. Inversion of the sky model into visibilities in the MODEL_DATA column (wsclean)

4. Subtraction of MODEL_DATA from CORRECTED_DATA to leave residual visibilities (probably a custom pyrap tool?)

5. Determine the RA and Dec of the Sun at the time centroid of each scan (pyrap / astropy / CASA measures?)

6. Split the Measurement Set into individual scans (CASA?)

7. Rephase each Measurement Set to the RA and Dec of the Sun specific to that scan (chgcentre)

8. Image the rephased Measurement Set for each scan to make solar images (wsclean)

The addition of the -shift parameter and the -subtract-model parameter mean that steps 4 and 7 are likely no longer necessary, and this should simplify the workflow considerably:

1. Reference calibration / flagging / selfcal (likely using an automated pipeline of some kind)

2. Deconvolution of the self-calibrated visibilities in CORRECTED_DATA to construct a sky model (wsclean)

3. Inversion of the sky model into visibilities in the MODEL_DATA column (wsclean)

4. Determine the RA and Dec of the Sun at the time centroid of each scan (pyrap / astropy / CASA measures?)

5. Split the Measurement Set into individual scans (CASA?)

6. Image the rephased Measurement Set for each scan to make solar images, with -shift set to the RA and Dec of the Sun for that scan, and -subtract-model enabled (wsclean)

(EDIT: Corrected some typos, formatting)

[IanHeywood]

On the subject of point 2, you can determine how 'complete' your model is by comparing the residual image (*-MFS-residual.fits) to the restored image (*-MFS-image.fits). The sky model is constructed by iterative deconvolution of the target field (the original pointing direction of the MS), and as such generally will not contain the Sun. But we do want to subtract the target field from the data, as when we rephase the data to image the Sun, strong sources in the target field might leak sidelobes into the solar image.

If there is significant structure or lots of bright sources that remain in your residual image, then those sources have not been deconvolved. Consequently they will not be in the model images, and therefore will not be in the visibility model that gets predicted from those models.

Ensuring that the model is as complete as possible is important not only for this use case, but also for self-calibration. Although the sky model will never be entirely complete, a model that is as complete as possible will greatly benefit the self-cal process.

Subtracting the residual image from the restored image will result in the model convolved with the restoring beam which is a good way to visually check the quality of the model. Another thing to watch out for in addition to incompleteness is the inclusion of spurious features in the model that are not part of the true sky (e.g. PSF-like artefacts that do not deconvolve properly).

[IanHeywood]

Here is an example where I am blinking the restored and residual images for a MeerKAT field. You can see most of the sources have been captured by the deconvolution, and do not feature in the residual image. The residual is dominated by PSF-like artefacts that are associated with the stronger sources in the field. Calibration deficiencies (primarily due to direction dependent effects) mean that these off axis sources no longer have the ideal PSF associated with them, and that fundamental assumption of the clean algorithm doesn't apply. The sidelobe structures do not deconvolve and thus remain in the residual.

The other thing one will typically find in the residual even when things have gone entirely to plan are low level sources at the ~few sigma level and below. These are hard to pick up in an automated way because of the risk of also including false positives due to high-sigma noise events in the image (source finding software has the same limitation).

[IanHeywood]

A closer look at the double-lobed radio galaxy at the top of the image shows that the mask used for the deconvolution here has generally performed well. Almost all of the extended structure is captured, and deconvolved leading to a fairly uniform flat residual within the mask. You can see from the profile plots on the right that the clean termination threshold here was set to 4 uJy/beam, which is probably about 1.5 sigma or so, bringing the typical residual plateau close to the general noise level.

Numerous peaks in the image remain, those are the low-level, few-sigma sources mentioned above. Some spokes also remain around the core of this radio galaxy, which is good, as these are calibration artefacts and we don't want them in the model.

[Victoria-Samboco]

Thank you Prof @IanHeywood for the notes I will very useful for me. I thing I don't need anymore to do a summary of the meeting because is basically everything summarised here. I will read it carefully and if I have some question by the way I will let you know.