RAP-918 Add step to calculate the apparent flux of the brightest sources

docker/pipeline/Dockerfile

@@ -3,6 +3,9 @@ FROM ubuntu:22.04 AS builder
 # This Docker image builds the dependencies for the Pre-Processing Pipeline.
 # Credits: this image is based on the Docker images for LINC and Rapthor.
 
+# Suppress warning from pip when installing packages as root
+ENV PIP_ROOT_USER_ACTION=ignore
+
 # Install all build-time dependencies
 RUN export DEBIAN_FRONTEND=noninteractive && \
     apt-get update && \
@@ -40,6 +43,26 @@ RUN git config --global alias.shallow-clone "!git clone --depth 1 --recurse-subm
 ARG PORTABLE=FALSE
 ARG TARGET_CPU=haswell
 
+# Do not use `pip` from the Debian repository, but fetch it from PyPA.
+# This way, we are sure that the latest versions of `pip`, `setuptools`, and
+# `wheel` are installed in /usr/local, the only directory we're going to copy
+# over to the next build stage.
+RUN wget -q https://bootstrap.pypa.io/get-pip.py && \
+    python3 get-pip.py
+
+# Install required python packages
+RUN python3 -m pip install \
+    'numpy<2'
+
+# Install the casacore measures data. We purposely do not install these from
+# the Ubuntu repository, but download the latest version directly from the
+# ASTRON ftp site.
+# Note: The file on the ftp site is updated daily. When warnings regarding
+# leap seconds appear, ignore them or regenerate the docker image.
+RUN mkdir -p /usr/local/share/casacore/data && \
+wget -qO - https://www.astron.nl/iers/WSRT_Measures.ztar | \
+    tar -C /usr/local/share/casacore/data -xzf -
+
 # Install casacore
 RUN git shallow-clone https://github.com/casacore/casacore.git
 RUN cmake \
@@ -77,17 +100,6 @@ RUN cmake \
     -G Ninja
 RUN ninja -C /src/EveryBeam/build install
 
-# Install Dysco 
-RUN git shallow-clone https://github.com/aroffringa/dysco.git
-RUN cmake \
-    -DCMAKE_BUILD_TYPE:STRING=Release \
-    -DPORTABLE=${PORTABLE} \
-    -DTARGET_CPU=${TARGET_CPU} \
-    -H/src/dysco \
-    -B/src/dysco/build \
-    -G Ninja
-RUN ninja -C /src/dysco/build install
-
 # Install AOFlagger
 RUN git shallow-clone https://gitlab.com/aroffringa/aoflagger.git
 RUN cmake \
@@ -123,17 +135,20 @@ RUN cmake \
     -G Ninja
 RUN ninja -C /src/DP3/build install
 
-# Do not use `pip` from the Debian repository, but fetch it from PyPA.
-# This way, we are sure that the latest versions of `pip`, `setuptools`, and
-# `wheel` are installed in /usr/local, the only directory we're going to copy
-# over to the next build stage.
-RUN wget https://bootstrap.pypa.io/get-pip.py && \
-    python3 get-pip.py
+# Install python-casacore
+RUN git clone https://github.com/casacore/python-casacore.git
+RUN CASACORE_DATA=/usr/local/share/casacore/data \
+    python3 -m pip install -v /src/python-casacore
 
-# Install all Python dependencies for Rapthor
-RUN python3 -m pip install \
-    cwltool \
-    toil[cwl]
+# We need to `pip install` EveryBeam from source as well, because the C++
+# build does not produce a proper Python package.
+RUN python3 -m pip install -v /src/EveryBeam
+
+# Install current version of the Pre-Processing Pipeline.
+COPY . /usr/local/share/prep
+
+# Install all Python dependencies as root.
+RUN python3 -m pip install -r /usr/local/share/prep/requirements.txt
 
 #---------------------------------------------------------------------------
 # The image will now be rebuilt without adding the sources, in order to
@@ -153,6 +168,7 @@ RUN export DEBIAN_FRONTEND=noninteractive && \
         libblas3 \
         libboost-filesystem1.74.0 \
         libboost-program-options1.74.0 \
+        libboost-python1.74.0 \
         libcairomm-1.0-1v5 \
         libcasa-casa6 \
         libcasa-fits6 \
@@ -184,15 +200,6 @@ RUN export DEBIAN_FRONTEND=noninteractive && \
 
 RUN rm -rf /var/lib/apt/lists/*
 
-# Install the casacore measures data. We purposely do not install these from
-# the Ubuntu repository, but download the latest version directly from the
-# ASTRON ftp site.
-# Note: The file on the ftp site is updated daily. When warnings regarding
-# leap seconds appear, ignore them or regenerate the docker image.
-RUN mkdir -p /usr/local/share/casacore/data && \
-wget -qO - ftp://ftp.astron.nl/outgoing/Measures/WSRT_Measures.ztar | \
-    tar -C /usr/local/share/casacore/data -xzf -
-
 # Try to run the compiled tools to make sure they run without
 # a problem (e.g. no missing libraries).
 RUN aoflagger --version && \
@@ -203,5 +210,5 @@ RUN aoflagger --version && \
 RUN cwltool --version && \
     toil-cwl-runner --version
 
-# Install current version of the Pre-Processing Pipeline.
-COPY . /usr/local/share/prep
+# Set environment variables needed at run-time.
+ENV EVERYBEAM_DATADIR=/usr/local/share/everybeam

requirements.txt

+cwltool
+toil[cwl]
+lsmtool
+pandas
+python-casacore

scripts/find_brightest_sources.py

+#!/usr/bin/env python3
+"""
+Script to find the apparent flux of the brightest sources in a skymodel or known catalogue, such as TGSS.
+"""
+import argparse
+
+import numpy as np
+import pandas as pd
+import lsmtool
+
+from casacore import tables
+
+
+def main(obs_ms, n_sources, skymodel_file=None, catalogue=None, search_radius=None, output_csv=None, apply_beam=True):
+    """Retrieves the apparent flux of the brightest sources in the direction of
+    a Measurement Set and writes these to disk or to standard output.
+
+    Parameters
+    ----------
+    obs_ms : str
+        Name of the Observation's MeasurementSet (MS), used for calculating
+        the beam and phase centre.
+    n_sources : int
+        The maximum number of sources to obtain.
+    skymodel_file : str
+        Name of the skymodel file; takes precedent over a catalogue.
+    catalogue : str
+        Query sources from a catalogue instead of skymodel file. Can be either
+        one of the following catalogues supported by LSMTool: TGSS, WENSS, GSM,
+        or NVSS.
+    search_radius : float
+        If using a catalogue , LSMTool will query sources in a cone with this
+        radius (in degrees) around the observation's phase centre.
+    output_csv : str
+        Name of the file to store the brightest sources in. If provided, these
+        are written to CSV, else the brightest sources will be displayed in
+        the terminal.
+    apply_beam : bool
+        Setting whether to apply the beam to get the apparent flux, or not.
+    """
+    skymodel = None
+    if skymodel_file:
+        skymodel = skymodel_file
+    else:
+        skymodel = catalogue
+
+        if not search_radius:
+            raise Exception('when querying a catalogue, please also specify --radius')
+
+    pointing_deg = read_measurementset_phase_dir(obs_ms, degrees=True)
+    lsm = lsmtool.load(skymodel, beamMS=obs_ms, VOPosition=pointing_deg, VORadius=search_radius)
+
+    read_skymodel_from_file = bool(skymodel_file)
+    brightest_sources = brightest_skymodel_sources(lsm, apply_beam=apply_beam, n_sources=n_sources, sum_patches=read_skymodel_from_file)
+
+    if output_csv:
+        brightest_sources.to_csv(output_csv, index=False)
+    else:
+        print(brightest_sources.to_string(index=False))
+
+
+def brightest_skymodel_sources(skymodel, apply_beam=True, n_sources=10, sum_patches=False):
+    """
+    Retrieves the brightest sources in the skymodel and returns the names,
+    coordinates (RA/Dec), and (apparent) flux of the sources.
+
+    Parameters
+    ----------
+    lofarskymodel : LSMTool.skymodel.SkyModel
+        Skymodel from which to extract the brightest sources
+    apply_beam : bool
+        Setting whether to apply the beam to get the apparent flux, or not.
+        By default we return the apparent flux, hence, apply_beam = True.
+    n_sources : int
+        The maximum number of sources to obtain.
+    sum_patches : bool
+        Whether to sum patch components. Used when the user provides their
+        own sky model file, in that case sources might be described by mutliple
+        components combined into a single patch.
+
+    Returns
+    -------
+    sources : pd.DataFrame
+        Pandas DataFrame with information about the brightest sources in the skymodel.
+    """
+    if sum_patches:
+        # Sum over the patch components to get the total flux.
+        fluxes = skymodel.getColValues('I', units='mJy', applyBeam=apply_beam, aggregate='sum')
+        source_names = skymodel.getPatchNames()
+        source_ra, source_dec = skymodel.getPatchPositions(patchName=source_names, method='wmean', asArray=True)
+    else:
+        fluxes = skymodel.getColValues('I', units='mJy', applyBeam=apply_beam)
+        source_ra = skymodel.getColValues('Ra')
+        source_dec = skymodel.getColValues('Dec')
+        source_names = skymodel.getColValues('Name')
+
+    sorted_indices = np.argsort(-fluxes) # Sort from highest to lowest flux
+    n_brightest_source_indices = sorted_indices[:n_sources]
+
+    brightest_sources = {
+        'name': source_names[n_brightest_source_indices],
+        'ra': source_ra[n_brightest_source_indices],
+        'dec': source_dec[n_brightest_source_indices],
+        'flux': fluxes[n_brightest_source_indices]
+    }
+
+    return pd.DataFrame(data=brightest_sources)
+
+
+def read_measurementset_phase_dir(ms, degrees=False):
+    """Read the phase centre from a MeasurementSet.
+
+    Parameters
+    ----------
+    ms : str
+        Name of the Observation's MeasurementSet (MS), used for calculating
+        the beam and phase centre.
+    degrees : bool
+        If True, return the phase centre in degrees, else return it in radians.
+
+    Returns
+    -------
+    phase_dir : list of floats
+        Observation's phase centre (RA/Dec) in degrees or radians.
+    """
+    field_table = tables.table(f'{ms}::FIELD')
+    phase_dir = field_table[0]['PHASE_DIR'][0]
+
+    return np.degrees(phase_dir) if degrees else phase_dir
+
+
+if __name__ == '__main__':
+    parser = argparse.ArgumentParser(description='Obtain a selection of the brightest sources in a skymodel or source catalogue')
+
+    parser.add_argument('measurementset', type=str, help='Observation MeasurementSet')
+    parser.add_argument('-c', '--catalogue', type=str, default='TGSS', help='Select catalogue to query. Uses LSMTool, which allows for querying either one of TGSS, WENSS, GSM, or NVSS. Default: TGSS.')
+    parser.add_argument('-s', '--skymodel', type=str, help='Skymodel read from disk.')
+    parser.add_argument('-r', '--radius', type=float, help='Catalalogue search radius in degrees. Uses the phase centre from the MeasurementSet.')
+    parser.add_argument('-n', '--sources', type=int, default=100, help='Number of brightest sources in the output CSV file.')
+    parser.add_argument('-o', '--output', type=str, help='Name of the output CSV file containing the brightest sources, with the following format: name, ra (deg), dec (deg), apparent flux (mJy).')
+    parser.add_argument('-b', '--no-beam', action='store_false', help='Don\'t apply the beam and return the absolute flux instead of the -- default -- apparent flux.')
+
+    cli_arguments = parser.parse_args()
+
+    main(cli_arguments.measurementset, cli_arguments.sources, skymodel_file=cli_arguments.skymodel,
+         catalogue=cli_arguments.catalogue, search_radius=cli_arguments.radius,
+         output_csv=cli_arguments.output, apply_beam=cli_arguments.no_beam)

steps/find_brightest_sources.cwl

+class: CommandLineTool
+cwlVersion: v1.2
+baseCommand:
+  - find_brightest_sources.py
+label: Find the brightest sources in the sky
+doc: |
+  This tool finds the brightest sources in the sky around the
+  phase centre and returns these sources in a CSV file.
+
+inputs:
+  - id: msin
+    type: Directory
+    inputBinding:
+      position: 6
+    doc: Observation MeasurementSet used to compute beam
+  - id: skymodel?
+    type: File?
+    inputBinding:
+      prefix: --skymodel
+    doc: Name of sky model in text format
+  - id: catalogue?
+    type: string
+    inputBinding:
+      prefix: --catalogue
+    default: TGSS
+    doc: Source catalogue
+  - id: radius
+    type: float?
+    inputBinding:
+      prefix: --radius
+    doc: Catalogue search radius in degrees
+  - id: sources?
+    type: int
+    inputBinding:
+      prefix: --sources
+    default: 10
+    doc: Maximum number of sources in the output
+  - id: no_beam
+    type: boolean?
+    inputBinding:
+      prefix: --no-beam
+    doc: Don't apply beam correction
+  - id: output_csv
+    type: string
+    inputBinding:
+      prefix: --output
+    doc: Name of output CSV source list
+
+outputs:
+  - id: csv
+    type: File
+    outputBinding:
+      glob: "$(inputs.output_csv)"