EveryBeam OSKAR SKA integration tests

Not to put a dampener on WODEN, but if you’re using the the OSKAR primary beam, I would just use OSKAR to make visibilities. Sure, WODEN does precession, which OSKAR doesn’t, but OSKAR has a fully GPU version of the beam which is so much faster. The EveryBeam OSKAR implementation takes a large, large amount of time to run.

These tests are all setup and explained in test_MWA.ipynb. We’ll look a the beam shape, test Stokes recovery, and make an image from visibilities.

Goes without saying, but these tests rely on everybeam being installed.

import os
from subprocess import call
from astropy.io import fits
import numpy as np
from astropy.table import Column, Table
from astropy.time import Time
from astropy.coordinates import SkyCoord
from astropy.coordinates import EarthLocation
from astropy import units as u
import matplotlib.pyplot as plt
import numpy.testing as npt
from astropy.constants import c
from astropy.wcs import WCS
import everybeam as eb
from wodenpy.primary_beam.use_everybeam import run_everybeam, run_everybeam_over_threads
import erfa
import mwa_hyperbeam
from wodenpy.use_libwoden.create_woden_struct_classes import Woden_Struct_Classes
from ctypes import c_double
from astropy.time import TimeDelta
from wodenpy.array_layout.precession import RTS_Precess_LST_Lat_to_J2000
from wodenpy.array_layout.create_array_layout import convert_ecef_to_enh

import sys
sys.path.append('../../scripts/')
from run_woden import main as run_woden

##A bunch of test and plotting code lives here so we can use it in multile notebooks
from eb_testing_code import create_WCS, plot_jones_on_sky, plot_everybeam_on_sky, make_sky_models, read_uvfits, convert_inst_to_stokes, test_stokes_recovery, getFDF, findpeaks, test_RM_recovery, make_RM_skymodel

# os.environ["MWA_BEAM_FILE"] = "/home/jline/software/MWA_beam_files/mwa_full_embedded_element_pattern.h5"

C = c.to('m/s').value

from casacore.tables import table

First up just quickly, let’s see what station layout I’m using

MWA_LAT=-26.703319405555554
MWA_LONG=116.67081523611111

with table('create_OSKAR-SKA_ms/OSKAR-SKA-layout.ms' + '/ANTENNA') as t: 
        
            num_ants = len(t)
            ant_locations = np.array([t.getcell('POSITION', ant) for ant in range(num_ants)])
            ##convert from ECEF to ENH, as WODEN starts with enh coords
            east, north, height = convert_ecef_to_enh(ant_locations[:,0],
                                        ant_locations[:,1], ant_locations[:,2],
                                        np.radians(MWA_LONG),
                                        np.radians(MWA_LAT))

east /= 1000
north /= 1000

fig, axs = plt.subplots(1, 3, figsize=(15, 5), layout='constrained')

for ax in axs:
    ax.plot(east, north, 'o',mfc='none', ms=2)

edge = 6
x_cent = 9.3
y_cent = -13.5

axs[1].set_xlim(x_cent-edge, x_cent+edge)
axs[1].set_ylim(y_cent-edge, y_cent+edge)

edge = 1.5
axs[2].set_xlim(x_cent-edge, x_cent+edge)
axs[2].set_ylim(y_cent-edge, y_cent+edge)

for ax in axs:
    ax.set_xlabel('East (km)')
    ax.set_aspect('equal', 'box')
axs[0].set_ylabel('North (km)')
plt.show()

Successful readonly open of default-locked table create_OSKAR-SKA_ms/OSKAR-SKA-layout.ms/ANTENNA: 8 columns, 512 rows

Image the beam on the sky

First up, let’s see if we can just plot an OSKAR SKA station beam to check it looks sensible.

##pick a time/date that sticks our phase centre overhead
date = "2024-07-21T20:13:00"
##Assume that the OSKAR telescope is near the MWA??
mwa_location = EarthLocation(lat=MWA_LAT*u.deg, 
                            lon=MWA_LONG*u.deg,
                            height=377.827)
observing_time = Time(date, scale='utc', location=mwa_location)

##Grab the LST
LST = observing_time.sidereal_time('mean')
LST_deg = LST.value*15
print(f"LST: {LST_deg} deg")

LST: 0.0045520088542088475 deg

ms_path="create_OSKAR-SKA_ms/OSKAR-SKA-layout.ms"

ra0 = LST_deg
dec0 = MWA_LAT

mjd_current = observing_time.mjd
lst_J2000, latitude_J2000 = RTS_Precess_LST_Lat_to_J2000(np.radians(LST_deg),
                                                         np.radians(MWA_LAT),
                                                         mjd_current)

nside=50
radec_reso=60/nside

header, wcs = create_WCS(ra0, dec0, nside, radec_reso)
x_mesh, y_mesh = np.meshgrid(np.arange(nside), np.arange(nside))
ras, decs = wcs.all_pix2world(x_mesh, y_mesh, 0)

ras = np.radians(ras.flatten())
decs = np.radians(decs.flatten())

freq = 100e+6

coeff_path=""
station_id = 0
station_ids = [0]
j2000_latitudes = [latitude_J2000]
j2000_lsts = [lst_J2000]
times = [observing_time]
freqs = [freq]
current_latitude = np.radians(MWA_LAT)
current_longitude = np.radians(MWA_LONG)


apply_beam_norms=True
parallactic_rotate=True
iau_order=True

element_response = "skala40_wave"

##How many cores to split the work over
num_threads=8

# telescope = load_OSKAR_telescope(ms_path)

all_jones = run_everybeam_over_threads(num_threads, ms_path,
                            coeff_path,
                            ras, decs,
                            np.radians(ra0), np.radians(dec0),
                            j2000_latitudes, j2000_lsts,
                            times, freqs, station_ids,
                            apply_beam_norms=apply_beam_norms,
                            iau_order=iau_order,
                            parallactic_rotate=parallactic_rotate,
                            element_response_model=element_response)

beam_ind, time_ind, freq_ind = 0, 0, 0

all_gx_eb = all_jones[beam_ind, time_ind, freq_ind, :, 0, 0]
all_Dx_eb = all_jones[beam_ind, time_ind, freq_ind, :, 0, 1]
all_Dy_eb = all_jones[beam_ind, time_ind, freq_ind, :, 1, 0]
all_gy_eb = all_jones[beam_ind, time_ind, freq_ind, :, 1, 1]

all_gx_eb.shape = (nside, nside)
all_Dx_eb.shape = (nside, nside)
all_Dy_eb.shape = (nside, nside)
all_gy_eb.shape = (nside, nside)

plot_jones_on_sky(all_gx_eb, all_Dx_eb, all_Dy_eb, all_gy_eb, wcs, 'Everybeam OSKAR SKA')

Thread 0 processing coords 0 to 313Thread 2 processing coords 626 to 939Thread 4 processing coords 1252 to 1565Thread 5 processing coords 1565 to 1878Thread 1 processing coords 313 to 626Thread 3 processing coords 939 to 1252


Thread 7 processing coords 2191 to 2504Thread 6 processing coords 1878 to 2191

2025-03-05 23:29:06	SEVERE	MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)	Leap second table TAI_UTC seems out-of-date.
2025-03-05 23:29:06	SEVERE	MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+	Until the table is updated (see the CASA documentation or your system admin),
2025-03-05 23:29:06	SEVERE	MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+	times and coordinates derived from UTC could be wrong by 1s or more.
2025-03-05 23:29:06	SEVERE	MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)	Leap second table TAI_UTC seems out-of-date.
2025-03-05 23:29:06	SEVERE	MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+	Until the table is updated (see the CASA documentation or your system admin),
2025-03-05 23:29:06	SEVERE	MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+	times and coordinates derived from UTC could be wrong by 1s or more.
2025-03-05 23:29:06	SEVERE	MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)	Leap second table TAI_UTC seems out-of-date.
2025-03-05 23:29:06	SEVERE	MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+	Until the table is updated (see the CASA documentation or your system admin),
2025-03-05 23:29:06	SEVERE	MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+	times and coordinates derived from UTC could be wrong by 1s or more.
2025-03-05 23:29:06	SEVERE	MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)	Leap second table TAI_UTC seems out-of-date.
2025-03-05 23:29:06	SEVERE	MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+	Until the table is updated (see the CASA documentation or your system admin),
2025-03-05 23:29:06	SEVERE	MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+	times and coordinates derived from UTC could be wrong by 1s or more.
2025-03-05 23:29:06	SEVERE	MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)	Leap second table TAI_UTC seems out-of-date.
2025-03-05 23:29:06	SEVERE	MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+	Until the table is updated (see the CASA documentation or your system admin),
2025-03-05 23:29:06	SEVERE	MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+	times and coordinates derived from UTC could be wrong by 1s or more.
2025-03-05 23:29:06	SEVERE	MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)	Leap second table TAI_UTC seems out-of-date.
2025-03-05 23:29:06	SEVERE	MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+	Until the table is updated (see the CASA documentation or your system admin),
2025-03-05 23:29:06	SEVERE	MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+	times and coordinates derived from UTC could be wrong by 1s or more.

Could not load dataset for frequency 100 MHz, using the nearest neighbor with frequency 110 MHz instead
Could not load dataset for frequency 100 MHz, using the nearest neighbor with frequency 110 MHz instead
Could not load dataset for frequency 100 MHz, using the nearest neighbor with frequency 110 MHz instead
Could not load dataset for frequency 100 MHz, using the nearest neighbor with frequency 110 MHz instead
Could not load dataset for frequency 100 MHz, using the nearest neighbor with frequency 110 MHz instead
Could not load dataset for frequency 100 MHz, using the nearest neighbor with frequency 110 MHz instead

2025-03-05 23:29:07	SEVERE	MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)	Leap second table TAI_UTC seems out-of-date.
2025-03-05 23:29:07	SEVERE	MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+	Until the table is updated (see the CASA documentation or your system admin),
2025-03-05 23:29:07	SEVERE	MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+	times and coordinates derived from UTC could be wrong by 1s or more.
2025-03-05 23:29:07	SEVERE	MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)	Leap second table TAI_UTC seems out-of-date.
2025-03-05 23:29:07	SEVERE	MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+	Until the table is updated (see the CASA documentation or your system admin),
2025-03-05 23:29:07	SEVERE	MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+	times and coordinates derived from UTC could be wrong by 1s or more.

Could not load dataset for frequency 100 MHz, using the nearest neighbor with frequency 110 MHz instead
Could not load dataset for frequency 100 MHz, using the nearest neighbor with frequency 110 MHz instead
Thread 5 finished
Thread 0 finishedThread 6 finished

Thread 7 finished
Thread 1 finished
Thread 4 finished
Thread 2 finished
Thread 3 finished

Stokes recovery

Now try to recover single point source of either Stokes I, Q, U or V.

make_sky_models(ra0, dec0)

freq_reso = 1e+6
# low_freq = 180e+6
# high_freq = 210e+6
low_freq = 100e+6
high_freq = 150e+6


num_freq_chans = int((high_freq - low_freq) / freq_reso)

primary_beam = "everybeam_OSKAR"

for pol in ['I', 'Q', 'U', 'V']:
# for pol in ['I']:

    uvfits_name = f"stokes{pol}_{primary_beam}"
    cat_name = f'{pol}_source.fits'
    
    ##The command to run WODEN
    command = f'run_woden.py --ra0={ra0} --dec0={dec0} '
    command += '--eb_point_to_phase ' ##this overrides whatever pointing is in beam_ms_path
                                      ##and sets it to ra0, dec0
    command += f'--latitude={MWA_LAT} --longitude={MWA_LONG} '
    command += f'--date={date} --output_uvfits_prepend={uvfits_name} '
    command += f'--cat_filename={cat_name} --primary_beam={primary_beam} '
    command += f'--lowest_channel_freq={low_freq} --freq_res={freq_reso} '
    command += f'--num_freq_channels={num_freq_chans} --band_nums=1 '
    command += f'--time_res=2 --num_time_steps=1 --IAU_order '
    command += f' --station_id=0 '
    command += f'--beam_ms_path=create_OSKAR-SKA_ms/OSKAR-SKA-layout.ms '
    command += f'--num_threads=1 '
    
    call(command, shell=True)

Note I’ve deleted the outputs of the cell above, because EveryBeam spits out a bunch of error messages with the OSKAR beam warning about frequencies mismatching, and not having a central direction set. The outputs seems fine. It’s a TODO to silence these warnings.

for pol in ['I', 'Q', 'U', 'V']:
    uvfits_name = f"stokes{pol}_everybeam_OSKAR"
    XX, XY, YX, YY = read_uvfits(f'{uvfits_name}_band01.uvfits')
    ##pick a random baseline to plot, they should all be the sam
    baseline = 0

    recover_I, recover_Q, recover_U, recover_V = convert_inst_to_stokes(XX[baseline], XY[baseline], YX[baseline], YY[baseline])
    
    print(f"Input {pol}, Recover I {recover_I[0].real:.2f}, Q {recover_Q[0].real:.2f}, U {recover_U[0].real:.2f}, V {recover_V[0].real:.2f}")
    # print(f"Input {pol}, Recover I {recover_I[-1].real:.2f}, Q {recover_Q[-1].real:.2f}, U {recover_U[-1].real:.2f}, V {recover_V[-1].real:.2f}")
    
    test_stokes_recovery(pol, 'everybeam_OSKAR', atol=9e-2)

Input I, Recover I 1.00, Q 0.00, U 0.00, V 0.00
Testing Stokes I
Stokes I passed
Input Q, Recover I 0.00, Q 1.00, U -0.00, V -0.00
Testing Stokes Q
Stokes Q passed
Input U, Recover I 0.00, Q 0.00, U 1.00, V 0.00
Testing Stokes U
Stokes U passed
Input V, Recover I 0.00, Q 0.00, U 0.00, V 1.00
Testing Stokes V
Stokes V passed

We get back what we put in so all is well.

RM recovery

phi_RM, pol_frac = make_RM_skymodel(ra0, dec0)

freq_reso = 0.1e+6
low_freq = 100e+6
high_freq = 130e+6
num_freq_chans = int((high_freq - low_freq) / freq_reso)


primary_beam = "everybeam_OSKAR"

uvfits_name = f"rm_source_{primary_beam}"
cat_name = 'RM_source.fits'

##The command to run WODEN
command = f'run_woden.py --ra0={ra0} --dec0={dec0} '
command += f'--date={date} --output_uvfits_prepend={uvfits_name} '
command += f'--latitude={MWA_LAT} --longitude={MWA_LONG} '
command += f'--cat_filename={cat_name} --primary_beam={primary_beam} '
command += f'--lowest_channel_freq={low_freq} --freq_res={freq_reso} '
command += f'--num_freq_channels={num_freq_chans} --band_nums=1 '
command += f'--time_res=2 --num_time_steps=1 --IAU_order '
command += f' --station_id=0 '
command += f'--beam_ms_path=create_OSKAR-SKA_ms/OSKAR-SKA-layout.ms '
command += f'--num_threads=1 '

call(command, shell=True)

Note I’ve deleted the outputs of the cell above, because EveryBeam spits out a bunch of error messages with the OSKAR beam warning about frequencies mismatching, and not having a central direction set. The outputs seems fine. It’s a TODO to silence these warnings.

XX, XY, YX, YY = read_uvfits('rm_source_everybeam_OSKAR_band01.uvfits')
baseline = 0    
recover_I, recover_Q, recover_U, recover_V = convert_inst_to_stokes(XX[baseline], XY[baseline], YX[baseline], YY[baseline])

    
freqs = np.arange(low_freq, high_freq, freq_reso)

fig, axs = plt.subplots(1, 1, figsize=(8, 6))

axs.plot(freqs / 1e+6, recover_I.real, 'k-', label='Recovered I')
axs.plot(freqs / 1e+6, recover_Q.real, 'r', label='Recovered Q')
axs.plot(freqs / 1e+6, recover_U.real, 'b', label='Recovered U')
axs.plot(freqs / 1e+6, recover_V.real, 'g', label='Recovered V')

axs.legend()

axs.set_xlabel('Frequency (MHz)')
axs.set_ylabel('Flux (Jy)')

plt.show()

test_RM_recovery(uvfits_name, phi_RM, pol_frac, freqs, atol=0.02)

Recovered RM: 50.0 Expected RM: 50
Recovered Pol. Fraction: 0.49999887 Expected Pol Fraction 0.5

We recover the correct RM which is nice.

Test unique primary beam per station

Righto, we’ll just run the first four stations as the skala40_wave model is slow and I don’t have all day. Run it for two time steps, two freqs to check things change sensibly.

After we’ll run a full WODEN simulation with an off-zenith point source to check that all visibilities are unique.

##Setup a grid of RA/Dec on the sky
##Setup a grid of RA/Dec on the sky
nside=64
radec_reso = 50/nside

header, wcs = create_WCS(ra0, dec0, nside, radec_reso)
x_mesh, y_mesh = np.meshgrid(np.arange(nside), np.arange(nside))
ras, decs = wcs.all_pix2world(x_mesh, y_mesh, 0)

ras = np.radians(ras.flatten())
decs = np.radians(decs.flatten())

num_comps = len(ras)

all_freqs = np.array([75e+6, 200e+6])
all_times = np.array([observing_time, observing_time + TimeDelta(3*3600.0, format='sec')])
num_beams = 1
##Do the first four stations
station_ids = np.arange(4)

lsts = []
latitudes = []

for obs_time in all_times:
    
    LST = obs_time.sidereal_time('mean')
    lst_current = np.radians(LST.value*15)
    
    lst_J2000, latitude_J2000 = RTS_Precess_LST_Lat_to_J2000(
                                lst_current,
                                np.radians(MWA_LAT),
                                obs_time.mjd)
    
    lsts.append(lst_J2000)
    latitudes.append(latitude_J2000)

coeff_path=""
ms_path="create_OSKAR-SKA_ms/OSKAR-SKA-layout.ms"
apply_beam_norms = True
iau_order = True
parallactic_rotate = True
element_response = "skala40_wave"

num_threads = 8

all_jones = run_everybeam_over_threads(num_threads, ms_path,
                            coeff_path,
                            ras, decs,
                            np.radians(ra0), np.radians(dec0),
                            latitudes, lsts,
                            all_times, all_freqs, station_ids,
                            apply_beam_norms=apply_beam_norms,
                            iau_order=iau_order,
                            parallactic_rotate=parallactic_rotate,
                            element_response_model=element_response)

Thread 0 processing coords 0 to 512Thread 1 processing coords 512 to 1024
Thread 4 processing coords 2048 to 2560Thread 3 processing coords 1536 to 2048Thread 2 processing coords 1024 to 1536Thread 5 processing coords 2560 to 3072

Thread 6 processing coords 3072 to 3584



Thread 7 processing coords 3584 to 4096

