MeerKAT Simulation

Here we show the stimela package can be used to write a simulation pipeline. The aim of the pipeline is to simulate a MeerKAT observation, and then image the simulated data.

Create and open a Python file using your preferred text editor (my file is named simulation_pipeline.py). First, import the stimela package.

import stimela

Next, we define the I/O flow of the pipeline. That is, the input and output directories and files. The pipeline will expect to find all input products in INPUT and all will dump all the output products to OUTPUT. Note that the I/O flow for each task (container) can be defined/modified at the stage when any given task is added to the Recipe instance, it is good practice to define I/O flow globally.

 INPUT = "input"
 OUTPUT = "output"
 PREFIX = "stimela-example"

The visibility data that we deal with in radio interferometry is usually stored as a {\it measurement set} (MS), and this MS is usually both the input and output product for many tasks. Therefore, it is treated differently from the standard input and output directories defined above. Now set the location of the MS file(s), as well as the name of the MS that will be simulated.

 MSDIR = "msdir"
 MSNAME = "meerkat_example_simulation.ms"

Next, we set the sky model that will be simulated. In this case, the sky model is catalogue of sources from the NVSS survey This sky model should be placed in the input folder.

 SKYMODEL = "nvss1deg.lsm.html"

We can now open a new Recipe instance.

 recipe = stimela.Recipe("Stimela simulation example", ms_dir=MSDIR)

Note that each Recipe is associated to a MS directory, and all MSs that will be simulated or processed in a given recipe should be place in this directory.

With our Recipe instance started, we can proceed to add tasks to the recipe. First, we create an simulated MeerKAT MS using the simms executor container.

 simms_params = {
    "msname"       :	MSNAME,
    "telescope"	   : 	"meerkat",
    "synthesis"	   :	4 # Synthesis time in Hours
    "dtime" :	5 # integration time in seconds
    "freq0"        :	"1400MHz" # Start frequency
    "dfreq"        :   	"1MHz",
    "nchan"        :	10,
 }
 
 recipe.add("cab/simms",	# this is the full name of the executor image
    "create_emmpty_ms",		# container label
    simms_params, 	# These are the options set above
    input=INPUT,
    output=OUTPUT,
    label="Create empty MS")	# Task label. For logging purposes

The MeerKAT antenna layout as well as the uv-coverage of the simulated simulated observation are shown below

The next step is to generate simulated visibilities from the sky model. We use the simulator executor image, which uses MeqTrees, to achieve this.

 simulator_params = {
    "msname"	:	MSNAME,
    "skymodel"	: 	SKYMODEL,
    "column"	:	"CORRECTED_DATA",
 }
 
 recipe.add("cab/simulator",
    "simulate_visibilities",
    simulator_params,
    input=INPUT,
    output=OUTPUT,
    label="simulate visibilities into MS")

Finally, we image the data. To make things more interesting, let us image the data using Briggs uv-weighting with three robustness parameters, -2,0,2

 imager_params = {
    "msname"	:	MSNAME,
    "npix"	: 	2048,
    "cellsize"	:	2 # in arcseconds,
    "column"	:	"CORRECTED_DATA",
    "clean_iterations"	:	3000,
 }
 
 for robust in -2,0,2:
 
    imager_params["weight"] = "briggs %f"%robust
    image_params["prefix"] = "%s-%d"%(prefix, robust)
    
    
    recipe.add("cab/wsclean",
	  "image_data_robust_%d"%robust,
	  image_params,
	  input=INPUT,
	  output=OUTPUT,
	  label="Image data. Robust=%d"%robust)

With our recipe fully described, we now execute it.

 recipe.run()

Also, one can also run manually giving those steps to the Recipe.run() function via the steps argument. For example, to only create an empty MS (step 1), and then generate simulated visibilities (step 2) one would run

 recipe.run(steps=[1, 2])

The script can then be executed on command line as

 stimela run <name of script>

The output images and the logfile can be found in the output directory.