diff --git a/CHANGELOG.md b/CHANGELOG.md index ea5d18c..75dd962 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 d59c996..ac66ed4 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 e0084b6..2a4f116 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 2ab62ee..483f703 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, )