From 8b5c6dab4715193e893c1a2f861f7b45b1873f24 Mon Sep 17 00:00:00 2001 From: Thomas Williams <32936518+thomaswilliamsastro@users.noreply.github.com> Date: Fri, 18 Sep 2026 09:54:54 +0100 Subject: [PATCH] Mask noise-only channels in cleaning This PR adds a function to test if channels in a cube are noise-only, and if so to mask them in the clean mask and skip through multiscale/singlescale cleaning. There are two checks that go on here: the first is that the sample-size-independent departure score (calculated from skew and kurtosis of the data) is above a sigma-threshold to a null test. The second is that the Jarque-Bera statistic is above a sigma-threshold. There is a little complication here that statistics need to be corrected for the effective number of independent measurements (number of beams), but this is all handled within the code. By default, this check is off. There's one exposed parameter, `channel_noise_sigma_threshold`, for which a value of 3 seems reasonable. Lower values are more stringent (i.e. will say there's signal in the channel at lower significance). Speed up is good! For an example 7m single field, processing time goes from 367s to 211s. For a 7m mosaic, it's TODO. The cubes are also basically identical. Plots attached for single field and mosaic, the blue show the median of the data and 1/99th percentiles, the black the same for the difference between the two cubes. Any deviations are well within the noise of the data. - Added functions to casaMaskingRoutines to calculate and mask noise-only channels - Add `task_mask_noise_channels` function to handlerImagingChunked to mask noise channels, post-dirty beam but pre-multiscale - Add `skip_multiscale_if_mask_empty` to handlerImagingChunked, which mimics the existing `skip_singlescale_if_mask_empty`. True by default - Add `task_mask_noise_channels` function to handlerImaging to mask noise channels, post-dirty beam but pre-multiscale - Update CHANGELOG.md --- CHANGELOG.md | 1 + phangsPipeline/casaMaskingRoutines.py | 414 +++++++++++++++++++++++- phangsPipeline/handlerImaging.py | 139 ++++++-- phangsPipeline/handlerImagingChunked.py | 231 +++++++++---- 4 files changed, 686 insertions(+), 99 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index ea5d18c3..75dd9622 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -27,6 +27,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Improved convolution with spectral-cube, including new analytic support for small kernels (#385). - Add tests for utilsLines, and tidy up formatting (#403). - Add optional file string in derived for nonstandard workflows (#406). +- Add optional tests to check if cube channels are noise only, and if so to mask them out (). ### Changed diff --git a/phangsPipeline/casaMaskingRoutines.py b/phangsPipeline/casaMaskingRoutines.py index d59c996f..ac66ed46 100644 --- a/phangsPipeline/casaMaskingRoutines.py +++ b/phangsPipeline/casaMaskingRoutines.py @@ -7,10 +7,15 @@ import os import analysisUtils as au +import astropy.units as u import numpy as np import scipy.ndimage as ndimage from astropy.io import fits -from scipy.special import erfc +from astropy.wcs.utils import proj_plane_pixel_scales +from radio_beam import Beam +from scipy.special import erfc, ndtri_exp +from scipy.stats import kurtosis, skew +from spectral_cube import Projection, SpectralCube from . import casaStuff from . import casaCubeRoutines as ccr @@ -352,6 +357,413 @@ def write_mask(infile, outfile, mask, huge_cube_workaround=True): return True +def get_beam_fft( + beam: Beam, + pixscale: u.Quantity | None = None, + shape: tuple[int, int] | None = None, +) -> tuple[np.ndarray, np.ndarray]: + """Get the normalised FFT of a beam kernel + + Args: + beam: Beam object + pixscale: u.Quantity: Pixel scale + shape: tuple[int, int]: Shape for the + output FFT + + Returns: + np.ndarray: Beam kernel, FFT of the beam kernel + """ + + if pixscale is None: + raise ValueError("pixscale must be defined!") + + if shape is None: + raise ValueError("shape must be defined!") + + # Create a kernel from the beam to convolve with the random field + k = beam.as_kernel(pixscale=pixscale, + x_size=shape[1], + y_size=shape[0], + ).array + + # Normalise the kernel + k = k / np.sqrt(np.sum(k ** 2)) + + # Get the FFT of the kernel + k_fft = np.fft.rfft2( + k, + ) + + return k, k_fft + + +def convolve_with_fft( + data: np.ndarray, + k_fft: np.ndarray, +) -> np.ndarray: + """Convolve data with a kernel using FFT + + Args: + data: np.ndarray: Input data to convolve + k_fft: np.ndarray: FFT of the kernel + + Returns: + np.ndarray: Convolved data + """ + + # Take FFT of the data + data_fft = np.fft.rfft2(data) + + # Multiply by the kernel FFT + convolved_fft = data_fft * k_fft + + # Take inverse FFT to get convolved data + convolved_data = np.fft.ifftshift(np.fft.irfft2(convolved_fft, s=data.shape)) + + return convolved_data + + +def calculate_departure( + data: np.ndarray, + n_eff: float | None = None, +): + """Calculate departure score for a dataset + + Args: + data (np.ndarray): Input data + n_eff (float | None): Effective number of independent measurements. + If None, will not correct for sample size. + + Returns: + dict: Dictionary containing skew, kurtosis, observed departure, and corrected departure + """ + + g1 = skew(data, bias=False) + g2 = kurtosis(data, fisher=True, bias=False) + + # Sample-size-independent departure score + observed_departure = np.sqrt(g1 ** 2 + g2 ** 2 / 4) + + corrected_departure = None + + # Approximations for Gaussian-beam correlation + if n_eff is not None: + n_eff_skew = 1.5 * n_eff + n_eff_kurt = 2.0 * n_eff + + expected_d2 = ( + 6 / n_eff_skew + + 6 / n_eff_kurt + ) + + corrected_departure = np.sqrt( + max(0.0, observed_departure ** 2 - expected_d2) + ) + + result = { + "g1": g1, + "g2": g2, + "observed_departure": observed_departure, + # "expected_gaussian_rms": np.sqrt(expected_d2), + "corrected_departure": corrected_departure, + } + + return result + + +def beam_corrected_gaussianity( + chan: Projection, +): + """Calculate beam-corrected statistics for a channel of a cube. + + Args: + chan (Projection): Channel to calculate statistics for. + + Returns: + dict: Dictionary containing skew, kurtosis, observed departure, corrected departure, + Jarque-Bera statistic, p-value, log p-value, log10 p-value, and sigma. + """ + + pix_per_beam = chan.pixels_per_beam + + # Pull out the data, use the underlying data + # since it's faster + data = chan._data + data = data[np.isfinite(data)] + + # If we have no valid data, return None + if data.size == 0: + return None + + # Calculate n_beams as the number of independent beams in the image + n_beams = data.size / pix_per_beam + + result = calculate_departure(data, n_eff=n_beams) + + if result["corrected_departure"] is not None: + departure = result["corrected_departure"] + else: + departure = result["observed_departure"] + + # Convert this to Jarque-Bera statistic + jb = n_beams * departure ** 2 / 6 + + # If we're below the noise floor, set + # to machine precision + if jb == 0: + jb = np.finfo(float).eps + + log_p = -jb / 2 + log10_p = log_p / np.log(10) + sigma = -ndtri_exp(log_p) + + # May still underflow, but sigma and log_p remain valid + p_value = np.exp(log_p) + + result.update( + { + "jarque_bera": jb, + "p_value": p_value, + "log_p_value": log_p, + "log10_p_value": log10_p, + "sigma": sigma, + } + ) + + return result + + +def calculate_departure_null( + beam: Beam, + pixscale: u.Quantity, + pix_per_beam: float | None = None, + shape: tuple[int, int] = (100, 100), + n_draws: int = 1000, +) -> np.ndarray: + """Calculate departure scores for a null test of Gaussianity, given a beam and pixel scale. + + Args: + beam (Beam): Beam object to convolve with pure noise + pixscale (u.Quantity): Pixel scale of the image + pix_per_beam (float | None): Number of pixels per beam. If None, + will not correct for the effective number of independent measurements. + shape (tuple[int, int]): Shape of the random noise field to generate. + Defaults to (100, 100). + n_draws (int): Number of random noise fields to generate. Defaults to 1000. + + Returns: + np.ndarray: Array of departure scores for each random noise field + """ + + _, k_fft = get_beam_fft( + beam, + pixscale=pixscale, + shape=shape, + ) + departure_null = np.full(n_draws, np.nan) + + rng = np.random.default_rng() + + for n_draw in range(n_draws): + + # Generate a random noise field, convolve and flatten + noise_field = rng.normal(size=shape) + noise_field = convolve_with_fft( + noise_field, + k_fft=k_fft, + ) + noise_field = noise_field.flatten() + + # Calculate the effective number of independent measurements + n_eff = None + if pix_per_beam is not None: + n_eff = noise_field.size / pix_per_beam + + result = calculate_departure(noise_field, + n_eff=n_eff, + ) + + if result["corrected_departure"] is not None: + departure = result["corrected_departure"] + else: + departure = result["observed_departure"] + + # Don't take the 0s, since they're meaningless + if departure == 0: + continue + + departure_null[n_draw] = departure + + return departure_null + + +def get_noise_only_channels( + f: str, + sigma_threshold: float | None = 2, +) -> list[bool]: + """Check where a cube only has noise channels + + There are two checks that go on here: the first is that the sample-size-independent + departure score (calculated from skew and kurtosis of the data) is above a sigma-threshold + to a null test. The second is that the Jarque-Bera statistic is above a sigma-threshold. + There is a little complication here that statistics need to be corrected for the effective + number of independent measurements (number of beams). + + Args: + f (str): Path to the cube + sigma_threshold (float | None): Threshold in sigma for departure from Gaussianity. + If None, will not check and just return False. + + Returns: + bool: True if all channels are consistent with being Gaussian noise, False otherwise + """ + + cube = SpectralCube.read(f) + cube.allow_huge_operations = True + + # If we're not checking, just return True + if sigma_threshold is None: + return [True] * cube.shape[0] + + # Get pixel scale in arcsec. Assume square pixels + pixscales = proj_plane_pixel_scales(cube.wcs.celestial) * u.deg + pixscale = [p.to(u.arcsec) for p in pixscales][0] + + # Take the first valid channel of the cube, assuming the beam stays relatively + # constant. This is a simplification, but should be fine for our purposes. + + mask = cube.get_mask_array() + valid_mask = list(np.sum(mask, axis=(1, 2)) > 0) + first_valid_chan = valid_mask.index(True) + + chan = cube[first_valid_chan] + beam = chan.beam + pix_per_beam = chan.pixels_per_beam + + departure_null = calculate_departure_null( + beam=beam, + pixscale=pixscale, + pix_per_beam=pix_per_beam, + ) + + # Calculate a "typical" departure from the Gaussian field + mean_departure_null = np.nanmean(departure_null) + std_departure_null = np.nanstd(departure_null) + + # Now calculate statistics for each channel in the cube + sigma = np.full(cube.shape[0], np.nan) + departure = np.full(cube.shape[0], np.nan) + + for chan_idx in range(cube.shape[0]): + + result = beam_corrected_gaussianity(cube[chan_idx]) + + if result is not None: + sigma[chan_idx] = result["sigma"] + if result["corrected_departure"] is not None: + d = result["corrected_departure"] + else: + d = result["observed_departure"] + departure[chan_idx] = d + + # Calculate noise-only channels + noise_only_channels = np.logical_or( + sigma < sigma_threshold, + departure < sigma_threshold * std_departure_null + mean_departure_null, + ) + + # Convert to a strict list of booleans + noise_only_channels = [bool(x) for x in noise_only_channels] + + return noise_only_channels + +def mask_noise_channels( + imaging_method='tclean', + cube_root=None, + suffix_in='', + suffix_out='', + operation='AND', + sigma_threshold=3, +): + """Mask out channels that are consistent with being noise-only.""" + + if imaging_method == 'sdintimaging': + cube_root += '.joint.cube' + + if not os.path.isdir(cube_root + '.image' + suffix_in): + logger.error('Data file not found: "' + cube_root + '.image' + suffix_in + '"') + logger.info('Need CUBE_ROOT.image to be an image file.') + logger.info('Returning. Generalize the code if you want different syntax.') + return + + header = casaStuff.imhead(cube_root + '.image' + suffix_in) + if header['axisnames'][2] == 'Frequency': + spec_axis = 2 + else: + spec_axis = 3 + + f = cube_root + '.image' + suffix_in + + logger.info('Reading cube.') + cube = read_cube(f, huge_cube_workaround=True) + + mask = np.ones(cube.shape, dtype=bool) + + logger.info("Finding noise-only channels") + noise_only_channels = get_noise_only_channels( + f, + sigma_threshold=sigma_threshold, + ) + total_noise_only = sum(noise_only_channels) + logger.info(f"{total_noise_only}/{len(noise_only_channels)} channels identified as noise-only.") + + for chan_idx in range(cube.shape[spec_axis]): + + if noise_only_channels[chan_idx]: + + # Slice out the channel + slc = [slice(None)] * len(mask.shape) + slc[spec_axis] = slice(chan_idx, chan_idx + 1) + slc = tuple(slc) + mask[slc] = False + + # Expect to be here with minimal memory footprint and mask + # created. + + if operation == 'AND' or operation == 'OR': + if os.path.isdir(cube_root + '.mask' + suffix_out): + old_mask = read_cube(cube_root + '.mask' + suffix_out, huge_cube_workaround=True) + else: + logger.info("Operation AND/OR requested but no previous mask found.") + logger.info("... will set operation=NEW.") + operation = 'NEW' + + logger.info('Joining with old mask.') + if operation == 'AND': + mask = mask * old_mask + if operation == 'OR': + mask = (mask + old_mask) > 0 + if operation == 'NEW': + mask = mask + else: + del old_mask + + logger.info('Recasting as an int.') + # this might be better: mask.astype(int, copy=False) + # mask = mask.astype(int) + mask = mask.astype(np.int32) + + # Export the image to fits, put in the mask and convert back to a CASA image + logger.info("Writing mask to disk") + + write_mask( + cube_root + ".image" + suffix_in, + cube_root + ".mask" + suffix_out, + mask, + huge_cube_workaround=True, + ) + + def signal_mask( imaging_method='tclean', cube_root=None, diff --git a/phangsPipeline/handlerImaging.py b/phangsPipeline/handlerImaging.py index e0084b69..2a4f116c 100644 --- a/phangsPipeline/handlerImaging.py +++ b/phangsPipeline/handlerImaging.py @@ -158,7 +158,7 @@ def _fname_dict( def loop_imaging( self, do_all=False, - imaging_method='tclean', + imaging_method="tclean", imaging_method_override=None, do_dirty_image=False, do_revert_to_dirty=False, @@ -174,11 +174,12 @@ def loop_imaging( do_export_to_fits=False, convergence_fracflux=0.01, convergence_noise_z_threshold=None, + channel_noise_sigma_threshold=None, singlescale_threshold_value=1.0, extra_ext_in=None, suffix_in=None, extra_ext_out=None, - recipe='phangsalma', + recipe="phangsalma", make_directories=True, dynamic_sizing=True, force_square=False, @@ -202,6 +203,10 @@ def loop_imaging( this is disabled by setting to None. We suggest setting this to 2.0 when enabled as a conservative choice. + We also have an optional check to mask out channels where there is no + signal. To use this, set channel_noise_sigma_threshold to a value (2 or 3 is + appropriate, lower values are stricter). + Parameters ---------- do_all : bool @@ -374,6 +379,7 @@ def loop_imaging( do_export_to_fits=do_export_to_fits, convergence_fracflux=convergence_fracflux, convergence_noise_z_threshold=convergence_noise_z_threshold, + channel_noise_sigma_threshold=channel_noise_sigma_threshold, singlescale_threshold_value=singlescale_threshold_value, dynamic_sizing=dynamic_sizing, force_square=force_square, @@ -916,6 +922,63 @@ def task_multiscale_clean( move_sd_psf(input_root=clean_call.get_param('imagename')) return () + + @CleanCallFunctionDecorator + def task_mask_noise_channels( + self, + clean_call=None, + product=None, + imaging_method='tclean', + channel_noise_sigma_threshold: float = 3.0, + ): + """ + Create a mask to remove noise-only channels + + Args: + channel_noise_sigma_threshold (float, optional): Channel noise threshold. + Defaults to 3. + """ + + if product is None: + logger.error("Require a product. Returning.") + raise Exception("Require a product. Returning.") + return () + + # Get fname dict + fname_dict = self._fname_dict(product=product, imagename=clean_call.get_param('imagename'), + imaging_method=imaging_method) + + # get imagename + + imagename = fname_dict['image'] + if not os.path.isdir(imagename): + logger.error("Image not found: " + imagename) + return () + + # print message + logger.info("") + logger.info("&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%") + logger.info("Creating noise channel mask for:") + logger.info(str(imagename)) + logger.info("&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%") + logger.info("") + + if self._dry_run: + return() + + # Noise channel mask + msr.mask_noise_channels( + imaging_method=imaging_method, + cube_root=fname_dict['root'], + sigma_threshold=channel_noise_sigma_threshold, + suffix_in="", + suffix_out="", + operation="AND", + ) + + clean_call.set_param('usemask', 'user') + + return () @CleanCallFunctionDecorator def task_singlescale_mask( @@ -1122,35 +1185,36 @@ def task_export_to_fits( ############################# def recipe_phangsalma_imaging( - self, - target=None, - product=None, - config=None, - extra_ext_in=None, - suffix_in=None, - extra_ext_out=None, - imaging_method='tclean', - imaging_method_override=None, - do_dirty_image=True, - do_revert_to_dirty=True, - do_read_clean_mask=True, - do_multiscale_clean=True, - do_revert_to_multiscale=True, - do_singlescale_mask=True, - singlescale_mask_high_snr=None, - singlescale_mask_low_snr=None, - singlescale_mask_absolute=False, - do_singlescale_clean=True, - do_revert_to_singlescale=True, - do_export_to_fits=True, - convergence_fracflux=0.01, - convergence_noise_z_threshold=None, - singlescale_threshold_value=1.0, - dynamic_sizing=True, - force_square=False, - export_dirty=False, - export_multiscale=False, - overwrite=False, + self, + target=None, + product=None, + config=None, + extra_ext_in=None, + suffix_in=None, + extra_ext_out=None, + imaging_method="tclean", + imaging_method_override=None, + do_dirty_image=True, + do_revert_to_dirty=True, + do_read_clean_mask=True, + do_multiscale_clean=True, + do_revert_to_multiscale=True, + do_singlescale_mask=True, + singlescale_mask_high_snr=None, + singlescale_mask_low_snr=None, + singlescale_mask_absolute=False, + do_singlescale_clean=True, + do_revert_to_singlescale=True, + do_export_to_fits=True, + convergence_fracflux=0.01, + convergence_noise_z_threshold=None, + channel_noise_sigma_threshold=None, + singlescale_threshold_value=1.0, + dynamic_sizing=True, + force_square=False, + export_dirty=False, + export_multiscale=False, + overwrite=False, ): """ PHANGS-ALMA basic imaging recipe. @@ -1325,6 +1389,15 @@ def recipe_phangsalma_imaging( # clean_call.set_param('usemask', 'pb') # clean_call.set_param('pbmask', 0.2) + # Mask noise channels if we have a sigma threshold defined + if channel_noise_sigma_threshold is not None: + self.task_mask_noise_channels( + clean_call=clean_call, + product=product, + imaging_method=imaging_method, + channel_noise_sigma_threshold=channel_noise_sigma_threshold, + ) + # Dynamic sizing if dynamic_sizing: @@ -1384,6 +1457,10 @@ def recipe_phangsalma_imaging( clean_call.set_param('cell', cell, nowarning=True) clean_call.set_param('imsize', imsize, nowarning=True) + # Since we've reset the clean call, need to reassign the noise channel mask if we have one + if channel_noise_sigma_threshold is not None: + clean_call.set_param('usemask', 'user') + # Make a signal-to-noise based mask for use in singlescale clean. if do_singlescale_mask: diff --git a/phangsPipeline/handlerImagingChunked.py b/phangsPipeline/handlerImagingChunked.py index 2ab62ee9..483f7030 100644 --- a/phangsPipeline/handlerImagingChunked.py +++ b/phangsPipeline/handlerImagingChunked.py @@ -310,6 +310,7 @@ def run_imaging( do_dirty_image=False, do_revert_to_dirty=False, do_read_clean_mask=False, + skip_multiscale_if_mask_empty=True, do_multiscale_clean=False, do_revert_to_multiscale=False, do_singlescale_mask=False, @@ -325,6 +326,7 @@ def run_imaging( convergence_fracflux=0.01, convergence_noise_z_threshold=None, singlescale_threshold_value=1.0, + channel_noise_sigma_threshold=None, extra_ext_in=None, suffix_in=None, extra_ext_out=None, @@ -355,6 +357,9 @@ def run_imaging( If True, revert to dirty image. do_read_clean_mask : bool If True, read clean mask. + skip_multiscale_if_mask_empty : bool + If True, will not perform multiscale cleaning if + clean mask is empty do_multiscale_clean : bool If True, multiscale clean. do_revert_to_multiscale : bool @@ -381,6 +386,9 @@ def run_imaging( greater than 2.0. singlescale_threshold_value : float Threshold value for singlescale clean. Default is 1.0. + channel_noise_sigma_threshold : float + Threshold value for defining a channel as purely noise. Default is None, + which will not perform this check. extra_ext_in : str Extra extension for input files. suffix_in : str @@ -450,6 +458,7 @@ def run_imaging( do_dirty_image=do_dirty_image, do_revert_to_dirty=do_revert_to_dirty, do_read_clean_mask=do_read_clean_mask, + skip_multiscale_if_mask_empty=skip_multiscale_if_mask_empty, do_multiscale_clean=do_multiscale_clean, do_revert_to_multiscale=do_revert_to_multiscale, do_singlescale_mask=do_singlescale_mask, @@ -463,6 +472,7 @@ def run_imaging( do_cleanup=do_cleanup, convergence_fracflux=convergence_fracflux, convergence_noise_z_threshold=convergence_noise_z_threshold, + channel_noise_sigma_threshold=channel_noise_sigma_threshold, singlescale_threshold_value=singlescale_threshold_value, dynamic_sizing=dynamic_sizing, force_square=force_square, @@ -1142,13 +1152,15 @@ def task_read_clean_mask( @CleanCallFunctionDecorator def task_multiscale_clean( - self, - chunk_num=None, - convergence_fracflux=0.01, - convergence_noise_z_threshold=None, - backup=True, - gather_chunks_into_cube=False, - remove_chunks=False, + self, + chunk_num=None, + convergence_fracflux=0.01, + convergence_noise_z_threshold=None, + channel_noise_sigma_threshold=None, + backup=True, + gather_chunks_into_cube=False, + remove_chunks=False, + skip_multiscale_if_mask_empty=True, ): """ Run a multiscale clean loop to convergence. This task @@ -1191,6 +1203,10 @@ def task_multiscale_clean( logger.warning("I expected a multiscale or mtmfs deconvolver but got " + str( this_clean_call.get_param('deconvolver')) + ".") raise Exception("Incorrect clean call! Should have a multiscale or mtmfs deconvolver.") + + # If we've masked noise channels, make sure we use the mask + if channel_noise_sigma_threshold is not None: + this_clean_call.set_param("usemask", "user") logger.info("") logger.info("&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%") @@ -1199,34 +1215,50 @@ def task_multiscale_clean( logger.info("This is {0} out of {1} to be imaged".format(ii+1, len(chunks_iter))) logger.info("&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%") logger.info("") + + skip_this_step = False + if skip_multiscale_if_mask_empty: + if self.imaging_method == 'sdintimaging': + mask_name = "{}.joint.cube.mask".format(self.chunk_params[this_chunk_num]['full_imagename']) + else: + mask_name = "{}.mask".format(self.chunk_params[this_chunk_num]['full_imagename']) - imr.clean_loop(clean_call=this_clean_call, - imaging_method=self.imaging_method, - record_file=this_clean_call.get_param('imagename') + '_multiscale_record.txt', - niter_base_perchan=10, - niter_growth_model='geometric', - niter_growth_factor=2.0, - niter_saturation_perchan=1000, - niter_other_input=None, - cycleniter_base=100, - cycleniter_growth_model='linear', - cycleniter_growth_factor=1.0, - cycleniter_saturation_value=1000, - cycleniter_other_input=None, - threshold_type='snr', - threshold_value=4.0, - min_loops=3, - max_loops=20, - max_total_niter=None, - convergence_fracflux=convergence_fracflux, - convergence_totalflux=None, - convergence_fluxperniter=None, - use_absolute_delta=True, - stop_at_negative=True, - remask_each_loop=False, - force_dirty_image=False, - convergence_noise_z_threshold=convergence_noise_z_threshold, - ) + if os.path.exists(mask_name): + mask_stats = msr.stat_cube(cube_file=mask_name) + if mask_stats['sum'] == 0: + skip_this_step = True + logger.info("") + logger.info("The clean mask is empty and SKIP_MULTISCALE_IF_MASK_EMPTY is True. Skipping the multiscale clean step.") + logger.info("") + + if not skip_this_step: + imr.clean_loop(clean_call=this_clean_call, + imaging_method=self.imaging_method, + record_file=this_clean_call.get_param('imagename') + '_multiscale_record.txt', + niter_base_perchan=10, + niter_growth_model='geometric', + niter_growth_factor=2.0, + niter_saturation_perchan=1000, + niter_other_input=None, + cycleniter_base=100, + cycleniter_growth_model='linear', + cycleniter_growth_factor=1.0, + cycleniter_saturation_value=1000, + cycleniter_other_input=None, + threshold_type='snr', + threshold_value=4.0, + min_loops=3, + max_loops=20, + max_total_niter=None, + convergence_fracflux=convergence_fracflux, + convergence_totalflux=None, + convergence_fluxperniter=None, + use_absolute_delta=True, + stop_at_negative=True, + remask_each_loop=False, + force_dirty_image=False, + convergence_noise_z_threshold=convergence_noise_z_threshold, + ) if backup: imr.copy_imaging( @@ -1245,6 +1277,59 @@ def task_multiscale_clean( os.chdir(cwd) + @CleanCallFunctionDecorator + def task_mask_noise_channels( + self, + chunk_num: int | None = None, + channel_noise_sigma_threshold: float = 3.0, + ): + """ + Create a mask to remove noise-only channels + + Args: + chunk_num (int, optional): Chunk number. Defaults + to None, which will loop over all chunks + channel_noise_sigma_threshold (float, optional): Channel noise threshold. + Defaults to 3. + """ + + # get imagename + + imagename = self._fname_dict(self.image_root, + imaging_method=self.imaging_method)['image'] + + # print message + logger.info("") + logger.info("&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%") + logger.info("Creating noise channel mask for:") + logger.info(str(imagename)) + logger.info("&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%") + logger.info("") + + if self._dry_run: + return True + + chunks_iter = self.return_valid_chunks(chunk_num=chunk_num) + + for ii, this_chunk_num in enumerate(chunks_iter): + logger.info("&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%") + logger.info("Noise channel masking chunk {0} of {1}".format(ii, len(chunks_iter))) + logger.info("&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%") + + image_root = self.chunk_params[this_chunk_num]["full_imagename"] + + # Noise channel mask + msr.mask_noise_channels( + imaging_method=self.imaging_method, + cube_root=image_root, + sigma_threshold=channel_noise_sigma_threshold, + suffix_in="", + suffix_out="", + operation="AND", + ) + + return True + @CleanCallFunctionDecorator def task_singlescale_mask( @@ -1422,12 +1507,13 @@ def task_singlescale_clean( else: mask_name = "{}.mask".format(self.chunk_params[this_chunk_num]['full_imagename']) - mask_stats = msr.stat_cube(cube_file=mask_name) - if mask_stats['sum'] == 0: - skip_this_step = True - logger.info("") - logger.info("The clean mask is empty and SKIP_SINGLESCALE_IF_MASK_EMPTY is True. Skipping the singlescale clean step.") - logger.info("") + if os.path.exists(mask_name): + mask_stats = msr.stat_cube(cube_file=mask_name) + if mask_stats['sum'] == 0: + skip_this_step = True + logger.info("") + logger.info("The clean mask is empty and SKIP_SINGLESCALE_IF_MASK_EMPTY is True. Skipping the singlescale clean step.") + logger.info("") if not skip_this_step: imr.clean_loop( @@ -1579,33 +1665,35 @@ def task_cleanup(self, tag=None, chunk_num=None): ############################# def recipe_phangsalma_imaging( - self, - chunk_num=None, - extra_ext_in=None, - suffix_in=None, - extra_ext_out=None, - do_dirty_image=True, - do_revert_to_dirty=True, - do_read_clean_mask=True, - do_multiscale_clean=True, - do_revert_to_multiscale=True, - do_singlescale_mask=True, - singlescale_mask_high_snr=None, - singlescale_mask_low_snr=None, - singlescale_mask_absolute=False, - skip_singlescale_if_mask_empty=True, - do_singlescale_clean=True, - do_revert_to_singlescale=True, - do_recombine_cubes=False, - do_export_to_fits=True, - do_cleanup=True, - convergence_fracflux=0.01, - convergence_noise_z_threshold=None, - singlescale_threshold_value=1.0, - dynamic_sizing=True, - force_square=False, - export_multiscale=False, - overwrite=False, + self, + chunk_num=None, + extra_ext_in=None, + suffix_in=None, + extra_ext_out=None, + do_dirty_image=True, + do_revert_to_dirty=True, + do_read_clean_mask=True, + skip_multiscale_if_mask_empty=True, + do_multiscale_clean=True, + do_revert_to_multiscale=True, + do_singlescale_mask=True, + singlescale_mask_high_snr=None, + singlescale_mask_low_snr=None, + singlescale_mask_absolute=False, + skip_singlescale_if_mask_empty=True, + do_singlescale_clean=True, + do_revert_to_singlescale=True, + do_recombine_cubes=False, + do_export_to_fits=True, + do_cleanup=True, + convergence_fracflux=0.01, + convergence_noise_z_threshold=None, + channel_noise_sigma_threshold=None, + singlescale_threshold_value=1.0, + dynamic_sizing=True, + force_square=False, + export_multiscale=False, + overwrite=False, ): """ PHANGS-ALMA basic imaging recipe. @@ -1673,11 +1761,20 @@ def recipe_phangsalma_imaging( # target=target, config=config, product=product, # imaging_method=imaging_method) + # Mask noise channels if we have a sigma threshold defined + if channel_noise_sigma_threshold is not None: + self.task_mask_noise_channels( + chunk_num=chunk_to_iter, + channel_noise_sigma_threshold=channel_noise_sigma_threshold, + ) + # Run a multiscale clean until it converges. if do_multiscale_clean: self.task_multiscale_clean(chunk_num=chunk_to_iter, convergence_fracflux=convergence_fracflux, convergence_noise_z_threshold=convergence_noise_z_threshold, + channel_noise_sigma_threshold=channel_noise_sigma_threshold, + skip_multiscale_if_mask_empty=skip_multiscale_if_mask_empty, gather_chunks_into_cube=False, )