Self-Calibration Demo with ASKAP RACS Data
This notebook demonstrates the self-calibration using ASKAPsoft. A short 5-minute scan from the ASKAP made for astrofest is used in this notebook. The data is available in the CASDA archive, under scheduling block id SB45305.
We will be using beam 0 visibility data of the scan and will be using the same parsets used by ASKAP SDP as reference, reducing the image size to be able to process the notebook on a simple laptop.
This notebook should take approximately 1+ hour and atleast a machine with 2 CPUs to run
Prerequisites
[10]:
# A function to plot FITS
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors
from astropy.io import fits
from astropy.wcs import WCS
%matplotlib inline
def plot_image(file, norm=None):
# Using one of the files to setup WCS projection
hdu = fits.open(file)[0]
wcs = WCS(hdu.header, naxis=2)
fig, ax = plt.subplots(1, 1, figsize=(6, 6), subplot_kw=dict(projection=wcs))
print(fits.info(file))
hdu = fits.open(file)[0]
if norm:
im = ax.imshow(hdu.data[0,0,:,:], origin='lower', cmap='viridis', norm=mcolors.PowerNorm(norm))
else:
im = ax.imshow(hdu.data[0,0,:,:], origin='lower', cmap='viridis')
ax.set_xlabel('RA')
ax.set_ylabel('Dec')
ax.set_autoscale_on(False)
ax.set_title(file)
# plt.colorbar(im, ax=ax)
plt.show()
Self calibration
[1]:
msname = "./scienceData.RACS_0012+00.SB45305.RACS_0012+00.beam00_averaged_cal.leakage.ms/"
obsname = "SB45305"
imagesize = [2048, 2048]
[2]:
from taitale import taitale_env
env = taitale_env(
image = "csirocass/askapsoft",
tag = "1.17.6-openmpi4"
)
The tools we would be using in this tutorial are 1. imager - The ASKAP imager 3. mslist - To get info from the measurement set 4. ccalibrator - To predict
[3]:
from taitale.askap import imager, mslist, ccalibrator, cmodel
from casacore.images import image
[4]:
mslist(env=env, args=f"--full {msname}", logfile=f"mslist_{obsname}.txt")
Starting mslist
mslist complete
[5]:
!cat mslist_SB45305.txt
Debug: registered context Global=0
2025-03-19 02:12:22 INFO Observer: Project:
2025-03-19 02:12:22 INFO + Observation: ASKAP(36 antennas)
2025-03-19 02:12:22 INFO MSMetaData::_computeScanAndSubScanProperties Computing scan and subscan properties...
2025-03-19 02:12:22 INFO Data records: 21312 Total elapsed time = 318.505 seconds
2025-03-19 02:12:22 INFO + Observed from 29-Oct-2022/14:11:38.4 to 29-Oct-2022/14:16:56.9 (UTC)
2025-03-19 02:12:22 INFO
2025-03-19 02:12:22 INFO Fields: 1
2025-03-19 02:12:22 INFO + ID Code Name RA Decl Epoch nRows
2025-03-19 02:12:22 INFO + 0 RACS_0012+00 00:12:24.827000 +00.00.00.00000 J2000 21312
2025-03-19 02:12:22 INFO Spectral Windows: (1 unique spectral windows and 1 unique polarization setups)
2025-03-19 02:12:22 INFO + SpwID Name #Chans Frame Ch0(MHz) ChanWid(kHz) TotBW(kHz) CtrFreq(MHz) Corrs
2025-03-19 02:12:22 INFO + 0 NO_NAME 288 TOPO 743.991 1000.000 288000.0 887.4907 XX XY YX YY
2025-03-19 02:12:22 INFO The SOURCE table is absent: see the FIELD table
2025-03-19 02:12:22 INFO + Antennas: 36 'name'='station'
2025-03-19 02:12:22 INFO + ID= 0-3: 'ak01'='ASKAP', 'ak02'='ASKAP', 'ak03'='ASKAP', 'ak04'='ASKAP',
2025-03-19 02:12:22 INFO + ID= 4-7: 'ak05'='ASKAP', 'ak06'='ASKAP', 'ak07'='ASKAP', 'ak08'='ASKAP',
2025-03-19 02:12:22 INFO + ID= 8-11: 'ak09'='ASKAP', 'ak10'='ASKAP', 'ak11'='ASKAP', 'ak12'='ASKAP',
2025-03-19 02:12:22 INFO + ID= 12-15: 'ak13'='ASKAP', 'ak14'='ASKAP', 'ak15'='ASKAP', 'ak16'='ASKAP',
2025-03-19 02:12:22 INFO + ID= 16-19: 'ak17'='ASKAP', 'ak18'='ASKAP', 'ak19'='ASKAP', 'ak20'='ASKAP',
2025-03-19 02:12:22 INFO + ID= 20-23: 'ak21'='ASKAP', 'ak22'='ASKAP', 'ak23'='ASKAP', 'ak24'='ASKAP',
2025-03-19 02:12:22 INFO + ID= 24-27: 'ak25'='ASKAP', 'ak26'='ASKAP', 'ak27'='ASKAP', 'ak28'='ASKAP',
2025-03-19 02:12:22 INFO + ID= 28-31: 'ak29'='ASKAP', 'ak30'='ASKAP', 'ak31'='ASKAP', 'ak32'='ASKAP',
2025-03-19 02:12:22 INFO + ID= 32-35: 'ak33'='ASKAP', 'ak34'='ASKAP', 'ak35'='ASKAP', 'ak36'='ASKAP'
2025-03-19 02:12:22 INFO The HISTORY table is empty
[6]:
frequency="887.4907MHz"
increment="1000.000kHz"
Loop 0
Let’s start with a template parset
[7]:
!cat ./SB45305/continuum_image_selfcal_F00_B00_LOOP0_14796240.in
##########
## Continuum imaging with cimager
##
#Standard Parameter set for Cimager
Cimager.dataset = scienceData.RACS_0012+00.SB45305.RACS_0012+00.beam00_averaged.ms
Cimager.datacolumn = DATA
Cimager.imagetype = fits
#
# No Channels selection performed
#
Cimager.Images.Names = [image.i.RACS_0012+00.SB45305.cont.RACS_0012+00.beam00]
Cimager.Images.shape = [6144, 6144]
Cimager.Images.cellsize = [2.5arcsec, 2.5arcsec]
Cimager.Images.image.i.RACS_0012+00.SB45305.cont.RACS_0012+00.beam00.direction = [00h10m18.828, +00.31.29.93, J2000]
# We will let the code derive the Nyquist-sampled griddingcellsize
Cimager.Images.nyquistgridding = true
# This is how many channels to write to the image - just a single one for continuum
Cimager.Images.image.i.RACS_0012+00.SB45305.cont.RACS_0012+00.beam00.nchan = 1
#
# The following are needed for MFS clean
# This one defines the number of Taylor terms
Cimager.Images.image.i.RACS_0012+00.SB45305.cont.RACS_0012+00.beam00.nterms = 2
# This one assigns one worker for each of the Taylor terms
Cimager.nworkergroups = 3
# Leave 'Cimager.visweights' to be determined by Cimager, based on nterms
# Leave 'Cimager.visweights.MFS.reffreq' to be determined by Cimager
#
# Options for the alternate imager are not required
#
# This defines the parameters for the gridding.
Cimager.gridder.snapshotimaging = false
Cimager.gridder.snapshotimaging.wtolerance = 2600
Cimager.gridder.snapshotimaging.longtrack = true
Cimager.gridder.snapshotimaging.clipping = 0.01
Cimager.gridder = WProject
Cimager.gridder.WProject.wmax = 26000
Cimager.gridder.WProject.nwplanes = 513
Cimager.gridder.WProject.oversample = 5
Cimager.gridder.WProject.maxsupport = 1024
Cimager.gridder.WProject.variablesupport = true
Cimager.gridder.WProject.offsetsupport = true
Cimager.gridder.WProject.sharecf = true
#
# These parameters define the clean algorithm
Cimager.solver = Clean
Cimager.solver.Clean.solutiontype = MAXBASE
Cimager.solver.Clean.verbose = false
Cimager.solver.Clean.decoupled = true
Cimager.solver.Clean.tolerance = 0.01
Cimager.solver.Clean.weightcutoff = zero
Cimager.solver.Clean.weightcutoff.clean = false
Cimager.threshold.masking = 0.9
Cimager.solver.Clean.logevery = 50
Cimager.solver.Clean.detectdivergence = true
Cimager.Images.writeAtMajorCycle = false
#
Cimager.restore = true
Cimager.restore.beam = fit
Cimager.restore.beam.cutoff = 0.5
#
# Cleaning parameters set for loop 0
Cimager.threshold.majorcycle = 0.75mJy
Cimager.ncycles = 3
Cimager.solver.Clean.algorithm = BasisfunctionMFS
Cimager.solver.Clean.niter = 200
Cimager.threshold.minorcycle = [40%,25.0mJy,0.5mJy]
Cimager.solver.Clean.gain = 0.3
Cimager.solver.Clean.psfwidth = 256
Cimager.solver.Clean.scales = [0]
#
# Preconditioning parameters set for loop 0
Cimager.preconditioner.Names = [Wiener]
Cimager.preconditioner.preservecf = true
Cimager.preconditioner.Wiener.robustness = 0.0
Cimager.MinUV = 100
# Output data selection parameters for continuum imaging
Cimager.write.residualimage = true
Cimager.write.psfrawimage = false
Cimager.write.psfimage = false
Cimager.write.weightsimage = false
Cimager.write.modelimage = true
Cimager.write.maskimage = false
Cimager.write.sensitivityim = false
Cimager.write.firstrestore = true
Cimager.write.weightslog = true
#
# No self-calibration as it is the first time around the loop
Cimager.calibrate = false
#
[8]:
source_name = "RACS_0012+00"
image_name = f"image.i.{source_name}"
nchan = 288
ra="00:12:24.827000"
dec="+00.00.00.00000"
direction=f"[{ra}, {dec}, J2000]"
imaging_workers = 1
nchanpercore = nchan // imaging_workers
print(f"channels: {nchan}")
print(f"nchanpercore: {nchanpercore}")
print(f"mpi_workers: {imaging_workers}")
channels: 288
nchanpercore: 288
mpi_workers: 1
[9]:
imager(
env=env,
workers=imaging_workers+1,
parset="./SB45305/continuum_image_selfcal_F00_B00_LOOP0_14796240.in",
args={
"dataset": f"./{msname}",
"Images.shape": imagesize,
"nchanpercore": nchanpercore,
"restore": True,
"nworkergroups": 1
}
)
Starting imager
Debug: registered context Global=0
Debug: registered context Global=0
imager complete
[33]:
# Let's have a look at our image
loop_0_image_name = "image.i.RACS_0012+00.SB45305.cont.RACS_0012+00.beam00"
plot_image(loop_0_image_name + ".taylor.0.restored.fits", norm=0.5)
image(loop_0_image_name + ".taylor.0.restored.fits").statistics()
Filename: image.i.RACS_0012+00.SB45305.cont.RACS_0012+00.beam00.taylor.0.restored.fits
No. Name Ver Type Cards Dimensions Format
0 PRIMARY 1 PrimaryHDU 64 (2048, 2048, 1, 1) float32
None
2025-03-19 05:04:50 INFO FITSCoordinateUtil::fromFITSHeader Neither SPECSYS nor VELREF keyword given, spectral reference frame not defined ...
[33]:
{'npts': array([4194304.]),
'sum': array([28.99557935]),
'sumsq': array([96.33087276]),
'min': array([-0.03034493]),
'max': array([0.58392441]),
'mean': array([6.91308483e-06]),
'median': array([-0.00016392]),
'medabsdevmed': array([0.00277174]),
'quartile': array([0.00554887]),
'sigma': array([0.00479239]),
'rms': array([0.0047924]),
'minpos': array([ 6, 1494, 0, 0], dtype=int32),
'maxpos': array([2038, 1510, 0, 0], dtype=int32)}
Loop 1
Calibration
[12]:
loop_1_image_name = "image.i.RACS_0012+00.SB45305.cont.RACS_0012+00.beam00.01"
[13]:
ccalibrator(
env=env,
workers=imaging_workers+1,
parset="./SB45305/continuum_image_selfcal_F00_B00_LOOP1_14796240.in",
args={
"dataset": msname,
"sources.lsm.model": loop_0_image_name
}
)
Starting ccalibrator
Debug: registered context Global=0
Debug: registered context Global=0
ccalibrator complete
Imaging
[17]:
imager(
env=env,
workers=imaging_workers+1,
parset="./SB45305/continuum_image_selfcal_F00_B00_LOOP1_14796240.in",
args={
"dataset": f"./{msname}",
"Images.Names": f"[{loop_1_image_name}]",
"Images.shape": imagesize,
"Images.cellsize": "[2.5arcsec, 2.5arcsec]",
f"Images.{loop_1_image_name}.direction": "[00h10m18.828, +00.31.29.93, J2000]",
"Images.nyquistgridding": "true",
f"fImages.{loop_1_image_name}.nchan": 1,
f"Images.{loop_1_image_name}.nterms": 2,
"nchanpercore": nchanpercore,
"restore": True,
"nworkergroups": 1,
"calibaccess.table": "caldata_loop1.tab"
}
)
Starting imager
Debug: registered context Global=0
Debug: registered context Global=0
imager complete
[32]:
plot_image(loop_1_image_name + ".taylor.0.restored.fits", norm=0.5)
image(loop_1_image_name + ".taylor.0.restored.fits").statistics()
Filename: image.i.RACS_0012+00.SB45305.cont.RACS_0012+00.beam00.01.taylor.0.restored.fits
No. Name Ver Type Cards Dimensions Format
0 PRIMARY 1 PrimaryHDU 64 (2048, 2048, 1, 1) float32
None
2025-03-19 05:04:46 INFO FITSCoordinateUtil::fromFITSHeader Neither SPECSYS nor VELREF keyword given, spectral reference frame not defined ...
[32]:
{'npts': array([4194304.]),
'sum': array([43.43023437]),
'sumsq': array([94.51874968]),
'min': array([-0.02792455]),
'max': array([0.64738053]),
'mean': array([1.03545748e-05]),
'median': array([-0.00014657]),
'medabsdevmed': array([0.00269761]),
'quartile': array([0.00539827]),
'sigma': array([0.0047471]),
'rms': array([0.00474711]),
'minpos': array([1975, 1620, 0, 0], dtype=int32),
'maxpos': array([2038, 1510, 0, 0], dtype=int32)}
Loop 2
Calibrator
[19]:
loop_2_image_name = "image.i.RACS_0012+00.SB45305.cont.RACS_0012+00.beam00.02"
[20]:
ccalibrator(
env=env,
workers=imaging_workers+1,
parset="./SB45305/continuum_image_selfcal_F00_B00_LOOP2_14796240.in",
args={
"dataset": msname,
"sources.lsm.model": loop_1_image_name
}
)
Starting ccalibrator
Debug: registered context Global=0
Debug: registered context Global=0
ccalibrator complete
[23]:
imager(
env=env,
workers=imaging_workers+1,
parset="./SB45305/continuum_image_selfcal_F00_B00_LOOP2_14796240.in",
args={
"dataset": f"./{msname}",
"Images.Names": f"[{loop_2_image_name}]",
"Images.shape": imagesize,
"Images.cellsize": "[2.5arcsec, 2.5arcsec]",
f"Images.{loop_2_image_name}.direction": "[00h10m18.828, +00.31.29.93, J2000]",
"Images.nyquistgridding": "true",
f"fImages.{loop_2_image_name}.nchan": 1,
f"Images.{loop_2_image_name}.nterms": 2,
"nchanpercore": nchanpercore,
"restore": True,
"nworkergroups": 1,
"calibaccess.table": "caldata_loop2.tab"
}
)
Starting imager
Debug: registered context Global=0
Debug: registered context Global=0
imager complete
[31]:
plot_image(loop_2_image_name + ".taylor.0.restored.fits", norm=0.5)
image(loop_2_image_name + ".taylor.0.restored.fits").statistics()
Filename: image.i.RACS_0012+00.SB45305.cont.RACS_0012+00.beam00.02.taylor.0.restored.fits
No. Name Ver Type Cards Dimensions Format
0 PRIMARY 1 PrimaryHDU 64 (2048, 2048, 1, 1) float32
None
2025-03-19 05:03:20 INFO FITSCoordinateUtil::fromFITSHeader Neither SPECSYS nor VELREF keyword given, spectral reference frame not defined ...
[31]:
{'npts': array([4194304.]),
'sum': array([339.46755438]),
'sumsq': array([69.80795965]),
'min': array([-0.0694611]),
'max': array([0.71428245]),
'mean': array([8.0935372e-05]),
'median': array([-2.62012381e-05]),
'medabsdevmed': array([0.0021464]),
'quartile': array([0.00429399]),
'sigma': array([0.00407884]),
'rms': array([0.00407965]),
'minpos': array([2029, 1509, 0, 0], dtype=int32),
'maxpos': array([2038, 1510, 0, 0], dtype=int32)}
Alternatively starting with an initial model
[44]:
from astropy.coordinates import SkyCoord
from astropy import units as u
[45]:
votablename = "selavy-image.i.RACS_0012+00.SB45305.cont.taylor.0.restored.conv.components.xml"
model_name = "model_name.model"
[52]:
cmodel(
env=env,
workers=2,
args={
"gsm.database": "votable",
"gsm.file": votablename,
"gsm.ref_freq": frequency,
"primarybeam": "GaussianPB",
"bunit": "Jy/pixel",
"frequency": frequency,
"increment": increment,
"flux_limit": "10uJy",
"shape": imagesize,
"cellsize": ['2.5arcsec','2.5arcsec'],
"direction": direction,
"nterms": 1,
"stokes": "[I]",
"output": "casa",
"filename": model_name,
}
)
Starting cmodel
Debug: registered context Global=0
Debug: registered context Global=0
cmodel complete
Let’s view our model
[48]:
import matplotlib.pyplot as plt
from astropy.io import fits
from astropy.wcs import WCS
import matplotlib.colors as mcolors
from casacore.tables import table
tbl = table(model_name)
data=tbl.getcol("map")
data = data[0,0,0,:,:]
# print(type(data), data.shape)
tbl.close()
im = plt.imshow(data, origin='lower', cmap='viridis', norm=mcolors.PowerNorm(0.2))
plt.title("Test")
plt.colorbar(im)
image(model_name).statistics()
Successful readonly open of default-locked table model_name.model: 1 columns, 1 rows
[48]:
{'npts': array([4194304.]),
'sum': array([2.06941897]),
'sumsq': array([0.01438303]),
'min': array([0.]),
'max': array([0.06215899]),
'mean': array([4.9338793e-07]),
'median': array([0.]),
'medabsdevmed': array([0.]),
'quartile': array([0.]),
'sigma': array([5.85571408e-05]),
'rms': array([5.85592123e-05]),
'minpos': array([0, 0, 0, 0], dtype=int32),
'maxpos': array([ 72, 183, 0, 0], dtype=int32)}
Using the model for the initial self calibration loop
TODO
[ ]: