Script for enforcing variable reliability cutoffs, and making clean masks

I've been applying SOURCERY to some ASKAP-12 data with the goal of automatically making some clean masks that avoid spurious features, both to aid deconvolution and to make more accurate sky models for self-cal.

If found that imposing a single reliability threshold was insufficient. Too low and artefacts around bright sources remained in the model, too high and many genuine faint sources were omitted, with the result being no better than a straightforward brightness cut.

To get around this I made the following script, which uses two reliability cut offs, depending on whether a component is within a certain angular distance of components above a certain brightness threshold. It requires a bit of tuning but it works very well.

Two other scrips which are useful in conjunction with this are this one (unity_fluxes.py) which makes all the source fluxes in a model unity, and this one (make_clean_mask.py) which takes a Tigger LSM and a reference image and outputs a clean mask based on the contents of the sky model (works best for unity fluxes, hence the first script).

import numpy
import Tigger
import glob
import os
from astLib import astCoords as ac

def r2d(x):
        return x*180.0/numpy.pi

def filter_lsm(mylsm,rel1,rel2,flux_cut,rad_cut):
        op_lsm = sourcery_lsm.replace('.lsm.html','_filtered.lsm.html')
        if os.path.isfile(op_lsm):
                print op_lsm,'already exists, skipping'
                op_lsm = False
        else:
                os.system('cp '+sourcery_lsm+' '+op_lsm)

                brights = []
                srcs = Tigger.load(mylsm).sources
                for src in srcs:
                        name = src.name
                        ra = r2d(src.pos.ra)
                        dec = r2d(src.pos.dec)
                        flux = src.flux.I
                        if flux > flux_cut:
                                brights.append((name,ra,dec,flux))

                model = Tigger.load(op_lsm)
                srcs = model.sources

                filtered_lsm = []

                for src in srcs:
                        near_bright = False

                        # Get source properties
                        name = src.name
                        ra = r2d(src.pos.ra)
                        dec = r2d(src.pos.dec)
                        flux = src.flux.I
                        rel = src.getTag('rel')

                        # Is it near a bright source?
                        for bright in brights:
                                b_ra = bright[1]
                                b_dec = bright[2]
                                sep = ac.calcAngSepDeg(ra,dec,b_ra,b_dec)
                                if sep < rad_cut:
                                        near_bright = True

                        # If it is then demand rel1
                        if near_bright and rel > rel1:
                                filtered_lsm.append(src)
                        # If it's not then relax and rel2 is fine
                        elif not near_bright and rel > rel2:
                                filtered_lsm.append(src)

                model.setSources(filtered_lsm)
                Tigger.save(model,op_lsm,verbose=True)
        return op_lsm

flux_cut = 0.065
rad_cut = 0.1
rel1 = 0.98
rel2 = 0.3

lsm_list = glob.glob('*_2016*/*image.lsm.html')

for sourcery_lsm in lsm_list:
        op_lsm = filter_lsm(sourcery_lsm,rel1,rel2,flux_cut,rad_cut)
        print op_lsm

radio-astro / sourcery

Script for enforcing variable reliability cutoffs, and making clean masks #15