Imaging solar system objects simultaneously with the main field

[landmanbester]

I have made a start towards this in the solarkat branch where I am trying to test this out on a slightly troublesome observation containing the sun. The idea is to do this "on the fly" without all the kludgy rephasing and copying that is currently required. You can already tell the init worker to apply QuartiCal gains and split your MS up into an arbitrary number of images per scan (controlled by the --integrations-per-image and --nband parameters). This will split your observation into a number of (possibly on the fly averaged) intermediary data products that can then be ingested by the grid worker. This worker will check the mean time of each dataset and locate the object (in this case the sun) using astropy (thanks to whoever put this together) and then, instead of rephasing, it finds the (l,m) coordinates of the object and simply instructs the gridder to image the patch around it. Here is the result of using the chgcentre tool and imaging the result

Here is the result doing it the on the fly way

They look pretty similar to me. If you blink between them it almost looks like there is a slight rotation between the images but I believe that is a projection effect due to the fact that they are reconstructed on different tangent planes (of course the PSF also needs to be computed differently when the center of the image does not correspond to (l=0, m=0)). There are a couple of oddities in getting to this result. In particular, I had to negate the signs of the lm coordinates computed with africanus for the images to match up but this is perhaps no more strange than the following message from chgcentre

-flipuvwsign
        Flips the UVW sign. Necessary for LOFAR, for unknown reasons.

@o-smirnov @IanHeywood @bennahugo @JSKenyon what do you think?

The next step is to try and restructure things so that you can deconvolve the object and the main field simultaneously. As a first attempt, borrowing ideas from SSD in DDF, it is probably easiest to run the deconvolution worker on the different "fields" separately and then just combine and subtract the visibilities correctly using separate specialised workers (which we can of course string together with stimela). However, I eventually want to create a consistent framework for this. I would appreciate comments, especially on what features would be useful and things that could possibly bite me when I try to implement them in the future

[landmanbester]

I have tried to deconvolve the above image but I am not getting very far. Here are the residual and model after cleaning

Cleaning is seriously slow with next to no progress being made in the later cycles

Here are a number of things to consider:

• negating uvw is not the same as negating lm because of the w-term. Negating uvw results in very smeared out images of the sun (well probably any source that is far enough from the phase center).
• the PSF corresponding to a "field" far away from the phase center should not be normalised to unity as dirty/residual images are normalised by wsum regardless of where in the sky they may be.
• the source can have both intrinsic as well as apparent time variability because of its relative motion with respect to the phase center. We may be able to do something about the latter by rephasing each unique time but it is probably better to add a time axis during deconvolution because that also takes care of the former.
• these sources will typically be in the sidelobes of the beam so, unless we specify a beam model, they can have essentially arbitrary spectra.
• We should not be overly concerned with a physically realistic reconstruction since we just want to remove these sources from the main field.

The last three points have given me an idea which I shall be testing shortly

[landmanbester]

screen-capture.webm

Animation over frequency axis of the dirty image

[IanHeywood]

these sources will typically be in the sidelobes of the beam so, unless we specify a beam model, they can have essentially arbitrary spectra.

I think any beam model that far out in the sidelobes will probably not be any good. I don't think it'll be worth the overheads.

We should not be overly concerned with a physically realistic reconstruction since we just want to remove these sources from the main field.

100% agree.

[IanHeywood]

Can't see the movie alas...

[landmanbester]

Ah, you may need the screen recorder plugin for chrome. Let me see if I can get it into a more friendly format

[IanHeywood]

Firefox user, sorry. 😅

[landmanbester]

freq_movie.mp4

How about now?

[IanHeywood]

Yeah that works, thanks.

[landmanbester]

I can get a deeper by increasing the frequency resolution to 16 bands and using the flux mop

but still not quite there (I suspect the remaining time variability might limit how deep we can go). It's also probable that some of the sidelobes from the main field have been absorbed into the model here. Next I will impose a mask and try to do the outer major cycle

[landmanbester]

Subtracting the above model and imaging the result gives the following

Clearly the subtraction did not go very well. Zooming in on the sun itself gives

which confirms that some level of subtraction did take place (note this is different from the previous residual because I used natural weighting). I'm not sure why the sidelobes in the ring around the sun are so pronounced. Possible reasons include:

1. The model of the sun obtained from deconvolving it separately has absorbed sidelobes from the main field and/or residual time variablity. In this case it is not clear why the residual of the sun itself looks ok but I now have the machinery to test this.
2. Calibration problems that are not properly accounted for during deconvolution and subsequent subtraction. Easy enough to test by peeling the source instead of doing simple subtraction. Leakage is probably quite pronounced this far away from the pointing direction.
3. There are low level structures outside the box I drew around the sun that hasn't been deconvolved. This is probably quite unlikely. Comments?
4. I haven't got proper tests for the new functionality so of course there can also be in bug in the machinery. At least my sanity tests check out.

I'm redoing the large image with the same weighting I used during deconvolution, hopefully that will narrow things down a bit

[landmanbester]

Ah, this is quite telling

Left is natural weighting, right is Briggs 0 (i.e. the same I used for deconvolution). Does this indicate problems with short spacings? I'm not sure. Gonna try and peel that source and see what happens

[IanHeywood]

It's probably just because MeerKAT's naturally-weighted PSF is a complete horrorshow. The robust 0.0 map looks much improved to me, although a closer-look before and after of the target field will be informative.

[o-smirnov]

Yeah I wouldn't worry about the naturally weighted residuals for now and just forge ahead... zoom in on the target before subtraction, after subtraction, and after subtraction and peeling...

[landmanbester]

Took a while to do the comparison because I accidentally deleted the original image. Here they are side by side (left=before, right=after)

And with a bit more zoom

It's not perfect but definitely an improvement. It should now be possible to clean the main field and peel the sun. I'm a bit weary of peeling the sun without a model for the main field, feels like things could go very wrong. Not sure how worried we should be about the Sgr A* complex. In principle we could make a separate model for that too, I would just need to add it as one of the target options (or allow --target option to take a fixed (ra,dec) since I'm guessing we can treat it as a fixed body?)

[IanHeywood]

Definitely an improvement. I think there should definitely be a model for the rest of the sky in place before peeling. The usual approach is to deconvolve everything, then partition that clean component model up into the problem source and the rest of the sky, then create two model visibility columns.

How did your original subtraction above work with pfb-clean's treatment of the Sun as a moving target? Does it generate residual visibilities on the fly? If the model doesn't have a time axis then predicting it will smear it back to the sidereal tracking case and we might lose the advantage of your solar system object tracking magic.

The easiest way to deal with Sgr A would be to put a small outlier facet on it and just deconvolve it as a fixed body along with everything else. For wsclean this would require either a huge image or some re-phasing and iterative imaging. Whether it also needs peeling or just deconvolving is TBD I guess...

[landmanbester]

How did your original subtraction above work with pfb-clean's treatment of the Sun as a moving target? Does it generate residual visibilities on the fly? If the model doesn't have a time axis then predicting it will smear it back to the sidereal tracking case and we might lose the advantage of your solar system object tracking magic.

Actually in this case I just did the full scan but in principle I can now do this (I think) for time variability within a scan. The model will have a time axis in that case (i.e. the average center of the image at that time step) which I can feed that into the wgridder so that it can degrid the image at every time step to the relavant rows in the MS. This column can be (approximately) correctly subtracted from data and then I just image the result. Does that sound sensible? The functionality I am currently lacking is to run clean on a model with a time axis but this shouldn't be hard if they can be deconvolved independently. I'm trying to get a recipe working for the meeting tomorrow, maybe run all this by you guys then?

[landmanbester]

I got slightly better results after predicting the model for the sun at a higher frequency resolution and using a bigger fraction of the "mopped" residual. Here is a before and after of the full field

and zooming in on the target

The recipe to do the above with pfb-clean and QuartiCal currently lives here (up to subtraction, peeling steps yet to be tested). Hopefully we can agree that it is a much simpler and less error prone procedure. Both the gridding and degridding steps now support time variable models. I think all that is needed is to add support for the clean worker to be able to cope with a time axis. As long as we are cleaning independently, this should not be difficult to add using the client interface of the distributed scheduler

[landmanbester]

And now it is possible to pass arbitrary coordinates to --target where they just need to be in a sensible format that SkyCoord can understand. As an example, passing in --target '17:45:40.04,-29:0:28.12' gives the following

I bet at least one person on this thread can tell me what that is ;-) One of the really cool aspects of doing things this way is that all the different fields, since they are stored in separate data sets and deconvolved separately, can have different resolutions. This should work as long as they are not too close together on the sky. I still have to figure out how this would work if you have significant sidelobe overlap between fields (I'm thinking about Tarantula of course)

[IanHeywood]

This is brilliant.