Source code for pipeline.hif.heuristics.cleanbox

import os.path

import numpy as np

import pipeline.infrastructure as infrastructure
import pipeline.infrastructure.utils as utils
from pipeline.infrastructure import casa_tools

LOG = infrastructure.get_logger(__name__)


[docs]def analyse_clean_result(multiterm, model, restored, residual, pb, cleanmask, pblimit_image=0.2, pblimit_cleanmask=0.3, cont_freq_ranges=None): qaTool = casa_tools.quanta if pb == '': pb = None if multiterm: extension = '.tt0' else: extension = '' # get the sum of the model image to find how much flux has been # cleaned model_sum = None if model is not None: with casa_tools.ImageReader(model + extension) as image: model_stats = image.statistics(robust=False) model_sum = model_stats['sum'][0] LOG.debug('Sum of model: %s' % model_sum) LOG.debug('Fixing coordsys of pb and cleanmask') with casa_tools.ImageReader(residual + extension) as image: csys = image.coordsys() if pb is not None: with casa_tools.ImageReader(pb + extension) as image: image.setcoordsys(csys.torecord()) if cleanmask is not None and os.path.exists(cleanmask): with casa_tools.ImageReader(cleanmask) as image: image.setcoordsys(csys.torecord()) csys.done() with casa_tools.ImageReader(residual + extension) as image: # get the rms of the residual image inside the cleaned area LOG.todo('Cannot use dirname in mask') residual_cleanmask_rms = None if cleanmask is not None and os.path.exists(cleanmask): # Area inside clean mask statsmask = '"%s" > 0.1' % (os.path.basename(cleanmask)) resid_clean_stats = image.statistics(mask=statsmask, robust=False) try: residual_cleanmask_rms = resid_clean_stats['rms'][0] LOG.info('Residual rms inside cleaned area: %s' % residual_cleanmask_rms) except: pass # and the rms of the residual image outside the cleaned area residual_non_cleanmask_rms = None if pb is not None and os.path.exists(pb+extension): have_mask = True statsmask = '("%s" > %f) && ("%s" < %f)' % (os.path.basename(pb)+extension, pblimit_image, os.path.basename(pb)+extension, pblimit_cleanmask) elif cleanmask is not None and os.path.exists(cleanmask): have_mask = True # Area outside clean mask statsmask = '"%s" < 0.1' % (os.path.basename(cleanmask)) else: have_mask = False statsmask = '' residual_stats = image.statistics(mask=statsmask, robust=False) try: residual_non_cleanmask_rms = residual_stats['rms'][0] if have_mask: LOG.info('Residual rms in the annulus: %s' % residual_non_cleanmask_rms) else: LOG.info('Residual rms across full area: %s' % residual_non_cleanmask_rms) except: pass # get the max, min of the residual image (avoiding the edges # where spikes can occur) if pb is not None and os.path.exists(pb+extension): residual_stats = image.statistics( mask='"%s" > %f' % (os.path.basename(pb)+extension, pblimit_image), robust=False) else: residual_stats = image.statistics(robust=False) try: residual_max = residual_stats['max'][0] residual_min = residual_stats['min'][0] except: residual_max = None residual_min = None LOG.info('Residual max: %s min: %s' % (residual_max, residual_min)) residual_stats = image.statistics(robust=True) residual_robust_rms = residual_stats['medabsdevmed'][0] * 1.4826 # see CAS-9631 LOG.debug('residual scaled MAD: %s' % residual_robust_rms) pbcor_image_min = None pbcor_image_max = None nonpbcor_imagename = None nonpbcor_image_non_cleanmask_rms = None nonpbcor_image_non_cleanmask_rms_min = None nonpbcor_image_non_cleanmask_rms_max = None nonpbcor_image_non_cleanmask_robust_rms = None nonpbcor_image_non_cleanmask_freq_ch1 = None nonpbcor_image_non_cleanmask_freq_chN = None nonpbcor_image_non_cleanmask_freq_frame = None nonpbcor_image_cleanmask_spectrum = None nonpbcor_image_cleanmask_spectrum_pblimit = None nonpbcor_image_cleanmask_npoints = None nonpbcor_image_statsmask = None if restored not in [None, '']: # get min and max of the pb-corrected cleaned result with casa_tools.ImageReader(restored.replace('.image', '.image%s' % extension)) as image: if pb is not None and os.path.exists(pb+extension): have_mask = True # Default is area pb > 0.3 statsmask = '"%s" > %f' % (os.path.basename(pb)+extension, pblimit_cleanmask) elif cleanmask is not None and os.path.exists(cleanmask): have_mask = True # Area inside clean mask statsmask = '"%s" > 0.1' % (os.path.basename(cleanmask)) else: have_mask = False statsmask = '' if 'TARGET' in image.miscinfo().get('intent', None): image_stats = image.statistics(mask=statsmask) else: # Restrict region to inner 25% x 25% of the image for calibrators to # avoid picking up sidelobes (PIPE-611) shape = image.shape() rgTool = casa_tools.regionmanager nPixels = max(shape[0], shape[1]) region = rgTool.box([nPixels*0.375-1, nPixels*0.375-1, 0, 0], [nPixels*0.625-1, nPixels*0.625-1, shape[1]-1, shape[2]-1]) image_stats = image.statistics(mask=statsmask, region=region) rgTool.done() pbcor_image_min = image_stats['min'][0] pbcor_image_max = image_stats['max'][0] if have_mask: LOG.debug('Clean pb-corrected image min in cleaned area: %s' % pbcor_image_min) LOG.debug('Clean pb-corrected image max in cleaned area: %s' % pbcor_image_max) else: LOG.debug('Clean pb-corrected image min in full area: %s' % pbcor_image_min) LOG.debug('Clean pb-corrected image max in full area: %s' % pbcor_image_max) # get RMS in non cleanmask area of non-pb-corrected cleaned result if restored.find('.image.pbcor') != -1: nonpbcor_imagename = restored.replace('.image.pbcor', '.image%s' % extension) else: nonpbcor_imagename = restored.replace('.image', '.image%s' % extension) # If possible use flattened clean mask for exclusion of areas for all spectral channels flattened_mask = None if cleanmask is not None and os.path.exists(cleanmask): if '.mask' in cleanmask: flattened_mask = cleanmask.replace('.mask', '.mask.flattened') elif '.cleanmask' in cleanmask: flattened_mask = cleanmask.replace('.cleanmask', '.cleanmask.flattened') else: raise 'Cannot handle clean mask name %s' % (os.path.basename(cleanmask)) with casa_tools.ImageReader(cleanmask) as image: flattened_mask_image = image.collapse(function='max', axes=[2, 3], outfile=flattened_mask) try: npoints_mask = flattened_mask_image.statistics(mask='"%s" > 0.1' % (os.path.basename(flattened_mask)), robust=False)['npts'] if npoints_mask.shape != (0,): nonpbcor_image_cleanmask_npoints = int(npoints_mask) else: nonpbcor_image_cleanmask_npoints = 0 except: nonpbcor_image_cleanmask_npoints = 0 flattened_mask_image.done() with casa_tools.ImageReader(nonpbcor_imagename) as image: # Get the image frequency axis for later plotting. imhead = image.summary(list=False) lcs = image.coordsys() try: csys = image.coordsys() freq_axis = csys.findaxisbyname('spectral') csys.done() except: num_axes = image.shape().shape[0] if num_axes > 3: LOG.warn("Can't find spectral axis. Assuming it is 3.") freq_axis = 3 elif num_axes > 2: LOG.warn("Can't find spectral axis. Assuming it is 2.") freq_axis = 2 elif num_axes == 2: LOG.error("No spectral axis found") freq_axis = -1 lcs.done() nonpbcor_image_non_cleanmask_freq_ch1 = qaTool.quantity(imhead['refval'][freq_axis] - imhead['refpix'][freq_axis] * imhead['incr'][freq_axis], imhead['axisunits'][freq_axis]) nonpbcor_image_non_cleanmask_freq_chN = qaTool.quantity(imhead['refval'][freq_axis] + (imhead['shape'][freq_axis] - imhead['refpix'][freq_axis]) * imhead['incr'][freq_axis], imhead['axisunits'][freq_axis]) # Get the spectral reference. Unfortunately this is coded in text # messages rather than a key/value pair. Hence the parsing code. try: for msg in imhead['messages'][1].split('\n'): msg_l = msg.lower() if 'spectral' in msg_l and 'reference' in msg_l: nonpbcor_image_non_cleanmask_freq_frame = msg.split(':')[1].strip() except: LOG.warn('Cannot determine spectral reference in %s. Assuming it is LSRK.' % (nonpbcor_imagename)) nonpbcor_image_non_cleanmask_freq_frame = 'LSRK' # define mask outside the cleaned area image_stats = None if pb is not None and os.path.exists(pb+extension) and cleanmask is not None and os.path.exists(cleanmask): pb_name = os.path.basename(pb)+extension have_mask = True # Annulus without clean mask statsmask = '("%s" < 0.1) && ("%s" > %f) && ("%s" < %f)' % \ (os.path.basename(flattened_mask), \ pb_name, pblimit_image, \ pb_name, pblimit_cleanmask) # Check for number of points per channel (PIPE-541): try: image_stats = image.statistics(mask=statsmask, robust=True, axes=[0, 1, 2], algorithm='chauvenet', maxiter=5) if image_stats['npts'].shape == (0,) or np.median(image_stats['npts']) < 10.0: # Switch to full annulus to avoid zero noise spectrum due to voluminous mask LOG.warn('Using full annulus for noise spectrum due to voluminous mask.') statsmask = '("%s" > %f) && ("%s" < %f)' % (pb_name, pblimit_image, pb_name, pblimit_cleanmask) image_stats = None except Exception as e: # Try full annulus as a fallback LOG.exception('Using full annulus for noise spectrum due to voluminous mask.', exc_info=e) statsmask = '("%s" > %f) && ("%s" < %f)' % (pb_name, pblimit_image, pb_name, pblimit_cleanmask) image_stats = None elif pb is not None and os.path.exists(pb+extension): pb_name = os.path.basename(pb)+extension have_mask = True # Full annulus statsmask = '("%s" > %f) && ("%s" < %f)' % (pb_name, pblimit_image, pb_name, pblimit_cleanmask) elif cleanmask is not None and os.path.exists(cleanmask): have_mask = True # Area outside clean mask statsmask = '"%s" < 0.1' % (os.path.basename(flattened_mask)) else: # Full area have_mask = False statsmask = '' try: # Get image RMS for all channels (this is for the weblog) # Avoid repeat if the check for npts was done and is OK. if image_stats is None: image_stats = image.statistics(mask=statsmask, robust=True, axes=[0, 1, 2], algorithm='chauvenet', maxiter=5) nonpbcor_image_statsmask = statsmask # Filter continuum frequency ranges if given if cont_freq_ranges not in (None, '', 'NONE', 'ALL'): # TODO: utils.freq_selection_to_channels uses casa_tools.image to get the frequency axis # and closes the global pipeline image tool. The context manager wrapped tool # used in this "with" statement is a different instance, so this is OK, but stacked # use of casa_tools.image might lead to unexpected results. cont_chan_ranges = utils.freq_selection_to_channels(nonpbcor_imagename, cont_freq_ranges) cont_chan_indices = np.hstack([np.arange(start, stop+1) for start, stop in cont_chan_ranges]) nonpbcor_image_non_cleanmask_rms_vs_chan = image_stats['rms'][cont_chan_indices] else: nonpbcor_image_non_cleanmask_rms_vs_chan = image_stats['rms'] nonpbcor_image_non_cleanmask_robust_rms = image_stats['medabsdevmed'] * 1.4826 nonpbcor_image_non_cleanmask_rms_median = np.median(nonpbcor_image_non_cleanmask_rms_vs_chan) nonpbcor_image_non_cleanmask_rms_mean = np.mean(nonpbcor_image_non_cleanmask_rms_vs_chan) nonpbcor_image_non_cleanmask_rms_min = np.min(nonpbcor_image_non_cleanmask_rms_vs_chan) nonpbcor_image_non_cleanmask_rms_max = np.max(nonpbcor_image_non_cleanmask_rms_vs_chan) nonpbcor_image_non_cleanmask_rms = nonpbcor_image_non_cleanmask_rms_median if have_mask: area_text = 'annulus' else: area_text = 'full image' LOG.info('Clean image statistics (%s) for %s: rmsmedian: %s Jy/bm rmsmean: %s Jy/bm rmsmin:' ' %s Jy/bm rmsmax: %s Jy/bm' % \ (area_text, os.path.basename(nonpbcor_imagename), nonpbcor_image_non_cleanmask_rms_median, nonpbcor_image_non_cleanmask_rms_mean, nonpbcor_image_non_cleanmask_rms_min, nonpbcor_image_non_cleanmask_rms_max)) except Exception as e: nonpbcor_image_non_cleanmask_rms_median, \ nonpbcor_image_non_cleanmask_rms_mean, \ nonpbcor_image_non_cleanmask_rms_min = \ nonpbcor_image_non_cleanmask_rms_max = \ nonpbcor_image_non_cleanmask_rms = \ -999.0 LOG.warn('Exception while determining image RMS for %s: %s' % (nonpbcor_imagename, e)) # Get the flux density spectrum in the clean mask area if available if nonpbcor_image_cleanmask_npoints not in (None, 0): # Area in flattened clean mask spectrum_mask = '"%s" > 0.1' % (os.path.basename(flattened_mask)) elif pb is not None: # Area of pb > pblimit_cleanmask spectrum_mask = '"%s" > %f' % (os.path.basename(pb)+extension, pblimit_cleanmask) nonpbcor_image_cleanmask_spectrum_pblimit = pblimit_cleanmask else: spectrum_mask = None nonpbcor_image_cleanmask_spectrum = image.getprofile(function='flux', mask=spectrum_mask, stretch=True, axis=freq_axis)['values'] return (residual_cleanmask_rms, residual_non_cleanmask_rms, residual_min, residual_max, nonpbcor_image_non_cleanmask_rms_min, nonpbcor_image_non_cleanmask_rms_max, nonpbcor_image_non_cleanmask_rms, pbcor_image_min, pbcor_image_max, residual_robust_rms, {'nonpbcor_imagename': nonpbcor_imagename, 'nonpbcor_image_non_cleanmask_freq_ch1': nonpbcor_image_non_cleanmask_freq_ch1, 'nonpbcor_image_non_cleanmask_freq_chN': nonpbcor_image_non_cleanmask_freq_chN, 'nonpbcor_image_non_cleanmask_freq_frame': nonpbcor_image_non_cleanmask_freq_frame, 'nonpbcor_image_non_cleanmask_robust_rms': nonpbcor_image_non_cleanmask_robust_rms, 'nonpbcor_image_cleanmask_spectrum': nonpbcor_image_cleanmask_spectrum, 'nonpbcor_image_cleanmask_spectrum_pblimit': nonpbcor_image_cleanmask_spectrum_pblimit, 'nonpbcor_image_cleanmask_npoints': nonpbcor_image_cleanmask_npoints, 'cont_freq_ranges': cont_freq_ranges, 'nonpbcor_image_statsmask': nonpbcor_image_statsmask})