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][https://skyview.gsfc.nasa.gov/current/cgi/query.pl], I’ve downloaded a WENSS image of the NCP at 325 MHz. 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.0029296875

I’ve downloaded the two catalogues from VizieR. The first thing I’ll do is crop them to a 3 degree radius around the NCP.

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)

<class 'numpy.ndarray'>

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)

axs.plot(ras_wenss, decs_wenss, 'C1o', mfc='none', transform=axs.get_transform('world'), label='WENSS')

axs.grid()
axs.set_xlabel('RA')
axs.set_ylabel('Dec')
plt.show()

Eh, what’s this? All of our positions are messed up! Time for an aside on plotting things centred on the NCP.

Plotting simulations centred at the NCP

Long story short, images centred at the NCP tend to render “upside down”. I’m going to do this from a fresh simulation and image it, so you can also run it and convince yourself. I’m going to use OSKAR in the following, to try and convince you this isn’t a WODEN bug. We’re just going to run a simulation centred at NCP, and check source positions. I should not this problem only happens at the NCP; everywhere else on the sky, the positions come out perfectly.

We’ll copy the OSKAR MWA telescope layout created using test_installation/everybeam/create_OSKAR-MWA_ms/make_oskar-mwa_beam.ipynb to run OSKAR.

cp -r ../../test_installation/everybeam/create_OSKAR-MWA_ms/MWA_phase1/ .

We’ll need to move the array location to be able to see the NCP though. Might as well chuck it at the North pole.

with open('MWA_phase1/position.txt', 'w') as outfile:
    outfile.write('0.0 90 0.0\n')

Next, we’ll make sky model with a number of point sources, close to the NCP. We’ll use those positions to test whether the WCS works with imaged simulations.

ras = [0.0, 90.0, 180.0, 270.0]
decs = [89.5, 89, 88.5, 88]

flux_I = 1
flux_Q = 0
flux_U = 0
flux_V = 0
ref_freq = 180e+6
SI = 0

with open('skymodel.osm', 'w') as outfile:
    for ra, dec in zip(ras, decs):
        outfile.write(f"{ra} {dec} {flux_I} {flux_Q} {flux_U} {flux_V} {ref_freq} {SI}\n")

Finally, we’ll write an OSKAR .ini script to simulate the observation centered on the NCP. We’ll use the sky model we just created.

with open('run_oskar_sim.ini', 'w') as outfile:
    outfile.write("[General]\n")
    outfile.write("app=oskar_sim_interferometer\n")
    outfile.write("[sky]\n")
    outfile.write("oskar_sky_model/file=skymodel.osm\n")
    outfile.write("[observation]\n")
    outfile.write("phase_centre_dec_deg=90.0\n")
    outfile.write("phase_centre_ra_deg=0.0\n")
    outfile.write("start_frequency_hz=1.8e+08\n")
    outfile.write("frequency_inc_hz=4.0e+04\n")
    outfile.write("start_time_utc=2024-07-21 20:13:00.0\n")
    outfile.write("length=2.0\n")
    outfile.write("[telescope]\n")
    outfile.write("input_directory=MWA_phase1\n")
    outfile.write("aperture_array/element_pattern/dipole_length_units=Metres\n")
    outfile.write("[interferometer]\n")
    outfile.write("channel_bandwidth_hz=40000.0\n")
    outfile.write("time_average_sec=2.0\n")
    outfile.write("ms_filename=OSKAR_NCP_obs.ms\n")

I then ran this on the command line as my local machine just has a container version of OSKAR.

%%bash
singularity exec --nv  /home/jline/software/OSKAR-2.8.3-Python3.sif \
    oskar_sim_interferometer run_oskar_sim.ini

W|                                                                   
W|== WARNING: No GPU capability available.
W|                                                                   
 |                                                                   
=|== Loaded settings file 'run_oskar_sim.ini'
 |                                                                   
 | + Simulator settings
 |   - Use double precision              : true
 |   - Use GPUs                          : true
 |   - CUDA device IDs to use            : all
 |   - Number of compute devices         : auto
 |   - Max. number of sources per chunk  : 16384
 | + Sky model settings
 |   - OSKAR sky model file settings
 |     * OSKAR sky model file(s)         : skymodel.osm
 | + Observation settings
 |   - Observation mode                  : Tracking
 |   - Phase centre RA [deg]             : 0.0
 |   - Phase centre Dec [deg]            : 90.0
 |   - Start frequency [Hz]              : 1.8e+08
 |   - Number of frequency channels      : 1
 |   - Frequency increment [Hz]          : 4.0e+04
 |   - Start time (UTC)                  : 2024-07-21 20:13:00.0
 |   - Observation length [sec, or hh:mm:ss]: 2.0
 |   - Number of time steps              : 1
 | + Telescope model settings
 |   - Input directory                   : MWA_phase1
 |   - Normalise beams at phase centre   : true
 |   - Allow station beam duplication    : false
 |   - Polarisation mode                 : Full
 |   - Station type                      : Aperture array
 |   - Aperture array settings
 |     * Element pattern settings
 |       + Dipole length units           : Metres
 | + Interferometer settings
 |   - Channel bandwidth [Hz]            : 40000.0
 |   - Time average [sec]                : 2.0
 |   - Max. time samples per block       : 8
 |   - Max. channels per block           : auto
 |   - Correlation type                  : Cross-correlations
 |   - Output OSKAR visibility file
 |   - Output Measurement Set            : OSKAR_NCP_obs.ms
 |   - Force polarised Measurement Set   : false
 |                                                                   
=|== Sky model set-up
 |                                                                   
 | + Loading OSKAR sky model file 'skymodel.osm' ...
 |   - done.
 |                                                                   
=|== Sky model summary
 |                                                                   
 | + Number of sources                   : 4
 | + Number of chunks                    : 1
 |                                                                   
=|== Telescope model summary
 |                                                                   
 | + Longitude [deg]                     : 0.000
 | + Latitude [deg]                      : 90.000
 | + Altitude [m]                        : 0
 | + Number of stations                  : 128
 | + Number of station models            : 128
 | + Max station size                    : 16
 | + Max station depth                   : 1
 |                                                                   
=|== Initial memory usage
 |                                                                   
 | + System memory is 47.4% (14.9 GB/31.3 GB) used.
 | + System memory used by current process: 231.0 MB.
 |                                                                   
=|== Starting simulation...
 |                                                                   
 |   - Time 1/1, Chunk 1/1, Channel 1/1 [Device 0, 4 sources]
 | + Block 1/1 (100%) complete. Simulation time elapsed: 0.018 s
 |                                                                   
=|== Final memory usage
 |                                                                   
 | + System memory is 47.2% (14.8 GB/31.3 GB) used.
 | + System memory used by current process: 127.0 MB.
 |                                                                   
=|== Simulation timing
 |                                                                   
 | + Total wall time      : 0.722 s
 | + Compute              : 0.012 s [Device 0]
 | + Compute              : 0.006 s [Device 1]
 | + Compute              : 0.006 s [Device 2]
 | + Compute              : 0.006 s [Device 3]
 | + Compute              : 0.005 s [Device 4]
 | + Compute              : 0.007 s [Device 5]
 | + Compute              : 0.007 s [Device 6]
 | + Compute              : 0.004 s [Device 7]
 | + Compute              : 0.001 s [Device 8]
 | + Compute              : 0.003 s [Device 9]
 | + Compute              : 0.002 s [Device 10]
 | + Compute              : 0.002 s [Device 11]
 | + Compute              : 0.002 s [Device 12]
 | + Compute              : 0.002 s [Device 13]
 | + Compute              : 0.002 s [Device 14]
 | + Compute              : 0.002 s [Device 15]
 | + Write                : 0.703 s
 | + Compute components:
 |   - Copy               : 10.3%
 |   - Horizon clip       :  0.2%
 |   - Jones E            :  7.4%
 |   - Jones K            :  0.1%
 |   - Jones join         :  0.1%
 |   - Jones correlate    :  2.2%
 |   - Other              : 79.7%
 |                                                                   
=|== Simulation complete
 |                                                                   
 | + Output(s):
 |   - Measurement Set    : OSKAR_NCP_obs.ms
 | + Run completed in 0.722 sec.
 |                                                                   
=|== OSKAR-2.8.3 ending at 2024-11-15, 10:02:04 (AEST).
 |                                                                   

OK, now let’s image things. Not actually going to bother cleaning as the dirty will be good enough here.

%%bash

###WARNING YOU PROBABLY DON'T NEED THIS ENV VARIABLE
###SOMETHING IS BORKED WITH MY WSCLEAN INSTALLATION
export OPENBLAS_NUM_THREADS=1

wsclean -name check_positions_NCP -size 256 256 -niter 0 \
  -auto-threshold 0.5 -auto-mask 3 \
  -pol I -multiscale -weight briggs 0 -scale 60asec -j 12 -mgain 0.85 \
  -no-update-model-required \
  OSKAR_NCP_obs.ms

WSClean version 3.0 (2021-08-26)
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  180-180 (1)
WARNING: This measurement set has no or an invalid WEIGHT_SPECTRUM column; will use less informative WEIGHT column.
Reordering OSKAR_NCP_obs.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 image ==
Precalculating weights for Briggs'(0) weighting...
Opening reordered part 0 spw 0 for OSKAR_NCP_obs.ms
Detected 31.3 GB of system memory, usage not limited.
Opening reordered part 0 spw 0 for OSKAR_NCP_obs.ms
Determining min and max w & theoretical beam size... DONE (w=[0:5.41033] lambdas, maxuvw=1577.28 lambda)
Theoretic beam = 2.18'
Minimal inversion size: 283 x 283, using optimal: 288 x 288
The theoretically suggested number of w-layers (1) is less than the number of availables
cores (12). Changing suggested number of w-layers to 12.
Suggested number of w-layers: 12
Will process 12/12 w-layers per pass.
Gridding pass 0... 
Rows that were required: 8128/0
Fourier transforms...
Gridded visibility count: 8128
Writing dirty image...
Rendering sources to restored image (beam is neither fitted nor estimated -- using delta scales!)... DONE
Writing restored image... DONE
Multi-scale cleaning summary:
Total: 0 components (0 Jy)
Inversion: 00:00:00.085352, prediction: 00:00:00, deconvolution: 00:00:00
Cleaning up temporary files...

OK, let’s read in our image, plot it, and overplot where our WCS says our sources should be.

with fits.open('check_positions_NCP-dirty.fits') as hdul:
    
    ##Images out of WSClean have 4 axes, as there is a time and freq axes
    ##Hence we use squeeze and celestial here to make everything 2D
    data = np.squeeze(hdul[0].data)
    header = hdul[0].header
    wcs = WCS(header).celestial
    
    
fig, axs = plt.subplots(1, 1, figsize=(5, 5), subplot_kw={'projection': wcs})
# fig, axs = plt.subplots(1, 1, figsize=(5, 5))

##Plot the data
axs.imshow(data, origin='lower', vmin=0, vmax=0.5)

##Now plot the sources as defined by the WCS
x, y = wcs.all_world2pix(ras, decs, 0)

axs.plot(x, y, 'C1o', mfc='none', label='Expected position')

axs.legend()

axs.grid()
plt.show()

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

OH NO, our sources aren’t where they are supposed to be! This comes down to slightly opaque definition of RA origin when we’re at the NCP. This isn’t an imshow “origin=lower’ mistake as well, as it turns out that is set automatically when you use the wcs in projection:

fig, axs = plt.subplots(1, 1, figsize=(5, 5), subplot_kw={'projection': wcs})
axs.imshow(data, vmin=0, vmax=0.5)
axs.plot(x, y, 'C1o', mfc='none', label='Expected position')
axs.legend()
axs.grid()
plt.show()

If you open this via DS9 or CARTA, you see the same positional differences. If you open things with kvis however, you actually get the WCS loaded up differently.

import matplotlib.image as mpimg
img = mpimg.imread('kvis_screenshot.png')
fig, axs = plt.subplots(1, 1, figsize=(6, 6))

axs.imshow(img)
plt.axis('off')  # Hide the axis
plt.show()

You can’t see my mouse, but you can see that Ra,Dec = 12h, 88d30m corresponds to x, y = 128,218. This is not the same as was astropy reports:

x, y = wcs.all_world2pix(180, 88.5, 0)
print(int(x), int(y))

128 38

Essentially, the image out of the simulation comes out with RA=180 pointed to the top, whereas the WCS sticks 180 pointed down. Somehow, kvis knows to put RA=180 at the top. Some ancient wisdom there that has been lost to time.

ANYWAYS, to plot anything that comes out of WODEN or OSKAR from the NCP, imaged through WSClean, you have to make one edit to the WCS, and then everything is fine:

from copy import deepcopy

header_edit = deepcopy(header)
# ##Just change this one thing
header_edit['CRVAL1'] = 180
wcs_edit = WCS(header_edit).celestial

x, y = wcs_edit.all_world2pix(ras, decs, 0)

fig, axs = plt.subplots(1, 1, figsize=(5, 5), subplot_kw={'projection': wcs_edit})
axs.imshow(data, vmin=0, vmax=0.5)
axs.plot(x, y, 'C1o', mfc='none', label='Expected position')
axs.legend()
axs.grid()
plt.show()

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

So, you just need to change [‘CRVAL1’] = 180, and then everything is fine. Dagnabbit.

Run the simulation

After all that, let’s check our sky model is legit by modifying our wcs and plotting our RA,Decs

header_wenss_edit = deepcopy(header_wenss)
header_wenss_edit['CRVAL1'] = 180
wcs_wenss_edit = WCS(header_wenss_edit)


fig, axs = plt.subplots(1, 1, figsize=(10, 10), subplot_kw={'projection': wcs_wenss_edit})

axs.imshow(data_wenss, origin='lower', cmap='Blues', vmin=0, vmax=0.05)

axs.plot(ras_wenss, decs_wenss, 'C1o', mfc='none', transform=axs.get_transform('world'), label='WENSS')

axs.grid()
axs.set_xlabel('RA')
axs.set_ylabel('Dec')
plt.show()

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} \
    --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/jline/software/anaconda3/envs/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.3.0')

Successful readonly open of default-locked table ../../test_installation/everybeam/lba.MS/ANTENNA: 10 columns, 37 rows
You are using WODEN commit: No git describe nor __version__ avaible
LOADING IN /home/jline/software/WODEN_dev/wodenpy/libwoden_double.so
Setting phase centre RA,DEC 0.00000deg 90.00000deg
Obs epoch initial LST was 294.9211333249 deg
Setting initial J2000 LST to 294.8284574349 deg
Setting initial mjd to 56946.7302777780
After precession initial latitude of the array is 52.8683555426 deg
We are precessing the array
Doing the initial reading/mapping of sky model into chunks
INFO: couldn't find second table containing shapelet information, so not attempting to load any shapelets.
Mapping took 0.0 mins
Have read in 602 components
After cropping there are 602 components
From Set 0 thread num 0 reading 0 points, 76 gauss, 0 shape, 0 shape coeffs using thread id 0
From Set 0 thread num 1 reading 0 points, 76 gauss, 0 shape, 0 shape coeffs using thread id 1
From Set 0 thread num 2 reading 0 points, 76 gauss, 0 shape, 0 shape coeffs using thread id 2
From Set 0 thread num 3 reading 0 points, 76 gauss, 0 shape, 0 shape coeffs using thread id 3
From Set 0 thread num 4 reading 0 points, 76 gauss, 0 shape, 0 shape coeffs using thread id 4
From Set 0 thread num 5 reading 0 points, 76 gauss, 0 shape, 0 shape coeffs using thread id 5
From Set 0 thread num 7 reading 0 points, 70 gauss, 0 shape, 0 shape coeffs using thread id 7From Set 0 thread num 6 reading 0 points, 76 gauss, 0 shape, 0 shape coeffs using thread id 6

Finshed thread 2 in 68.3 seconds
Finshed thread 1 in 70.9 seconds
Finshed thread 3 in 72.2 seconds
Finshed thread 7 in 72.5 seconds
Finshed thread 0 in 73.1 seconds
Finshed thread 5 in 73.5 seconds
Finshed thread 6 in 73.6 seconds
Finshed thread 4 in 76.2 seconds
Sending set 0 to GPU
Set 0 has returned from the GPU after 1.1 seconds

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.7410455463 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.7410455463 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 2048 2048 -scale 0.0029296875 \
   -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.0 (2021-08-26)
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 31.3 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:10286.7] lambdas, maxuvw=17187.2 lambda)
Theoretic beam = 12''
Small inversion enabled, but inversion resolution already smaller than beam size: not using optimization.
Suggested number of w-layers: 178
Will process 178/178 w-layers per pass.
Gridding pass 0... 
Rows that were required: 15792/0
Fourier transforms...
Gridded visibility count: 269273
Writing psf image... DONE
 == Constructing image ==
Opening reordered part 0 spw 0 for woden_LBA-NCP_band01.ms
Gridding pass 0... Fourier transforms...
Writing dirty image...
 == Deconvolving (1) ==
Estimated standard deviation of background noise: 15.47 mJy
Scale info:
- Scale 0, bias factor=1, psfpeak=1, gain=0.1, kernel peak=0.19888
- Scale 5, bias factor=1.7, psfpeak=0.395592, gain=0.252786, kernel peak=0.0815272
- Scale 9, bias factor=2.8, psfpeak=0.259315, gain=0.385631, kernel peak=0.0330771
- Scale 18, bias factor=4.6, psfpeak=0.132431, gain=0.755112, kernel peak=0.00966831
- Scale 36, bias factor=7.7, psfpeak=0.0687231, gain=1.45512, kernel peak=0.00287831
- Scale 73, bias factor=12.9, psfpeak=0.0234852, gain=4.258, kernel peak=0.000791489
- Scale 146, bias factor=21.4, psfpeak=0.00621147, gain=16.0992, kernel peak=0.000207982
- Scale 291, bias factor=35.7, psfpeak=0.00154831, gain=64.5865, kernel peak=5.26599e-05
- Scale 583, bias factor=59.5, psfpeak=0.000351752, gain=284.291, kernel peak=1.32493e-05
RMS per scale: {0: 21.3 mJy, 5: 16.52 mJy, 9: 14.12 mJy, 18: 11.02 mJy, 36: 7.91 mJy, 73: 4.97 mJy, 146: 2.81 mJy, 291: 1.49 mJy, 583: 709.59 µJy}
Starting multi-scale cleaning. Start peak=8.56 Jy, major iteration threshold=1.28 Jy
Iteration 3, scale 9 px : 6.84 Jy at 965,857
Iteration 6, scale 0 px : 6.1 Jy at 965,857
Iteration 9, scale 0 px : 4.45 Jy at 965,857
Iteration 12, scale 0 px : 3.53 Jy at 966,857
Iteration 16, scale 0 px : 2.74 Jy at 965,857
Iteration 20, scale 9 px : 2.34 Jy at 930,1360
Iteration 24, scale 9 px : 1.84 Jy at 965,856
Iteration 31, scale 0 px : 1.91 Jy at 966,856
Iteration 36, scale 0 px : 1.49 Jy at 930,1359
Subminor loop is near minor loop threshold. Initiating countdown.
(12) Iteration 41, scale 9 px : 1.46 Jy at 990,854
(11) Iteration 43, scale 9 px : 1.26 Jy at 991,854
Minor loop finished, continuing cleaning after inversion/prediction round.
 == Converting model image to visibilities ==
Opening reordered part 0 spw 0 for woden_LBA-NCP_band01.ms
Fourier transforms for pass 0... Predicting...
 == Constructing image ==
Opening reordered part 0 spw 0 for woden_LBA-NCP_band01.ms
Gridding pass 0... Fourier transforms...
 == Deconvolving (2) ==
Estimated standard deviation of background noise: 6.77 mJy
Scale info:
- Scale 0, bias factor=1, psfpeak=1, gain=0.1, kernel peak=0.19888
- Scale 5, bias factor=1.7, psfpeak=0.395592, gain=0.252786, kernel peak=0.0815272
- Scale 9, bias factor=2.8, psfpeak=0.259315, gain=0.385631, kernel peak=0.0330771
- Scale 18, bias factor=4.6, psfpeak=0.132431, gain=0.755112, kernel peak=0.00966831
- Scale 36, bias factor=7.7, psfpeak=0.0687231, gain=1.45512, kernel peak=0.00287831
- Scale 73, bias factor=12.9, psfpeak=0.0234852, gain=4.258, kernel peak=0.000791489
- Scale 146, bias factor=21.4, psfpeak=0.00621147, gain=16.0992, kernel peak=0.000207982
- Scale 291, bias factor=35.7, psfpeak=0.00154831, gain=64.5865, kernel peak=5.26599e-05
- Scale 583, bias factor=59.5, psfpeak=0.000351752, gain=284.291, kernel peak=1.32493e-05
RMS per scale: {0: 8.53 mJy, 5: 6.12 mJy, 9: 5.02 mJy, 18: 3.82 mJy, 36: 2.7 mJy, 73: 1.61 mJy, 146: 884.48 µJy, 291: 481.34 µJy, 583: 246.54 µJy}
Starting multi-scale cleaning. Start peak=1.27 Jy, major iteration threshold=190.53 mJy
Iteration 51, scale 9 px : 1.18 Jy at 991,854
Iteration 56, scale 9 px : 945.15 mJy at 991,854
Iteration 69, scale 0 px : 951.87 mJy at 965,856
Iteration 76, scale 0 px : 753.13 mJy at 965,857
Iteration 89, scale 18 px : 769.74 mJy at 1145,759
Iteration 96, scale 9 px : 682.5 mJy at 1237,467
Iteration 109, scale 36 px : -626.44 mJy at 964,857
Iteration 113, scale 0 px : 639.81 mJy at 964,857
Iteration 122, scale 0 px : 507.11 mJy at 965,856
Iteration 141, scale 18 px : -705.35 mJy at 965,858
Iteration 144, scale 18 px : -561.36 mJy at 965,858
Iteration 151, scale 9 px : 512.62 mJy at 1237,467
Iteration 163, scale 9 px : 409.16 mJy at 1237,467
Iteration 204, scale 0 px : 500.67 mJy at 964,857
Iteration 209, scale 0 px : 389.64 mJy at 930,1359
Iteration 228, scale 18 px : -598.18 mJy at 965,857
Iteration 231, scale 18 px : -476.06 mJy at 965,857
Iteration 234, scale 18 px : -378.88 mJy at 965,857
Iteration 240, scale 5 px : 304.55 mJy at 1237,467
Iteration 274, scale 0 px : 366.61 mJy at 965,856
Iteration 281, scale 0 px : 282.97 mJy at 871,836
Iteration 332, scale 18 px : -441.03 mJy at 965,857
Iteration 335, scale 18 px : -351.25 mJy at 965,857
Iteration 341, scale 36 px : 295.76 mJy at 554,965
Iteration 347, scale 9 px : 281.88 mJy at 836,1039
Iteration 381, scale 18 px : -280.48 mJy at 966,857
Iteration 385, scale 18 px : 224.26 mJy at 906,999
Subminor loop is near minor loop threshold. Initiating countdown.
(12) Iteration 399, scale 0 px : 301.58 mJy at 964,857
Iteration 404, scale 0 px : 237.82 mJy at 965,858
(11) Iteration 447, scale 9 px : -259.76 mJy at 963,854
Iteration 461, scale 9 px : 206.79 mJy at 1016,885
(10) Iteration 490, scale 18 px : -237.68 mJy at 965,857
(9) Iteration 495, scale 36 px : -205.39 mJy at 930,1359
(8) Iteration 497, scale 0 px : 222.75 mJy at 966,856
(7) Iteration 502, scale 0 px : 187.59 mJy at 977,877
Minor loop finished, continuing cleaning after inversion/prediction round.
 == Converting model image to visibilities ==
Opening reordered part 0 spw 0 for woden_LBA-NCP_band01.ms
Fourier transforms for pass 0... Predicting...
 == Constructing image ==
Opening reordered part 0 spw 0 for woden_LBA-NCP_band01.ms
Gridding pass 0... Fourier transforms...
 == Deconvolving (3) ==
Estimated standard deviation of background noise: 3.8 mJy
Scale info:
- Scale 0, bias factor=1, psfpeak=1, gain=0.1, kernel peak=0.19888
- Scale 5, bias factor=1.7, psfpeak=0.395592, gain=0.252786, kernel peak=0.0815272
- Scale 9, bias factor=2.8, psfpeak=0.259315, gain=0.385631, kernel peak=0.0330771
- Scale 18, bias factor=4.6, psfpeak=0.132431, gain=0.755112, kernel peak=0.00966831
- Scale 36, bias factor=7.7, psfpeak=0.0687231, gain=1.45512, kernel peak=0.00287831
- Scale 73, bias factor=12.9, psfpeak=0.0234852, gain=4.258, kernel peak=0.000791489
- Scale 146, bias factor=21.4, psfpeak=0.00621147, gain=16.0992, kernel peak=0.000207982
- Scale 291, bias factor=35.7, psfpeak=0.00154831, gain=64.5865, kernel peak=5.26599e-05
- Scale 583, bias factor=59.5, psfpeak=0.000351752, gain=284.291, kernel peak=1.32493e-05
RMS per scale: {0: 4.49 mJy, 5: 3.14 mJy, 9: 2.56 mJy, 18: 1.9 mJy, 36: 1.31 mJy, 73: 756.61 µJy, 146: 400.52 µJy, 291: 197.14 µJy, 583: 98.79 µJy}
Starting multi-scale cleaning. Start peak=-246.02 mJy, major iteration threshold=36.9 mJy
Iteration 505, scale 18 px : -195.8 mJy at 965,857
Iteration 526, scale 0 px : 211.81 mJy at 965,858
Iteration 560, scale 9 px : 186.16 mJy at 1145,759
Iteration 659, scale 0 px : 176.87 mJy at 964,856
Iteration 726, scale 36 px : -205.96 mJy at 930,1359
Iteration 730, scale 36 px : -158.54 mJy at 930,1359
Iteration 742, scale 18 px : -172.75 mJy at 965,857
Iteration 748, scale 18 px : -137.57 mJy at 965,857
Iteration 817, scale 0 px : 163.14 mJy at 930,1359
Iteration 861, scale 0 px : -140.18 mJy at 962,857
Iteration 969, scale 9 px : -145.9 mJy at 995,855
Iteration 1003, scale 9 px : 116.27 mJy at 1306,1187
Iteration 1187, scale 0 px : 124.26 mJy at 929,1359
Iteration 1240, scale 18 px : -146.18 mJy at 930,1359
Iteration 1243, scale 18 px : -116.33 mJy at 930,1359
Iteration 1250, scale 18 px : 92.9 mJy at 1326,654
Iteration 1306, scale 0 px : 119.35 mJy at 931,1360
Iteration 1341, scale 0 px : -96.93 mJy at 964,855
Iteration 1546, scale 9 px : -104.73 mJy at 987,853
Iteration 1574, scale 9 px : 83.45 mJy at 972,852
Iteration 1779, scale 0 px : 89.14 mJy at 990,853
Iteration 1844, scale 36 px : 85.46 mJy at 959,847
Iteration 1850, scale 73 px : 68.82 mJy at 961,853
Iteration 1854, scale 18 px : -103.04 mJy at 930,1359
Iteration 1858, scale 18 px : -82.02 mJy at 930,1359
Iteration 1868, scale 9 px : -84.69 mJy at 965,856
Iteration 1873, scale 0 px : 81.21 mJy at 929,1360
Iteration 1980, scale 0 px : -70.24 mJy at 933,1360
Iteration 2215, scale 36 px : 85.67 mJy at 969,857
Iteration 2227, scale 18 px : -75.8 mJy at 978,880
Iteration 2245, scale 9 px : -75.95 mJy at 1141,761
Iteration 2277, scale 9 px : 63.08 mJy at 971,851
Iteration 2464, scale 18 px : -58.47 mJy at 930,1359
Iteration 2482, scale 0 px : 64.29 mJy at 978,877
Iteration 2604, scale 0 px : -56.08 mJy at 932,1361
Iteration 2855, scale 36 px : 69.37 mJy at 965,852
Iteration 2865, scale 36 px : 55.03 mJy at 969,858
Iteration 2891, scale 9 px : -59.02 mJy at 1149,757
Iteration 2930, scale 9 px : 50.55 mJy at 973,858
Iteration 3121, scale 0 px : 50.58 mJy at 1145,760
Iteration 3299, scale 0 px : -43.2 mJy at 988,855
Subminor loop is near minor loop threshold. Initiating countdown.
(12) Iteration 3521, scale 36 px : 55.71 mJy at 969,858
Iteration 3530, scale 18 px : -55.42 mJy at 978,879
Iteration 3539, scale 18 px : -44.13 mJy at 978,879
(11) Iteration 3577, scale 9 px : -44.92 mJy at 1149,757
(10) Iteration 3647, scale 5 px : 41.38 mJy at 169,1030
(9) Iteration 3653, scale 0 px : 42.17 mJy at 978,877
(8) Iteration 3678, scale 36 px : 42.09 mJy at 971,859
(7) Iteration 3688, scale 9 px : -40.28 mJy at 973,879
(6) Iteration 3697, scale 18 px : -41.24 mJy at 872,835
(5) Iteration 3701, scale 0 px : 38.25 mJy at 1124,814
(4) Iteration 3706, scale 18 px : -39.48 mJy at 1125,815
(3) Iteration 3708, scale 36 px : 38.12 mJy at 971,859
(2) Iteration 3709, scale 9 px : -38.24 mJy at 965,857
(1) Iteration 3710, scale 36 px : 38 mJy at 970,859
(0) Iteration 3711, scale 5 px : -37.46 mJy at 987,853
Minor loop finished, continuing cleaning after inversion/prediction round.
 == Converting model image to visibilities ==
Opening reordered part 0 spw 0 for woden_LBA-NCP_band01.ms
Fourier transforms for pass 0... Predicting...
 == Constructing image ==
Opening reordered part 0 spw 0 for woden_LBA-NCP_band01.ms
Gridding pass 0... Fourier transforms...
 == Deconvolving (4) ==
Estimated standard deviation of background noise: 2.34 mJy
Scale info:
- Scale 0, bias factor=1, psfpeak=1, gain=0.1, kernel peak=0.19888
- Scale 5, bias factor=1.7, psfpeak=0.395592, gain=0.252786, kernel peak=0.0815272
- Scale 9, bias factor=2.8, psfpeak=0.259315, gain=0.385631, kernel peak=0.0330771
- Scale 18, bias factor=4.6, psfpeak=0.132431, gain=0.755112, kernel peak=0.00966831
- Scale 36, bias factor=7.7, psfpeak=0.0687231, gain=1.45512, kernel peak=0.00287831
- Scale 73, bias factor=12.9, psfpeak=0.0234852, gain=4.258, kernel peak=0.000791489
- Scale 146, bias factor=21.4, psfpeak=0.00621147, gain=16.0992, kernel peak=0.000207982
- Scale 291, bias factor=35.7, psfpeak=0.00154831, gain=64.5865, kernel peak=5.26599e-05
- Scale 583, bias factor=59.5, psfpeak=0.000351752, gain=284.291, kernel peak=1.32493e-05
RMS per scale: {0: 2.51 mJy, 5: 1.72 mJy, 9: 1.35 mJy, 18: 933.09 µJy, 36: 589.85 µJy, 73: 326.14 µJy, 146: 182.91 µJy, 291: 95.89 µJy, 583: 46.95 µJy}
Starting multi-scale cleaning. Start peak=44.59 mJy, major iteration threshold=7.02 mJy (final)
Iteration 3807, scale 18 px : -39.8 mJy at 1125,816
Iteration 3870, scale 9 px : 39.07 mJy at 1752,1288
Iteration 4071, scale 0 px : 40.08 mJy at 964,856
Iteration 4244, scale 0 px : -34.32 mJy at 963,857
Iteration 4692, scale 18 px : -43.8 mJy at 1125,816
Iteration 4699, scale 36 px : -39.05 mJy at 555,967
Iteration 4712, scale 9 px : 33.85 mJy at 922,1358
Iteration 4890, scale 5 px : 31.78 mJy at 169,1029
Iteration 4945, scale 18 px : -32.44 mJy at 1125,816
Iteration 4978, scale 9 px : -29.7 mJy at 177,1030
Iteration 5240, scale 0 px : 31.76 mJy at 1125,813
Iteration 5398, scale 0 px : 26.83 mJy at 965,863
Iteration 6036, scale 18 px : -41.03 mJy at 1125,815
Iteration 6039, scale 36 px : -34.5 mJy at 555,966
Iteration 6052, scale 18 px : -31.21 mJy at 1125,815
Iteration 6078, scale 73 px : 25.03 mJy at 927,1359
Iteration 6088, scale 36 px : -26.7 mJy at 989,783
Iteration 6164, scale 9 px : -26.28 mJy at 835,1035
Iteration 6387, scale 0 px : 25.47 mJy at 1124,813
Iteration 6579, scale 18 px : -32.02 mJy at 1125,815
Iteration 6583, scale 18 px : -25.51 mJy at 1126,815
Iteration 6648, scale 9 px : -23.06 mJy at 1148,756
Iteration 7030, scale 0 px : 23.47 mJy at 1125,815
Iteration 7363, scale 36 px : -23.45 mJy at 998,1311
Iteration 7410, scale 5 px : -21.74 mJy at 987,853
Iteration 7515, scale 18 px : -24.41 mJy at 1125,815
Iteration 7529, scale 18 px : 19.51 mJy at 516,991
Iteration 7851, scale 0 px : 21.77 mJy at 1126,814
Iteration 8176, scale 9 px : -20.11 mJy at 1109,1363
Iteration 8504, scale 0 px : 18.64 mJy at 966,858
Iteration 9408, scale 36 px : -23.01 mJy at 554,965
Iteration 9426, scale 36 px : -18.37 mJy at 868,496
Iteration 9556, scale 9 px : -17.88 mJy at 718,1654
Iteration 10011, scale 5 px : 16 mJy at 779,253
Iteration 10354, scale 18 px : -17.17 mJy at 172,1048
Iteration 10558, scale 0 px : 17.17 mJy at 1124,814
Iteration 10881, scale 36 px : -18.97 mJy at 555,965
Iteration 10897, scale 9 px : -15.95 mJy at 718,1655
Iteration 11489, scale 36 px : 14.96 mJy at 1151,1701
Iteration 11772, scale 0 px : 14.77 mJy at 552,970
Iteration 12921, scale 18 px : -17.4 mJy at 808,1332
Iteration 12955, scale 18 px : -13.86 mJy at 808,1332
Iteration 13935, scale 9 px : 14.69 mJy at 998,856
Iteration 14365, scale 0 px : 13.45 mJy at 1232,771
Iteration 15255, scale 36 px : -15.48 mJy at 555,964
Iteration 15281, scale 36 px : -12.21 mJy at 1046,1423
Iteration 15961, scale 18 px : -13.87 mJy at 581,1084
Iteration 16028, scale 9 px : -12.93 mJy at 719,1664
Iteration 17610, scale 0 px : 11.95 mJy at 556,961
Iteration 19143, scale 18 px : -13.56 mJy at 872,496
Iteration 19232, scale 9 px : -12.12 mJy at 1109,1362
Iteration 20000, scale 73 px : -11.6 mJy at 306,481
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.
Writing model image...
 == Converting model image to visibilities ==
Opening reordered part 0 spw 0 for woden_LBA-NCP_band01.ms
Fourier transforms for pass 0... Predicting...
 == Constructing image ==
Opening reordered part 0 spw 0 for woden_LBA-NCP_band01.ms
Gridding pass 0... Fourier transforms...
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: 9133 (28.51 Jy)
- Scale 5 px, nr of components cleaned: 543 (1.69 Jy)
- Scale 9 px, nr of components cleaned: 6665 (27.76 Jy)
- Scale 18 px, nr of components cleaned: 2240 (392.94 mJy)
- Scale 36 px, nr of components cleaned: 1405 (-802.41 mJy)
- Scale 73 px, nr of components cleaned: 14 (99.63 mJy)
- Scale 146 px, nr of components cleaned: 0 (0 Jy)
- Scale 291 px, nr of components cleaned: 0 (0 Jy)
- Scale 583 px, nr of components cleaned: 0 (0 Jy)
Total: 20000 components (57.64 Jy)
Inversion: 00:01:55.010394, prediction: 00:01:05.742549, deconvolution: 00:01:31.751814
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))


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=600
cent_x = 1024
cent_y = 1024

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} \
    --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/jline/software/anaconda3/envs/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.3.0')

Successful readonly open of default-locked table ../../test_installation/everybeam/lba.MS/ANTENNA: 10 columns, 37 rows
You are using WODEN commit: No git describe nor __version__ avaible
LOADING IN /home/jline/software/WODEN_dev/wodenpy/libwoden_double.so
Setting phase centre RA,DEC 0.00000deg 90.00000deg
Obs epoch initial LST was 294.9211333249 deg
Setting initial J2000 LST to 294.8284574349 deg
Setting initial mjd to 56946.7302777780
After precession initial latitude of the array is 52.8683555426 deg
We are precessing the array
Doing the initial reading/mapping of sky model into chunks
INFO: couldn't find second table containing shapelet information, so not attempting to load any shapelets.
Mapping took 0.0 mins
Have read in 1272 components
After cropping there are 1272 components
From Set 0 thread num 0 reading 0 points, 159 gauss, 0 shape, 0 shape coeffs using thread id 0
From Set 0 thread num 1 reading 0 points, 159 gauss, 0 shape, 0 shape coeffs using thread id 1
From Set 0 thread num 3 reading 0 points, 159 gauss, 0 shape, 0 shape coeffs using thread id 3
From Set 0 thread num 4 reading 0 points, 159 gauss, 0 shape, 0 shape coeffs using thread id 4
From Set 0 thread num 2 reading 0 points, 159 gauss, 0 shape, 0 shape coeffs using thread id 2
From Set 0 thread num 5 reading 0 points, 159 gauss, 0 shape, 0 shape coeffs using thread id 5
From Set 0 thread num 6 reading 0 points, 159 gauss, 0 shape, 0 shape coeffs using thread id 6
From Set 0 thread num 7 reading 0 points, 159 gauss, 0 shape, 0 shape coeffs using thread id 7
Finshed thread 1 in 128.1 seconds
Finshed thread 0 in 128.9 seconds
Finshed thread 7 in 130.1 seconds
Finshed thread 6 in 132.7 seconds
Finshed thread 5 in 140.1 seconds
Finshed thread 2 in 144.2 seconds
Finshed thread 4 in 146.5 seconds
Finshed thread 3 in 146.7 seconds
Sending set 0 to GPU
Set 0 has returned from the GPU after 1.5 seconds

make_ms('woden_LBA-NCP_nvss_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.7410455463 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.7410455463 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

##YOU PROBABLY DON'T NEED THIS ENV VARIABLE
export OPENBLAS_NUM_THREADS=1

wsclean -name woden_LBA-NCP_nvss -niter 20000 \
   -size 2048 2048 -scale 0.0029296875 \
   -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.0 (2021-08-26)
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 31.3 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:10286.7] lambdas, maxuvw=17187.2 lambda)
Theoretic beam = 12''
Small inversion enabled, but inversion resolution already smaller than beam size: not using optimization.
Suggested number of w-layers: 178
Will process 178/178 w-layers per pass.
Gridding pass 0... 
Rows that were required: 15792/0
Fourier transforms...
Gridded visibility count: 269273
Writing psf image... DONE
 == Constructing image ==
Opening reordered part 0 spw 0 for woden_LBA-NCP_nvss_band01.ms
Gridding pass 0... Fourier transforms...
Writing dirty image...
 == Deconvolving (1) ==
Estimated standard deviation of background noise: 20.79 mJy
Scale info:
- Scale 0, bias factor=1, psfpeak=1, gain=0.1, kernel peak=0.19888
- Scale 5, bias factor=1.7, psfpeak=0.395592, gain=0.252786, kernel peak=0.0815272
- Scale 9, bias factor=2.8, psfpeak=0.259315, gain=0.385631, kernel peak=0.0330771
- Scale 18, bias factor=4.6, psfpeak=0.132431, gain=0.755112, kernel peak=0.00966831
- Scale 36, bias factor=7.7, psfpeak=0.0687231, gain=1.45512, kernel peak=0.00287831
- Scale 73, bias factor=12.9, psfpeak=0.0234852, gain=4.258, kernel peak=0.000791489
- Scale 146, bias factor=21.4, psfpeak=0.00621147, gain=16.0992, kernel peak=0.000207982
- Scale 291, bias factor=35.7, psfpeak=0.00154831, gain=64.5865, kernel peak=5.26599e-05
- Scale 583, bias factor=59.5, psfpeak=0.000351752, gain=284.291, kernel peak=1.32493e-05
RMS per scale: {0: 28.64 mJy, 5: 22.14 mJy, 9: 18.89 mJy, 18: 14.74 mJy, 36: 10.48 mJy, 73: 6.28 mJy, 146: 3.4 mJy, 291: 1.7 mJy, 583: 765.98 µJy}
Starting multi-scale cleaning. Start peak=12.14 Jy, major iteration threshold=1.82 Jy
Iteration 3, scale 9 px : 9.7 Jy at 965,857
Iteration 6, scale 0 px : 8.59 Jy at 965,857
Iteration 9, scale 0 px : 6.59 Jy at 965,857
Iteration 12, scale 0 px : 5.05 Jy at 965,856
Iteration 16, scale 0 px : 3.82 Jy at 966,857
Iteration 20, scale 0 px : 2.91 Jy at 966,856
Iteration 24, scale 0 px : 2.25 Jy at 966,856
Subminor loop is near minor loop threshold. Initiating countdown.
(12) Iteration 28, scale 0 px : 1.8 Jy at 964,857
Minor loop finished, continuing cleaning after inversion/prediction round.
 == Converting model image to visibilities ==
Opening reordered part 0 spw 0 for woden_LBA-NCP_nvss_band01.ms
Fourier transforms for pass 0... Predicting...
 == Constructing image ==
Opening reordered part 0 spw 0 for woden_LBA-NCP_nvss_band01.ms
Gridding pass 0... Fourier transforms...
 == Deconvolving (2) ==
Estimated standard deviation of background noise: 7.87 mJy
Scale info:
- Scale 0, bias factor=1, psfpeak=1, gain=0.1, kernel peak=0.19888
- Scale 5, bias factor=1.7, psfpeak=0.395592, gain=0.252786, kernel peak=0.0815272
- Scale 9, bias factor=2.8, psfpeak=0.259315, gain=0.385631, kernel peak=0.0330771
- Scale 18, bias factor=4.6, psfpeak=0.132431, gain=0.755112, kernel peak=0.00966831
- Scale 36, bias factor=7.7, psfpeak=0.0687231, gain=1.45512, kernel peak=0.00287831
- Scale 73, bias factor=12.9, psfpeak=0.0234852, gain=4.258, kernel peak=0.000791489
- Scale 146, bias factor=21.4, psfpeak=0.00621147, gain=16.0992, kernel peak=0.000207982
- Scale 291, bias factor=35.7, psfpeak=0.00154831, gain=64.5865, kernel peak=5.26599e-05
- Scale 583, bias factor=59.5, psfpeak=0.000351752, gain=284.291, kernel peak=1.32493e-05
RMS per scale: {0: 9.78 mJy, 5: 6.85 mJy, 9: 5.55 mJy, 18: 4.17 mJy, 36: 2.93 mJy, 73: 1.7 mJy, 146: 886.69 µJy, 291: 422.75 µJy, 583: 179.92 µJy}
Starting multi-scale cleaning. Start peak=1.8 Jy, major iteration threshold=269.66 mJy
Iteration 33, scale 0 px : 1.4 Jy at 964,857
Iteration 40, scale 9 px : 1.44 Jy at 990,854
Iteration 46, scale 9 px : 1.09 Jy at 930,1360
Iteration 56, scale 0 px : 1.06 Jy at 966,856
Iteration 63, scale 18 px : -812.82 mJy at 965,860
Iteration 69, scale 0 px : 851.89 mJy at 964,857
Iteration 81, scale 9 px : -743.89 mJy at 963,861
Iteration 96, scale 9 px : 594.56 mJy at 930,1360
Iteration 122, scale 0 px : 738.71 mJy at 964,857
Iteration 127, scale 0 px : 574.66 mJy at 930,1359
Iteration 140, scale 36 px : -821.32 mJy at 965,857
Iteration 144, scale 36 px : -632.2 mJy at 965,857
Iteration 148, scale 36 px : -486.74 mJy at 965,857
Iteration 159, scale 0 px : 550.28 mJy at 964,857
Iteration 167, scale 0 px : 439.06 mJy at 1540,689
Iteration 199, scale 18 px : -664.73 mJy at 965,857
Iteration 202, scale 18 px : -529.03 mJy at 965,857
Iteration 207, scale 9 px : 435.53 mJy at 836,1039
Iteration 227, scale 18 px : -422.67 mJy at 965,857
Iteration 232, scale 18 px : -336.71 mJy at 965,857
Subminor loop is near minor loop threshold. Initiating countdown.
(12) Iteration 260, scale 0 px : 448.45 mJy at 966,856
Iteration 265, scale 0 px : 348.82 mJy at 966,856
Iteration 294, scale 9 px : -368.3 mJy at 969,857
Iteration 304, scale 9 px : 294.04 mJy at 836,1039
(11) Iteration 323, scale 18 px : -330.05 mJy at 965,857
(10) Iteration 326, scale 0 px : 329.05 mJy at 966,856
(9) Iteration 333, scale 0 px : 267.06 mJy at 630,921
Minor loop finished, continuing cleaning after inversion/prediction round.
 == Converting model image to visibilities ==
Opening reordered part 0 spw 0 for woden_LBA-NCP_nvss_band01.ms
Fourier transforms for pass 0... Predicting...
 == Constructing image ==
Opening reordered part 0 spw 0 for woden_LBA-NCP_nvss_band01.ms
Gridding pass 0... Fourier transforms...
 == Deconvolving (3) ==
Estimated standard deviation of background noise: 4.53 mJy
Scale info:
- Scale 0, bias factor=1, psfpeak=1, gain=0.1, kernel peak=0.19888
- Scale 5, bias factor=1.7, psfpeak=0.395592, gain=0.252786, kernel peak=0.0815272
- Scale 9, bias factor=2.8, psfpeak=0.259315, gain=0.385631, kernel peak=0.0330771
- Scale 18, bias factor=4.6, psfpeak=0.132431, gain=0.755112, kernel peak=0.00966831
- Scale 36, bias factor=7.7, psfpeak=0.0687231, gain=1.45512, kernel peak=0.00287831
- Scale 73, bias factor=12.9, psfpeak=0.0234852, gain=4.258, kernel peak=0.000791489
- Scale 146, bias factor=21.4, psfpeak=0.00621147, gain=16.0992, kernel peak=0.000207982
- Scale 291, bias factor=35.7, psfpeak=0.00154831, gain=64.5865, kernel peak=5.26599e-05
- Scale 583, bias factor=59.5, psfpeak=0.000351752, gain=284.291, kernel peak=1.32493e-05
RMS per scale: {0: 5.41 mJy, 5: 3.83 mJy, 9: 3.15 mJy, 18: 2.38 mJy, 36: 1.64 mJy, 73: 945.68 µJy, 146: 493.39 µJy, 291: 232.98 µJy, 583: 104.89 µJy}
Starting multi-scale cleaning. Start peak=-359.17 mJy, major iteration threshold=53.88 mJy
Iteration 336, scale 18 px : -286.04 mJy at 965,857
Iteration 352, scale 0 px : 298.22 mJy at 965,858
Iteration 377, scale 9 px : -258.22 mJy at 962,857
Iteration 442, scale 0 px : 256.93 mJy at 965,858
Iteration 476, scale 36 px : 240.94 mJy at 906,1002
Iteration 484, scale 18 px : -267.15 mJy at 965,856
Iteration 487, scale 18 px : -212.61 mJy at 965,856
Iteration 510, scale 0 px : 226.15 mJy at 964,856
Iteration 554, scale 9 px : 204.19 mJy at 837,1039
Iteration 615, scale 0 px : 190.36 mJy at 965,856
Iteration 683, scale 36 px : 184.85 mJy at 972,869
Iteration 691, scale 18 px : 156.06 mJy at 1000,1312
Iteration 749, scale 9 px : 158.05 mJy at 1145,759
Iteration 851, scale 0 px : 159.79 mJy at 966,856
Iteration 926, scale 0 px : 127.69 mJy at 1637,1248
Iteration 1099, scale 73 px : 189.74 mJy at 961,857
Iteration 1103, scale 73 px : 149.87 mJy at 961,857
Iteration 1107, scale 36 px : -156.53 mJy at 930,1359
Iteration 1113, scale 18 px : -131.71 mJy at 990,854
Iteration 1136, scale 9 px : -130.79 mJy at 964,856
Iteration 1185, scale 73 px : 119.92 mJy at 963,855
Iteration 1189, scale 5 px : -100.07 mJy at 962,856
Iteration 1342, scale 73 px : 121.28 mJy at 964,855
Iteration 1346, scale 0 px : 116.86 mJy at 990,853
Iteration 1375, scale 0 px : -95.93 mJy at 963,854
Iteration 1541, scale 18 px : -137.91 mJy at 930,1359
Iteration 1548, scale 36 px : 120.81 mJy at 969,860
Iteration 1559, scale 18 px : -102.83 mJy at 991,855
Iteration 1590, scale 9 px : -100.25 mJy at 965,856
Iteration 1658, scale 36 px : -90.48 mJy at 1099,676
Iteration 1672, scale 18 px : -78.84 mJy at 930,1359
Iteration 1758, scale 0 px : 98.2 mJy at 930,1360
Iteration 1773, scale 0 px : 79.55 mJy at 964,856
Iteration 2008, scale 9 px : 87.73 mJy at 971,851
Iteration 2048, scale 9 px : -76.94 mJy at 965,856
Iteration 2203, scale 36 px : -67.59 mJy at 1098,676
Iteration 2226, scale 0 px : 72.32 mJy at 991,853
Iteration 2345, scale 18 px : -83.65 mJy at 930,1359
Iteration 2352, scale 18 px : -66.6 mJy at 930,1359
Subminor loop is near minor loop threshold. Initiating countdown.
(12) Iteration 2387, scale 9 px : -64.77 mJy at 719,1666
(11) Iteration 2558, scale 9 px : 59.54 mJy at 722,1660
(10) Iteration 2560, scale 0 px : 65.79 mJy at 990,853
(9) Iteration 2646, scale 0 px : -55.88 mJy at 964,855
(8) Iteration 2649, scale 18 px : -66.88 mJy at 930,1359
(7) Iteration 2659, scale 9 px : 57.12 mJy at 959,862
(6) Iteration 2665, scale 0 px : 57.24 mJy at 990,855
(5) Iteration 2670, scale 18 px : -57.15 mJy at 1100,676
(4) Iteration 2673, scale 5 px : -55.22 mJy at 987,852
(3) Iteration 2675, scale 0 px : 55.93 mJy at 990,853
(2) Iteration 2678, scale 18 px : -56.23 mJy at 1100,676
(1) Iteration 2680, scale 0 px : 54.24 mJy at 1099,677
(0) Iteration 2681, scale 18 px : -55.36 mJy at 1100,676
Minor loop finished, continuing cleaning after inversion/prediction round.
 == Converting model image to visibilities ==
Opening reordered part 0 spw 0 for woden_LBA-NCP_nvss_band01.ms
Fourier transforms for pass 0... Predicting...
 == Constructing image ==
Opening reordered part 0 spw 0 for woden_LBA-NCP_nvss_band01.ms
Gridding pass 0... Fourier transforms...
 == Deconvolving (4) ==
Estimated standard deviation of background noise: 2.89 mJy
Scale info:
- Scale 0, bias factor=1, psfpeak=1, gain=0.1, kernel peak=0.19888
- Scale 5, bias factor=1.7, psfpeak=0.395592, gain=0.252786, kernel peak=0.0815272
- Scale 9, bias factor=2.8, psfpeak=0.259315, gain=0.385631, kernel peak=0.0330771
- Scale 18, bias factor=4.6, psfpeak=0.132431, gain=0.755112, kernel peak=0.00966831
- Scale 36, bias factor=7.7, psfpeak=0.0687231, gain=1.45512, kernel peak=0.00287831
- Scale 73, bias factor=12.9, psfpeak=0.0234852, gain=4.258, kernel peak=0.000791489
- Scale 146, bias factor=21.4, psfpeak=0.00621147, gain=16.0992, kernel peak=0.000207982
- Scale 291, bias factor=35.7, psfpeak=0.00154831, gain=64.5865, kernel peak=5.26599e-05
- Scale 583, bias factor=59.5, psfpeak=0.000351752, gain=284.291, kernel peak=1.32493e-05
RMS per scale: {0: 3.17 mJy, 5: 2.17 mJy, 9: 1.72 mJy, 18: 1.21 mJy, 36: 790.8 µJy, 73: 437.97 µJy, 146: 237.28 µJy, 291: 123.16 µJy, 583: 66.32 µJy}
Starting multi-scale cleaning. Start peak=60.21 mJy, major iteration threshold=9.03 mJy
Iteration 2864, scale 0 px : 59.15 mJy at 1238,468
Iteration 3050, scale 36 px : 56.85 mJy at 963,852
Iteration 3067, scale 18 px : -55.8 mJy at 930,1360
Iteration 3088, scale 0 px : -50.33 mJy at 964,855
Iteration 3472, scale 9 px : -56.38 mJy at 1149,756
Iteration 3532, scale 36 px : -49.64 mJy at 1305,1187
Iteration 3554, scale 9 px : -48.94 mJy at 965,856
Iteration 3760, scale 18 px : -44.9 mJy at 930,1360
Iteration 3800, scale 0 px : 45.08 mJy at 931,1360
Iteration 4061, scale 9 px : -44.08 mJy at 1141,762
Iteration 4271, scale 36 px : -41.4 mJy at 1305,1188
Iteration 4287, scale 0 px : 38 mJy at 1144,757
Iteration 4848, scale 18 px : -43.74 mJy at 1305,1187
Iteration 4871, scale 9 px : -40.34 mJy at 1149,756
Iteration 4993, scale 9 px : -32.92 mJy at 965,856
Iteration 5633, scale 0 px : 34.32 mJy at 966,856
Iteration 5887, scale 18 px : -38.51 mJy at 1100,677
Iteration 5904, scale 36 px : -33.87 mJy at 1325,1026
Iteration 5945, scale 5 px : -31.17 mJy at 836,1043
Iteration 6055, scale 0 px : 30.63 mJy at 989,854
Iteration 6468, scale 18 px : -35.87 mJy at 1125,815
Iteration 6485, scale 18 px : -28.57 mJy at 629,921
Iteration 6670, scale 36 px : -28.85 mJy at 1145,759
Iteration 6696, scale 9 px : -28.16 mJy at 943,850
Iteration 7020, scale 0 px : 27.55 mJy at 1307,1188
Iteration 7539, scale 5 px : -24.68 mJy at 1148,757
Iteration 7754, scale 18 px : -29.15 mJy at 1126,815
Iteration 7768, scale 36 px : -25.48 mJy at 868,835
Iteration 7853, scale 9 px : -25.66 mJy at 722,1654
Iteration 8054, scale 0 px : 24.19 mJy at 1144,758
Iteration 8699, scale 18 px : -28.85 mJy at 1126,815
Iteration 8714, scale 18 px : 23.06 mJy at 1032,1102
Iteration 8982, scale 9 px : -22.55 mJy at 1149,757
Iteration 9430, scale 0 px : 22.1 mJy at 1125,815
Iteration 10008, scale 73 px : -20.43 mJy at 1147,761
Iteration 10022, scale 18 px : -24.78 mJy at 1126,814
Iteration 10042, scale 36 px : -21.63 mJy at 867,835
Iteration 10119, scale 0 px : -20.39 mJy at 965,857
Iteration 10751, scale 18 px : -23.08 mJy at 1125,814
Iteration 10779, scale 9 px : -20.1 mJy at 722,1654
Iteration 11379, scale 18 px : 18.76 mJy at 1075,961
Iteration 11791, scale 0 px : 18.76 mJy at 1125,815
Iteration 12476, scale 36 px : -19.36 mJy at 868,835
Iteration 12542, scale 9 px : -18 mJy at 717,1665
Iteration 13370, scale 5 px : 16.61 mJy at 779,260
Iteration 13719, scale 18 px : -16.77 mJy at 1125,815
Iteration 14241, scale 0 px : 16.66 mJy at 709,1024
Iteration 15155, scale 36 px : -16.84 mJy at 868,836
Iteration 15300, scale 9 px : -16.08 mJy at 981,821
Iteration 16515, scale 9 px : -14.71 mJy at 929,1362
Iteration 18685, scale 18 px : 15.03 mJy at 1074,961
Iteration 19055, scale 0 px : 14.76 mJy at 710,1025
Iteration 20000, scale 73 px : -14.71 mJy at 1598,765
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.
Writing model image...
 == Converting model image to visibilities ==
Opening reordered part 0 spw 0 for woden_LBA-NCP_nvss_band01.ms
Fourier transforms for pass 0... Predicting...
 == Constructing image ==
Opening reordered part 0 spw 0 for woden_LBA-NCP_nvss_band01.ms
Gridding pass 0... Fourier transforms...
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: 8210 (33.63 Jy)
- Scale 5 px, nr of components cleaned: 829 (2.28 Jy)
- Scale 9 px, nr of components cleaned: 8038 (29.67 Jy)
- Scale 18 px, nr of components cleaned: 2309 (2.92 Jy)
- Scale 36 px, nr of components cleaned: 584 (-313.67 mJy)
- Scale 73 px, nr of components cleaned: 30 (717.29 mJy)
- Scale 146 px, nr of components cleaned: 0 (0 Jy)
- Scale 291 px, nr of components cleaned: 0 (0 Jy)
- Scale 583 px, nr of components cleaned: 0 (0 Jy)
Total: 20000 components (68.91 Jy)
Inversion: 00:01:57.988430, prediction: 00:01:08.047451, deconvolution: 00:01:27.787669
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=600
cent_x = 1024
cent_y = 1024

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.