Initial continuum-subtraction and spectral line prefect flow (#192)

* initial beginnings of cube imaging script

* Initial toying around with subtracted cube

* added subtract cube operaiton

* fixed installed flint flow entrypoint

* relocated SubtractFieldOptions/fixed argparser / test

* calling the flow

* using ummaped

* patching in channel range into naming

* rejigging the wsclean names that are generated

* filter corrected#

* adding linmos filter

* some logging output on valueerror

* added a NaN beam shape

* setting a zero'd beam when NaN beam detected

* bumped the default number

* Added copyfile for bad conv beams

* copy zero beams

* no-reorder added to WSCleanOptions

* corected local variable unbound

* rejigged the wsclean delete cleanups

* prefix log

* updated prefix creation

* tweaked logs

* unlinking preconvolved images after convolve

* typo fix

* added to the chanelog

* batching the results

* batched limit added to SubtractFieldOptions

* added else to for loop

* fixed the check on list

* treat as rate limit

* single item waiting

* silly list index missing

* added a more aggressive loop for complete

* sleep

* rejigged again

* inverted state test

* no sleep

* sleeping

* 10 to 5

* added a stagger_delay_seconds

* removed comments

* fixed subtract column output / test

* initial subtraction

* commited options to SubtractFieldOptions

* added mandatory argument

* a silly

* initiall addmodel prediction

* options fix

* unmapped

* fixede mode input

* changelog

* optional attempt_subtract

* typo fix

* silly typo pirate

* added subflow for addmodel

* missing keyword argument

* another keyword typo

* silly

* more silly pirate business

* typing and binding error

* corrected column

* updated log message

* added a flint_ skip options in the wsclean logger

* added some notes to the changelog

* fixed type errors raised on pre-commit

---------

Co-authored-by: tgalvin <[email protected]>

CHANGELOG.md

@@ -1,5 +1,22 @@
 # Change log
 
+# dev
+
+- Created `subtract_cube_pipeline.py`. Associated changes include:
+  - wsclean imaging will delete files while still in temporary storage (for
+    instance is a ramdisk is used)
+  - added option to prefect convol task to delete files once they have be
+    convoled
+  - predicting model produced by `wsclean -save-source-list` into measurement
+    set
+  - using `taql` to subtract model from nominated data column
+  - added a 'flow_addmodel_to_mss` to allow different task runner to be
+    specified
+  - added `flint_no_log_wsclean_output` to `WSCleanOptions` to disable wsclean
+    logging
+  - added flag to disable logging singularity command output to the
+    `flint.logging.logger`
+
 # 0.2.8
 
 - added `wrapper_options_from_strategy` decorator helper function

README.md

@@ -112,12 +112,8 @@ tasks together (outlined above) into a single data-processing pipeline.
   observation sequence.
 - `flint_flow_continuum_pipeline`: Performs bandpass calibration, solution
   copying, imaging, self-calibration and mosaicing.
-- `flint_flow_cointinuum_mask_pipeline`: Performs bandpass calibration, solution
-  copying, imaging, self-calibration and mosaicing. In this flow a process to
-  construct a robust clean mask is performed by exploiting an initial imaging
-  round. The field image is constructed across all beams, S/N clipping is
-  performed, then guard masks on a per-beam basis are extracted. This pipeline
-  has fallen out of use and could be removed.
+- `flint_flow_subtract_cube_pipeline`: Subtract a continuum model and image the
+  residual data.
 
 ## Sky-model catalogues
 

flint/coadd/linmos.py

@@ -526,6 +526,7 @@ def linmos_images(
     container: Path = Path("yandasoft.sif"),
     cutoff: float = 0.001,
     pol_axis: Optional[float] = None,
+    trim_linmos_fits: bool = True,
 ) -> LinmosCommand:
     """Create a linmos parset file and execute it.
 
@@ -538,6 +539,7 @@ def linmos_images(
         container (Path, optional): Path to the singularity container that has the yandasoft tools. Defaults to Path('yandasoft.sif').
         cutoff (float, optional): Pixels whose primary beam attenuation is below this cutoff value are blanked. Defaults to 0.001.
         pol_axis (Optional[float], optional): The physical oritentation of the ASKAP third-axis in radians. Defaults to None.
+        trim_linmos_fits (bool, optional): Attempt to trim the output linmos files of as much empty space as possible. Defaults to True.
 
     Returns:
         LinmosCommand: The linmos command executed and the associated parset file
@@ -577,11 +579,12 @@ def linmos_images(
     )
 
     # Trim the fits image to remove empty pixels
-    image_trim_results = trim_fits_image(image_path=linmos_names.image_fits)
-    trim_fits_image(
-        image_path=linmos_names.weight_fits,
-        bounding_box=image_trim_results.bounding_box,
-    )
+    if trim_linmos_fits:
+        image_trim_results = trim_fits_image(image_path=linmos_names.image_fits)
+        trim_fits_image(
+            image_path=linmos_names.weight_fits,
+            bounding_box=image_trim_results.bounding_box,
+        )
 
     return linmos_cmd
 

flint/configuration.py

@@ -30,7 +30,7 @@
 # into there
 
 KNOWN_HEADERS = ("defaults", "initial", "selfcal", "version")
-KNOWN_OPERATIONS = ("stokesv",)
+KNOWN_OPERATIONS = ("stokesv", "subtractcube")
 FORMAT_VERSION = 0.1
 MODE_OPTIONS_MAPPING = {
     "wsclean": WSCleanOptions,
@@ -81,6 +81,32 @@ def copy_and_timestamp_strategy_file(output_dir: Path, input_yaml: Path) -> Path
     return Path(stamped_imaging_strategy)
 
 
+def _load_and_copy_strategy(
+    output_split_science_path: Path, imaging_strategy: Optional[Path] = None
+) -> Union[Strategy, None]:
+    """Load a strategy file and copy a timestamped version into the output directory
+    that would contain the science processing.
+
+    Args:
+        output_split_science_path (Path): Where the strategy file should be copied to (where the data would be processed)
+        imaging_strategy (Optional[Path], optional): Location of the strategy file. Defaults to None.
+
+    Returns:
+        Union[Strategy, None]: The loadded strategy file if provided, `None` otherwise
+    """
+    return (
+        load_strategy_yaml(
+            input_yaml=copy_and_timestamp_strategy_file(
+                output_dir=output_split_science_path,
+                input_yaml=imaging_strategy,
+            ),
+            verify=True,
+        )
+        if imaging_strategy
+        else None
+    )
+
+
 def get_selfcal_options_from_yaml(input_yaml: Optional[Path] = None) -> Dict:
     """Stub to represent interaction with a configurationf ile
 

flint/convol.py

@@ -7,6 +7,7 @@
 import warnings
 from argparse import ArgumentParser
 from pathlib import Path
+from shutil import copyfile
 from typing import Collection, List, Literal, NamedTuple, Optional
 
 import astropy.units as u
@@ -198,10 +199,14 @@ def get_common_beam(
     if cutoff:
         logger.info(f"Setting beam cutoff to {cutoff} arcseconds. ")
 
-    beam, beams = beamcon_2D.get_common_beam(files=list(image_paths), cutoff=cutoff)
+    try:
+        beam, beams = beamcon_2D.get_common_beam(files=list(image_paths), cutoff=cutoff)
 
-    beam_shape = BeamShape.from_radio_beam(beam)
-    logger.info(f"Constructed {beam_shape=}")
+        beam_shape = BeamShape.from_radio_beam(beam)
+        logger.info(f"Constructed {beam_shape=}")
+    except ValueError:
+        logger.info("The beam was not constrained. Setting to NaNs")
+        beam_shape = BeamShape(bmaj_arcsec=np.nan, bmin_arcsec=np.nan, bpa_deg=np.nan)
 
     return beam_shape
 
@@ -231,29 +236,50 @@ def convolve_images(
     if cutoff:
         logger.info(f"Supplied cutoff of {cutoff} arcsecond")
 
+    if not np.isfinite(beam_shape.bmaj_arcsec):
+        logger.info("Beam shape is not defined. Copying files into place. ")
+
+        conv_image_paths = [
+            Path(str(image_path).replace(".fits", f".{convol_suffix}.fits"))
+            for image_path in image_paths
+        ]
+        # If the beam is not defined, simply copy the file into place. Although
+        # this takes up more space, it is not more than otherwise
+        for original_path, copy_path in zip(image_paths, conv_image_paths):
+            logger.info(f"Copying {original_path=} {copy_path=}")
+            copyfile(original_path, copy_path)
+
+        return conv_image_paths
+
     radio_beam = Beam(
         major=beam_shape.bmaj_arcsec * u.arcsecond,
         minor=beam_shape.bmin_arcsec * u.arcsecond,
         pa=beam_shape.bpa_deg * u.deg,
     )
 
-    conv_image_paths: List[Path] = []
+    return_conv_image_paths: List[Path] = []
 
     for image_path in image_paths:
-        logger.info(f"Convolving {str(image_path.name)}")
-        beamcon_2D.beamcon_2d_on_fits(
-            file=image_path,
-            outdir=None,
-            new_beam=radio_beam,
-            conv_mode="robust",
-            suffix=convol_suffix,
-            cutoff=cutoff,
-        )
-        conv_image_paths.append(
-            Path(str(image_path).replace(".fits", f".{convol_suffix}.fits"))
+        convol_output_path = Path(
+            str(image_path).replace(".fits", f".{convol_suffix}.fits")
         )
+        header = fits.getheader(image_path)
+        if header["BMAJ"] == 0.0:
+            logger.info(f"Copying {image_path} to {convol_output_path=} for empty beam")
+            copyfile(image_path, convol_output_path)
+        else:
+            logger.info(f"Convolving {str(image_path.name)}")
+            beamcon_2D.beamcon_2d_on_fits(
+                file=image_path,
+                outdir=None,
+                new_beam=radio_beam,
+                conv_mode="robust",
+                suffix=convol_suffix,
+                cutoff=cutoff,
+            )
+        return_conv_image_paths.append(convol_output_path)
 
-    return conv_image_paths
+    return return_conv_image_paths
 
 
 def get_parser() -> ArgumentParser:

flint/exceptions.py

@@ -4,6 +4,10 @@ class MSError(Exception):
     pass
 
 
+class FrequencyMismatchError(Exception):
+    """Raised when there are differences in frequencies"""
+
+
 class PhaseOutlierFitError(Exception):
     """Raised when the phase outlier fit routine fails."""