WSClean

Usage

Although the documentation for WSClean is extensive, I find it difficult to parse at times. I will show some examples of how I would use WSClean, explaining some of its parameters as I go along.

Access Examples on Oscar

For the remainder of this section, I will assume that you have access to the jpober condo on Oscar. If you navigate to the group's data/shared folder, you should find a directory called Tutorials.

$ cd /oscar/data/jpober/shared/Tutorials
                    

This is where all the tutorial folders I use as demos are located. For now, cd into the wsclean folder.

[Tutorials]$ cd wsclean
                    

You should find a directory structure that looks like this:

• wsclean/
  • out/
  • satellite_sample.ms/
  • plot_wsclean_output.ipynb
  • wsclean_advanced.sh
  • wsclean_basic.sh

out is where we will send the outputs of WSClean runs. satellite_sample.ms is the sample observation I prepared for this tutorial. It is composed of 37 frequency channels and 1 time-step of this observation, available on MWA ASVO. As its name suggests, this observation features radio-frequency interference from a satellite, which we can image using WSClean! plot_wsclean_output.ipynb contains some Python that will allow you to interpret the various WSClean output files and actually plot some images. Finally, wsclean_advanced.sh and wsclean_basic.sh contain some bash code for an advanced and a basic run of WSClean, respectively.

Basic WSClean run

From inside the wsclean tutorial folder, launch the WSClean container:

[wsclean]$ apptainer shell /path/to/containers/wsclean.sif
                    

(Make sure to replace /path/to with the actual path to your containers folder.) This should launch the WSClean container directly in your wsclean folder, which you can verify via pwd:

Apptainer> pwd
/oscar/data/jpober/shared/Tutorials/wsclean
                    

An example for a very basic run of WSClean is shown in wsclean_basic.sh. If you opened this file, you would find:

#!/bin/bash

wsclean \
-name out/satellite_sample_basic \
-size 2048 2048 \
-scale 2.25amin \
satellite_sample.ms
                    

This launches WSClean and tells it to prepend all output files with out/satellite_sample_basic. From the satellite_sample.ms data, it will generate a 2048x2048 image at an angular resolution of 2.25 arcminutes.

You can launch this run as follows:

Apptainer> bash wsclean_basic.sh
                    

Once it's done, your out folder should look like this:

• out/
  • satellite_sample_basic-dirty.fits
  • satellite_sample_basic-image.fits

As you can see, WSClean creates a "dirty" and a "clean" image. Note however that in this very basic run, WSClean is not actually performing any iterative cleaning. This is enabled via the niter flag, which is included in the advanced example below.

fits files are somewhat inconvenient because you cannot simply open it and look at it. To do that, you need specific software such as DS9 or CARTA. Instead of using either of these, however, I will show you how to plot these images in Python. The plot_wsclean_output.ipynb Jupyter Notebook file contains the following lines of code:

import numpy as np
import matplotlib.pyplot as plt
from astropy.io import fits
from astropy.wcs import WCS
from astropy.coordinates import SkyCoord
import pandas as pd
import astropy.units as u
plt.style.use("seaborn-v0_8")

file = "./out/satellite_sample_basic-image.fits"

hdul = fits.open(file)
image = hdul[0].data[0,0]
wcs = WCS(header=hdul[0].header)

fig, ax = plt.subplots(figsize=(6, 3),
                 subplot_kw={'projection': wcs, 'slices': ['x', 'y', 0, 0]},
                 )
im = ax.imshow(np.abs(image),
               cmap='inferno',
               norm='log',
               vmin=1e-1,
               aspect='auto',
              )
    
plt.colorbar(im, label="Intensity (Jy/beam)")
plt.title(f"Imaging WSClean's output",y=1.02)
plt.xlabel("Right Ascension")
plt.ylabel("Declination")
ax.coords.grid(color='white', alpha=0.2, linestyle='solid')

plt.show()
                    

In short, you open the appropriate fits file using astropy, read in its World Coordinate System (WCS), and plot the image using matplotlib. The WCS projection is what lets matplotlib understand that the coordinate grid is curved.

This code should produce an image that looks like this:

You can clearly see a bright spot near the bottom-right corner — this is the satellite!

Advanced WSClean run

Let's now take a look at the wsclean_advanced.sh file:

#!/bin/bash

wsclean \
-name out/satellite_sample_advanced \
-no-update-model-required \
-parallel-gridding 12 \
-size 2048 2048 \
-scale 2.25amin \
-niter 10000 \
-mgain 0.6 \
-super-weight 4 \
-weight briggs -1.0 \
-auto-threshold 1 \
-auto-mask 5 \
-multiscale \
-save-source-list \
-multiscale-scale-bias 0.85 \
-join-channels -channels-out 2  \
-fit-spectral-pol 4 \
satellite_sample.ms
                

This is a completely over-the-top run, but let's take a look at some of the new parameters:

• parallel-gridding: "Tells WSClean to use N threads during gridding and degridding." (source)
• no-update-model-required: Does not write model directly to the input Measurement Set. (This means you cannot continue a previous run.)
• niter: "Sets the maximum number of minor cleaning iterations that are allowed to be used. Once the specified number of iterations are reached, cleaning will stop." (source)
• mgain: "Sets the major iteration gain: during every major iteration, the peak is reduced by the given factor. With an mgain of 0.8, the peak is reduced by 80%." (source)
• super-weight: "Super weighting refers to performing the counting of visibilities on larger grid cells, such that neighbouring uv-samples share their weight." (source)
• weight: Determine the weight to apply to each individual baseline when forming the image. Some options include uniform, natural, or briggs.
• auto-threshold: "With this option, WSClean will calculate the standard deviation of the residual image before the start of every major deconvolution iteration, and clean until the given auto-threshold factor times the standard deviation of the image." (source)
• auto-mask: "Auto-masking allows automated deep cleaning and solves the two problems mentioned above: Only one run of wsclean is required; It maintains scale-dependent masks, which improves multi-scale cleaning." (source)
• multiscale: "WSClean's multi-scale deconvolution selects the highest peak and subtracts it with the best fitting 'scale'" (source)
• save-source-list: Make WSClean return the txt file containing the coordinates of all sources it finds during the cleaning process.
• multiscale-scale-bias: "This parameter balances between how sensitive the algorithm is towards large scales compared to smaller scales. Lower values will clean larger scales earlier and deeper. Its default is 0.6, which means something like 'if a peak is 0.6 times larger at a 2x larger scale, select the larger scale'" (source)
• join-channels channels-out: "The multi-scale algorithm works well in combination with the parameters '-join-channels -channels-out'. This woul perform peak finding and scale selection on the integrated image, and fit the found scale to each output frequency." (source)
• fit-spectral-pol: Fits the WSClean generated model to a polynomial (useful if it known that the sources have smooth behavior).

In practice you never need THAT many WSClean flags. I'm just trying to be comprehensive. You can launch that run from inside the wsclean container via:

Apptainer> bash wsclean_advanced.sh
                    

Your out folder should now contain the following:

• out/
  • satellite_sample_advanced-0000-dirty.fits
  • satellite_sample_advanced-0000-image.fits
  • satellite_sample_advanced-0000-model.fits
  • satellite_sample_advanced-0000-psf.fits
  • satellite_sample_advanced-0000-residual.fits
  • satellite_sample_advanced-0001-dirty.fits
  • satellite_sample_advanced-0001-image.fits
  • satellite_sample_advanced-0001-model.fits
  • satellite_sample_advanced-0001-psf.fits
  • satellite_sample_advanced-0001-residual.fits
  • satellite_sample_advanced-MFS-dirty.fits
  • satellite_sample_advanced-MFS-image.fits
  • satellite_sample_advanced-MFS-model.fits
  • satellite_sample_advanced-MFS-psf.fits
  • satellite_sample_advanced-MFS-residual.fits
  • satellite_sample_advanced-sources.txt

There are a lot more files now! All the 0000 and 0001 files correspond to the two frequency bands that we told WSClean to image via -channels-out 2. The MFS files correspond to the synthesized files (combined across frequencies) that we asked for via -join-channels. Finally, the sources file contains the source list we asked for via -save-source-list.

During its cleaning procedure, WSClean uses the point-spread function (psf) to produce a model, which it then compares with the dirty data to generate a residual. The end-goal product is the clean image.

Let's take a look at the rest of the plot_wsclean_output.ipynb notebook:

file = "./out/satellite_sample_advanced-MFS-image.fits"
sources = "./out/satellite_sample_advanced-sources.txt"

# Load the component list using pandas
df = pd.read_csv(sources, delimiter=',', 
                    index_col=False,
                    on_bad_lines='skip',
                )

# Remove unparsable entries
df = df[~(df == 'd(:d(:d(.\x1c').any(axis=1)]

ra = df[" Ra"].values  # Right Ascension in hh:mm:ss.sss format
bad_dec = df[" Dec"].values  # Declination in dd.mm.ss.sss format (bad format for astropy)

dec = []
for i,c in enumerate(bad_dec):
    dec.append(c.replace(".",":",2))  # Convert declination to the correct dd:mm:ss.sss format

# Convert using SkyCoord
coord = SkyCoord(ra, dec, unit=(u.hourangle, u.deg))

# Also grab the intensity of each source
intensity = np.float32(df[" I"].values)

# Load the image fits file and grab its WCS
hdul = fits.open(file)
image = hdul[0].data[0,0]
wcs = WCS(header=hdul[0].header)

# Convert RA/Dec to pixel coordinates
pix_coords = wcs.wcs_world2pix(np.column_stack(
    (coord.ra.degree, 
        coord.dec.degree, 
        np.zeros(len(coord.ra.degree)), 
        np.zeros(len(coord.ra.degree)),
    ),
), 0)

# Extract pixel coordinates
ra_pix, dec_pix = pix_coords[:, 0], pix_coords[:, 1]

# Plot the image...
fig, ax = plt.subplots(figsize=(6, 3), subplot_kw={'projection': wcs, 'slices': ['x', 'y', 0, 0]})
im = ax.imshow(np.abs(image),
                cmap='viridis',
                norm='log',
                vmin=1e-1,
                aspect='auto',
                )

# ... with the sources on top
ax.scatter(ra_pix, 
            dec_pix, 
            s=np.abs(intensity)*4, 
            color=sns.color_palette()[2],
            alpha=0.25
            )

plt.colorbar(im, label="Intensity (Jy/beam)")
plt.title(f"Imaging WSClean's output (advanced)",y=1.01)

plt.xlabel("Right Ascension")
plt.ylabel("Declination")
ax.coords.grid(color='white', alpha=0.2, linestyle='solid')

# Save x and y limits
xlim = ax.get_xlim()
ylim = ax.get_ylim()

plt.show()
                    

Without going into too much detail, this notebook loads the image fits file and the source list, converts the source list Ra/Dec coordinates into pixel coordinates, and plots the sources on top of the image. If you run this, you should get something like this:

Immediately, we see that the image is a lot "smoother" than the one made from the basic run. The sources (plotted in red) also seem to line up pretty well with where we know the satellite to be.

This concludes my WSClean tutorial. Hope this was helpful!