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).
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.h5Always 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/.cacheBind 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.
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.msCheck 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.
Rapthor reads makesourcedb-format sky models (true-sky fluxes; patches define calibration directions). Three ways to get one for an MWA field:
- 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. - From GGSM / a hyperdrive list:
hyperdrive srclist-by-beamto the top-N sources, then convert to makesourcedb (hyperdrive srclist-convertto a text format and a short script writingName, Type, Patch, Ra, Dec, I, SpectralIndex, ReferenceFrequency). Group into patches with lsmtool (group('threshold')orgroup('voronoi')). - 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.
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 = Falsestrategy.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 cycleOne 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.parsetOutputs 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).
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.
| 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.
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.
- Sun or a single dominant source:
max_directions = 1,regroup_model = True,target_fluxwell 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; withdde_method = singlerapthor 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_flux5 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-beamis removed), so the primary beam is applied once per sector; keep sectors small enough (grid_nsectors_ra> 0 or explicitsector_*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 (setgrid_width_*tighter than the beam), but both trigger extra predict steps; leave off for solar runs.
| 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 |
do_slowgain_solvefrom cycle 2, never cycle 1.do_check = Truefrom 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.
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.
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.
- DP3
applycalaborts withSolTab has no element <tile> in antfor fully flagged tiles →missingantennabehavior=flagin the prepare/predict steps. process_gains.pyreferenced 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/threshpixabove. - 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.
| 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 |
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).
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.metafitsTest 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:
-
UVW storage manager. DP3 6.6 writes new MeasurementSets with casacore's
UvwStMan(msout.uvwcompression=trueby 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 withInconsistent UVW value written for antenna 23on both variants. The patch setsmsout.uvwcompression=Falsein those steps. If you write your own DP3 parsets for these MSs, add the same key. -
Far-sidelobe sources in the sky model. hyperdrive
srclist-by-beamranks 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 DP3Error 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.
-
DP3 FastPredict + MWA beam. The
rapthor-mwaimage ships DP3 6.6.20260819 built with FastPredict, which is used by default byddecalandpredict. 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 thesp5505image do not. The patch setsusefastpredict=Falsein rapthor'sddecalandpredictsteps; in your own DP3 parsets usepredict.usefastpredict=false(oronebeamperpatch=true, which also avoids it and is much faster for large models). -
beammode=array_factoris wrong for the MWA. rapthor'sddecal/predictsteps usebeammode=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 (fullandelementgive the expected 0.91-0.93 of the no-beam model at 167-198 MHz;array_factorgives 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 thesp5505image. The patch now usesbeammode=fullin the MWA steps; use the same in your own DP3 parsets.
| 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 | - |
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).
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).