Gain calibration demo on simulated data
This notebook demonstrates gain calibration with ASKAPsoft. We create a simulated data with some corruption and then solve for the gains with ccalibrator and image the data after applying the gain solution
Tools used in this notebook 1. askapsoft * cmodel * csimulator * mslist * imager * ccalibrator 3. bptool 4. python-casacore
[1]:
source_name = 'PKS 1934-638'
source_ms = 'PKS1934-638'
ms_filename = f'{source_ms}.ms'
model_name = f'{source_ms}.model'
nchan = 48
frequency="887MHz"
increment="1MHz"
[2]:
# Clean up
import os
os.system(f"rm *.log *.i{source_ms}.fits")
os.system(f"rm -rf {ms_filename} {model_name}")
rm: *.iPKS1934-638.fits: No such file or directory
[2]:
0
[51]:
# 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
import numpy as np
from casacore.images import image
import matplotlib.patches as patches
%matplotlib inline
def plot_image(file, color_norm=None, norm=None, color_bar=False, box=None, print_stats=False):
# 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 color_norm:
im = ax.imshow(hdu.data[0,0,:,:], origin='lower', cmap='viridis', norm=mcolors.PowerNorm(color_norm))
elif norm:
sigma = np.std(hdu.data[0,0,:,:])
mean = np.std(hdu.data[0,0,:,:])
im = ax.imshow(hdu.data[0,0,:,:], origin='lower', cmap='viridis', vmin=-norm*sigma, vmax=norm*sigma)
else:
im = ax.imshow(hdu.data[0,0,:,:], origin='lower', cmap='viridis')
if print_stats:
if box:
stats = image(file).subimage(
(0, 0, box["f0"] - box["width"] // 2, box["f1"] - box["height"] // 2),
(0, 0, box["f0"] + box["width"] // 2, box["f1"] + box["height"] // 2)
).statistics()
else:
stats = image(file).statistics()
print("------------------------------------------")
print("------------ Mean and RMS region ---------")
print("------------------------------------------")
print(f"MEAN: {stats["mean"]}")
print(f"RMS: {stats["rms"]}")
print("------------------------------------------")
if box:
ax.add_patch(
patches.Rectangle(
xy=(box["f0"], box["f1"]), # point of origin.
width=box["width"], height=box["height"], linewidth=1,
color='red', fill=False))
ax.set_xlabel('RA')
ax.set_ylabel('Dec')
ax.set_autoscale_on(False)
ax.set_title(file)
plt.colorbar(im, ax=ax)
plt.show()
[4]:
from astropy.coordinates import SkyCoord
# Get the coordinates from Sesame
coord = SkyCoord.from_name(source_name)
ra ,dec = coord.to_string('hmsdms').split(" ")
# If you are offline, uncomment this
# ra ,dec = "08h16m36.24s -70d39m55.08s"
print(ra, dec)
direction=f"[{ra}, {dec}, J2000]"
19h39m25.0245312s -63d42m45.640368s
[5]:
from taitale import taitale_env
env = taitale_env(
image = "csirocass/askapsoft",
tag = "1.17.6-openmpi4",
)
[6]:
from taitale.askap import cmodel
votablename = "AS110/RACS1951-62A/image.i.SB8676.cont.RACS_1951-62A.linmos.taylor.0.restored.cres.fcor.components.xml"
cmodel(
env=env,
workers=3,
args={
"gsm.database": "votable",
"gsm.file": votablename,
"gsm.ref_freq": frequency,
"primarybeam": "GaussianPB",
"bunit": "Jy/pixel",
"frequency": frequency,
"increment": increment,
"flux_limit": "10uJy",
"shape": [2048, 2048],
"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
Debug: registered context Global=0
cmodel complete
Simulating the MS with some corruption
Creating a gain corruption table
[7]:
import numpy as np
nantennas = 36
def get_noisy_flux(gain = 1 + 0j, noise_level = 0.1):
magnitude = np.abs(gain)
noise_std = 0.1 * magnitude
noisy_real = np.random.normal(np.real(gain), noise_std)
noisy_imag = np.random.normal(np.imag(gain), noise_std)
return noisy_real, noisy_imag
gain_corruption_file = "gains.core.true"
with open(gain_corruption_file, "w+") as f:
for antenna_idx in range(nantennas):
g11 = get_noisy_flux()
g22 = get_noisy_flux()
# d12 = get_noisy_flux()
# d21 = get_noisy_flux()
f.write(f"gain.g11.{antenna_idx}.0 = [{g11[0]}, {g11[1]}]\n")
f.write(f"gain.g22.{antenna_idx}.0 = [{g22[0]}, {g22[1]}]\n")
Simulating data with gain corruption
[8]:
from taitale.askap import csimulator
import os
os.system(f"rm -rf {ms_filename}")
csimulator(
env=env,
parset="./parsets/csimulator.in",
args={
"dataset": ms_filename,
"sources.names": f"[{source_name}]",
f"sources.{source_name}.direction": direction,
f"sources.{source_name}.model": model_name,
"spws.names": "[TestBand]",
"spws.TestBand": f'[48, {frequency}, 1.0MHz, "XX XY YX YY"]',
"observe.number": 1,
"observe.scan0": f"[{source_name}, TestBand, -2h, 2h]",
"noise": "false",
"noise.rms": 1,
"corrupt": "true",
"calibaccess": "parset",
"calibaccess.parset": "gains.core.true"
},
)
Starting csimulator
Debug: registered context Global=0
csimulator complete
[9]:
from taitale.askap import mslist
mslist(env=env, args=f"--full {ms_filename}", logfile="mslist.txt")
Starting mslist
mslist complete
[10]:
!cat mslist.txt
Debug: registered context Global=0
2025-04-11 06:52:20 INFO Observer: ASKAP simulator Project:
2025-04-11 06:52:20 INFO + Observation: ASKAP36(36 antennas)
2025-04-11 06:52:20 INFO MSMetaData::_computeScanAndSubScanProperties Computing scan and subscan properties...
2025-04-11 06:52:20 INFO Data records: 60480 Total elapsed time = 14400 seconds
2025-04-11 06:52:20 INFO + Observed from 06-Mar-2007/22:56:44.6 to 07-Mar-2007/02:56:44.6 (UTC)
2025-04-11 06:52:20 INFO
2025-04-11 06:52:20 INFO Fields: 1
2025-04-11 06:52:20 INFO + ID Code Name RA Decl Epoch nRows
2025-04-11 06:52:20 INFO + 0 PKS 1934-638 19:39:25.024531 -63.42.45.64037 J2000 60480
2025-04-11 06:52:20 INFO Spectral Windows: (1 unique spectral windows and 1 unique polarization setups)
2025-04-11 06:52:20 INFO + SpwID Name #Chans Frame Ch0(MHz) ChanWid(kHz) TotBW(kHz) CtrFreq(MHz) Corrs
2025-04-11 06:52:20 INFO + 0 TestBand 48 TOPO 887.000 1000.000 48000.0 910.5000 XX XY YX YY
2025-04-11 06:52:20 INFO The SOURCE table is absent: see the FIELD table
2025-04-11 06:52:20 INFO + Antennas: 36 'name'='station'
2025-04-11 06:52:20 INFO + ID= 0-4: 'Pad01'='', 'Pad02'='', 'Pad03'='', 'Pad04'='', 'Pad05'='',
2025-04-11 06:52:20 INFO + ID= 5-9: 'Pad06'='', 'Pad07'='', 'Pad08'='', 'Pad09'='', 'Pad10'='',
2025-04-11 06:52:20 INFO + ID= 10-14: 'Pad11'='', 'Pad12'='', 'Pad13'='', 'Pad14'='', 'Pad15'='',
2025-04-11 06:52:20 INFO + ID= 15-19: 'Pad16'='', 'Pad17'='', 'Pad18'='', 'Pad19'='', 'Pad20'='',
2025-04-11 06:52:20 INFO + ID= 20-24: 'Pad21'='', 'Pad22'='', 'Pad23'='', 'Pad24'='', 'Pad25'='',
2025-04-11 06:52:20 INFO + ID= 25-29: 'Pad26'='', 'Pad27'='', 'Pad28'='', 'Pad29'='', 'Pad30'='',
2025-04-11 06:52:20 INFO + ID= 30-34: 'Pad31'='', 'Pad32'='', 'Pad33'='', 'Pad34'='', 'Pad35'='',
2025-04-11 06:52:20 INFO + ID= 35-35: 'Pad36'=''
2025-04-11 06:52:20 INFO The HISTORY table is empty
Imaging without calibration
[11]:
from taitale.askap import imager
image_name = f"image.i.{source_ms}"
imaging_workers = 1
nchanpercore = nchan // imaging_workers
print(f"channels: {nchan}")
print(f"nchanpercore: {nchanpercore}")
print(f"mpi_workers: {imaging_workers}")
# Default parset for Clean used
imager(
env=env,
workers=imaging_workers+1,
parset="./parsets/cimager.in",
args={
"dataset": f"./{ms_filename}",
"Images.Names": image_name,
f"Images.{image_name}.direction": direction,
f"Images.{image_name}.nchan": 1,
f"Images.{image_name}.nterms": 2,
"Images.shape": "[512, 512]",
"Images.cellsize": "['2.5arcsec', '2.5arcsec']",
"nchanpercore": nchanpercore,
"restore": True,
}
)
channels: 48
nchanpercore: 48
mpi_workers: 1
Starting imager
Debug: registered context Global=0
Debug: registered context Global=0
imager complete
[52]:
plot_image(image_name + ".taylor.0.restored.fits", norm=3, color_bar=True, box=dict(f0=206, f1=206, width=100, height=100), print_stats=True)
Filename: image.i.PKS1934-638.taylor.0.restored.fits
No. Name Ver Type Cards Dimensions Format
0 PRIMARY 1 PrimaryHDU 64 (512, 512, 1, 1) float32
None
------------------------------------------
------------ Mean and RMS region ---------
------------------------------------------
MEAN: [0.01395344]
RMS: [0.32471413]
------------------------------------------
2025-04-11 07:16:25 INFO FITSCoordinateUtil::fromFITSHeader Neither SPECSYS nor VELREF keyword given, spectral reference frame not defined ...
Calibration
As a simple demonstration we are going to use the model we used for simulation for the calibration. ccalibrator will solve the gains and create a gain table
[21]:
from taitale.askap import ccalibrator
ccalibrator(
env=env,
workers=2,
args={
"dataset": ms_filename,
"refgain": "gain.g11.0.0",
"sources.names": "[PKS1934638]",
"sources.PKS1934638.direction": direction,
"sources.PKS1934638.model": model_name,
"gridder.snapshotimaging": "false",
"gridder": "WProject",
"gridder.WProject.wmax": 35000,
"gridder.WProject.nwplanes": 257,
"gridder.WProject.oversample": 7,
"gridder.WProject.maxsupport": 512,
"gridder.WProject.variablesupport": "true",
"gridder.WProject.offsetsupport": "true",
"gridder.WProject.sharecf": "true",
"gridder.WProject.frequencydependent": "true",
"ncycles": "5",
"solver": "SVD",
# "calibaccess": "parset",
# "calibaccess.parset": "caldata.parset",
"calibaccess": "table",
"calibaccess.table": "caldata.tab",
"calibaccess.table.maxant": 36,
"calibaccess.table.maxbeam": 36,
"calibaccess.table.maxchan": 48,
"calibaccess.table.reuse": "false"
}
)
Starting ccalibrator
Debug: registered context Global=0
Debug: registered context Global=0
ccalibrator complete
Viewing our calibration solution
[22]:
bptool_env = taitale_env(
image = "wasimraja81/askappy-ubuntu",
tag = "bptool-2.4",
)
[23]:
from taitale.utils import CliArgs
from taitale.utils.wrapper import cli_app
def plot_gains(**kwargs):
cli_args = CliArgs()
bptool_app = cli_app(
name="plot_gains",
cmd="plot_gains.py {option}",
cli_args=cli_args,
)
bptool_app(**kwargs)
[24]:
plot_gains(env=bptool_env, args=f"-t caldata.tab -b1 0 -b2 0 -s -m {ms_filename}", logfile="bptool.txt")
Starting plot_gains
/usr/local/lib/python3.6/dist-packages/astropy/config/configuration.py:581: ConfigurationMissingWarning: Configuration defaults will be used due to PermissionError:13 on None
warn(ConfigurationMissingWarning(msg))
Matplotlib created a temporary config/cache directory at /tmp/matplotlib-gqtmsnla because the default path (/.config/matplotlib) is not a writable directory; it is highly recommended to set the MPLCONFIGDIR environment variable to a writable directory, in particular to speed up the import of Matplotlib and to better support multiprocessing.
WARNING: CacheMissingWarning: Remote data cache could not be accessed due to PermissionError: [Errno 13] Permission denied: '/.astropy' [astropy.utils.data]
WARNING: CacheMissingWarning: Not clearing data cache - cache inaccessible due to PermissionError: [Errno 13] Permission denied: '/.astropy' [astropy.utils.data]
WARNING: CacheMissingWarning: ("Cache directory cannot be read or created ([Errno 13] Permission denied: '/.astropy'), providing data in temporary file instead.", '/tmp/astropy-download-1-ohxmqp3l') [astropy.utils.data]
/_build/bptool/bptool/bin/plot_gains.py:240: RuntimeWarning: invalid value encountered in double_scalars
dT = (time_arr[nTime-1] - time_arr[0])/(nTime-1)
plot_gains complete
[25]:
!cat bptool.txt
Successful readonly open of default-locked table caldata.tab: 3 columns, 1 rows
Plotting gain solutions for Input table: caldata.tab
For:
Beam Nums: 0 - 0
Ante Nums: 0 - 0
Ignore Beam: False
[26]:
from IPython.display import Image
Image(filename='caldata.beam00.ante00.png')
[26]:
Applying the calibration solution
Finally we’ll use the calibration table for imaging
[27]:
image_name_cal = f"image.i.{source_ms}.cal"
imager(
env=env,
workers=imaging_workers+1,
parset="./parsets/cimager.in",
args={
"dataset": f"./{ms_filename}",
"Images.Names": image_name_cal,
f"Images.{image_name_cal}.direction": direction,
f"Images.{image_name_cal}.nchan": 1,
f"Images.{image_name_cal}.nterms": 2,
"Images.shape": "[512, 512]",
"Images.cellsize": "['2.5arcsec', '2.5arcsec']",
"nchanpercore": nchanpercore,
"restore": True,
# Calibrate with the generated gain table
"calibrate": "true",
"calibrate.ignorebeam": "true",
"calibrate.allowflag": "true",
"calibrate.scalenoise": "true",
"calibrate.interpolatetime": "true",
"calibaccess": "table",
"calibaccess.table": "caldata.tab",
"calibaccess.table.maxant": 36,
"calibaccess.table.maxbeam": 36,
"calibaccess.table.maxchan": 48,
"calibaccess.table.reuse": "false"
}
)
Starting imager
Debug: registered context Global=0
Debug: registered context Global=0
imager complete
Image after calibration
[53]:
plot_image(image_name_cal + ".taylor.0.restored.fits", norm=3, color_bar=True, box=dict(f0=206, f1=206, width=100, height=100), print_stats=True)
Filename: image.i.PKS1934-638.cal.taylor.0.restored.fits
No. Name Ver Type Cards Dimensions Format
0 PRIMARY 1 PrimaryHDU 64 (512, 512, 1, 1) float32
None
------------------------------------------
------------ Mean and RMS region ---------
------------------------------------------
MEAN: [0.01338182]
RMS: [0.3210266]
------------------------------------------
2025-04-11 07:16:38 INFO FITSCoordinateUtil::fromFITSHeader Neither SPECSYS nor VELREF keyword given, spectral reference frame not defined ...
[ ]: