LOFAR LBA NCP¶
In this example, we’ll make an observation of the North Celestial Pole (NCP) using the Low Band Antennas (LBA) of the LOFAR radio telescope. I think the LOFAR EoR group often observe towards the NCP, so it’s a good example to try. I couldn’t easily find a LOFAR catalogue that covers the NCP however, so create a sky model using the WENSS and NVSS catalogues, to try and get realistic spectral indicies.
Then we’ll decided the morphology of the WENSS cat isn’t good enough, and we’ll just use NVSS. But we’ll leave all the cross-matching stuff in there cos it’s a useful reference, and otherwise I wasted my time writing it.
Sky model WENSS x NVSS¶
The WENSS catalogue is a 325 MHz survey of the NCP, so sits at a frequency close to where the LBA operates. It’s relatively low resolution compared to what LBA is capable of at 54 arcsecs, but that will make the example simulation less computationally intensive which is sweet. We want to include spectra indicies if we can, so we’ll cross-match it to the prolific NVSS catalogue (NVSS paper has a whopping 5000+ citations). NVSS is at 1.4 GHz, and we’ll just assume everything is a power law between the two frequencies.
Using the postage stamp server SkyView, I’ve downloaded a WENSS image of the NCP at 325 MHz (setting desired coords to 0,90 deg and asking for a 5 deg image with 1024 pixels). This is what we want to try and emulate:
from astropy.io import fits
from astropy.wcs import WCS
import matplotlib.pyplot as plt
with fits.open('wenss_sky.fits') as hdu:
data_wenss = hdu[0].data
header_wenss = hdu[0].header
wcs_wenss = WCS(hdu[0].header)
img_reso = hdu[0].header['CDELT2']
print(img_reso)
fig, axs = plt.subplots(1, 1, figsize=(10, 10), subplot_kw={'projection': wcs_wenss})
axs.imshow(data_wenss, origin='lower', cmap='Blues', vmin=0, vmax=0.1)
axs.grid()
axs.set_xlabel('RA')
axs.set_ylabel('Dec')
plt.show()
0.005859375
I’ve downloaded the two catalogues from VizieR. For those unfamiliar, type WENSS into the free text search. Click on the ViZieR button next to the The Westerbork Northern Sky Survey (Leiden, 1998) result that pops up. You’ll be given a set of columns to choose to include. On the left, in the Preferences toolbar, change ‘max’ to unlimited (which will download all sources) and change HTML table into a “FITS binary” table, so download the FITS catalogue. Finally, click “Submit”. Do something similar for the NVSS catalogue. Note that the PA columns in NVSS isn’t included by default, and I use it, so make sure to include that.
The first thing I’ll do is crop them to a 3 degree radius around the NCP. Note that you can actually crop the catalogues on Vizier itself if you don’t want to download the entire tables. I find it useful to have the complete catalogues lying around though. Up to you.
from astropy.table import Table
def crop_table(name):
table = Table.read(name)
decs = table['DEJ2000']
crop_table = table[decs > 87.0]
crop_table.write(name.replace('.fits', '_crop.fits'), overwrite=True)
return crop_table
crop_wenss = crop_table('vizier_wenss.fits')
crop_nvss = crop_table('vizier_nvss.fits')
print(len(crop_wenss), len(crop_nvss))
642 1272
OK, so we’ve ended up with 642 WENSS sources. Let’s cross match them with the NVSS sources. Normally I’d just use TOPCAT to do this graphically. Check out this tutorial YouTube video if you’re unfamiliar. The video seems to have happened during the Covid-19 years, so if you lived through them you can have lovely flashbacks to endless online events. I’ll use the astropy package to do this instead as I can do it in this notebook.
from astropy.coordinates import SkyCoord
from astropy import units as u
import numpy as np
ra_wenss = crop_wenss['RAJ2000'].data
dec_wenss = crop_wenss['DEJ2000'].data
ra_nvss = crop_nvss['RAJ2000'].data
dec_nvss = crop_nvss['DEJ2000'].data
coords_wenss = SkyCoord(ra=ra_wenss*u.degree, dec=dec_wenss*u.degree)
coords_nvss = SkyCoord(ra=ra_nvss*u.degree, dec=dec_nvss*u.degree)
nvss_indexes, d2d, d3d = coords_wenss.match_to_catalog_sky(coords_nvss)
OK, so nvss_indexes returns the indexes of NVSS sources that are the closest match to each WENSS source. d2d is the angular separation between the two sources, so if you want to chuck away anything that isnit nearby on the sky you can do that. Normally something like the resolution of the base catalogue you are matching to is a good start. This is why TOPCAT is sweet, because you can quickly try different values and plot them up quickly. Anyways, we’ll just check that a) all of the cross-matches are unique, and b) none of the matches are too far away.
print("Number of WENSS sources: ", len(crop_wenss))
print("Number of unique matches: ", len(np.unique(nvss_indexes)))
Number of WENSS sources: 642
Number of unique matches: 603
Huh OK so we don’t have matches to everything in WENSS. First thing to do is to discard any repeated matches. We’ll just go for the closest match in the case of repeats.
##Sort everything by angular distance
argorder = np.argsort(d2d)
wenss_inds_sorted = argorder
nvss_inds_sorted = nvss_indexes[argorder]
nvss_match = np.empty(len(wenss_inds_sorted))
for wenss_ind, nvss_ind in zip(wenss_inds_sorted, nvss_inds_sorted):
if nvss_ind in nvss_match:
nvss_match[wenss_ind] = np.nan
else:
nvss_match[wenss_ind] = int(nvss_ind)
Let’s plot up matches to see what’s going on. We’ll use a WCS so we can plot them up on the sky.
from astropy.wcs import WCS
from astropy.io import fits
# Define the FITS header template
header = fits.Header()
nside = 256
reso = 6 / nside
# Basic WCS parameters
header['NAXIS'] = 2 # Number of axes
header['NAXIS1'] = nside # Reference pixel along X-axis
header['NAXIS2'] = nside # Reference pixel along Y-axis
header['CRPIX1'] = nside // 2 # Reference pixel along X-axis
header['CRPIX2'] = nside // 2 # Reference pixel along Y-axis
header['CRVAL1'] = 0.0 # Reference RA (degrees)
header['CRVAL2'] = 90.0 # Reference Dec (degrees)
header['CTYPE1'] = 'RA---SIN' # RA projection type
header['CTYPE2'] = 'DEC--SIN' # Dec projection type
header['CDELT1'] = -reso # Pixel scale in RA (degrees/pixel)
header['CDELT2'] = reso # Pixel scale in Dec (degrees/pixel)
header['CUNIT1'] = 'deg' # Units for RA
header['CUNIT2'] = 'deg' # Units for Dec
# Optionally, set the equinox and epoch
header['EQUINOX'] = 2000.0 # Equatorial coordinates J2000
header['RADESYS'] = 'ICRS' # Coordinate reference system
# Create a WCS object from the header
wcs = WCS(header)
fig, axs = plt.subplots(1, 1, figsize=(7, 7), subplot_kw={'projection': wcs})
nvss_match_without_nan = nvss_match[~np.isnan(nvss_match)].astype(int)
axs.plot(ra_wenss, dec_wenss, 'C1o', mfc='none', transform=axs.get_transform('world'), label='WENSS')
axs.plot(ra_nvss[nvss_match_without_nan], dec_nvss[nvss_match_without_nan], 'C1o', alpha=0.5,
transform=axs.get_transform('world'), label='NVSS matched')
axs.plot(ra_nvss, dec_nvss, 'k.', ms=2, transform=axs.get_transform('world'), label='NVSS')
axs.legend()
plt.show()
Eh, we have a couple of places where two NVSS sources are close to a WENSS source, and just few instances where there is no match in NVSS. As this is just a simple example, we’ll just use everything that has a match and call it good enough. If you were doing this for real, you could always use try and use (PUMA)[https://github.com/JLBLine/PUMA] to do a more sophisticated cross-match SHAMELESS PLUG OVER.
Right, let’s take our matches and make a WODEN style sky model. First, grab all the relevant information, calculate spectral indices, check they look sensible.
from astropy.table import Column
##subsets we need
wenss_indexes = ~np.isnan(nvss_match)
use_nvss = nvss_match[~np.isnan(nvss_match)].astype(int)
# nvss_match_without_nan = nvss_match[~np.isnan(nvss_match)].astype(int)
unq_source_ids = np.array(crop_wenss['Name'][wenss_indexes], dtype='U32')
##order everything by source name to make adding component names easier
ras_wenss = crop_wenss['RAJ2000'][wenss_indexes].data
decs_wenss = crop_wenss['DEJ2000'][wenss_indexes].data
fluxes_330 = crop_wenss['Sint'][wenss_indexes].data/1000.0
ras_nvss = crop_nvss['RAJ2000'][use_nvss].data
decs_nvss = crop_nvss['DEJ2000'][use_nvss].data
ras_nvss_full = crop_nvss['RAJ2000'].data
decs_nvss_full = crop_nvss['DEJ2000'].data
##use the gaussians from NVSS as it has better resolution
fluxes_1400 = crop_nvss['S1_4'][use_nvss].data/1000.0
majors = crop_nvss['MajAxis'][use_nvss].data/3600.0
minors = crop_nvss['MinAxis'][use_nvss].data/3600.0
pas = np.array(crop_nvss['PA'][use_nvss], dtype=np.float64)
pas[np.isnan(pas)] = 0.0
alphas = np.log10(fluxes_1400/fluxes_330)/np.log10(1400.0/330.0)
##WODEN references all power law fluxes to 200MHz
fluxes_extrap = fluxes_330*(200.0/330.0)**alphas
plt.hist(alphas, bins=20, histtype='step')
plt.xlabel('Spectral index')
plt.ylabel('Number of sources')
plt.show()
fig, axs = plt.subplots(1, 1, figsize=(7, 7), subplot_kw={'projection': wcs})
axs.plot(ras_nvss, decs_nvss, 'C1o', alpha=0.5, transform=axs.get_transform('world'), label='NVSS matched')
axs.plot(ras_wenss, decs_wenss, 'C0o', mfc='none', transform=axs.get_transform('world'), label='WENSS')
axs.plot(ras_nvss_full, decs_nvss_full, 'k.', ms=2, alpha=0.5, transform=axs.get_transform('world'), label='NVSS full')
edge = 60
axs.set_xlim(nside//2-edge, nside//2+edge)
axs.set_ylim(nside//2-edge, nside//2+edge)
axs.legend()
plt.show()
Yeah looks like we have a couple of mismatches but the histogram of the spectral index is good enough. Let’s just crack on.
num_comps = len(unq_source_ids)
unq_source_id_col = Column(unq_source_ids, name='UNQ_SOURCE_ID')
names = np.array([name+'_C000' for name in unq_source_ids], dtype='U32')
names_col = Column(names, name='NAME')
ras_col = Column(ras_wenss, name='RA', unit='deg')
decs_col = Column(decs_wenss, name='DEC', unit='deg')
fluxes_col = Column(fluxes_extrap, name='NORM_COMP_PL', unit='Jy')
si_col = Column(alphas, name='ALPHA_PL', unit='Jy')
majors_col = Column(majors, name='MAJOR_DC', unit='deg')
minors_col = Column(minors, name='MINOR_DC', unit='deg')
pas_col = Column(pas, name='PA_DC', unit='deg')
##Anything that has a major axis is a Gaussian, other make it a point source
comp_types = np.empty(unq_source_ids.shape, dtype='U3')
comp_types[majors == 0.0] = 'P'
comp_types[majors != 0.0] = 'G'
comp_types_col = Column(comp_types, name='COMP_TYPE')
##We're inputting everything as a power law
mod_types_col = Column(np.full(num_comps, 'pl'), name='MOD_TYPE')
cols = [unq_source_id_col, names_col,
ras_col, decs_col, comp_types_col, mod_types_col,
fluxes_col, si_col,
majors_col, minors_col, pas_col]
new_table = Table(cols)
new_table.write('srclist_wenss_x_nvss.fits', overwrite=True)
Finally, let’s check out selected subset positions match the image we’re trying to reproduce.
fig, axs = plt.subplots(1, 1, figsize=(10, 10), subplot_kw={'projection': wcs_wenss})
axs.imshow(data_wenss, origin='lower', cmap='Blues', vmin=0, vmax=0.1)
x, y = wcs_wenss.all_world2pix(ras_wenss, decs_wenss, 0)
# axs.plot(ras_wenss, decs_wenss, 'C1o', mfc='none', transform=axs.get_transform('world'), label='WENSS')
axs.plot(x, y, 'C1o', mfc='none', label='WENSS')
axs.grid()
axs.set_xlabel('RA')
axs.set_ylabel('Dec')
plt.show()
Plotting this at the NCP¶
Eh, what’s this? All of our positions are messed up! It turns out you need to specify what counts as “up” in the WCS when the NCP is at the field centre. You can fix that by editing the header of the image. This is one way to do it:
from astropy.io import fits
# Load the FITS file
with fits.open("wenss_sky.fits", mode="update") as hdul:
header = hdul[0].header
# Add or update the LONPOLE keyword
header["LONPOLE"] = 180 # Set to appropriate value (e.g., 180 for common celestial cases)
hdul[0].header = header
# Save changes
hdul.writeto("wenss_sky_edit.fits", overwrite=True)
with fits.open('wenss_sky_edit.fits') as hdu:
data_wenss = hdu[0].data
header = hdu[0].header
wcs_wenss = WCS(header)
fig, axs = plt.subplots(1, 1, figsize=(10, 10), subplot_kw={'projection': wcs_wenss})
axs.imshow(data_wenss, origin='lower', cmap='Blues', vmin=0, vmax=0.1)
x, y = wcs_wenss.all_world2pix(ras_wenss, decs_wenss, 0)
# axs.plot(ras_wenss, decs_wenss, 'C1o', mfc='none', transform=axs.get_transform('world'), label='WENSS')
axs.plot(x, y, 'C1o', mfc='none', label='WENSS')
axs.grid()
axs.set_xlabel('RA')
axs.set_ylabel('Dec')
plt.show()
and now it’s all good. Some FITS files are missing this keyword, some aren’t, so just be vigilant and sanity check your WCS.
Run the simulation¶
Ah, there go, source positions match up now. Let’s run the simulation. Seeing as though this is the NCP, which is always visibile from LOFAR, we can do so mega \(uv\)-synthesis. We’ll do a 24 hour obs, sampling once every hour. We’ll do a 20MHz bandwidth, with 1MHz resoultion.
%%bash
primary_beam="everybeam_LOFAR"
ra0=0.0
dec0=90.0
lofar_lat=52.905329712
lofar_long=6.867996528
date="2014-10-16T17:01:36.000"
ms_path=../../test_installation/everybeam/lba.MS
uvfits_name=woden_LBA-NCP
cat_name=srclist_wenss_x_nvss.fits
freq_reso=1e+6
low_freq=75e+6
time_res=3600.0
num_freq_chans=20
num_time_steps=24
run_woden.py \
--ra0=${ra0} \
--dec0=${dec0} \
--eb_point_to_phase \
--latitude=${lofar_lat} \
--longitude=${lofar_long} \
--date=${date} \
--output_uvfits_prepend=${uvfits_name} \
--cat_filename=${cat_name} \
--primary_beam=${primary_beam} \
--lowest_channel_freq=${low_freq} \
--freq_res=${freq_reso} \
--num_freq_channels=${num_freq_chans} \
--band_nums=1 \
--time_res=${time_res} \
--num_time_steps=${num_time_steps} \
--beam_ms_path=${ms_path}
/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 ../../test_installation/everybeam/lba.MS: 24 columns, 4218 rows
Successful readonly open of default-locked table ../../test_installation/everybeam/lba.MS/SPECTRAL_WINDOW: 14 columns, 1 rows
Successful read/write open of default-locked table ../../test_installation/everybeam/lba.MS::FIELD: 10 columns, 1 rows
Successful readonly open of default-locked table ../../test_installation/everybeam/lba.MS: 24 columns, 4218 rows
Successful read/write open of default-locked table /home/jack-line/software/WODEN_dev/examples/LOFAR_LBA_NCP/pointed_woden_LBA-NCP_band01.ms::FIELD: 10 columns, 1 rows
Successful readonly open of default-locked table ../../test_installation/everybeam/lba.MS/ANTENNA: 10 columns, 37 rows
2025-03-17 10:19:42 - INFO -
) ( )
( ( ( /( )\ ) ( /(
)\))( ')\())(()/( ( )\())
((_)()\ )((_)\ /(_)) )\ ((_)\
_(())\_)() ((_)(_))_ ((_) _((_)
\ \((_)/ // _ \ | \ | __|| \| |
\ \/\/ /| (_) || |) || _| | .` |
\_/\_/ \___/ |___/ |___||_|\_|
You are using wodenpy version/git hash: a43c383
2025-03-17 10:19:42 - INFO - Input arguments after parsing:
Phase centre: 0.00000, 90.00000 deg
Array central latitude: 52.905 deg
Array central longitude: 6.868 deg
Array central height: 0.000 m
Lowest channel frequency: 75.000 MHz
Channel frequency resolution: 1000.000 kHz
Coarse band bandwidth: 0.195 MHz
Num channels per coarse band: 20
Start date: 2014-10-16T17:01:36.000
Time resolution: 3600.0 (s)
Num time steps: 24
Have read 37 antenna positions from measurement set: ../../test_installation/everybeam/lba.MS
Will write outputs to: /home/jack-line/software/WODEN_dev/examples/LOFAR_LBA_NCP
2025-03-17 10:19:43 - INFO - Obs epoch initial LST was 294.9211333249 deg
2025-03-17 10:19:43 - INFO - Setting initial J2000 LST to 294.8284574349 deg
2025-03-17 10:19:43 - INFO - Setting initial mjd to 56946.7302777780
2025-03-17 10:19:43 - INFO - After precession initial latitude of the array is 52.8683555426 deg
2025-03-17 10:19:43 - INFO - Precessing array layout to J2000
2025-03-17 10:19:43 - INFO - Will run with EveryBeam LOFAR primary beam, based on this measurement set:
../../test_installation/everybeam/lba.MS
Created the following minimal MS to point the beam:
/home/jack-line/software/WODEN_dev/examples/LOFAR_LBA_NCP/pointed_woden_LBA-NCP_band01.ms
Primary beam is pointed at RA,Dec = 0.0, 90.0 deg
2025-03-17 10:19:43 - 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-17 10:19:43 - INFO - Sky model mapping took 0.0 mins
2025-03-17 10:19:43 - INFO - Have read in 602 components
2025-03-17 10:19:43 - INFO - After cropping there are 602 components
2025-03-17 10:19:43 - INFO - Will load libwoden from /home/jack-line/software/WODEN_dev/wodenpy/libwoden_double.so
2025-03-17 10:19:43 - INFO - Running in parallel on GPU.
Will read sky model using 8 threads
There are 1 sets of sky models to run
2025-03-17 10:19:43 - INFO - Reading set 0 sky models
2025-03-17 10:19:43 - INFO - From sky set 0 thread num 0 reading 0 points, 76 gauss, 0 shape, 0 shape coeffs
2025-03-17 10:19:43 - INFO - From sky set 0 thread num 1 reading 0 points, 76 gauss, 0 shape, 0 shape coeffs
2025-03-17 10:19:43 - INFO - Finshed sky set 0 reading thread num 0 in 0.0 seconds
2025-03-17 10:19:43 - INFO - Finshed sky set 0 reading thread num 1 in 0.0 seconds
2025-03-17 10:19:43 - INFO - From sky set 0 thread num 2 reading 0 points, 76 gauss, 0 shape, 0 shape coeffs
2025-03-17 10:19:43 - INFO - From sky set 0 thread num 3 reading 0 points, 76 gauss, 0 shape, 0 shape coeffs
2025-03-17 10:19:43 - INFO - Finshed sky set 0 reading thread num 2 in 0.0 seconds
2025-03-17 10:19:43 - INFO - From sky set 0 thread num 4 reading 0 points, 76 gauss, 0 shape, 0 shape coeffs
2025-03-17 10:19:43 - INFO - Finshed sky set 0 reading thread num 3 in 0.0 seconds
2025-03-17 10:19:43 - INFO - From sky set 0 thread num 5 reading 0 points, 76 gauss, 0 shape, 0 shape coeffs
2025-03-17 10:19:43 - INFO - Finshed sky set 0 reading thread num 4 in 0.0 seconds
2025-03-17 10:19:43 - INFO - From sky set 0 thread num 6 reading 0 points, 76 gauss, 0 shape, 0 shape coeffs
2025-03-17 10:19:43 - INFO - Finshed sky set 0 reading thread num 5 in 0.0 seconds
2025-03-17 10:19:43 - INFO - Finshed sky set 0 reading thread num 6 in 0.0 seconds
2025-03-17 10:19:43 - INFO - From sky set 0 thread num 7 reading 0 points, 70 gauss, 0 shape, 0 shape coeffs
2025-03-17 10:19:43 - INFO - Finshed sky set 0 reading thread num 7 in 0.0 seconds
2025-03-17 10:19:43 - INFO - Sending Sky set 0 to the GPU
2025-03-17 00:19:44 SEVERE MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292) Leap second table TAI_UTC seems out-of-date.
2025-03-17 00:19:44 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-17 00:19:44 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-17 10:21:19 - INFO - Sky set 0 has returned from the GPU after 95.4 seconds
2025-03-17 10:21:19 - INFO - Have completed 602 of 602 components calcs (100.0%)
2025-03-17 10:21:19 - INFO - Finished all rounds of processing
2025-03-17 10:21:19 - INFO - Deleting /home/jack-line/software/WODEN_dev/examples/LOFAR_LBA_NCP/pointed_woden_LBA-NCP_band01.ms...
2025-03-17 10:21:19 - INFO - Done
2025-03-17 10:21:19 - INFO - Full run took 0:01:36.999311
2025-03-17 10:21:19 - INFO - WODEN is done. Closing the log. S'later
Let’s see what uv coverage we got:
from astropy.constants import c
c = c.to('km/s').value
with fits.open('woden_LBA-NCP_band01.uvfits') as hdu:
uu = (hdu[0].data['UU']*c)
vv = (hdu[0].data['VV']*c)
fig, axs = plt.subplots(1, 1, figsize=(8, 8))
axs.plot(uu, vv, 'k.', ms=1)
axs.plot(-uu, -vv, 'k.', ms=1)
axs.set_xlabel('U (km)')
axs.set_ylabel('V (km)')
plt.show()
ALL HAIL THE COSMIC EYE. Convert the UVFITS to an MS. Below is the function that’s in woden_uv2ms.py. Just copied it here to show what’s going on.
from pyuvdata import UVData
from subprocess import call
import os
def make_ms(uvfits_file, no_delete=False):
"""Takes a path to a uvfits file, checks it exists, and converts it to
a measurement set of the same name at the same location"""
##If file exists, attempt to tranform it
if os.path.isfile(uvfits_file):
##Get rid of '.uvfits' off the end of the file
name = uvfits_file[:-7]
##Delete old measurement set if requested
if no_delete:
pass
else:
call("rm -r %s.ms" %name,shell=True)
UV = UVData()
UV.read(uvfits_file)
UV.write_ms("{:s}.ms".format(name))
else:
##print warning that uvfits file doesn't exist
print("Could not find the uvfits specified by the path: "
"\t{:s}. Skipping this tranformation".format(uvfits_file))
make_ms('woden_LBA-NCP_band01.uvfits')
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 191664.74533642438 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 191664.74533642438 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).
Righto, let’s image it up! We’ll use the pixel resolution of the WENSS image, and also use a restoring beam that matches the WNESS resolutions. This should give us a decent match image wise.
%%bash
##YOU PROBABLY DON'T NEED THIS ENV VARIABLE
# export OPENBLAS_NUM_THREADS=1
wsclean -name woden_LBA-NCP -niter 20000 \
-size 1024 1024 -scale 0.005859375 \
-weight briggs 0 \
-auto-threshold 0.5 -auto-mask 3 -pol I -multiscale \
-j 12 -mgain 0.85 -no-update-model-required \
-beam-size 56 \
woden_LBA-NCP_band01.ms
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 75-95 (20)
Reordering woden_LBA-NCP_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 woden_LBA-NCP_band01.ms
Detected 62.7 GB of system memory, usage not limited.
Opening reordered part 0 spw 0 for woden_LBA-NCP_band01.ms
Determining min and max w & theoretical beam size... DONE (w=[0.0061668:5135.65] lambdas, maxuvw=8519.12 lambda)
Theoretic beam = 24.21''
Small inversion enabled, but inversion resolution already smaller than beam size: not using optimization.
Loading data in memory...
Gridding 15984 rows...
Gridded visibility count: 221909
Writing psf image... DONE
== Constructing image ==
Opening reordered part 0 spw 0 for woden_LBA-NCP_band01.ms
Loading data in memory...
Gridding 15984 rows...
Gridded visibility count: 221909
Writing dirty image...
== Deconvolving (1) ==
Estimated standard deviation of background noise: 21.35 mJy
Scale info:
- Scale 0, bias factor=1, psfpeak=1, gain=0.1, kernel peak=0.19661
- Scale 5, bias factor=1.7, psfpeak=0.291008, gain=0.343633, kernel peak=0.0812188
- Scale 9, bias factor=2.8, psfpeak=0.180118, gain=0.555191, kernel peak=0.0330129
- Scale 18, bias factor=4.6, psfpeak=0.09095, gain=1.0995, kernel peak=0.00965284
- Scale 37, bias factor=7.7, psfpeak=0.0324371, gain=3.08289, kernel peak=0.00287457
- Scale 73, bias factor=12.9, psfpeak=0.00893642, gain=11.1902, kernel peak=0.000790563
- Scale 147, bias factor=21.4, psfpeak=0.00224912, gain=44.4619, kernel peak=0.000202596
- Scale 294, bias factor=35.7, psfpeak=0.000518892, gain=192.718, kernel peak=5.1938e-05
RMS per scale: {0: 30.23 mJy, 5: 18.86 mJy, 9: 15.8 mJy, 18: 11.21 mJy, 37: 7.21 mJy, 73: 4.13 mJy, 147: 2.21 mJy, 294: 1.05 mJy}
Starting multi-scale cleaning. Start peak=8.29 Jy, major iteration threshold=1.24 Jy
Iteration 3, scale 0 px : 6.37 Jy at 482,428
Iteration 7, scale 0 px : 4.84 Jy at 483,429
Iteration 12, scale 0 px : 3.8 Jy at 482,428
Iteration 17, scale 0 px : 2.87 Jy at 483,428
Iteration 24, scale 0 px : 2.24 Jy at 482,429
Iteration 31, scale 0 px : 1.79 Jy at 465,680
Iteration 40, scale 0 px : 1.32 Jy at 483,429
Subminor loop is near minor loop threshold. Initiating countdown.
(11) Iteration 45, scale 0 px : 1.21 Jy at 483,428
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 woden_LBA-NCP_band01.ms
Predicting 15984 rows...
Writing...
== Constructing image ==
Opening reordered part 0 spw 0 for woden_LBA-NCP_band01.ms
Loading data in memory...
Gridding 15984 rows...
Gridded visibility count: 221909
== Deconvolving (2) ==
Estimated standard deviation of background noise: 9.58 mJy
Scale info:
- Scale 0, bias factor=1, psfpeak=1, gain=0.1, kernel peak=0.19661
- Scale 5, bias factor=1.7, psfpeak=0.291008, gain=0.343633, kernel peak=0.0812188
- Scale 9, bias factor=2.8, psfpeak=0.180118, gain=0.555191, kernel peak=0.0330129
- Scale 18, bias factor=4.6, psfpeak=0.09095, gain=1.0995, kernel peak=0.00965284
- Scale 37, bias factor=7.7, psfpeak=0.0324371, gain=3.08289, kernel peak=0.00287457
- Scale 73, bias factor=12.9, psfpeak=0.00893642, gain=11.1902, kernel peak=0.000790563
- Scale 147, bias factor=21.4, psfpeak=0.00224912, gain=44.4619, kernel peak=0.000202596
- Scale 294, bias factor=35.7, psfpeak=0.000518892, gain=192.718, kernel peak=5.1938e-05
RMS per scale: {0: 12.33 mJy, 5: 6.91 mJy, 9: 5.79 mJy, 18: 4.1 mJy, 37: 2.64 mJy, 73: 1.53 mJy, 147: 850.18 µJy, 294: 429.23 µJy}
Starting multi-scale cleaning. Start peak=1.21 Jy, major iteration threshold=181.46 mJy
Iteration 54, scale 0 px : 967.78 mJy at 465,680
Iteration 71, scale 0 px : 772.56 mJy at 488,439
Iteration 89, scale 9 px : 651.87 mJy at 572,379
Iteration 96, scale 0 px : 608.09 mJy at 483,428
Iteration 121, scale 0 px : -534.96 mJy at 481,428
Iteration 151, scale 9 px : 449.15 mJy at 572,379
Iteration 161, scale 0 px : -427.13 mJy at 481,428
Iteration 221, scale 0 px : 341.61 mJy at 495,426
Iteration 311, scale 9 px : 325.35 mJy at 453,500
Iteration 319, scale 9 px : 259.38 mJy at 453,499
Iteration 345, scale 0 px : 284.7 mJy at 483,426
Iteration 442, scale 0 px : 247.79 mJy at 485,428
Iteration 557, scale 0 px : -200.74 mJy at 479,428
Subminor loop is near minor loop threshold. Initiating countdown.
(11) Iteration 640, scale 0 px : 184.71 mJy at 483,433
(10) Iteration 641, scale 0 px : 181.21 mJy at 435,247
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 woden_LBA-NCP_band01.ms
Predicting 15984 rows...
Writing...
== Constructing image ==
Opening reordered part 0 spw 0 for woden_LBA-NCP_band01.ms
Loading data in memory...
Gridding 15984 rows...
Gridded visibility count: 221909
== Deconvolving (3) ==
Estimated standard deviation of background noise: 5.81 mJy
Scale info:
- Scale 0, bias factor=1, psfpeak=1, gain=0.1, kernel peak=0.19661
- Scale 5, bias factor=1.7, psfpeak=0.291008, gain=0.343633, kernel peak=0.0812188
- Scale 9, bias factor=2.8, psfpeak=0.180118, gain=0.555191, kernel peak=0.0330129
- Scale 18, bias factor=4.6, psfpeak=0.09095, gain=1.0995, kernel peak=0.00965284
- Scale 37, bias factor=7.7, psfpeak=0.0324371, gain=3.08289, kernel peak=0.00287457
- Scale 73, bias factor=12.9, psfpeak=0.00893642, gain=11.1902, kernel peak=0.000790563
- Scale 147, bias factor=21.4, psfpeak=0.00224912, gain=44.4619, kernel peak=0.000202596
- Scale 294, bias factor=35.7, psfpeak=0.000518892, gain=192.718, kernel peak=5.1938e-05
RMS per scale: {0: 6.88 mJy, 5: 3.41 mJy, 9: 2.81 mJy, 18: 1.92 mJy, 37: 1.16 mJy, 73: 627.79 µJy, 147: 319.49 µJy, 294: 158.98 µJy}
Starting multi-scale cleaning. Start peak=187.9 mJy, major iteration threshold=28.19 mJy
Iteration 861, scale 9 px : 157.11 mJy at 453,499
Iteration 877, scale 0 px : 150.2 mJy at 495,428
Iteration 1181, scale 0 px : 124.6 mJy at 482,422
Iteration 1516, scale 9 px : 109.54 mJy at 663,327
Iteration 1532, scale 0 px : 99.76 mJy at 693,715
Iteration 1992, scale 0 px : 79.76 mJy at 661,327
Iteration 2551, scale 9 px : 81.41 mJy at 475,480
Iteration 2557, scale 9 px : 64.95 mJy at 466,491
Iteration 2614, scale 0 px : 68.58 mJy at 465,677
Iteration 3035, scale 0 px : 54.77 mJy at 623,685
Iteration 3823, scale 0 px : -48.72 mJy at 617,233
Iteration 4375, scale 18 px : -48.41 mJy at 569,378
Iteration 4389, scale 9 px : 42.59 mJy at 475,479
Iteration 4426, scale 0 px : 39.78 mJy at 572,378
Iteration 5489, scale 18 px : -39.51 mJy at 570,378
Iteration 5494, scale 18 px : 31.1 mJy at 820,870
Subminor loop is near minor loop threshold. Initiating countdown.
(11) Iteration 5506, scale 0 px : -34.84 mJy at 806,390
(10) Iteration 6393, scale 0 px : 31.81 mJy at 428,606
(9) Iteration 6415, scale 9 px : 30.99 mJy at 258,495
(8) Iteration 6429, scale 18 px : -30.74 mJy at 571,378
(7) Iteration 6431, scale 0 px : 30.67 mJy at 572,379
(6) Iteration 6442, scale 0 px : -28.46 mJy at 483,411
(5) Iteration 6444, scale 18 px : -30.21 mJy at 571,378
(4) Iteration 6446, scale 0 px : 29.3 mJy at 571,377
(3) Iteration 6450, scale 18 px : -28.87 mJy at 571,378
(2) Iteration 6451, scale 0 px : 28.35 mJy at 563,406
(1) Iteration 6453, scale 0 px : -28.82 mJy at 475,428
(0) Iteration 6454, scale 0 px : 28.18 mJy at 483,412
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 woden_LBA-NCP_band01.ms
Predicting 15984 rows...
Writing...
== Constructing image ==
Opening reordered part 0 spw 0 for woden_LBA-NCP_band01.ms
Loading data in memory...
Gridding 15984 rows...
Gridded visibility count: 221909
== Deconvolving (4) ==
Estimated standard deviation of background noise: 3.34 mJy
Scale info:
- Scale 0, bias factor=1, psfpeak=1, gain=0.1, kernel peak=0.19661
- Scale 5, bias factor=1.7, psfpeak=0.291008, gain=0.343633, kernel peak=0.0812188
- Scale 9, bias factor=2.8, psfpeak=0.180118, gain=0.555191, kernel peak=0.0330129
- Scale 18, bias factor=4.6, psfpeak=0.09095, gain=1.0995, kernel peak=0.00965284
- Scale 37, bias factor=7.7, psfpeak=0.0324371, gain=3.08289, kernel peak=0.00287457
- Scale 73, bias factor=12.9, psfpeak=0.00893642, gain=11.1902, kernel peak=0.000790563
- Scale 147, bias factor=21.4, psfpeak=0.00224912, gain=44.4619, kernel peak=0.000202596
- Scale 294, bias factor=35.7, psfpeak=0.000518892, gain=192.718, kernel peak=5.1938e-05
RMS per scale: {0: 3.56 mJy, 5: 1.65 mJy, 9: 1.29 mJy, 18: 808.86 µJy, 37: 471.21 µJy, 73: 267.46 µJy, 147: 138.49 µJy, 294: 70.04 µJy}
Starting multi-scale cleaning. Start peak=38.34 mJy, major iteration threshold=10.03 mJy (final)
Iteration 6534, scale 0 px : -31.49 mJy at 482,411
Iteration 7220, scale 9 px : 29.75 mJy at 94,511
Iteration 7243, scale 0 px : -27.79 mJy at 486,428
Iteration 8198, scale 18 px : -25.03 mJy at 571,379
Iteration 8235, scale 0 px : 24.63 mJy at 474,428
Iteration 9371, scale 9 px : 24.53 mJy at 93,512
Iteration 9399, scale 5 px : 20.18 mJy at 92,512
Iteration 9473, scale 0 px : 21.65 mJy at 573,379
Iteration 10895, scale 9 px : -20.2 mJy at 81,505
Iteration 10964, scale 0 px : 18.68 mJy at 458,428
Iteration 13075, scale 18 px : -17.42 mJy at 452,497
Iteration 13116, scale 9 px : -17.46 mJy at 82,505
Iteration 13200, scale 0 px : 16.37 mJy at 481,425
Iteration 15523, scale 5 px : 14.83 mJy at 92,512
Iteration 15552, scale 18 px : -15.19 mJy at 452,498
Iteration 15611, scale 9 px : -14.54 mJy at 87,524
Iteration 15855, scale 0 px : 14.32 mJy at 84,509
Iteration 19288, scale 73 px : 12.74 mJy at 558,813
Iteration 19298, scale 18 px : -13.46 mJy at 454,498
Iteration 19361, scale 0 px : 12.46 mJy at 481,423
Subminor loop is near minor loop threshold. Initiating countdown.
(11) Iteration 20000, scale 9 px : -11.94 mJy at 78,506
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 woden_LBA-NCP_band01.ms
Predicting 15984 rows...
Writing...
== Constructing image ==
Opening reordered part 0 spw 0 for woden_LBA-NCP_band01.ms
Loading data in memory...
Gridding 15984 rows...
Gridded visibility count: 221909
4 major iterations were performed.
Rendering sources to restored image (beam=56''-56'', PA=0 masec)... DONE
Writing restored image... DONE
Multi-scale cleaning summary:
- Scale 0 px, nr of components cleaned: 19006 (48.59 Jy)
- Scale 5 px, nr of components cleaned: 103 (141.07 mJy)
- Scale 9 px, nr of components cleaned: 645 (5.43 Jy)
- Scale 18 px, nr of components cleaned: 236 (169.18 mJy)
- Scale 37 px, nr of components cleaned: 0 (0 Jy)
- Scale 73 px, nr of components cleaned: 10 (61.7 mJy)
- Scale 147 px, nr of components cleaned: 0 (0 Jy)
- Scale 294 px, nr of components cleaned: 0 (0 Jy)
Total: 20000 components (54.39 Jy)
Inversion: 00:00:02.721378, prediction: 00:00:05.854667, deconvolution: 00:00:07.043551
Cleaning up temporary files...
I’m sure this isn’t an optimal way to do the CLEANing, and I might be able to go deeper, but this is good enough for a comparison. Let’s see what we get!
fig = plt.figure(figsize=(9, 11))
with fits.open('wenss_sky_edit.fits') as hdu:
data_wenss = hdu[0].data
header = hdu[0].header
wcs_wenss = WCS(header)
vmax = 5e-2
ax1 = fig.add_subplot(211, projection=wcs_wenss)
im = ax1.imshow(data_wenss, origin='lower', cmap='Blues', vmin=0, vmax=vmax)
ax1.grid()
plt.colorbar(im, ax=ax1, label='Flux (Jy/beam)')
ax1.set_title('WENSS image')
with fits.open('woden_LBA-NCP-image.fits') as hdu:
wcs_woden = WCS(hdu[0].header).celestial
data_woden = np.squeeze(hdu[0].data)
ax2 = fig.add_subplot(212, projection=wcs_woden)
vmax=1e-1
im = ax2.imshow(data_woden, origin='lower', cmap='Blues', vmin=0, vmax=vmax)
plt.colorbar(im, ax=ax2, label='Flux (Jy/beam)')
ax2.grid()
ax2.set_title('Imaged WENSSxNVSS catalogue sim')
ax1.set_xlabel(' ')
ax1.set_ylabel('Dec')
ax1.set_xticklabels([])
ax2.set_xlabel('RA')
ax2.set_ylabel('Dec')
half_width=300
cent_x = 512
cent_y = 512
for ax in [ax1, ax2]:
ax.set_xlim(cent_x - half_width, cent_x + half_width)
ax.set_ylim(cent_y - half_width, cent_y + half_width)
# plt.tight_layout()
plt.show()
WARNING: FITSFixedWarning: 'datfix' made the change 'Set MJD-OBS to 56946.730277 from DATE-OBS'. [astropy.wcs.wcs]
Righto, so the LBA primary beam is still in the image below, so the sources are dropping off at the edges. The units are different, because we have a different beam for LBA in the Jy/Beam units. Finally, the WENSS catalogue isn’t perfect, and we did a lazy cross-match, so the sky model isn’t perfect. But overall, this looks pretty good!
Do it all again, but with just NVSS¶
Out of interest, do we get better morphology if we just use the NVSS catalogue, and assume a spectral index of -0.8?
unq_source_ids = np.array(crop_nvss['NVSS'], dtype='U32')
##order everything by source name to make adding component names easier
ras_nvss_full = crop_nvss['RAJ2000'].data
decs_nvss_full = crop_nvss['DEJ2000'].data
##use the gaussians from NVSS as it has better resolution
fluxes_1400 = crop_nvss['S1_4'].data/1000.0
majors = crop_nvss['MajAxis'].data/3600.0
minors = crop_nvss['MinAxis'].data/3600.0
pas = np.array(crop_nvss['PA'], dtype=np.float64)
pas[np.isnan(pas)] = 0.0
fluxes_extrap = fluxes_1400*(200.0/1400.0)**-0.8
num_comps = len(unq_source_ids)
unq_source_id_col = Column(unq_source_ids, name='UNQ_SOURCE_ID')
names = np.array([name+'_C000' for name in unq_source_ids], dtype='U32')
names_col = Column(names, name='NAME')
ras_col = Column(ras_nvss_full, name='RA', unit='deg')
decs_col = Column(decs_nvss_full, name='DEC', unit='deg')
fluxes_col = Column(fluxes_extrap, name='NORM_COMP_PL', unit='Jy')
si_col = Column(np.full_like(ras_nvss_full, -0.8), name='ALPHA_PL', unit='Jy')
majors_col = Column(majors, name='MAJOR_DC', unit='deg')
minors_col = Column(minors, name='MINOR_DC', unit='deg')
pas_col = Column(pas, name='PA_DC', unit='deg')
##Anything that has a major axis is a Gaussian, other make it a point source
comp_types = np.empty(unq_source_ids.shape, dtype='U3')
comp_types[majors == 0.0] = 'P'
comp_types[majors != 0.0] = 'G'
comp_types_col = Column(comp_types, name='COMP_TYPE')
##We're inputting everything as a power law
mod_types_col = Column(np.full(num_comps, 'pl'), name='MOD_TYPE')
cols = [unq_source_id_col, names_col,
ras_col, decs_col, comp_types_col, mod_types_col,
fluxes_col, si_col,
majors_col, minors_col, pas_col]
new_table = Table(cols)
new_table.write('srclist_nvss_NCP.fits', overwrite=True)
%%bash
primary_beam="everybeam_LOFAR"
ra0=0.0
dec0=90.0
lofar_lat=52.905329712
lofar_long=6.867996528
date="2014-10-16T17:01:36.000"
ms_path=../../test_installation/everybeam/lba.MS
uvfits_name=woden_LBA-NCP_nvss
cat_name=srclist_nvss_NCP.fits
freq_reso=1e+6
low_freq=75e+6
num_freq_chans=20
time_res=3600.0
num_time_steps=24
run_woden.py \
--ra0=${ra0} \
--dec0=${dec0} \
--eb_point_to_phase \
--latitude=${lofar_lat} \
--longitude=${lofar_long} \
--date=${date} \
--output_uvfits_prepend=${uvfits_name} \
--cat_filename=${cat_name} \
--primary_beam=${primary_beam} \
--lowest_channel_freq=${low_freq} \
--freq_res=${freq_reso} \
--num_freq_channels=${num_freq_chans} \
--band_nums=1 \
--time_res=${time_res} \
--num_time_steps=${num_time_steps} \
--beam_ms_path=${ms_path}
/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 ../../test_installation/everybeam/lba.MS: 24 columns, 4218 rows
Successful readonly open of default-locked table ../../test_installation/everybeam/lba.MS/SPECTRAL_WINDOW: 14 columns, 1 rows
Successful read/write open of default-locked table ../../test_installation/everybeam/lba.MS::FIELD: 10 columns, 1 rows
Successful readonly open of default-locked table ../../test_installation/everybeam/lba.MS: 24 columns, 4218 rows
Successful read/write open of default-locked table /home/jack-line/software/WODEN_dev/examples/LOFAR_LBA_NCP/pointed_woden_LBA-NCP_nvss_band01.ms::FIELD: 10 columns, 1 rows
Successful readonly open of default-locked table ../../test_installation/everybeam/lba.MS/ANTENNA: 10 columns, 37 rows
2025-03-17 10:29:57 - INFO -
) ( )
( ( ( /( )\ ) ( /(
)\))( ')\())(()/( ( )\())
((_)()\ )((_)\ /(_)) )\ ((_)\
_(())\_)() ((_)(_))_ ((_) _((_)
\ \((_)/ // _ \ | \ | __|| \| |
\ \/\/ /| (_) || |) || _| | .` |
\_/\_/ \___/ |___/ |___||_|\_|
You are using wodenpy version/git hash: a43c383
2025-03-17 10:29:57 - INFO - Input arguments after parsing:
Phase centre: 0.00000, 90.00000 deg
Array central latitude: 52.905 deg
Array central longitude: 6.868 deg
Array central height: 0.000 m
Lowest channel frequency: 75.000 MHz
Channel frequency resolution: 1000.000 kHz
Coarse band bandwidth: 0.195 MHz
Num channels per coarse band: 20
Start date: 2014-10-16T17:01:36.000
Time resolution: 3600.0 (s)
Num time steps: 24
Have read 37 antenna positions from measurement set: ../../test_installation/everybeam/lba.MS
Will write outputs to: /home/jack-line/software/WODEN_dev/examples/LOFAR_LBA_NCP
2025-03-17 10:29:57 - INFO - Obs epoch initial LST was 294.9211333249 deg
2025-03-17 10:29:57 - INFO - Setting initial J2000 LST to 294.8284574349 deg
2025-03-17 10:29:57 - INFO - Setting initial mjd to 56946.7302777780
2025-03-17 10:29:57 - INFO - After precession initial latitude of the array is 52.8683555426 deg
2025-03-17 10:29:57 - INFO - Precessing array layout to J2000
2025-03-17 10:29:57 - INFO - Will run with EveryBeam LOFAR primary beam, based on this measurement set:
../../test_installation/everybeam/lba.MS
Created the following minimal MS to point the beam:
/home/jack-line/software/WODEN_dev/examples/LOFAR_LBA_NCP/pointed_woden_LBA-NCP_nvss_band01.ms
Primary beam is pointed at RA,Dec = 0.0, 90.0 deg
2025-03-17 10:29:57 - 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-17 10:29:57 - INFO - Sky model mapping took 0.0 mins
2025-03-17 10:29:57 - INFO - Have read in 1272 components
2025-03-17 10:29:57 - INFO - After cropping there are 1272 components
2025-03-17 10:29:57 - INFO - Will load libwoden from /home/jack-line/software/WODEN_dev/wodenpy/libwoden_double.so
2025-03-17 10:29:58 - INFO - Running in parallel on GPU.
Will read sky model using 8 threads
There are 1 sets of sky models to run
2025-03-17 10:29:58 - INFO - Reading set 0 sky models
2025-03-17 10:29:58 - INFO - From sky set 0 thread num 0 reading 0 points, 159 gauss, 0 shape, 0 shape coeffs
2025-03-17 10:29:58 - INFO - Finshed sky set 0 reading thread num 0 in 0.0 seconds
2025-03-17 10:29:58 - INFO - From sky set 0 thread num 1 reading 0 points, 159 gauss, 0 shape, 0 shape coeffs
2025-03-17 10:29:58 - INFO - From sky set 0 thread num 3 reading 0 points, 159 gauss, 0 shape, 0 shape coeffs
2025-03-17 10:29:58 - INFO - From sky set 0 thread num 2 reading 0 points, 159 gauss, 0 shape, 0 shape coeffs
2025-03-17 10:29:58 - INFO - Finshed sky set 0 reading thread num 1 in 0.0 seconds
2025-03-17 10:29:58 - INFO - Finshed sky set 0 reading thread num 3 in 0.0 seconds
2025-03-17 10:29:58 - INFO - From sky set 0 thread num 4 reading 0 points, 159 gauss, 0 shape, 0 shape coeffs
2025-03-17 10:29:58 - INFO - From sky set 0 thread num 5 reading 0 points, 159 gauss, 0 shape, 0 shape coeffs
2025-03-17 10:29:58 - INFO - Finshed sky set 0 reading thread num 2 in 0.0 seconds
2025-03-17 10:29:58 - INFO - Finshed sky set 0 reading thread num 4 in 0.0 seconds
2025-03-17 10:29:58 - INFO - Finshed sky set 0 reading thread num 5 in 0.0 seconds
2025-03-17 10:29:58 - INFO - From sky set 0 thread num 6 reading 0 points, 159 gauss, 0 shape, 0 shape coeffs
2025-03-17 10:29:58 - INFO - Finshed sky set 0 reading thread num 6 in 0.0 seconds
2025-03-17 10:29:58 - INFO - From sky set 0 thread num 7 reading 0 points, 159 gauss, 0 shape, 0 shape coeffs
2025-03-17 10:29:58 - INFO - Finshed sky set 0 reading thread num 7 in 0.0 seconds
2025-03-17 10:29:58 - INFO - Sending Sky set 0 to the GPU
2025-03-17 00:29:59 SEVERE MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292) Leap second table TAI_UTC seems out-of-date.
2025-03-17 00:29:59 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-17 00:29:59 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-17 10:33:04 - INFO - Sky set 0 has returned from the GPU after 185.7 seconds
2025-03-17 10:33:04 - INFO - Have completed 1272 of 1272 components calcs (100.0%)
2025-03-17 10:33:04 - INFO - Finished all rounds of processing
2025-03-17 10:33:04 - INFO - Deleting /home/jack-line/software/WODEN_dev/examples/LOFAR_LBA_NCP/pointed_woden_LBA-NCP_nvss_band01.ms...
2025-03-17 10:33:04 - INFO - Done
2025-03-17 10:33:04 - INFO - Full run took 0:03:07.319572
2025-03-17 10:33:04 - INFO - WODEN is done. Closing the log. S'later
make_ms('woden_LBA-NCP_nvss_band01.uvfits')
rm: cannot remove 'woden_LBA-NCP_nvss_band01.ms': No such file or directory
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 191664.74533642438 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 191664.74533642438 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).
%%bash
wsclean -name woden_LBA-NCP_nvss -niter 20000 \
-size 1024 1024 -scale 0.005859375 \
-weight briggs 0 \
-auto-threshold 0.5 -auto-mask 3 -pol I -multiscale \
-j 12 -mgain 0.85 -no-update-model-required \
-beam-size 56 \
woden_LBA-NCP_nvss_band01.ms
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 75-95 (20)
Reordering woden_LBA-NCP_nvss_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 woden_LBA-NCP_nvss_band01.ms
Detected 62.7 GB of system memory, usage not limited.
Opening reordered part 0 spw 0 for woden_LBA-NCP_nvss_band01.ms
Determining min and max w & theoretical beam size... DONE (w=[0.0061668:5135.65] lambdas, maxuvw=8519.12 lambda)
Theoretic beam = 24.21''
Small inversion enabled, but inversion resolution already smaller than beam size: not using optimization.
Loading data in memory...
Gridding 15984 rows...
Gridded visibility count: 221909
Writing psf image... DONE
== Constructing image ==
Opening reordered part 0 spw 0 for woden_LBA-NCP_nvss_band01.ms
Loading data in memory...
Gridding 15984 rows...
Gridded visibility count: 221909
Writing dirty image...
== Deconvolving (1) ==
Estimated standard deviation of background noise: 28.65 mJy
Scale info:
- Scale 0, bias factor=1, psfpeak=1, gain=0.1, kernel peak=0.19661
- Scale 5, bias factor=1.7, psfpeak=0.291008, gain=0.343633, kernel peak=0.0812188
- Scale 9, bias factor=2.8, psfpeak=0.180118, gain=0.555191, kernel peak=0.0330129
- Scale 18, bias factor=4.6, psfpeak=0.09095, gain=1.0995, kernel peak=0.00965284
- Scale 37, bias factor=7.7, psfpeak=0.0324371, gain=3.08289, kernel peak=0.00287457
- Scale 73, bias factor=12.9, psfpeak=0.00893642, gain=11.1902, kernel peak=0.000790563
- Scale 147, bias factor=21.4, psfpeak=0.00224912, gain=44.4619, kernel peak=0.000202596
- Scale 294, bias factor=35.7, psfpeak=0.000518892, gain=192.718, kernel peak=5.1938e-05
RMS per scale: {0: 40.63 mJy, 5: 25.23 mJy, 9: 21.13 mJy, 18: 14.78 mJy, 37: 9.08 mJy, 73: 5 mJy, 147: 2.51 mJy, 294: 1.13 mJy}
Starting multi-scale cleaning. Start peak=11.72 Jy, major iteration threshold=1.76 Jy
Iteration 3, scale 0 px : 9.1 Jy at 483,428
Iteration 7, scale 0 px : 6.92 Jy at 483,429
Iteration 12, scale 0 px : 5.28 Jy at 482,428
Iteration 17, scale 0 px : 4.07 Jy at 483,428
Iteration 22, scale 0 px : 3.11 Jy at 482,429
Iteration 27, scale 0 px : 2.37 Jy at 483,429
Iteration 32, scale 0 px : 1.8 Jy at 482,428
Subminor loop is near minor loop threshold. Initiating countdown.
(11) Iteration 33, scale 0 px : 1.73 Jy at 483,428
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 woden_LBA-NCP_nvss_band01.ms
Predicting 15984 rows...
Writing...
== Constructing image ==
Opening reordered part 0 spw 0 for woden_LBA-NCP_nvss_band01.ms
Loading data in memory...
Gridding 15984 rows...
Gridded visibility count: 221909
== Deconvolving (2) ==
Estimated standard deviation of background noise: 11.76 mJy
Scale info:
- Scale 0, bias factor=1, psfpeak=1, gain=0.1, kernel peak=0.19661
- Scale 5, bias factor=1.7, psfpeak=0.291008, gain=0.343633, kernel peak=0.0812188
- Scale 9, bias factor=2.8, psfpeak=0.180118, gain=0.555191, kernel peak=0.0330129
- Scale 18, bias factor=4.6, psfpeak=0.09095, gain=1.0995, kernel peak=0.00965284
- Scale 37, bias factor=7.7, psfpeak=0.0324371, gain=3.08289, kernel peak=0.00287457
- Scale 73, bias factor=12.9, psfpeak=0.00893642, gain=11.1902, kernel peak=0.000790563
- Scale 147, bias factor=21.4, psfpeak=0.00224912, gain=44.4619, kernel peak=0.000202596
- Scale 294, bias factor=35.7, psfpeak=0.000518892, gain=192.718, kernel peak=5.1938e-05
RMS per scale: {0: 14.6 mJy, 5: 7.34 mJy, 9: 6.12 mJy, 18: 4.29 mJy, 37: 2.64 mJy, 73: 1.44 mJy, 147: 697.94 µJy, 294: 299.09 µJy}
Starting multi-scale cleaning. Start peak=1.73 Jy, major iteration threshold=258.77 mJy
Iteration 42, scale 0 px : 1.32 Jy at 495,427
Iteration 54, scale 0 px : 1.04 Jy at 465,680
Iteration 66, scale 0 px : 821.5 mJy at 495,427
Iteration 89, scale 0 px : -658.31 mJy at 483,430
Iteration 129, scale 0 px : 523.58 mJy at 619,233
Iteration 183, scale 9 px : 424.5 mJy at 495,427
Iteration 197, scale 0 px : 418.72 mJy at 483,429
Iteration 277, scale 0 px : 360.09 mJy at 483,426
Iteration 364, scale 9 px : 309.63 mJy at 360,830
Subminor loop is near minor loop threshold. Initiating countdown.
(11) Iteration 372, scale 0 px : -306.86 mJy at 483,432
(10) Iteration 446, scale 0 px : 262.34 mJy at 487,428
(9) Iteration 448, scale 0 px : 258.3 mJy at 316,460
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 woden_LBA-NCP_nvss_band01.ms
Predicting 15984 rows...
Writing...
== Constructing image ==
Opening reordered part 0 spw 0 for woden_LBA-NCP_nvss_band01.ms
Loading data in memory...
Gridding 15984 rows...
Gridded visibility count: 221909
== Deconvolving (3) ==
Estimated standard deviation of background noise: 7.29 mJy
Scale info:
- Scale 0, bias factor=1, psfpeak=1, gain=0.1, kernel peak=0.19661
- Scale 5, bias factor=1.7, psfpeak=0.291008, gain=0.343633, kernel peak=0.0812188
- Scale 9, bias factor=2.8, psfpeak=0.180118, gain=0.555191, kernel peak=0.0330129
- Scale 18, bias factor=4.6, psfpeak=0.09095, gain=1.0995, kernel peak=0.00965284
- Scale 37, bias factor=7.7, psfpeak=0.0324371, gain=3.08289, kernel peak=0.00287457
- Scale 73, bias factor=12.9, psfpeak=0.00893642, gain=11.1902, kernel peak=0.000790563
- Scale 147, bias factor=21.4, psfpeak=0.00224912, gain=44.4619, kernel peak=0.000202596
- Scale 294, bias factor=35.7, psfpeak=0.000518892, gain=192.718, kernel peak=5.1938e-05
RMS per scale: {0: 8.69 mJy, 5: 4.29 mJy, 9: 3.58 mJy, 18: 2.47 mJy, 37: 1.48 mJy, 73: 786.39 µJy, 147: 370.54 µJy, 294: 162.68 µJy}
Starting multi-scale cleaning. Start peak=270 mJy, major iteration threshold=40.5 mJy
Iteration 594, scale 9 px : 245.59 mJy at 453,500
Iteration 608, scale 0 px : 226.18 mJy at 482,424
Iteration 809, scale 0 px : -185.27 mJy at 482,423
Iteration 1081, scale 9 px : 197.35 mJy at 453,501
Iteration 1085, scale 9 px : 157.79 mJy at 453,500
Iteration 1109, scale 0 px : 147.81 mJy at 495,428
Iteration 1438, scale 0 px : -119.38 mJy at 482,434
Iteration 1829, scale 9 px : 124.11 mJy at 453,501
Iteration 1834, scale 9 px : 98.78 mJy at 453,501
Iteration 1875, scale 0 px : -105.42 mJy at 490,428
Iteration 2146, scale 0 px : -88.61 mJy at 483,436
Iteration 2601, scale 18 px : 81.78 mJy at 336,248
Iteration 2608, scale 0 px : 77.98 mJy at 482,437
Iteration 3027, scale 9 px : 68.39 mJy at 452,501
Iteration 3054, scale 0 px : 62.49 mJy at 843,786
Iteration 3916, scale 0 px : -53.35 mJy at 391,128
Iteration 4763, scale 9 px : 49.54 mJy at 437,408
Subminor loop is near minor loop threshold. Initiating countdown.
(11) Iteration 4787, scale 18 px : 42.99 mJy at 518,479
(10) Iteration 4790, scale 0 px : -44.17 mJy at 836,512
(9) Iteration 5078, scale 0 px : 42.8 mJy at 481,430
(8) Iteration 5108, scale 0 px : -42.15 mJy at 482,434
(7) Iteration 5127, scale 0 px : 42.41 mJy at 482,435
(6) Iteration 5135, scale 0 px : 42.06 mJy at 495,428
(5) Iteration 5142, scale 0 px : -41.66 mJy at 467,428
(4) Iteration 5143, scale 0 px : 40.49 mJy at 234,538
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 woden_LBA-NCP_nvss_band01.ms
Predicting 15984 rows...
Writing...
== Constructing image ==
Opening reordered part 0 spw 0 for woden_LBA-NCP_nvss_band01.ms
Loading data in memory...
Gridding 15984 rows...
Gridded visibility count: 221909
== Deconvolving (4) ==
Estimated standard deviation of background noise: 4.12 mJy
Scale info:
- Scale 0, bias factor=1, psfpeak=1, gain=0.1, kernel peak=0.19661
- Scale 5, bias factor=1.7, psfpeak=0.291008, gain=0.343633, kernel peak=0.0812188
- Scale 9, bias factor=2.8, psfpeak=0.180118, gain=0.555191, kernel peak=0.0330129
- Scale 18, bias factor=4.6, psfpeak=0.09095, gain=1.0995, kernel peak=0.00965284
- Scale 37, bias factor=7.7, psfpeak=0.0324371, gain=3.08289, kernel peak=0.00287457
- Scale 73, bias factor=12.9, psfpeak=0.00893642, gain=11.1902, kernel peak=0.000790563
- Scale 147, bias factor=21.4, psfpeak=0.00224912, gain=44.4619, kernel peak=0.000202596
- Scale 294, bias factor=35.7, psfpeak=0.000518892, gain=192.718, kernel peak=5.1938e-05
RMS per scale: {0: 4.47 mJy, 5: 2.08 mJy, 9: 1.65 mJy, 18: 1.07 mJy, 37: 626.78 µJy, 73: 343.89 µJy, 147: 169.74 µJy, 294: 87.96 µJy}
Starting multi-scale cleaning. Start peak=56.67 mJy, major iteration threshold=12.36 mJy (final)
Iteration 5184, scale 0 px : -46.37 mJy at 483,446
Iteration 5708, scale 0 px : 40.23 mJy at 466,428
Iteration 6657, scale 18 px : 38.75 mJy at 425,528
Iteration 6670, scale 0 px : 36.51 mJy at 507,428
Iteration 7421, scale 9 px : 34.45 mJy at 682,370
Iteration 7466, scale 0 px : -32.88 mJy at 482,419
Iteration 8412, scale 18 px : -29.56 mJy at 568,378
Iteration 8452, scale 0 px : 29.28 mJy at 481,430
Iteration 9676, scale 9 px : 26.68 mJy at 447,513
Iteration 9778, scale 0 px : 26.43 mJy at 493,428
Iteration 11001, scale 5 px : -22.95 mJy at 360,833
Iteration 11030, scale 18 px : -24.15 mJy at 569,378
Iteration 11061, scale 0 px : 22.87 mJy at 360,830
Iteration 13218, scale 9 px : 21.7 mJy at 447,513
Iteration 13324, scale 0 px : -20.17 mJy at 571,380
Iteration 15704, scale 18 px : -19.88 mJy at 570,379
Iteration 15740, scale 18 px : 15.91 mJy at 441,861
Iteration 16094, scale 0 px : 18.66 mJy at 572,380
Iteration 17920, scale 9 px : -17.34 mJy at 511,380
Iteration 18102, scale 5 px : -16.6 mJy at 360,833
Iteration 18122, scale 0 px : 17 mJy at 361,831
Iteration 20000, scale 0 px : 15.25 mJy at 529,428
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 woden_LBA-NCP_nvss_band01.ms
Predicting 15984 rows...
Writing...
== Constructing image ==
Opening reordered part 0 spw 0 for woden_LBA-NCP_nvss_band01.ms
Loading data in memory...
Gridding 15984 rows...
Gridded visibility count: 221909
4 major iterations were performed.
Rendering sources to restored image (beam=56''-56'', PA=0 masec)... DONE
Writing restored image... DONE
Multi-scale cleaning summary:
- Scale 0 px, nr of components cleaned: 18871 (56.26 Jy)
- Scale 5 px, nr of components cleaned: 49 (95.06 mJy)
- Scale 9 px, nr of components cleaned: 596 (5.32 Jy)
- Scale 18 px, nr of components cleaned: 484 (1.21 Jy)
- Scale 37 px, nr of components cleaned: 0 (0 Jy)
- Scale 73 px, nr of components cleaned: 0 (0 Jy)
- Scale 147 px, nr of components cleaned: 0 (0 Jy)
- Scale 294 px, nr of components cleaned: 0 (0 Jy)
Total: 20000 components (62.89 Jy)
Inversion: 00:00:02.776069, prediction: 00:00:06.041787, deconvolution: 00:00:06.225899
Cleaning up temporary files...
fig = plt.figure(figsize=(9, 11))
vmax = 5e-2
ax1 = fig.add_subplot(211, projection=wcs_wenss)
im = ax1.imshow(data_wenss, origin='lower', cmap='Blues', vmin=0, vmax=vmax)
ax1.grid()
plt.colorbar(im, ax=ax1, label='Flux (Jy/beam)')
ax1.set_title('WENSS image')
with fits.open('woden_LBA-NCP_nvss-image.fits') as hdu:
wcs_woden = WCS(hdu[0].header).celestial
data_woden = np.squeeze(hdu[0].data)
ax2 = fig.add_subplot(212, projection=wcs_woden)
vmax=1e-1
im = ax2.imshow(data_woden, origin='lower', cmap='Blues', vmin=0, vmax=vmax)
plt.colorbar(im, ax=ax2, label='Flux (Jy/beam)')
ax2.grid()
ax2.set_title('Imaged NVSS catalogue sim')
ax1.set_xlabel(' ')
ax1.set_ylabel('Dec')
ax1.set_xticklabels([])
ax2.set_xlabel('RA')
ax2.set_ylabel('Dec')
half_width=300
cent_x = 512
cent_y = 512
for ax in [ax1, ax2]:
ax.set_xlim(cent_x - half_width, cent_x + half_width)
ax.set_ylim(cent_y - half_width, cent_y + half_width)
# plt.tight_layout()
plt.show()
WARNING: FITSFixedWarning: 'datfix' made the change 'Set MJD-OBS to 56946.730277 from DATE-OBS'. [astropy.wcs.wcs]
Huh so in terms of morphology, the NVSS sky model is better. Couple of spectral indexes are clearly off, where we have sources that are relatively too bright. But this just goes to highlight how difficult it is to make a good sky model.