Skip to content

Instantly share code, notes, and snippets.

@d3v-null
Last active October 1, 2026 09:12
Show Gist options
  • Select an option

  • Save d3v-null/63eb2255d8823783866da9e77debfce0 to your computer and use it in GitHub Desktop.

Select an option

Save d3v-null/63eb2255d8823783866da9e77debfce0 to your computer and use it in GitHub Desktop.
Rapthor on MWA data at Pawsey Setonix with the rapthor-mwa Singularity image (ghcr.io/d3v-null/rapthor-mwa), incl. a guide to choosing rapthor settings for MWA

Rapthor on MWA data at Pawsey Setonix with the rapthor-mwa container

ghcr.io/d3v-null/rapthor-mwa is a lean (no Jupyter/conda) Spack build of the LOFAR Rapthor direction-dependent calibration and imaging pipeline with MWA support compiled in:

component version MWA notes
rapthor 2.1.20260630 (py-rapthor +mwa) patched: DP3/WSClean pointed at the MWA beam, facet-beam options dropped, flagged-tile and solar-model robustness fixes (list below)
DP3 6.6.20260819 applybeam.coefficients_path / solve.coefficients_path = /opt/mwa_full_embedded_element_pattern.h5
EveryBeam 0.8.3 (+python) MWA full-embedded-element tile beam
WSClean 3.6.20260109 -mwa-path /opt, -apply-primary-beam
Python 3.12 (lsmtool, losoto, PyBDSF, python-casacore)

The MWA beam file is shipped at /opt/mwa_full_embedded_element_pattern.h5; you do not need to download it or set MWA_BEAM_FILE. The image is CPU-only.

Tag used here: docker://ghcr.io/d3v-null/rapthor-mwa:sha-48f973a-pass (-pretest tags are pushed as soon as the build finishes; the matching sha-48f973a-pass tag exists too, published after the image's own smoke tests passed, and numpy2-pass tracks the branch). Multi-arch: linux/amd64 and linux/arm64. Source and Dockerfile: https://github.com/d3v-null/Karabo-Pipeline (rapthor-lean.Dockerfile, RAPTHOR_MWA=1; patch in spack-overlay/packages/py-rapthor/mwa-beam-support-20260630.patch).

1. Pull the image on Setonix

Singularity's cache defaults to $HOME, whose quota on Setonix is tiny; put it on scratch. Use a plain (-nohost) Singularity module for CPU work.

ssh setonix
module load singularity/4.1.0-nohost
export SINGULARITY_CACHEDIR=$MYSCRATCH/.singularity
export SINGULARITY_TMPDIR=$SINGULARITY_CACHEDIR/tmp
mkdir -p $SINGULARITY_TMPDIR $MYSCRATCH/containers

SIF=$MYSCRATCH/containers/rapthor-mwa_sha-48f973a-pass.sif
singularity pull $SIF docker://ghcr.io/d3v-null/rapthor-mwa:sha-48f973a-pass   # 590 MB SIF; ~11 min on a login node

# smoke test
export OPENBLAS_NUM_THREADS=1      # WSClean refuses to start otherwise
singularity exec $SIF bash -c 'rapthor --version; DP3 --version | head -1; wsclean --version | sed -n 2p; ls -l /opt/mwa_full_embedded_element_pattern.h5'

Expected (verified on Setonix, 2026-10-01):

rapthor v2.2.dev117+g01a81e11.d20261001
DP3 6.6.0
WSClean version 3.6 (2025-02-07)
-rw-r--r-- 1 root root 133209592 ... /opt/mwa_full_embedded_element_pattern.h5

A quick end-to-end beam check uses EveryBeam's own MWA test MS:

cd $MYSCRATCH && mkdir -p mwa-smoke && cd mwa-smoke
wget -q https://support.astron.nl/software/ci_data/EveryBeam/MWA-single-timeslot.tar.bz2
mkdir -p MWA_MOCK.ms && tar -xf MWA-single-timeslot.tar.bz2 -C MWA_MOCK.ms --strip-components=1
singularity exec -B $MYSCRATCH $SIF python3 -c "
import everybeam; t = everybeam.load_telescope('MWA_MOCK.ms', coeff_path='/opt/mwa_full_embedded_element_pattern.h5')
print(type(t).__name__, t.nr_stations)"                                   # -> MWA 128
singularity exec -B $MYSCRATCH $SIF DP3 msin=MWA_MOCK.ms msout=. msout.datacolumn=CORRECTED_DATA \
    steps=[applybeam] applybeam.coefficients_path=/opt/mwa_full_embedded_element_pattern.h5

Always run with these environment variables (put them in your job script):

export OPENBLAS_NUM_THREADS=1          # mandatory for WSClean and DP3
export TOIL_HISTORY=False TOIL_JOB_HISTORY=False   # rapthor's CWL runner bookkeeping does not need a writable DB
export HOME=$WORKDIR/.home             # rapthor/astropy/matplotlib write under $HOME; keep it on scratch
export XDG_CACHE_HOME=$HOME/.cache

Bind your scratch tree explicitly (-B $MYSCRATCH); the image's own /opt/xdg already holds an offline astropy IERS table, so nothing needs internet access at run time.

2. Prepare MWA data for rapthor

Rapthor expects a MeasurementSet with a single field, a DATA column, and the MWA-specific subtables MWA_TILE_POINTING (tile delays, used by EveryBeam) and MWA_SUBBAND. Birli, Cotter and ASVO conversion jobs all write these.

One MS per coarse channel. MWA observations are often "picket fence" (non-contiguous coarse channels); Birli writes one MS per contiguous range and rapthor is happiest with one run per 1.28 MHz coarse channel (frequency-dependent solution intervals and image cell sizes then fall out naturally; see section 6). For a contiguous 30.72 MHz band a single MS/run is fine.

Calibrator transfer first. Rapthor has no "transfer a calibrator's bandpass" step: it self-calibrates against a sky model. On MWA data you therefore start from visibilities that already carry a direction-independent calibration from a calibrator observation, e.g. with hyperdrive (available as a module on Setonix):

module use /software/projects/pawsey1094/setonix/2024.05/modules/zen3/gcc/12.2.0
module load hyperdrive/default birli/default
export MWA_BEAM_FILE=/software/projects/$PAWSEY_PROJECT/mwa_full_embedded_element_pattern.h5   # or any copy

# per coarse channel: calibrate the calibrator, apply to the target
hyperdrive srclist-by-beam -m cal.metafits -n 300 GGSM_updated.fits -- cal_top300.yaml
hyperdrive di-calibrate --data cal.metafits cal_ch113.ms --source-list cal_top300.yaml \
    --uvw-min 20l --max-iterations 300 --outputs soln_ch113.fits
hyperdrive solutions-apply --data target.metafits target_ch113.ms --solutions soln_ch113.fits \
    --tile-flags HexE21 HexE22 ... --outputs target_ch113_cal.ms

Check the solutions (hyperdrive solutions-plot) and flag tiles with random phases or wild amplitudes before applying; rapthor cannot recover a tile whose calibrator solution is garbage. Tiles that are flagged in the calibrator but not the target are lost at this step.

Where to write MeasurementSets. casacore writes MS one cell at a time through Marlu/rubbl (Birli, hyperdrive); on networked filesystems this is very slow (hundreds of GB of read-modify-write per MS). On Setonix write to /scratch (Lustre, striped) rather than /software or $HOME, average as far as the science allows (birli --avg-time-res, --avg-freq-res), and prefer ASVO-converted MS where possible. A bulk-write fix for Marlu/Birli is in progress (pkgw/rubbl#601, MWATelescope/Marlu#38).

Flagging. Do not run AOFlagger on solar data (bursts look like RFI). For other fields the Birli default MWA strategy is fine. Pre-flag known bad tiles in the MS FLAG column; rapthor's DDECal leaves fully flagged tiles out of the h5parm and the patched image applies missingantennabehavior=flag for them.

3. Sky model

Rapthor reads makesourcedb-format sky models (true-sky fluxes; patches define calibration directions). Three ways to get one for an MWA field:

  1. Let rapthor make it (generate_initial_skymodel = True): rapthor images the field and extracts a model with PyBDSF. Works for ordinary fields once the calibrator transfer is applied.
  2. From GGSM / a hyperdrive list: hyperdrive srclist-by-beam to the top-N sources, then convert to makesourcedb (hyperdrive srclist-convert to a text format and a short script writing Name, Type, Patch, Ra, Dec, I, SpectralIndex, ReferenceFrequency). Group into patches with lsmtool (group('threshold') or group('voronoi')).
  3. The Sun: a single Gaussian at the phase centre is enough to start (the same starting point P-AIRCARS uses):
FORMAT = Name, Type, Patch, Ra, Dec, I, SpectralIndex, LogarithmicSI, ReferenceFrequency='144560000.0', MajorAxis, MinorAxis, Orientation
 , , Sun, 341.106457, -7.989477
QuietSun, GAUSSIAN, Sun, 341.106457, -7.989477, 42000.0, [0.0], false, 144560000.0, 2250.0, 2250.0, 0.0

with FWHM = 32' + 2.2'·(ν/GHz)^-0.6 and the quiet-Sun flux in Jy (≈1.4×10⁴ at 74 MHz, 4×10⁴ at 145 MHz, 1×10⁵ at 240 MHz). Pass the same file as input_skymodel and apparent_skymodel; DP3 applies the MWA beam itself.

4. Parset and strategy (solar, direction-independent)

These are the settings used for the MWA solar self-calibration study (49 live tiles, Sun at the phase centre, one coarse channel per run); see section 6 for how to adapt them.

rapthor.parset:

[global]
dir_working = /work/rapthor
input_ms = /work/input.ms
data_colname = DATA
generate_initial_skymodel = False
download_initial_skymodel = False
input_skymodel = /work/sun.skymodel
apparent_skymodel = /work/sun.skymodel
regroup_input_skymodel = False
strategy = /work/strategy.py
selfcal_data_fraction = 1.0
final_data_fraction = 1.0
dde_mode = faceting

[calibration]
use_included_skymodels = False
use_image_based_predict = False
dd_interval_factor = 1
dd_smoothness_factor = 1
llssolver = qr
maxiter = 150
propagatesolutions = True
solveralgorithm = directioniterative
onebeamperpatch = False
stepsize = 0.02
stepsigma = 2.0
tolerance = 5e-3
fast_freqstep_hz = 160e3
fast_smoothnessconstraint = 1e6
fast_datause = single
medium_freqstep_hz = 160e3
medium_smoothnessconstraint = 1e6
medium_datause = single
slow_freqstep_hz = 320e3
slow_smoothnessconstraint = 1.5e6
slow_datause = dual
fulljones_freqstep_hz = 1.28e6
fulljones_smoothnessconstraint = 0.0
correct_time_frequency_smearing = False

[imaging]
cellsize_arcsec = 20.0        ; ~ (lambda / B_max) / 4 ; see section 6
robust = 0.0
min_uv_lambda = 10
max_uv_lambda = 1e7
mgain = 0.8
do_multiscale_clean = True
dde_method = single           ; DI: one solution applied per sector
filter_skymodel = True
source_finder = bdsf
save_visibilities = True
save_supplementary_images = True
compress_selfcal_images = False
compress_final_images = False
idg_mode = cpu
mem_gb = 100
apply_diagonal_solutions = True
make_quv_images = False
use_mpi = False
reweight = False
grid_width_ra_deg = 4.0
grid_width_dec_deg = 4.0
grid_nsectors_ra = 0
skip_final_major_iteration = False

[cluster]
batch_system = single_machine
max_nodes = 1
cpus_per_task = 0
max_cores = 32
max_threads = 32
deconvolution_threads = 16
mem_per_node_gb = 100
cwl_runner = toil
debug_workflow = False
keep_temporary_files = False

strategy.py (single direction, amplitude self-cal from cycle 2, rapthor's own convergence check on so it stops if the image degrades):

strategy_steps = []
for i in range(3):
    s = {}
    s["do_calibrate"] = True
    s["do_slowgain_solve"] = i >= 1
    s["do_fulljones_solve"] = False
    s["peel_outliers"] = False
    s["peel_bright_sources"] = False
    s["max_normalization_delta"] = 0.3
    s["scale_normalization_delta"] = True
    s["solve_min_uv_lambda"] = 10
    s["fast_timestep_sec"] = 2.0
    s["medium_timestep_sec"] = 8.0
    s["slow_timestep_sec"] = 32.0
    s["fulljones_timestep_sec"] = 176.0
    s["do_normalize"] = False
    s["do_image"] = True
    s["auto_mask"] = 4.0
    s["auto_mask_nmiter"] = 2
    s["max_nmiter"] = 8
    s["channel_width_hz"] = 1.28e6
    s["threshisl"] = 4.0
    s["threshpix"] = 5.0
    s["regroup_model"] = True
    s["max_directions"] = 1
    s["max_distance"] = None
    s["target_flux"] = 100.0
    s["do_check"] = i > 0
    s["convergence_ratio"] = 0.95
    s["divergence_ratio"] = 1.1
    s["failure_ratio"] = 10.0
    strategy_steps.append(s)
strategy_steps.append(dict(strategy_steps[-1]))   # final pass = last cycle

5. Slurm job script

One coarse channel per task; several tasks share a 256-CPU / 245 GB work node. --cpus-per-task 32 with max_cores = 32 in the parset gives 8 channels per node.

#!/bin/bash
#SBATCH --job-name=rapthor-mwa
#SBATCH --account=${PAWSEY_PROJECT}
#SBATCH --partition=work
#SBATCH --array=0-23                 # 24 coarse channels
#SBATCH --cpus-per-task=32
#SBATCH --mem=100G
#SBATCH --time=06:00:00
#SBATCH --output=%x-%A_%a.out

set -euo pipefail
module load singularity/4.1.0-nohost
CHANS=(58 61 65 69 73 77 81 86 91 96 101 107 113 120 127 134 142 150 158 167 177 187 210 226)
CH=${CHANS[$SLURM_ARRAY_TASK_ID]}; CHS=$(printf '%03d' $CH)

SIF=$MYSCRATCH/containers/rapthor-mwa_sha-48f973a-pass.sif
MS=$MYSCRATCH/1424757616/hypcal/hyp_cal_1424757616_ch${CHS}.ms      # calibrator-transferred MS
RUN=$MYSCRATCH/rapthor/ch${CHS}; mkdir -p $RUN/.home

# per-channel parset: cell size ~ (lambda/B_max)/4 for a 5.3 km array
FMHZ=$(singularity exec -B $MYSCRATCH $SIF python3 -c "
from casacore.tables import table; print('%.3f' % (table('$MS/SPECTRAL_WINDOW', ack=False).getcol('CHAN_FREQ')[0].mean()/1e6))")
CELL=$(python3 -c "print('%.1f' % (2900.0/$FMHZ))")
sed -e "s/^cellsize_arcsec = .*/cellsize_arcsec = $CELL/" rapthor.parset.template > $RUN/rapthor.parset
cp strategy.py $RUN/strategy.py
singularity exec -B $MYSCRATCH $SIF python3 make_sun_skymodel.py $MS $RUN/sun.skymodel   # or your own model

export OPENBLAS_NUM_THREADS=1 TOIL_HISTORY=False TOIL_JOB_HISTORY=False
export HOME=/work/.home XDG_CACHE_HOME=/work/.home/.cache
# rapthor resumes completed operations; clear stale Toil job stores from an interrupted run
find $RUN/rapthor/pipelines -maxdepth 2 -name jobstore -type d -exec rm -rf {} + 2>/dev/null || true

singularity exec \
  -B $MYSCRATCH \
  -B $RUN:/work \
  -B $MS:/work/input.ms:ro \
  $SIF rapthor /work/rapthor.parset

Outputs land in $RUN/rapthor/: solutions/calibrate_N/field-solutions.h5 (losoto h5parm; plotrapthor field-solutions.h5 phase), images/image_N/ (field-MFS-image.fits, -pb = MWA primary-beam corrected), skymodels/, plots/, logs/<operation>/ (CWL step logs; failed_* on errors). With save_visibilities = True, visibilities/image_N/ holds the imaging MS, but note it is averaged to the smearing limit; to re-image at full resolution apply the h5parm to the input MS yourself:

singularity exec -B $MYSCRATCH $SIF DP3 msin=$MS msout=$RUN/cal.ms steps=[applycal] \
  applycal.parmdb=$RUN/rapthor/solutions/calibrate_1/field-solutions.h5 applycal.solset=sol000 \
  applycal.correction=phase000 applycal.missingantennabehavior=flag

(for cycles with amplitudes add applycal.steps=[ph,amp] applycal.ph.correction=phase000 applycal.amp.correction=amplitude000, and applycal.direction=[Sun] if the h5parm has several directions).

6. Choosing rapthor settings for MWA

Rapthor's defaults are tuned for LOFAR HBA (48 MHz of contiguous band, 1.5" cells, thousands of sources, 8 h tracks). MWA data differ in almost every parameter that matters, so start from the table below rather than the defaults.

Frequency structure

situation *_freqstep_hz *_smoothnessconstraint runs
picket fence, one 1.28 MHz MS per coarse channel (8×160 kHz or 32×40 kHz) 160e3 (fast/medium), 320e3 (slow) 1e6 to 1.5e6 (never larger than the band) one rapthor run per channel
contiguous 30.72 MHz 1e6 (fast), 2e6 (slow) 3e6 to 6e6 one run; channel_width_hz 4e6 to 8e6 for cube/model images

The smoothness constraint is a kernel in Hz; set it to a fraction of the bandwidth you actually have or every channel collapses to one solution.

Solution intervals (SNR)

Per-baseline noise is σ ≈ SEFD / √(2 Δν Δt). With an MWA tile SEFD of roughly 20 to 50 kJy (worse at 70 to 90 MHz), 160 kHz and 2 s give σ ≈ 25 to 60 Jy per visibility. Scalar-phase solves need the model flux in the solve interval to exceed a few σ per baseline:

target model flux fast_timestep_sec slow_timestep_sec
Sun (10⁴ to 10⁵ Jy) huge 1 to 4 s (ionosphere and bursts change on seconds) 16 to 32 s
strong field (A-team, >100 Jy patch) 100 Jy+ 8 to 20 s with full 1.28 MHz 60 to 120 s
EoR-style field (brightest patch 10 to 50 Jy) tens of Jy 20 to 120 s, full band, fast_freqstep_hz = whole MS 120 to 600 s, or skip

medium_timestep_sec sits between the two (rapthor solves fast→medium→slow). Keep fulljones_* off unless you have a polarised calibrator model.

Directions

  • Sun or a single dominant source: max_directions = 1, regroup_model = True, target_flux well below the source flux (100 Jy), dde_method = single. Without regrouping, cycle 2 splits the model into the main source plus a leftover "Patch" direction; with dde_method = single rapthor may then pre-apply the wrong direction and image almost no data.
  • Wide fields: MWA's ~25° (150 MHz) field of view contains hundreds of sources; group with target_flux 5 to 20 Jy into 4 to 12 directions, regroup_model = True, dde_method = full (facets). Facet imaging works with the MWA patch, but WSClean's per-facet beam is not available for MWA (-apply-facet-beam is removed), so the primary beam is applied once per sector; keep sectors small enough (grid_nsectors_ra > 0 or explicit sector_* lists) that the beam does not vary much across a sector.
  • peel_outliers / peel_bright_sources: useful for A-team sources outside the imaged area (set grid_width_* tighter than the beam), but both trigger extra predict steps; leave off for solar runs.

Imaging

parameter guidance
cellsize_arcsec ≈ (λ/B_max)/4: for the compact/extended 5.3 km configuration that is 2900″/ν(MHz) (39″ at 74 MHz, 20″ at 145 MHz, 10″ at 290 MHz); for the full 5.3 km with long baselines flagged check B_max from your ANTENNA table
grid_width_*_deg Sun: 3 to 4° (just the disc and sidelobes). Fields: the beam FWHM (≈ 25° × 150 MHz/ν) or less; MWA images at 20″ over 25° are 4500 px, so budget mem_gb and threads accordingly or use several sectors
robust 0 (Sun; keeps short-baseline disc flux), −0.5 to −1 (fields, point sources)
min_uv_lambda 10 (solar; keep the disc), 30 to 50 (fields; suppress Galactic emission)
do_multiscale_clean True for the Sun / extended emission, False for point-source fields
threshisl/threshpix 4/5 is fine; on a bright resolved Sun PyBDSF may still find nothing at >120 MHz, in which case the patched filter_skymodel retries with global rms boxes and falls back to a 5σ threshold mask and the unfiltered WSClean model
max_nmiter 8 for the Sun; 6 to 10 for fields
skip_final_major_iteration False for short observations where imaging is cheap

Cycles and safety nets

  • do_slowgain_solve from cycle 2, never cycle 1.
  • do_check = True from cycle 2 (convergence_ratio 0.95, divergence_ratio 1.1): on MWA solar data amplitude self-cal reduced the dynamic range in both this pipeline and P-AIRCARS; rapthor's check stops the run at the best cycle instead of continuing. Phase-only (cycle 1) is often already the best image.
  • do_normalize = False: there is no external flux scale to normalise to; be aware that unconstrained amplitude solutions can drift the absolute flux scale (factors of 2 to 4 were seen). Check the peak flux cycle to cycle.
  • selfcal_data_fraction = 1.0: MWA scans are short (2 to 5 min); there is no benefit in subsampling.

Beam

The patch sets usebeammodel=True, beammode=array_factor, beam_interval=120 and the coefficient path for every DP3 step; the MWA tile beam is a slowly varying analytic beamformer response so this is sufficient. Images are PB-corrected with WSClean's -apply-primary-beam -mwa-path /opt. Nothing to configure, but the MS must retain MWA_TILE_POINTING (EveryBeam reads the delays from it) and the delays must match the observation.

Resources on Setonix

Per 1.28 MHz channel of a 3 min, 128-tile, 1 s / 160 kHz observation a three-cycle run takes 20 to 30 minutes on 10 to 12 cores and well under 100 GB. work nodes have 256 logical CPUs and 245 GB, so 8 concurrent channel runs per node (--cpus-per-task=32, max_cores = 32) is a sensible array layout; Toil adds a fixed overhead of a minute or two per CWL step, so very small max_cores does not help.

Known MWA-specific behaviour fixed in this image

  • DP3 applycal aborts with SolTab has no element <tile> in ant for fully flagged tiles → missingantennabehavior=flag in the prepare/predict steps.
  • process_gains.py referenced phases to station 0, which on MWA is often a flagged tile → all slow-gain phases NaN → blank images; it now uses the first unflagged station.
  • PyBDSF finds zero sources on a bright resolved Sun → no mask → crash in the mosaic step; see threshisl/threshpix above.
  • The LOFAR SEFD table used for the "expected noise" diagnostic ends at 240 MHz; it is clamped instead of raising above that.
  • WSClean facet-beam options (-apply-facet-beam, -facet-beam-update) are unsupported for MWA in EveryBeam and are removed.

7. Troubleshooting

symptom cause / fix
OpenBLAS multi-threading error at WSClean start export OPENBLAS_NUM_THREADS=1
Permission denied under /home/jovyan or ~/.astropy set HOME to a bound, writable scratch path
toil ... --restart failure on rerun delete rapthor/pipelines/*/jobstore and rerun; completed operations are skipped
Gridded visibility count: 0, blank images all data flagged after applycal: check the h5parm has solutions for your tiles (losoto/h5py), and that the direction being applied is the one with solutions
insufficient flux density in the model to meet the target flux density target_flux larger than any patch; lower it
Antenna type not recognized (only LBA and HBA data are supported) harmless warning for MWA
slow MS writes on /scratch see section 2; average, or use ASVO-produced MS

Appendix: make_sun_skymodel.py

Writes the quiet-Sun Gaussian starting model used above from an MS (phase centre and mean frequency are read from the MS; flux from the P-AIRCARS quiet-Sun polynomial, size from White 2016).

#!/usr/bin/env python3
import sys
import numpy as np
from casacore.tables import table

def quiet_sun_flux_jy(f_mhz):
    p = np.poly1d([-1.93715165e-06, 7.84627718e-04, -3.15744433e-02, 2.32834400e-01])  # SFU
    return float(p(f_mhz) * 1e4)

def sun_dia_arcmin(f_mhz):
    return 32 + 2.2 * (f_mhz / 1e3) ** (-0.6)

ms, out = sys.argv[1], sys.argv[2]
freqs = table(ms + "/SPECTRAL_WINDOW", ack=False).getcol("CHAN_FREQ")[0]
ra, dec = np.degrees(table(ms + "/FIELD", ack=False).getcol("PHASE_DIR")[0][0])
ra %= 360.0
f_mhz = freqs.mean() / 1e6
fwhm = sun_dia_arcmin(f_mhz) * 60.0
with open(out, "w") as fh:
    fh.write("FORMAT = Name, Type, Patch, Ra, Dec, I, SpectralIndex, LogarithmicSI, "
             f"ReferenceFrequency='{freqs.mean():.1f}', MajorAxis, MinorAxis, Orientation\n")
    fh.write(f" , , Sun, {ra:.6f}, {dec:.6f}\n")
    fh.write(f"QuietSun, GAUSSIAN, Sun, {ra:.6f}, {dec:.6f}, {quiet_sun_flux_jy(f_mhz):.1f}, [0.0], false, "
             f"{freqs.mean():.1f}, {fwhm:.1f}, {fwhm:.1f}, 0.0\n")

Built and validated on MWA observation 1424757616 (Sun, 24 picket-fence coarse channels, PicA calibrator 1424775768), September/October 2026. Image source: https://github.com/d3v-null/Karabo-Pipeline (branch numpy2).

Measurement sets converted from uvfits

Hyperdrive writes uvfits; WSClean/DP3 need a MeasurementSet. The route tested here is CASA importuvfits followed by cotter's fixmwams (adds the MWA_TILE_POINTING / MWA_SUBBAND subtables and MWA keywords that EveryBeam needs for the tile delays):

curl -O https://projects.pawsey.org.au/high0.uvfits/hyp_1224424640_ssins_30l_src8k_300it_8s_80kHz.uvfits
docker run --rm -v $PWD:/w -w /w d3vnull0/casa casa --nologger --nogui -c \
  'importuvfits(fitsfile="hyp_1224424640_ssins_30l_src8k_300it_8s_80kHz.uvfits", vis="hyp_1224424640_ssins_30l_src8k_300it_8s_80kHz.ms", antnamescheme="old")'
docker run --rm -v $PWD:/w -w /w mwatelescope/cotter fixmwams hyp_1224424640_ssins_30l_src8k_300it_8s_80kHz.ms 1224424640.metafits

Test data: obs 1224424640 (EoR0, phase centre RA 0h Dec -27), 14 x 8 s, 384 x 80 kHz (167-198 MHz), 128 tiles. Two variants were run through rapthor:

variant MAIN rows / timestep tiles in ANTENNA tiles with rows in MAIN note
asis (straight from importuvfits + fixmwams) 8256 128 128 (4 fully flagged) hyperdrive kept the flagged tiles' rows
norows (rows of 4 flagged tiles deleted) 7750 128 124 emulates a uvfits whose flagged tiles have no baselines

The second variant reproduces the report that "flagged tiles remain in ANTENNA but their rows are absent from the MAIN table". DP3 reads, applybeams and predicts both variants without complaint: missing baselines are simply missing, and the +mwa patch already runs applycal with missingantennabehavior=flag, so tiles without solutions do not make DP3 abort when the model is subtracted. The problems that did show up were elsewhere:

What broke, and what the +mwa image now does about it

  1. UVW storage manager. DP3 6.6 writes new MeasurementSets with casacore's UvwStMan (msout.uvwcompression=true by default), which re-derives per-baseline UVWs from per-antenna values and refuses UVWs that are not self-consistent to its tolerance. uvfits UVWs are not, so every rapthor step that writes an MS (prepare_imaging_data, predict_model_data, concat_ms) aborted with Inconsistent UVW value written for antenna 23 on both variants. The patch sets msout.uvwcompression=False in those steps. If you write your own DP3 parsets for these MSs, add the same key.

  2. Far-sidelobe sources in the sky model. hyperdrive srclist-by-beam ranks by apparent flux, so a model can contain e.g. Cas A at 86 deg from the phase centre. lsmtool (used by rapthor for grouping and patch positions) projects the sources with a TAN projection centred on one of them, which gives NaN pixel coordinates for such sources, a NaN patch position, and finally DP3 Error in reading position nan. Filter the model to the main lobe before giving it to rapthor (the model used here keeps 299 of 300 sources):

    import lsmtool
    s = lsmtool.load('ggsm_top300.skymodel')
    s.select(s.getDistance(0.0, -27.0) < 45.0)   # phase centre RA, Dec in deg
    s.group('single'); s.setPatchPositions(method='mid')
    s.write('ggsm_r45.skymodel', clobber=True)

    The patched rapthor now logs a warning when sources are more than 45 deg from the phase centre, keeps the input patches if lsmtool's threshold grouping fails, and skips the facet outlines in the field-overview plot instead of crashing when the facet bounding box leaves the SIN projection.

  3. DP3 FastPredict + MWA beam. The rapthor-mwa image ships DP3 6.6.20260819 built with FastPredict, which is used by default by ddecal and predict. With the MWA beam and more than one thread it races on the HDF5 coefficients file and dies at start-up (H5Ovisit1(): invalid location identifier, free(): invalid pointer). This is not uvfits-specific: a 300-source model crashes identically on a Birli-written MS, while a one-source (solar) model and the older DP3 in the sp5505 image do not. The patch sets usefastpredict=False in rapthor's ddecal and predict steps; in your own DP3 parsets use predict.usefastpredict=false (or onebeamperpatch=true, which also avoids it and is much faster for large models).

  4. beammode=array_factor is wrong for the MWA. rapthor's ddecal/predict steps use beammode=array_factor (correct for LOFAR, where the element beam is handled separately). DP3 implements that mode by taking element (0,0) of the EveryBeam response as a scalar; EveryBeam's MWA tile beam ignores the mode and returns the full Jones matrix, and the resulting scalar is ~6e-4 instead of ~0.9, so the predicted model is ~1700x too faint (full and element give the expected 0.91-0.93 of the no-beam model at 167-198 MHz; array_factor gives 0.0006). Scalar-phase solutions are insensitive to the model scale, but amplitude solutions and model subtraction are not. The same happens with the older DP3 6.6 in the sp5505 image. The patch now uses beammode=full in the MWA steps; use the same in your own DP3 parsets.

Checks run

check asis norows
python-casacore structure (rows, tiles, subtables, UVW dtype) ok ok
DP3 msin read + applybeam (MWA beam) ok ok
DP3 predict, no beam, 16 threads, 300 sources ok (34 s) ok
DP3 predict, MWA beam, 16 threads, FastPredict (default) HDF5 crash HDF5 crash (also at 4 threads)
DP3 predict, MWA beam, 16 threads, usefastpredict=false - ok (238 s)
DP3 predict, MWA beam, 16 threads, onebeamperpatch=true - ok (190 s)
DP3 predict model scale: full / element / array_factor vs no beam - 0.92 / 0.92 / 0.0006
DP3 msout new MS, default UvwStMan Inconsistent UVW Inconsistent UVW
DP3 msout new MS, msout.uvwcompression=false ok ok
rapthor one DI cycle (calibrate_1 -> image_1 -> mosaic_1), 299-source GGSM model, 16 cores ok (5 min) ok (5 min)
same, with the published rapthor-mwa:sha-48f973a-pass image and no local overrides ok -

Result of the rapthor runs

Both variants ran the same one-cycle DI strategy (scalar phase, 16 s solution interval, 1 MHz channel blocks, 3 MHz smoothness; 480 x 480 pixel image at 60", robust -0.5, min_uv_lambda = 30) with the fixes above:

asis norows
tiles in h5parm 128 (4 fully flagged: 67, 74, 88, 109) 124 (1 fully flagged: 74; 26, 67, 88, 109 absent)
fast phase, median per-tile rms over time and frequency 0.3 deg 0.3 deg
image_1 peak / edge rms (Jy/beam) 7.74 / 0.287 7.74 / 0.287
PyBDSF sources in image_1 53 52
wall time (calibrate + image + mosaic) 4.5 min 4.5 min

The data were already calibrated by hyperdrive, so near-zero residual phases are the expected answer; the two variants give the same image and the same solutions for the tiles they share. Tiles without rows simply have no entry in the h5parm (and are flagged by missingantennabehavior=flag when solutions are applied), so an MS whose flagged tiles have no baselines is handled the same way as one that keeps them as flagged rows. Plot: plots/uvfits_rapthor_check.png in the comparison work directory (image_1 and fast phases per tile for both variants).

Which image has the fixes

Items 1-4 are in the +mwa rapthor patch from Karabo-Pipeline commit 48f973a (branch numpy2, fork d3v-null), i.e. ghcr.io/d3v-null/rapthor-mwa:sha-48f973a-pass and later; sha-989f724-* predates them. That image was also run end-to-end on the as-converted MS with no local overrides (calibrate_1, image_1, mosaic_1 all completed). With the older image you can avoid item 3 by setting onebeamperpatch = True in the [calibration] section of the parset and item 2 by filtering the sky model, but items 1 and 4 need the newer image (no parset switch exists for msout.uvwcompression or beammode).

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment