Add qu imaging

scripts/make_pol_cubes.py

@@ -0,0 +1,142 @@
+#!/usr/bin/env python3
+"""Create Stokes Q and U cubes from WSClean channel images"""
+
+
+import numpy as np
+from astropy.io import fits
+from typing import List
+import os
+import argparse
+
+from validate_lofar_images import get_rms
+
+
+def cube_maker(q_images: List[str], u_images: List[str], nchan: int, imsize: List[int]) -> None:
+    """
+    Constructs Stokes Q and U data cubes from individual FITS channel images from WSClean,
+    filtering out noisy channels based on an RMS threshold.
+    This function reads a sequence of Stokes Q and U FITS files, extracts their
+    data and frequencies, and calculates the average RMS noise per channel. Channels
+    with an average noise greater than 5 times the median noise across all channels
+    are flagged and excluded. The filtered results are written out as new 4D FITS 
+    cubes alongside metadata text files.
+
+    Parameters
+    ----------
+    q_images : list of str
+        File paths to the Stokes Q FITS images, ordered by channel.
+    u_images : list of str
+        File paths to the Stokes U FITS images, ordered by channel.
+    nchan : int
+        The total number of frequency channels to process.
+    imsize : tuple of int or list of int
+        The spatial dimensions of the input images in pixels.
+
+    Returns
+    -------
+    None
+        Outputs are written directly to disk as FITS files 
+        ('<imagename>-polcube-Q.fits', '<imagename>-polcube-U.fits')
+        and text data files.
+
+    Notes
+    -----
+    - The output files use a naming prefix extracted from the first Stokes Q filename.
+    - Temporary logs of frequencies and noise values are generated and then 
+        overwritten with the final filtered subsets.
+    """
+
+    # Extract the name for the frequency and rms noise files
+    firstq = os.path.basename(q_images[0])
+    imagename = firstq.split('-')[0]
+
+    cube_q = np.zeros((1,nchan,imsize[1],imsize[0]))
+    cube_u = np.zeros((1,nchan,imsize[1],imsize[0]))
+
+    header_q = None
+    header_u = None
+
+    with open('' + str(imagename) + '_frequency_list.dat','w') as lfr:
+        with open('' + str(imagename) + '_avg_qunoise_list.dat','w') as avgqunoise:
+            for i in range(0,nchan):
+                hdu_q = fits.open(q_images[i])
+                hdu_u = fits.open(u_images[i])
+                data_q = hdu_q[0].data[:,:]
+                data_u = hdu_u[0].data[:,:]
+                if np.isnan(data_q).any() or np.isnan(data_u).any():
+                    print("Channel " + str(i) + "is NaN, move to the next one...")
+                    continue
+
+                if header_q is None:      # This is required to store the first frequency channel in the header of the final cube, useful for visualization
+                    header_q = hdu_q[0].header
+                    header_u = hdu_u[0].header
+
+                avgnoise = 0.5 * ( get_rms(q_images[i]) + get_rms(u_images[i]) )
+                frequ = hdu_q[0].header['CRVAL3']
+                lfr.write(str(frequ)+'\n')
+                avgqunoise.write(str(avgnoise)+'\n')
+                cube_q[0,i,:,:] = data_q
+                cube_u[0,i,:,:] = data_u
+
+
+    # Check the average QU noise for each channel. If it is more than 5 times the
+    # median noise, then it is required to exclude that channel.
+
+    with open('' + str(imagename) + '_avg_qunoise_list.dat','r') as avgqunoise:
+        avgqunoise_values = np.array([float(line.strip()) for line in avgqunoise])
+
+    with open('' + str(imagename) + '_frequency_list.dat','r') as lfr:
+        frequencies = np.array([float(line.strip()) for line in lfr])
+
+    mask = avgqunoise_values <= 5. * np.median(avgqunoise_values)
+
+    correct_frequencies = frequencies[mask]
+    correct_avgqunoise = avgqunoise_values[mask]
+    correct_cube_q = cube_q[:,mask,:,:]
+    correct_cube_u = cube_u[:,mask,:,:]
+
+    with open('' + str(imagename) + '_frequency_list.dat','w') as lfr:
+        for freq in correct_frequencies:
+            lfr.write(str(freq)+'\n')
+
+    with open('' + str(imagename) + '_avg_qunoise_list.dat','w') as avgqunoise:
+        for noise in correct_avgqunoise:
+            avgqunoise.write(str(noise)+'\n')
+
+
+    # Writing the cubes
+    hdu_cube_q = fits.PrimaryHDU(correct_cube_q,header_q)
+    hdu_cube_q.writeto('' + str(imagename) + '-polcube-Q.fits', overwrite=True)
+    print("Stokes Q cube written")
+
+    hdu_cube_u = fits.PrimaryHDU(correct_cube_u,header_u)
+    hdu_cube_u.writeto('' + str(imagename) + '-polcube-U.fits', overwrite=True)
+    print("Stokes U cube written")
+
+
+def main():
+    parser = argparse.ArgumentParser(description='Cubes maker from WSClean images for RM-synthesis with RMtools')
+    parser.add_argument('--qimages', help='List of input Stokes Q channel images', type=str, required=True)
+    parser.add_argument('--uimages', help='List of input Stokes U channel images', type=str, required=True)
+    parser.add_argument('--nchan', help='Channels out as in WSClean', default=480, type=int)
+    parser.add_argument('--imsize', help='Image size in pixels as in WSClean (width,height or single value)', default='1024,1024', type=str)
+    args = parser.parse_args()
+
+    q_images = args.qimages.split(',')
+    u_images = args.uimages.split(',')
+
+    q_images = sorted(q_images)
+    u_images = sorted(u_images)
+
+    imsize_parts = args.imsize.split(',')
+    if len(imsize_parts) == 2:
+        imsize = [int(imsize_parts[0]), int(imsize_parts[1])]
+    else:
+        single_size = int(imsize_parts[0])
+        imsize = [single_size, single_size]
+
+    cube_maker(q_images, u_images, args.nchan, imsize)
+
+
+if __name__ == '__main__':
+    main()

