Skip to content

Creating a MeerKAT DD corrected intrinsic flux image with DDF kMS

Jump to bottom
bennahugo edited this page Jul 9, 2021 · 1 revision

Prologue

DDF/kMS are toolboxes that you can use to do many things around DD/DI imaging and calibration... This tutorial is just one approach for one dataset already DI precalibrated and presents a way to conduct the DD calibration and imaging with DDFacet and kMS. We use a MeerKAT dataset with a 1283.8 MHz central frequency and 855.7 MHz bandwidth and calibrate in 13 direction, modeling flux spectral behaviours using 3rd order polynomial (in log space) and taking into account he MeerKAT's primary beam (described in FITS files - but other beam models are available). The DD calibration is made on top of this DD modeling, and corrects for additional effects (that are likely to be dominated by pointings errors in that case).

Major steps can go as follows:

  1. Create a first image/skymodel (this step includes defining the clustering, making mask, and deconvolving) that takes the primary beam model into account, but not the DD solutions from kMS that are to be estimated in the next step.
  2. Calibrate the Direction Dependent Jones in different directions
  3. Apply these in a final imaging step to create a final DD corrected image

Do not copy/paste the following blindly, adapt both the imaging parameters as well and the data chunking, NCPU etc to match your imaging/system requirements...

Branches that have been used for this processing:

DDF: TestAroundDeconv_ChangeParallel
kMS: InWeights

1- Making the clustered DI image/skymodel

This section name suggests that we don't apply DD Jones matrices, altough we will by using a time-frequency-direction dependent primary beam model. We refer to this image as a Direction Independent because on top of the DD beam model we do not use any kMS-estimated DD Jones matrices. These will be estimated in Sec. 2 to produce the DD image in Sec. 3.

Make a dirty image that we'll use to define a manual clustering (that will be defined once for all).

If you have multiple Measurement Sets, you can put then in an mslist.txt as

ms1.avg.ms
ms2.avg.ms
etc...

and after defining your image size, cell size and additional parameters (beam, chunking, etc), you can create a dirty image using

DDF.py --Data-MS <YourData.ms/txt> --Data-ColName CORRECTED_DATA --Output-Name image_DI --Image-Cell 1.5 --Image-NPix 6000 --Output-Mode Dirty --Facets-DiamMin 0.03 --Facets-DiamMax 0.3 --Parallel-NCPU 96 --Data-ChunkHours 2 --Freq-NBand 3 --Freq-NDegridBand 10 --Beam-Model FITS --Beam-NBand 20 --Beam-FITSParAngleIncDeg 0.5 --Beam-FITSFile='/SomePath/meerkat_pb_jones_cube_95channels_$(xy)_$(reim).fits' --Beam-CenterNorm 1 --Beam-FITSFeed xy --Beam-FITSFeedSwap 1 --Beam-ApplyPJones 1 --Beam-FlipVisibilityHands 1 --RIME-DecorrMode FT

Check the individual parameters... if you need help on a particular one, you can type DDF.py -h

This will create different images, the one you're interested in at that stage is image_DI.dirty.fits. Read the DDF.py -h to see what is what but the names should speak for themselves.

Create the ds9 ClusterFile manually

Here you will set the direction in which the calibration will be done in the later stage. We do this before any DI deconvolution, since the piecewise-constant beam applied in both kMS/DDF is better consistently defined over the recursion. This is mostly to save processing time. (We do this manually in this example but there are automatic ways to do this using MakeModel.py or MakeCatalog.py+ClusterCat.py.)

  • Open your image_DI.dirty.fits with ds9
  • Make sure edit > region is enabled
  • Put green circles (default region in ds9) where you want your tessels center (preferably where there are bright sources). Only the position of the circles matter here. The number of tessels/regions you can calibrate against depends on you telescope uv-coverage, FOV, sky content, etc. For meerKAT, we'll pick maybe 6-10 directions (for LOFAR/LoTSS we do 45)...
  • Save your region file as a .reg file with regions > save regions > ds9.reg (using default values for saving wcs etc)

Convert the .reg file to .npy

MakeModel.py --ds9PreClusterFile ds9.reg --NCluster=0 --BaseImageName image_DI

The --NCluster=0 will give you as many directions to solve for as ds9 regions you've saved

Do an initial deconvolution with SSD2/automasking/clustering

Now you can run the a preliminary deconvolution using SSD2 and automasking. SSD is magical and converges very quickly. SSD needs islands to be definable, so either --Mask-External or --Mask-Auto (or both) should be defined. The automasking is a little conservative but you can use at this stage.

DDF.py image_DI.parset --Output-Name image_DI_Clustered --Output-Mode Clean --Facets-CatNodes ds9.reg.ClusterCat.npy --Deconv-Mode SSD2 --SSD2-PolyFreqOrder 3 --Deconv-MaxMajorIter 2 --Mask-Auto 1

You can also check weither you're happy with how the sky is tesseled opening your image with ds9 and superposing the tessels to it (region > load regions > image_DI_Clustered.tessel.reg) and you will get something like

Make a deeper deconvolution with external masking

Having an extra step making an external mask is in general a good idea (even when pipelining) since the brightest source are deconvolved at this stage.

MakeMask.py --RestoredIm image_DI_Clustered.app.restored.fits --Th 5

The lower Th the more the physical emission will be decovonvolved/modeled by SSD in the later stages, but the more there will be artifacts too. And the less artifacts the better for the calibration of course (since it will think you mean these artifacts are real astrophysical objects)

So check that you're happy with the mask:

dsm.py image_DI_Clustered.app.restored.fits.mask.fits image_DI_Clustered.app.restored.fits

(dsm.py is simply a wrapper to ds9 that matching the wcs, colorscale, colorbar, and setting the scale time something sensible based on the measured rms in the image)

Finalise the deconvolution of the DI image

Do a final major cycle (explicitly saying DDFacet to continue from the residual cache/image), using that ideal complete and robust mask and disabling automasking:

DDF.py image_DI_Clustered.parset --Output-Name image_DI_Clustered.DeeperDeconv --Predict-InitDicoModel image_DI_Clustered.DicoModel --Cache-Reset 0 --Cache-Dirty forceresidual --Cache-PSF force --Deconv-MaxMajorIter 1 --Mask-Auto 0 --Mask-External image_DI_Clustered.app.restored.fits.mask.fits

You image image_DI_Clustered.DeeperDeconv.int.restored.fits is the restored image with intrinsic fluxes. The skymodel itself is described in the .DicoModel file.

2- DDE calibration

Now tell kMS to derive the jones matrices. Type kMS.py -h if you want to know more about the kMS paramaters. Many of the things it needs are contained inside image_DI_Clustered.DeeperDeconv products (inclusing the skymodel and clustering) but you can override the tesselation, maximum facet size, decorrelation mode etc... Here we go simple, we just configure the beam (mirroring what we've set in DDF) and the solver to KAFCA (the Kalman filter), and the solution time and freq intervals (set by dt and NChanSols).

kMS.py --MSName <YourData.ms/txt> --FieldID 0 --SolverType KAFCA --PolMode Scalar --BaseImageName image_DI_Clustered.DeeperDeconv --dt 5 --NCPU 96 --OutSolsName DD0 --NChanSols 20 --InCol CORRECTED_DATA --TChunk 2 --BeamModel FITS --FITSParAngleIncDeg 0.5 --FITSFile='/SomePath/meerkat_pb_jones_cube_95channels_$(xy)_$(reim).fits' --CenterNorm 1 --FITSFeed xy --FITSFeedSwap 1 --ApplyPJones 1 --FlipVisibilityHands 1 --NChanBeamPerMS 20

You can inspect your solutions if you need to using

PlotSols.py --SolsFile <Your>.ms/killMS.DD0.sols.npz

3- Final DDE imaging

Now you have your DDE solutions, just run a single major iteration initialising the skymodel with the one you've used for calibration. The dirty image will the residual image you obtained in the last imaging step (but DD-corrected). So to update that residual and the associated skymodel and restored images (adding .AP to the name for "Amplitude" and "Phase" since we apply both):

DDF.py image_DI_Clustered.DeeperDeconv.parset --Output-Name image_DI_Clustered.DeeperDeconv.AP --Cache-Reset 1 --Cache-PSF auto --Cache-Dirty auto --Predict-InitDicoModel image_DI_Clustered.DeeperDeconv.DicoModel --DDESolutions-DDSols DD0 --Beam-Smooth 1 --Weight-ColName [WEIGHT_SPECTRUM,IMAGING_WEIGHT]

The artifacts in the .dirty.fits (a residual since you used a non-zero --Predict-InitDicoModel) should look better than the .residual.fits of the previous run (and this is the first thing you should check - use dsm.py call as above to compare these images). Since you apply amplitude and phase of the DD Jones estimated in the previous step, sources can become more focused and or morphologically different, so positive and negatives can be present in this dirty image. The SSD minor cycle(s) will update/correct the sky model to better estimates, and your DD app/int should look better than the DI ones... You can now compare the DI to the DD images such as the intrinsic restored image: dsm.py image_DI_Clustered.DeeperDeconv.int.restored.fits image_DI_Clustered.DeeperDeconv.AP.int.restored.fits

Many things can be explored to improve this image: DD-selfcalibration (calling kMS/DDF once more and using the current sky model), the sky clustering, the polarisation mode of the DD solver, the solution intervals time and frequency intervals, running further steps of high resolution full polarisation DI calibration based on the current DD-visibility modeling, etc etc...

Here the comparison between redisual images at some key steps above

and the comparison between the DI and DD apparant flux restored images

Clone this wiki locally