2025-03-05 23:41:02	SEVERE	MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)	Leap second table TAI_UTC seems out-of-date.
2025-03-05 23:41:02	SEVERE	MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+	Until the table is updated (see the CASA documentation or your system admin),
2025-03-05 23:41:02	SEVERE	MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+	times and coordinates derived from UTC could be wrong by 1s or more.
2025-03-05 23:41:02	SEVERE	MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)	Leap second table TAI_UTC seems out-of-date.
2025-03-05 23:41:02	SEVERE	MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+	Until the table is updated (see the CASA documentation or your system admin),
2025-03-05 23:41:02	SEVERE	MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+	times and coordinates derived from UTC could be wrong by 1s or more.
2025-03-05 23:41:02	SEVERE	MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)	Leap second table TAI_UTC seems out-of-date.
2025-03-05 23:41:02	SEVERE	MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+	Until the table is updated (see the CASA documentation or your system admin),
2025-03-05 23:41:02	SEVERE	MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+	times and coordinates derived from UTC could be wrong by 1s or more.
2025-03-05 23:41:02	SEVERE	MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)	Leap second table TAI_UTC seems out-of-date.
2025-03-05 23:41:02	SEVERE	MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+	Until the table is updated (see the CASA documentation or your system admin),
2025-03-05 23:41:02	SEVERE	MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+	times and coordinates derived from UTC could be wrong by 1s or more.
2025-03-05 23:41:02	SEVERE	MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)	Leap second table TAI_UTC seems out-of-date.
2025-03-05 23:41:02	SEVERE	MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+	Until the table is updated (see the CASA documentation or your system admin),
2025-03-05 23:41:02	SEVERE	MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+	times and coordinates derived from UTC could be wrong by 1s or more.
2025-03-05 23:41:02	SEVERE	MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)	Leap second table TAI_UTC seems out-of-date.
2025-03-05 23:41:02	SEVERE	MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+	Until the table is updated (see the CASA documentation or your system admin),
2025-03-05 23:41:02	SEVERE	MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+	times and coordinates derived from UTC could be wrong by 1s or more.
2025-03-05 23:41:02	SEVERE	MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)	Leap second table TAI_UTC seems out-of-date.
2025-03-05 23:41:02	SEVERE	MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+	Until the table is updated (see the CASA documentation or your system admin),
2025-03-05 23:41:02	SEVERE	MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+	times and coordinates derived from UTC could be wrong by 1s or more.
2025-03-05 23:41:02	SEVERE	MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)	Leap second table TAI_UTC seems out-of-date.
2025-03-05 23:41:02	SEVERE	MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+	Until the table is updated (see the CASA documentation or your system admin),
2025-03-05 23:41:02	SEVERE	MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+	times and coordinates derived from UTC could be wrong by 1s or more.

Could not load dataset for frequency 75 MHz, using the nearest neighbor with frequency 70 MHz instead
Could not load dataset for frequency 75 MHz, using the nearest neighbor with frequency 70 MHz instead
Could not load dataset for frequency 75 MHz, using the nearest neighbor with frequency 70 MHz instead
Could not load dataset for frequency 75 MHz, using the nearest neighbor with frequency 70 MHz instead
Could not load dataset for frequency 75 MHz, using the nearest neighbor with frequency 70 MHz instead
Could not load dataset for frequency 75 MHz, using the nearest neighbor with frequency 70 MHz instead
Could not load dataset for frequency 75 MHz, using the nearest neighbor with frequency 70 MHz instead
Could not load dataset for frequency 75 MHz, using the nearest neighbor with frequency 70 MHz instead
Could not load dataset for frequency 200 MHz, using the nearest neighbor with frequency 230 MHz instead
Could not load dataset for frequency 200 MHz, using the nearest neighbor with frequency 230 MHz instead
Could not load dataset for frequency 200 MHz, using the nearest neighbor with frequency 230 MHz instead
Could not load dataset for frequency 200 MHz, using the nearest neighbor with frequency 230 MHz instead
Could not load dataset for frequency 200 MHz, using the nearest neighbor with frequency 230 MHz instead
Could not load dataset for frequency 200 MHz, using the nearest neighbor with frequency 230 MHz instead
Could not load dataset for frequency 200 MHz, using the nearest neighbor with frequency 230 MHz instead
Could not load dataset for frequency 200 MHz, using the nearest neighbor with frequency 230 MHz instead
Thread 5 finished
Thread 1 finished
Thread 2 finished
Thread 7 finished
Thread 0 finished
Thread 4 finished
Thread 6 finished
Thread 3 finished

def plot_beam_selection(all_jones, time_ind, freq_ind, nside, wcs):
    fig, axs = plt.subplots(2, 2, figsize=(12, 12), layout='constrained',
                            subplot_kw={'projection': wcs})
    
    for col in range(2):
        for row in range(2):
            beam_ind = col*2 + row
            
            gx = all_jones[beam_ind, time_ind, freq_ind, :, 0,0]
            gx.shape = (nside, nside)
            
            im = axs[row, col].imshow(np.log10(np.abs(gx)), origin='lower')
            
            axs[row, col].grid(alpha=0.5)
            
            lon = axs[row, col].coords[0]
            lat = axs[row, col].coords[1]
            
            lon.set_ticks_visible(False)
            lat.set_ticks_visible(False)
            lon.set_axislabel(' ')
            lat.set_axislabel(' ')
            
            if col > 0:
                lat.set_ticklabel_visible(False)
            if row < 2:
                lon.set_ticklabel_visible(False)
                
            if col == 2:
                plt.colorbar(im, ax=axs[row, col], label='log10(abs(gx))')
            

    fig.suptitle(f"Time: {all_times[time_ind].isot}, Freq: {all_freqs[freq_ind]/1e+6} MHz")

    plt.show()

time_ind = 0
freq_ind = 0
plot_beam_selection(all_jones, time_ind, freq_ind, nside, wcs)

time_ind = 0
freq_ind = 1
plot_beam_selection(all_jones, time_ind, freq_ind, nside, wcs)

time_ind = 1
freq_ind = 0
plot_beam_selection(all_jones, time_ind, freq_ind, nside, wcs)

time_ind = 1
freq_ind = 1
plot_beam_selection(all_jones, time_ind, freq_ind, nside, wcs)

OK, so 3.6s to calculate the LOFAR Harmaker beam for a single station (but for the whole observation)

Let’s make an image

Just create a random sky image of some point sources for fun

ra_coord_hw = 15
dec_coord_hw = 10

num_coords = 500

ras = np.random.uniform(ra0 - ra_coord_hw, ra0 + ra_coord_hw, num_coords)
decs = np.random.uniform(dec0 - dec_coord_hw, dec0 + dec_coord_hw, num_coords)

ras[ras < 0] += 360

num_comps = len(ras)

c_ids = Column(data=np.array(['source']*num_comps), name='UNQ_SOURCE_ID', dtype='|S20')
c_names = Column(data=np.array([f'source_C{i:04d}' for i in range(num_comps)]), name='NAME', dtype='|S20')

##Component position
c_ras = Column(data=ras, name='RA')
c_decs = Column(data=decs, name='DEC')

##This says we have a point source
c_comp_types = Column(data=np.array(['P']*num_comps, dtype='|S1'), name="COMP_TYPE", dtype='|S1')
##This says we have a Stokes I power-law SED
c_mod_types = Column(data=np.array(['pl']*num_comps, dtype='|S3'), name="MOD_TYPE", dtype='|S3')