scripts/run_rmtools.py

@@ -0,0 +1,96 @@
+#!/usr/bin/env python3
+"""Run RM-Tools rmsynth3d locally."""
+
+import argparse
+import glob
+import os
+import shutil
+import subprocess
+import sys
+from typing import List
+from pathlib import Path
+
+
+def parse_args() -> argparse.Namespace:
+    parser = argparse.ArgumentParser(
+        description="Run RM-Tools rmsynth3d on Stokes Q/U cubes and a frequency list.",
+    )
+    parser.add_argument("--stokes-q", required=True, help="Stokes Q cube (per-channel).")
+    parser.add_argument("--stokes-u", required=True, help="Stokes U cube (per-channel).")
+    parser.add_argument("--freqs", required=True, help="Frequency list in Hz (one per channel).")
+    parser.add_argument("--l", dest="max_lam2", default="150", help="Maximum lambda-squared (-l).")
+    parser.add_argument("--d", dest="dlam2", default="0.3", help="Lambda-squared channel width (-d).")
+    parser.add_argument(
+        "--output-prefix",
+        dest="output_prefix",
+        default="",
+        help="Prefix to prepend to output files (passed to rmsynth3d -o).",
+    )
+    parser.add_argument(
+        "--extra-args",
+        default="",
+        help="Extra arguments passed to rmsynth3d (as a single string).",
+    )
+    return parser.parse_args()
+
+def stage_inputs(args: argparse.Namespace):
+    workdir = Path.cwd()
+
+    q_local = workdir / "image-polcube-Q.fits"
+    u_local = workdir / "image-polcube-U.fits"
+    freq_local = workdir / "image_frequency_list.dat"
+
+    shutil.copy2(args.stokes_q, q_local)
+    shutil.copy2(args.stokes_u, u_local)
+    shutil.copy2(args.freqs, freq_local)
+
+    return str(q_local), str(u_local), str(freq_local)
+
+def build_rmsynth_cmd(args: argparse.Namespace, qfile: str, ufile: str, freqfile: str,) -> List[str]:
+    cmd = [
+        "rmsynth3d",
+        qfile,
+        ufile,
+        freqfile,
+        "-l",
+        str(args.max_lam2),
+        "-d",
+        str(args.dlam2),
+        "-v",
+        "-R",
+    ]
+    if args.output_prefix:
+        cmd += ["-o", args.output_prefix]
+    if args.extra_args:
+        cmd += args.extra_args.split()
+    return cmd
+
+
+def move_outputs(stokes_q_path: str) -> None:
+    # RM-Tools writes outputs into the same directory as the staged input FITS.
+    input_dir = os.path.dirname(stokes_q_path)
+    moved = []
+
+    for path in glob.glob(os.path.join(input_dir, "*FDF_*.fits")):
+        dest = os.path.basename(path)
+        shutil.move(path, dest)
+        moved.append(dest)
+
+    if moved:
+        print("Moved outputs to workdir:")
+        for name in moved:
+            print(f"  {name}")
+    else:
+        print("No staged outputs found to move.")
+
+
+def main() -> int:
+    args = parse_args()
+    qfile,ufile,freqfile = stage_inputs(args)
+    rmsynth_cmd = build_rmsynth_cmd(args,qfile,ufile,freqfile,)
+    print("Running:", " ".join(rmsynth_cmd))
+    print("Working directory: ", Path.cwd())
+    return subprocess.call(rmsynth_cmd)
+
+if __name__ == "__main__":
+    sys.exit(main())

steps/concat_pol.cwl

@@ -0,0 +1,62 @@
+class: CommandLineTool
+cwlVersion: v1.2
+id: make_cubes
+doc: Concatenate Q and U images into cubes for RM-synthesis
+
+baseCommand:
+  - make_pol_cubes.py
+
+inputs:
+  - id: Q_images
+    type: File[]
+    doc: List of Stokes Q images from /steps/wsclean_pol.cwl.
+    inputBinding:
+      position: 0
+      prefix: '--qimages'
+      itemSeparator: ','
+  - id: U_images
+    type: File[]
+    doc: List of Stokes U images from /steps/wsclean_pol.cwl.
+    inputBinding:
+      position: 0
+      prefix: '--uimages'
+      itemSeparator: ','
+  - id: image_size
+    type: int[]
+    doc: Image size in pixels as [x, y], matching WSClean -size.
+    inputBinding:
+      position: 0
+      prefix: '--imsize'
+      itemSeparator: ','
+  - id: nchannels
+    type: int
+    doc: Number of channels as provided in WSClean as -channels-out.
+    default: 480
+    inputBinding:
+      position: 0
+      prefix: '--nchan'
+
+outputs:
+  - id: stokesQcube
+    type: File
+    doc: Name of the Stokes Q that will be created.
+    outputBinding:
+      glob: "*-polcube-Q.fits"
+  - id: stokesUcube
+    type: File
+    doc: Name of the Stokes U that will be created.
+    outputBinding:
+      glob: "*-polcube-U.fits"
+  - id: frequencies_list
+    type: File
+    doc: List of channels frequencies.
+    outputBinding:
+      glob: "*frequency_list.dat"
+  - id: rms_list
+    type: File
+    doc: List of average Q and U rms noise per channel.
+    outputBinding:
+      glob: "*avg_qunoise_list.dat"
+
+requirements:
+    - class: InlineJavascriptRequirement

steps/run_rmtools.cwl

@@ -0,0 +1,93 @@
+class: CommandLineTool
+cwlVersion: v1.2
+id: run_rmtools
+label: Run RM-Tools rmsynth3d
+doc: |
+  Run RM-Tools rmsynth3d on Stokes Q/U cubes and a frequency list.
+
+baseCommand: run_rmtools.py
+inputs:
+  - id: stokes_q
+    type: File
+    doc: Stokes Q cube (per-channel).
+    inputBinding:
+      position: 1
+      prefix: --stokes-q
+  - id: stokes_u
+    type: File
+    doc: Stokes U cube (per-channel).
+    inputBinding:
+      position: 1
+      prefix: --stokes-u
+  - id: freqs_hz
+    type: File
+    doc: Frequency list in Hz (one per channel).
+    inputBinding:
+      position: 1
+      prefix: --freqs
+  - id: max_lam2
+    type: float?
+    default: 150
+    doc: Maximum lambda-squared value for rmsynth3d (-l).
+    inputBinding:
+      position: 1
+      prefix: --l
+  - id: dlam2
+    type: float?
+    default: 0.3
+    doc: Lambda-squared channel width for rmsynth3d (-d).
+    inputBinding:
+      position: 1
+      prefix: --d
+  - id: output_prefix
+    type: string?
+    doc: Prefix for output products passed to rmsynth3d (-o). Defaults to Stokes Q basename if not provided.
+    inputBinding:
+      position: 1
+      prefix: --output-prefix
+  - id: extra_args
+    type: string?
+    doc: Extra arguments passed to rmsynth3d.
+    inputBinding:
+      position: 1
+      prefix: --extra-args
+outputs:
+  - id: fdf_im_dirty
+    type: File
+    outputBinding:
+      glob: '*FDF_im_dirty.fits'
+  - id: fdf_real_dirty
+    type: File
+    outputBinding:
+      glob: '*FDF_real_dirty.fits'
+  - id: fdf_tot_dirty
+    type: File
+    outputBinding:
+      glob: '*FDF_tot_dirty.fits'
+  - id: fdf_maxpi
+    type: File
+    outputBinding:
+      glob: '*FDF_maxPI.fits'
+  - id: fdf_peakrm
+    type: File
+    outputBinding:
+      glob: '*FDF_peakRM.fits'
+  - id: rmsynth_stdout
+    type: File
+    outputBinding:
+      glob: rmsynth3d_stdout.log
+  - id: rmsynth_stderr
+    type: File
+    outputBinding:
+      glob: rmsynth3d_stderr.log
+
+requirements:
+  - class: InlineJavascriptRequirement
+hints:
+  - class: DockerRequirement
+    dockerPull: vlbi-cwl
+
+
+stdout: rmsynth3d_stdout.log
+stderr: rmsynth3d_stderr.log
+