##Set everything as a flat spectrum 1 Jy stokes I source. That way we can image it and see the beam pattern
c_stokes_I_ref = Column(data=np.random.uniform(1, 10, num_comps), name='NORM_COMP_PL')
c_stokes_I_SI = Column(data=np.zeros(num_comps), name='ALPHA_PL')


cols = [c_ids, c_names, c_ras, c_decs, c_comp_types, c_mod_types, c_stokes_I_ref, c_stokes_I_SI]

main_table = Table(cols)

profile_cat = 'profiling_source.fits'

main_table.write(profile_cat, format='fits', overwrite=True)

cat_name = "profiling_source.fits"
low_freq = 50e+6
freq_reso = 100e+3
num_freq_chans = 1
uvfits_name = 'ska_image_test'
primary_beam = 'everybeam_OSKAR'

##The command to run WODEN
command = f'run_woden.py --ra0={ra0} --dec0={dec0} '
command += '--eb_point_to_phase ' ##this overrides whatever pointing is in beam_ms_path
                                  ##and sets it to ra0, dec0
command += f'--latitude={MWA_LAT} --longitude={MWA_LONG} '
command += f'--date={date} --output_uvfits_prepend={uvfits_name} '
command += f'--cat_filename={cat_name} --primary_beam={primary_beam} '
command += f'--lowest_channel_freq={low_freq} --freq_res={freq_reso} '
command += f'--num_freq_channels={num_freq_chans} --band_nums=1 '
command += f'--time_res=2 --num_time_steps=2 --IAU_order '
command += ' --cpu_mode '
# command += f' --station_id=0 '
command += f'--beam_ms_path=create_OSKAR-SKA_ms/OSKAR-SKA-layout.ms'

call(command, shell=True)

/home/jack-line/software/WODEN_dev/woden_dev/bin/run_woden.py:4: DeprecationWarning: pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html
  __import__('pkg_resources').require('wodenpy==2.5.0')

Successful readonly open of default-locked table create_OSKAR-SKA_ms/OSKAR-SKA-layout.ms: 22 columns, 130816 rows
Successful readonly open of default-locked table create_OSKAR-SKA_ms/OSKAR-SKA-layout.ms/SPECTRAL_WINDOW: 14 columns, 1 rows
Successful read/write open of default-locked table create_OSKAR-SKA_ms/OSKAR-SKA-layout.ms::FIELD: 9 columns, 1 rows
Successful readonly open of default-locked table create_OSKAR-SKA_ms/OSKAR-SKA-layout.ms: 22 columns, 130816 rows
Successful read/write open of default-locked table /home/jack-line/software/WODEN_dev/test_installation/everybeam/pointed_ska_image_test_band01.ms::FIELD: 9 columns, 1 rows
Successful readonly open of default-locked table create_OSKAR-SKA_ms/OSKAR-SKA-layout.ms/ANTENNA: 8 columns, 512 rows
2025-03-06 10:07:29 - INFO - 
                                              )  (              )  
                                  (  (     ( /(  )\ )        ( /(  
                                  )\))(   ')\())(()/(   (    )\()) 
                                 ((_)()\ )((_)\  /(_))  )\  ((_)\  
                                 _(())\_)() ((_)(_))_  ((_)  _((_) 
                                 \ \((_)/ // _ \ |   \ | __|| \| | 
                                  \ \/\/ /| (_) || |) || _| | .` | 
                                   \_/\_/  \___/ |___/ |___||_|\_| 
                                   
                                 You are using wodenpy version/git hash: a43c383
                                 
2025-03-06 10:07:29 - INFO - Input arguments after parsing:
                             	Phase centre: 0.00455, -26.70332 deg
                             	Array central latitude: -26.703 deg
                             	Array central longitude: 116.671 deg
                             	Array central height: 377.827 m
                             	Lowest channel frequency: 50.000 MHz
                             	Channel frequency resolution: 100.000 kHz
                             	Coarse band bandwidth: 0.040 MHz
                             	Num channels per coarse band: 1
                             	Start date: 2024-07-21T20:13:00
                             	Time resolution: 2.0 (s)
                             	Num time steps: 2
                             	Have read 512 antenna positions from measurement set: create_OSKAR-SKA_ms/OSKAR-SKA-layout.ms
                             	Will write outputs to: /home/jack-line/software/WODEN_dev/test_installation/everybeam
2025-03-06 10:07:30 - INFO - Obs epoch initial LST was 0.0087300922 deg
2025-03-06 10:07:30 - INFO - Setting initial J2000 LST to 359.6931723607 deg
2025-03-06 10:07:30 - INFO - Setting initial mjd to 60512.8423726853
2025-03-06 10:07:30 - INFO - After precession initial latitude of the array is -26.8398178265 deg
2025-03-06 10:07:30 - INFO - Precessing array layout to J2000
2025-03-06 10:07:30 - INFO - Will run with EveryBeam OSKAR primary beam, based on this measurement set:
                             	create_OSKAR-SKA_ms/OSKAR-SKA-layout.ms
                             Created the following minimal MS to point the beam:
                             	/home/jack-line/software/WODEN_dev/test_installation/everybeam/pointed_ska_image_test_band01.ms
                             Primary beam is pointed at RA,Dec = 0.0045520088542088475, -26.703319405555554 deg

2025-03-06 10:07:30 - INFO - Doing the initial mapping of sky model
INFO: couldn't find second table containing shapelet information, so not attempting to load any shapelets.
2025-03-06 10:07:30 - INFO - Sky model mapping took 0.0 mins
2025-03-06 10:07:30 - INFO - Have read in 500 components
2025-03-06 10:07:30 - INFO - After cropping there are 500 components
2025-03-06 10:07:30 - INFO - Will load libwoden from /home/jack-line/software/WODEN_dev/wodenpy/libwoden_double.so
2025-03-06 10:07:30 - INFO - Running in parallel on CPU mode with 8 threads
                             There are 1 sets of sky models to run
2025-03-06 10:07:30 - INFO - Reading set 0 sky models
2025-03-06 10:07:30 - INFO - From sky set 0 thread num 0 reading 63 points, 0 gauss, 0 shape, 0 shape coeffs
2025-03-06 10:07:30 - INFO - Finshed sky set 0 reading thread num 0 in 0.0 seconds
2025-03-06 10:07:30 - INFO - From sky set 0 thread num 1 reading 63 points, 0 gauss, 0 shape, 0 shape coeffs
2025-03-06 10:07:30 - INFO - Finshed sky set 0 reading thread num 1 in 0.0 seconds
2025-03-06 10:07:30 - INFO - From sky set 0 thread num 2 reading 63 points, 0 gauss, 0 shape, 0 shape coeffs
2025-03-06 10:07:30 - INFO - From sky set 0 thread num 4 reading 63 points, 0 gauss, 0 shape, 0 shape coeffs
2025-03-06 10:07:30 - INFO - Finshed sky set 0 reading thread num 2 in 0.0 seconds
2025-03-06 10:07:30 - INFO - Finshed sky set 0 reading thread num 4 in 0.0 seconds
2025-03-06 10:07:30 - INFO - From sky set 0 thread num 3 reading 63 points, 0 gauss, 0 shape, 0 shape coeffs
2025-03-06 10:07:30 - INFO - Finshed sky set 0 reading thread num 3 in 0.0 seconds
2025-03-06 10:07:30 - INFO - From sky set 0 thread num 5 reading 63 points, 0 gauss, 0 shape, 0 shape coeffs
2025-03-06 10:07:30 - INFO - From sky set 0 thread num 6 reading 63 points, 0 gauss, 0 shape, 0 shape coeffs
2025-03-06 10:07:30 - INFO - Finshed sky set 0 reading thread num 5 in 0.0 seconds
2025-03-06 10:07:30 - INFO - Finshed sky set 0 reading thread num 6 in 0.0 seconds
2025-03-06 10:07:30 - INFO - From sky set 0 thread num 7 reading 59 points, 0 gauss, 0 shape, 0 shape coeffs
2025-03-06 10:07:30 - INFO - Finshed sky set 0 reading thread num 7 in 0.0 seconds
2025-03-06 10:07:31 - INFO - Sending Sky set 0 chunk 0 to the CPU
2025-03-06 10:07:31 - INFO - Sending Sky set 0 chunk 1 to the CPU
2025-03-06 10:07:32 - INFO - Sending Sky set 0 chunk 2 to the CPU
2025-03-06 10:07:32 - INFO - Sending Sky set 0 chunk 3 to the CPU
POINTING LOFAR beam 0.0 -26.7
2025-03-06 10:07:32 - INFO - Sending Sky set 0 chunk 4 to the CPU
POINTING LOFAR beam 0.0 -26.7
2025-03-06 10:07:32 - INFO - Sending Sky set 0 chunk 5 to the CPU

2025-03-06 00:07:32	SEVERE	MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)	Leap second table TAI_UTC seems out-of-date.
2025-03-06 00:07:32	SEVERE	MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+	Until the table is updated (see the CASA documentation or your system admin),
2025-03-06 00:07:32	SEVERE	MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+	times and coordinates derived from UTC could be wrong by 1s or more.
2025-03-06 00:07:33	SEVERE	MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)	Leap second table TAI_UTC seems out-of-date.
2025-03-06 00:07:33	SEVERE	MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+	Until the table is updated (see the CASA documentation or your system admin),
2025-03-06 00:07:33	SEVERE	MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+	times and coordinates derived from UTC could be wrong by 1s or more.

2025-03-06 10:07:32 - INFO - Sending Sky set 0 chunk 6 to the CPU
POINTING LOFAR beam 0.0 -26.7
POINTING LOFAR beam 0.0 -26.7
2025-03-06 10:07:32 - INFO - Sending Sky set 0 chunk 7 to the CPU

2025-03-06 00:07:33	SEVERE	MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)	Leap second table TAI_UTC seems out-of-date.
2025-03-06 00:07:33	SEVERE	MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+	Until the table is updated (see the CASA documentation or your system admin),
2025-03-06 00:07:33	SEVERE	MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+	times and coordinates derived from UTC could be wrong by 1s or more.
2025-03-06 00:07:33	SEVERE	MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)	Leap second table TAI_UTC seems out-of-date.
2025-03-06 00:07:33	SEVERE	MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+	Until the table is updated (see the CASA documentation or your system admin),
2025-03-06 00:07:33	SEVERE	MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+	times and coordinates derived from UTC could be wrong by 1s or more.

POINTING LOFAR beam 0.0 -26.7
POINTING LOFAR beam 0.0 -26.7

2025-03-06 00:07:33	SEVERE	MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)	Leap second table TAI_UTC seems out-of-date.
2025-03-06 00:07:33	SEVERE	MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+	Until the table is updated (see the CASA documentation or your system admin),
2025-03-06 00:07:33	SEVERE	MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+	times and coordinates derived from UTC could be wrong by 1s or more.
2025-03-06 00:07:33	SEVERE	MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)	Leap second table TAI_UTC seems out-of-date.
2025-03-06 00:07:33	SEVERE	MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+	Until the table is updated (see the CASA documentation or your system admin),
2025-03-06 00:07:33	SEVERE	MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+	times and coordinates derived from UTC could be wrong by 1s or more.

POINTING LOFAR beam 0.0 -26.7
POINTING LOFAR beam 0.0 -26.7

2025-03-06 00:07:34	SEVERE	MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)	Leap second table TAI_UTC seems out-of-date.
2025-03-06 00:07:34	SEVERE	MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+	Until the table is updated (see the CASA documentation or your system admin),
2025-03-06 00:07:34	SEVERE	MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+	times and coordinates derived from UTC could be wrong by 1s or more.
2025-03-06 00:07:34	SEVERE	MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)	Leap second table TAI_UTC seems out-of-date.
2025-03-06 00:07:34	SEVERE	MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+	Until the table is updated (see the CASA documentation or your system admin),
2025-03-06 00:07:34	SEVERE	MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+	times and coordinates derived from UTC could be wrong by 1s or more.

2025-03-06 10:18:42 - INFO - Sky set 0 chunk 1 has returned from the CPU after 671.1 seconds
2025-03-06 10:18:50 - INFO - Sky set 0 chunk 6 has returned from the CPU after 677.2 seconds
2025-03-06 10:18:50 - INFO - Sky set 0 chunk 7 has returned from the CPU after 677.4 seconds
2025-03-06 10:19:02 - INFO - Sky set 0 chunk 5 has returned from the CPU after 689.5 seconds
2025-03-06 10:19:09 - INFO - Sky set 0 chunk 0 has returned from the CPU after 698.0 seconds
2025-03-06 10:19:10 - INFO - Sky set 0 chunk 2 has returned from the CPU after 698.0 seconds
2025-03-06 10:19:10 - INFO - Sky set 0 chunk 3 has returned from the CPU after 698.3 seconds
2025-03-06 10:19:14 - INFO - Sky set 0 chunk 4 has returned from the CPU after 701.7 seconds
2025-03-06 10:19:14 - INFO - Have completed 500 of 500 components calcs (100.0%)
2025-03-06 10:19:14 - INFO - Finished all rounds of processing
2025-03-06 10:19:15 - INFO - Deleting /home/jack-line/software/WODEN_dev/test_installation/everybeam/pointed_ska_image_test_band01.ms...
2025-03-06 10:19:15 - INFO - Done
2025-03-06 10:19:15 - INFO - Full run took 0:11:47.560312
2025-03-06 10:19:15 - INFO - WODEN is done. Closing the log. S'later

0

Note I’ve deleted the outputs of the cell above, because EveryBeam spits out a bunch of error messages with the OSKAR beam warning about frequencies mismatching, and not having a central direction set. The outputs seems fine. It’s a TODO to silence these warnings.

Also note that this simulation took 281 minutes on my desktop.

cmd = "woden_uv2ms.py "
cmd += "  --uvfits_prepend=ska_image_test_band "
cmd += "  --band_nums=1  "

call(cmd, shell=True)

cmd = "wsclean -name ska_image_test -size 4096 4096 -niter 2000 "
cmd += "  -auto-threshold 0.5 -auto-mask 3 "
cmd += "  -pol I -multiscale -weight briggs 0 -scale 0.005 -j 12 -mgain 0.85 "
cmd += "  -no-update-model-required "
cmd += "  ska_image_test_band*.ms "

call(cmd, shell=True)

/home/jack-line/software/WODEN_dev/woden_dev/bin/woden_uv2ms.py:4: DeprecationWarning: pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html
  __import__('pkg_resources').require('wodenpy==2.5.0')
The telescope frame is set to '????', which generally indicates ignorance. Defaulting the frame to 'itrs', but this may lead to other warnings or errors.
The uvw_array does not match the expected values given the antenna positions. The largest discrepancy is 4.7435968157024035 meters. This is a fairly common situation but might indicate an error in the antenna positions, the uvws or the phasing.
The uvw_array does not match the expected values given the antenna positions. The largest discrepancy is 4.7435968157024035 meters. This is a fairly common situation but might indicate an error in the antenna positions, the uvws or the phasing.
Writing in the MS file that the units of the data are uncalib, although some CASA process will ignore this and assume the units are all in Jy (or may not know how to handle data in these units).

WSClean version 3.4 (2023-10-11)
This software package is released under the GPL version 3.
Author: André Offringa (offringa@gmail.com).

No corrected data in first measurement set: tasks will be applied on the data column.
=== IMAGING TABLE ===
       # Pol Ch JG ²G FG FI In Freq(MHz)
| Independent group:
+-+-J- 0  I   0  0  0  0  0  0  50-50 (1)
Reordering ska_image_test_band01.ms into 1 x 1 parts.
Reordering: 0%....10%....20%....30%....40%....50%....60%....70%....80%....90%....100%
Initializing model visibilities: 0%....10%....20%....30%....40%....50%....60%....70%....80%....90%....100%
 == Constructing PSF ==
Precalculating weights for Briggs'(0) weighting...
Opening reordered part 0 spw 0 for ska_image_test_band01.ms
Detected 62.7 GB of system memory, usage not limited.
Opening reordered part 0 spw 0 for ska_image_test_band01.ms
Determining min and max w & theoretical beam size... DONE (w=[3.24584e-06:47.8641] lambdas, maxuvw=7253.73 lambda)
Theoretic beam = 28.44''
Small inversion enabled, but inversion resolution already smaller than beam size: not using optimization.
Loading data in memory...
Gridding 261632 rows...
Gridded visibility count: 256664
Fitting beam... major=110.32'', minor=73.99'', PA=136.76 deg, theoretical=28.44''.
Writing psf image... DONE
 == Constructing image ==
Opening reordered part 0 spw 0 for ska_image_test_band01.ms
Loading data in memory...
Gridding 261632 rows...
Gridded visibility count: 256664
Writing dirty image...
 == Deconvolving (1) ==
Estimated standard deviation of background noise: 100.34 mJy
Scale info:
- Scale 0, bias factor=1, psfpeak=1, gain=0.1, kernel peak=0.150298
- Scale 6, bias factor=1.7, psfpeak=0.445397, gain=0.224519, kernel peak=0.0501087
- Scale 13, bias factor=2.8, psfpeak=0.332387, gain=0.300854, kernel peak=0.0184458
- Scale 25, bias factor=4.6, psfpeak=0.219276, gain=0.456047, kernel peak=0.00588155
- Scale 51, bias factor=7.7, psfpeak=0.118502, gain=0.843867, kernel peak=0.00157286
- Scale 101, bias factor=12.9, psfpeak=0.0520027, gain=1.92298, kernel peak=0.00042186
- Scale 202, bias factor=21.4, psfpeak=0.0159038, gain=6.28782, kernel peak=0.000107389
- Scale 404, bias factor=35.7, psfpeak=0.00402302, gain=24.857, kernel peak=2.73429e-05
- Scale 809, bias factor=59.5, psfpeak=0.00096395, gain=103.74, kernel peak=6.90052e-06
- Scale 1618, bias factor=99.2, psfpeak=0.000191813, gain=521.342, kernel peak=1.7359e-06
RMS per scale: {0: 114.08 mJy, 6: 69.79 mJy, 13: 60.1 mJy, 25: 48.03 mJy, 51: 34.8 mJy, 101: 22.33 mJy, 202: 11.68 mJy, 404: 5.57 mJy, 809: 2.71 mJy, 1618: 1.11 mJy}
Starting multi-scale cleaning. Start peak=9.11 Jy, major iteration threshold=1.37 Jy
Iteration 7, scale 0 px : 8.1 Jy at 2031,1872
Iteration 14, scale 0 px : 6.42 Jy at 2114,1598
Iteration 35, scale 25 px : 5.79 Jy at 2502,2397
Iteration 79, scale 25 px : 4.59 Jy at 2617,1278
Iteration 189, scale 0 px : 4.71 Jy at 2031,1872
Iteration 239, scale 0 px : 3.76 Jy at 2197,2656
Iteration 355, scale 25 px : 3.3 Jy at 1278,1915
Iteration 459, scale 0 px : 3 Jy at 2503,2397
Iteration 614, scale 0 px : 2.4 Jy at 1801,1527
Iteration 826, scale 25 px : 2.27 Jy at 1741,2938
Iteration 878, scale 25 px : 1.82 Jy at 3020,2906
Iteration 1077, scale 0 px : 1.9 Jy at 2330,1890
Iteration 1294, scale 0 px : 1.52 Jy at 2502,1621
Subminor loop is near minor loop threshold. Initiating countdown.
(14) Iteration 1430, scale 25 px : 1.46 Jy at 2351,1022
(13) Iteration 1442, scale 0 px : 1.37 Jy at 3016,1773
Minor loop finished, continuing cleaning after inversion/prediction round.
Assigning from 1 to 1 channels...
 == Converting model image to visibilities ==
Opening reordered part 0 spw 0 for ska_image_test_band01.ms
Predicting 261632 rows...
Writing...
 == Constructing image ==
Opening reordered part 0 spw 0 for ska_image_test_band01.ms
Loading data in memory...
Gridding 261632 rows...
Gridded visibility count: 256664
 == Deconvolving (2) ==
Estimated standard deviation of background noise: 54.34 mJy
Scale info:
- Scale 0, bias factor=1, psfpeak=1, gain=0.1, kernel peak=0.150298
- Scale 6, bias factor=1.7, psfpeak=0.445397, gain=0.224519, kernel peak=0.0501087
- Scale 13, bias factor=2.8, psfpeak=0.332387, gain=0.300854, kernel peak=0.0184458
- Scale 25, bias factor=4.6, psfpeak=0.219276, gain=0.456047, kernel peak=0.00588155
- Scale 51, bias factor=7.7, psfpeak=0.118502, gain=0.843867, kernel peak=0.00157286
- Scale 101, bias factor=12.9, psfpeak=0.0520027, gain=1.92298, kernel peak=0.00042186
- Scale 202, bias factor=21.4, psfpeak=0.0159038, gain=6.28782, kernel peak=0.000107389
- Scale 404, bias factor=35.7, psfpeak=0.00402302, gain=24.857, kernel peak=2.73429e-05
- Scale 809, bias factor=59.5, psfpeak=0.00096395, gain=103.74, kernel peak=6.90052e-06
- Scale 1618, bias factor=99.2, psfpeak=0.000191813, gain=521.342, kernel peak=1.7359e-06
RMS per scale: {0: 58.55 mJy, 6: 25.61 mJy, 13: 20.03 mJy, 25: 14.67 mJy, 51: 10.12 mJy, 101: 6.41 mJy, 202: 3.47 mJy, 404: 1.92 mJy, 809: 1.15 mJy, 1618: 673.43 µJy}
Starting multi-scale cleaning. Start peak=1.39 Jy, major iteration threshold=208.45 mJy
Iteration 1733, scale 25 px : 1.34 Jy at 2982,1995
Iteration 1765, scale 25 px : 1.07 Jy at 1868,712
Iteration 1914, scale 0 px : 1.11 Jy at 1278,1914
Iteration 2000, scale 0 px : 1.04 Jy at 2619,3227
Cleaning finished because maximum number of iterations was reached.
Auto-masking threshold reached; continuing next major iteration with deeper threshold and mask.
Maximum number of minor deconvolution iterations was reached: not continuing deconvolution.
Assigning from 1 to 1 channels...
Writing model image...
 == Converting model image to visibilities ==
Opening reordered part 0 spw 0 for ska_image_test_band01.ms
Predicting 261632 rows...
Writing...
 == Constructing image ==
Opening reordered part 0 spw 0 for ska_image_test_band01.ms
Loading data in memory...
Gridding 261632 rows...
Gridded visibility count: 256664
2 major iterations were performed.
Rendering sources to restored image (beam=73.99''-110.32'', PA=136.76 deg)... DONE
Writing restored image... DONE
Multi-scale cleaning summary:
- Scale 0 px, nr of components cleaned: 1291 (263.53 Jy)
- Scale 6 px, nr of components cleaned: 0 (0 Jy)
- Scale 13 px, nr of components cleaned: 0 (0 Jy)
- Scale 25 px, nr of components cleaned: 709 (162.74 Jy)
- Scale 51 px, nr of components cleaned: 0 (0 Jy)
- Scale 101 px, nr of components cleaned: 0 (0 Jy)
- Scale 202 px, nr of components cleaned: 0 (0 Jy)
- Scale 404 px, nr of components cleaned: 0 (0 Jy)
- Scale 809 px, nr of components cleaned: 0 (0 Jy)
- Scale 1618 px, nr of components cleaned: 0 (0 Jy)
Total: 2000 components (426.27 Jy)
Inversion: 00:00:18.661843, prediction: 00:00:07.969461, deconvolution: 00:01:24.374294
Cleaning up temporary files...

0

with fits.open('ska_image_test-image.fits') as hdu:
    image = np.squeeze(hdu[0].data)
    wcs = WCS(hdu[0].header).celestial
    
fig, axs = plt.subplots(1, 1, figsize=(12, 8), subplot_kw={'projection': wcs})

im = axs.imshow(image, origin='lower', cmap='Blues_r', vmin=0, vmax=1.0)
plt.colorbar(im, ax=axs, label='Jy/beam')

half_width = 2000
mid = 2048

axs.set_xlim(mid-half_width, mid+half_width)
axs.set_ylim(mid-half_width, mid+half_width)

plt.show()

WARNING: FITSFixedWarning: 'datfix' made the change 'Set MJD-OBS to 60512.842372 from DATE-OBS'. [astropy.wcs.wcs]

Pretty incredible image quality, given this is a single frequency channel, and only 2 time steps. We’ll overplot the source posisitons to show that the beam is properly down-weighting the sources at image edge.

fig, axs = plt.subplots(1, 1, figsize=(12, 8), subplot_kw={'projection': wcs})


im = axs.imshow(image, origin='lower', cmap='Blues_r', vmin=0, vmax=1.0)
plt.colorbar(im, ax=axs, label='Jy/beam')

x, y = wcs.all_world2pix(ras, decs, 0)
axs.scatter(x, y, s=2, c='C1', marker='x', alpha=0.3)

half_width = 2000
mid = 2048

axs.set_xlim(mid-half_width, mid+half_width)
axs.set_ylim(mid-half_width, mid+half_width)

plt.show()

Out of interest, what does the PSF look like?

with fits.open('ska_image_test-psf.fits') as hdu:
    image = np.squeeze(hdu[0].data)
    wcs = WCS(hdu[0].header).celestial
    
fig, axs = plt.subplots(1, 1, figsize=(12, 8), subplot_kw={'projection': wcs})

im = axs.imshow(image, origin='lower', cmap='Blues_r', vmin=0, vmax=1.0)
plt.colorbar(im, ax=axs, label='Jy/beam')

half_width = 100
mid = 2048

axs.set_xlim(mid-half_width, mid+half_width)
axs.set_ylim(mid-half_width, mid+half_width)

plt.show()

WARNING: FITSFixedWarning: 'datfix' made the change 'Set MJD-OBS to 0.000000 from DATE-OBS'. [astropy.wcs.wcs]

Like a snowflake, pretty cool.