some fixes, add preflag and facet steps

Dockerfile

@@ -185,7 +185,6 @@ FROM ubuntu:20.04 as runner
 RUN mkdir /src
 COPY --from=builder /usr/local /usr/local
 RUN chmod +rx /usr/local/bin/*
-COPY --from=builder /src/wsclean /opt/
 
 SHELL ["/bin/bash", "-c"]
 

cluster.py

@@ -605,21 +605,12 @@ def voronoi_clustering(fig, ax, df, wcs, resid_data, nbright, nclusters,
 
     # logging.info('Checking {} brightest model components'.format(nbright))
 
-    bright_df = df.sort_values('I')[::-1][:nbright][['ra', 'dec', 'I']]
-    # csnrs = []
-    # cfluxes = []
-    # cmeasures = [] # the main value to clusterize by
-    # cellipse_params = []
-
-    # fmin, fmax = min(a.I), max(a.I)
-    rms = mad(resid_data)
-    # resid_mean = np.mean(resid_data)
-
+    bright_df = df.sort_values('I')[::-1][['ra', 'dec', 'I']]
+    # rms = mad(resid_data)
     logging.info('Getting measures for the potential clusters...')
     clusters = []
     clusters_centers = [] #
 
-
     for ra, dec, flux in bright_df.values:
         c = SkyCoord(ra, dec, unit='deg')
         px, py = np.round(wcs.all_world2pix(c.ra, c.dec, 0)).astype(int)
@@ -627,17 +618,15 @@ def voronoi_clustering(fig, ax, df, wcs, resid_data, nbright, nclusters,
         if (abs(px-resid_data.shape[1]) < boxsize) or (abs(py-resid_data.shape[0]) < boxsize):
             # logging.debug('Skipping the edge source')
             continue
-# Check if the component already in a cluster
+# Check if the component is nearby
         if clusters and any([c.separation(_).arcmin<same_source_radius for _ in clusters]):
             continue
 
         small_resid = resid_data[py-boxsize:py+boxsize, px-boxsize:px+boxsize]
         ellipse_mean, ecc, amaj, numpix = ellipses_coh(small_resid, amin=20, amax=boxsize-1, dr=1.0)
 
-        if abs(ellipse_mean/rms) > 1.4:
-             # rect = plt.Rectangle((px-boxsize, py-boxsize), 2*boxsize, 2*boxsize,
-             #                      lw=1, color='k', fc='none', alpha=0.3)
-             # ax.add_artist(rect)
+# Uncomment the condition to search for artifacts:
+        # if abs(ellipse_mean/rms) > 1.4:
         clusters_centers.append([ra, dec])
         clusters.append(c)
         print(ra, dec)
@@ -650,12 +639,6 @@ def voronoi_clustering(fig, ax, df, wcs, resid_data, nbright, nclusters,
         logging.warning('Decreasing number of clusters')
         nclusters = len(clusters_centers)
 
-    # X = np.vstack([df.ra, df.dec]).T
-    # kmeans = KMeans(n_clusters=nclusters, init=clusters_centers, n_init=1, random_state=0, max_iter=1)
-    # kmeans.fit(X)
-    # y_kmeans = kmeans.predict(X) # cluster index for each observation
-    # ax.scatter(clusters_centers[:, 0], clusters_centers[:, 1], c='black', s=100, alpha=0.5, transform=ax.get_transform("world"))
-
     vor = Voronoi(np.array(clusters_centers))
     voronoi_plot_2d_world(vor, ax=ax, show_vertices=False)
 

imcal.py

@@ -40,7 +40,6 @@ _MAX_TIME = 1 * 3600 # SECONDS
 _MAX_POOL = _MAX_TIME // _POOL_TIME
 
 def execute_binary(binary, args, ):
-    
     command = [f'{binary}'] + args
     logging.debug('executing %s', ','.join(command))
     dppp_process = subprocess.Popen(command)
@@ -65,6 +64,7 @@ def fft_psf(bmaj, bmin, bpa, size=3073):
     fkern = fkern.array
     return fkern
 
+
 def reconvolve_gaussian_kernel(img, old_maj, old_min, old_pa, new_maj, new_min, new_pa):
     """
     convolve image with a gaussian kernel without FFTing it
@@ -178,12 +178,13 @@ def smoothImage(imgfits, psf=30, out=None) :
     return fits_reconvolve_psf(imgfits, Beam(psf*u.arcsec), out=out)
 
 
-def create_mask(imgfits, residfits, clipval, outname='mask.fits', ):
+def create_mask(imgfits, residfits, clipval, outname='mask.fits'):
     """
     Create mask using Tom's code (e-mail on 1 Jul 2021)
     """
-    cmd1 = f'makeNoiseMapFitsLow {imgfits} {residfits} noise.fits noiseMap.fits'
-    cmd2 = f'makeMaskFits noiseMap.fits {outname} {clipval}'
+    outbase = os.path.splitext(imgfits)[0]
+    cmd1 = f'makeNoiseMapFitsLow {imgfits} {residfits} {outbase}_noise.fits {outbase}_noiseMap.fits'
+    cmd2 = f'makeMaskFits {outbase}_noiseMap.fits {outname} {clipval}'
     logging.debug("Running command: %s", cmd1)
     subprocess.call(cmd1, shell=True)
     logging.debug("Running command: %s", cmd2)
@@ -191,42 +192,50 @@ def create_mask(imgfits, residfits, clipval, outname='mask.fits', ):
     return outname
 
 
-def makeNoiseImage(imgfits, residfits, low=False, ) :
+def makeNoiseImage(imgfits, residfits, low=False) :
     """
     Create mask using Tom's code (e-mail on 1 Jul 2021)
     """
+    # if outbase is None:
+    outbase = os.path.splitext(imgfits)[0]
     if low:
-      cmd = f'makeNoiseMapFitsLow {imgfits} {residfits} noiseLow.fits noiseMapLow.fits'
+        img1, img2 = f'{outbase}_noiseLow.fits', f'{outbase}_noiseMapLow.fits'
+        cmd = f'makeNoiseMapFitsLow {imgfits} {residfits} {img1} {img2}'
     else :
-      cmd = f'makeNoiseMapFits {imgfits} {residfits} noise.fits noiseMap.fits'
-      
+        img1, img2 = f'{outbase}_noise.fits', f'{outbase}_noiseMap.fits'
+        cmd = f'makeNoiseMapFits {imgfits} {residfits} {img1} {img2}'
     logging.debug("Running command: %s", cmd)
     subprocess.call(cmd, shell=True)
-    return
+    return img2
 
 
-def makeCombMask(ima1='noiseMap.fits', ima2='noiseMapLow.fits',
-		 clip1=5, clip2=7, outname='mask.fits', ) :
+def makeCombMask(img1, img2, clip1=5, clip2=7, outname=None) :
     """
     Create mask using Tom's code (e-mail on 1 Jul 2021)
     """
-    cmd = f'makeCombMaskFits {ima1} {ima2} {outname} {clip1} {clip2}'
-    
+    if outname is None:
+        outname = os.path.splitext(img1)[0] + '_mask.fits'
+    cmd = f'makeCombMaskFits {img1} {img2} {outname} {clip1} {clip2}'
     logging.debug("Running command: %s", cmd)
     subprocess.call(cmd, shell=True)
-    return
+    return outname
 
 
 def get_image_ra_dec_min_max(msin):
     """
     Determine image center coords, min and max values for msin
     """
-    cmd = f'wsclean -niter 0 -size 3072 3072 -scale 3arcsec -use-wgridder {msin}'
+    outbase = os.path.splitext(msin.rstrip('/'))[0]+'-wsclean'
+    outname = outbase+'-image.fits'
+    cmd = f'wsclean -name {outbase} -niter 0 -size 3072 3072 -scale 3arcsec -use-wgridder {msin}'
+    if os.path.exists(outname):
+        logging.debug('Image exists. Skipping cleaning...')
+    else:
         logging.debug('Running command: %s', cmd)
         subprocess.call(cmd, shell=True)
-    data = fits.getdata('wsclean-image.fits')
-    header = fits.getheader('wsclean-image.fits')
-    return header['CRVAL1'], header['CRVAL2'], np.nanmin(data), np.nanmax(data)
+    data = fits.getdata(outname)
+    header = fits.getheader(outname)
+    return outname, header['CRVAL1'], header['CRVAL2'], np.nanmin(data), np.nanmax(data)
 
 
 def makesourcedb(modelfile, out=None, ):
@@ -255,9 +264,7 @@ def render(bkgr, model, out=None, ):
     return out
 
 
-        
 def execute_dppp(args, ):
-    
     command = ['DP3'] + args
     logging.debug('executing %s', ','.join(command))
     dppp_process = subprocess.Popen(command)
@@ -298,6 +305,25 @@ def split_ms(msin_path, startchan, nchan=0, msout_path='', ):
     return msout_path
 
 
+def preflag(msin, msout=None, **kwargs):
+    """
+    preflag data using DP3 preflag module
+    """
+    if not any(kwargs.values()):
+        logging.debug('No preflag options specified. Skipping...')
+        return msin
+    msout = msout or '.'
+    command_args = ['steps=[preflag]',
+                    f'msin={msin}',
+                    f'msout={msout}'] + ['preflag.'+'='.join(_) for _ in kwargs.items() if _[1] is not None]
+    logging.info('Flagging data (%s)', command_args)
+    return_code = execute_dppp(command_args)
+    logging.debug('Preflag of %s returned status code %s', msin, return_code)
+    check_return_code(return_code)
+    if msout == '.': msout = msin
+    return msout
+
+
 def dical(msin, srcdb, msout=None, h5out=None, solint=1, startchan=0, split_nchan=0,
           mode='phaseonly', cal_nchan=0, uvlambdamin=500, ):
     """ direction independent calibration with DPPP """
@@ -467,7 +493,7 @@ def main(msin, steps='all', outbase=None, cfgfile='imcal.yml', force=False):
 
 
     if steps == 'all':
-        steps = ['nvss', 'mask', 'dical', 'ddcal']
+        steps = ['nvss', 'preflag', 'mask', 'dical', 'ddcal']
     else:
         steps = steps.split(',')
 
@@ -480,6 +506,7 @@ def main(msin, steps='all', outbase=None, cfgfile='imcal.yml', force=False):
 
     ms_split = msbase + '_splt.MS'
 
+
     img0 = outbase + '_0'
     img1 = outbase + '_1'
     img2 = outbase + '_2'
@@ -487,8 +514,8 @@ def main(msin, steps='all', outbase=None, cfgfile='imcal.yml', force=False):
     img_dical = outbase + '-dical'
     img_ddsub_1 = outbase + '-ddsub-1'
     img_ddsub_2 = outbase + '-ddsub-2'
-    img_ddcal_1 = outbase + '-ddcal'
-    img_ddcal_2 = outbase + '-ddcal'
+    img_ddcal_1 = outbase + '-ddcal-1'
+    img_ddcal_2 = outbase + '-ddcal-2'
 
     mask0 = outbase + '-mask0.fits'
     mask1 = outbase + '-mask1.fits'
@@ -520,10 +547,15 @@ def main(msin, steps='all', outbase=None, cfgfile='imcal.yml', force=False):
         return 0
 
 # get image parameters
-    img_ra, img_dec, img_min, img_max = get_image_ra_dec_min_max(msin)
+    # if 'init' in steps:
+    # if not force and os.path.exists(outbase+'-wsclean-image.fits'):
+    initial_img, img_ra, img_dec, img_min, img_max = get_image_ra_dec_min_max(msin)
+    logging.info('Image: %s', initial_img)
+    logging.info('Image RA, DEC: %s, %s', img_ra, img_dec)
+    logging.info('Image Min, Max: %s, %s', img_min, img_max)
 
     if 'nvss' in steps and cfg['nvss']:
-       nvss_model = nvss_cutout('wsclean-image.fits', nvsscat='/opt/nvss.csv.zip', clip=cfg['nvsscal']['clip_model'])
+       nvss_model = nvss_cutout(initial_img, nvsscat='/opt/nvss.csv.zip', clip=cfg['nvsscal']['clip_model'])
        if (not os.path.exists(ms_split)) and (cfg['split']['startchan'] or cfg['split']['nchan']):
            ms_split = split_ms(msin, msout_path=ms_split, **cfg['split'])
 
@@ -532,8 +564,13 @@ def main(msin, steps='all', outbase=None, cfgfile='imcal.yml', force=False):
        dical(ms_split, nvssMod, msout=dical0, h5out=h5_0, **cfg['nvsscal'])
        view_sols(h5_0, outname=msbase+'_sols_dical0')
     else :
-       # if (not os.path.exists(ms_split)) and (cfg['split']['startchan'] or cfg['split']['nchan']):
-       ms_split = split_ms(msin, msout_path=dical0, **cfg['split'])
+       if (not os.path.exists(ms_split)) and (cfg['split']['startchan'] or cfg['split']['nchan']):
+           dical0 = split_ms(msin, msout_path=dical0, **cfg['split'])
+       else:
+           dical0 = msin
+
+    if 'preflag' in steps:
+        dical0 = preflag(ms_split, msout=dical0, **cfg['preflag'])
 
     if 'mask' in steps:
         if not force and (os.path.exists(img0 +'-image.fits') or (os.path.exists(img0 +'-MFS-image.fits'))):
@@ -565,9 +602,9 @@ def main(msin, steps='all', outbase=None, cfgfile='imcal.yml', force=False):
         else:
             wsclean(dical1, fitsmask=mask0, outname=img2, **cfg['clean2'])
             smoothImage(img2+'-residual.fits')
-            makeNoiseImage(img2 +'-image.fits', img2 +'-residual.fits', )
-            makeNoiseImage(img2 +'-residual-smooth.fits', img2 +'-residual.fits', low=True, )
-            makeCombMask(outname=mask1, clip1=7, clip2=15, )
+            i1 = makeNoiseImage(img2 +'-image.fits', img2 +'-residual.fits', )
+            i2 = makeNoiseImage(img2 +'-residual-smooth.fits', img2 +'-residual.fits', low=True, )
+            makeCombMask(i1, i2, clip1=7, clip2=15, outname=mask1, )
 
             makesourcedb(img2+'-sources.txt', out=model2, )
 
@@ -583,9 +620,9 @@ def main(msin, steps='all', outbase=None, cfgfile='imcal.yml', force=False):
         else:
             wsclean(dical2, fitsmask=mask1, outname=img3, **cfg['clean3'])
             smoothImage(img3+'-residual.fits')
-            makeNoiseImage(img3 +'-image.fits', img3 +'-residual.fits', )
-            makeNoiseImage(img3 +'-residual-smooth.fits', img3 +'-residual.fits', low=True, )
-            makeCombMask(outname=mask2, clip1=5, clip2=10)
+            i1 = makeNoiseImage(img3 +'-image.fits', img3 +'-residual.fits', )
+            i2 = makeNoiseImage(img3 +'-residual-smooth.fits', img3 +'-residual.fits', low=True, )
+            makeCombMask(i1, i2, clip1=5, clip2=10, outname=mask2,)
 
             makesourcedb(img3+'-sources.txt', out=model3)
 
@@ -602,9 +639,9 @@ def main(msin, steps='all', outbase=None, cfgfile='imcal.yml', force=False):
         else:
             wsclean(dical3, fitsmask=mask2, outname=img_dical,  **cfg['clean4'])
             smoothImage(img_dical+'-residual.fits')
-            makeNoiseImage(img_dical +'-image.fits', img_dical +'-residual.fits', )
-            makeNoiseImage(img_dical +'-residual-smooth.fits', img_dical +'-residual.fits',low=True, )
-            makeCombMask(outname=mask3, clip1=5, clip2=7, )
+            i1 = makeNoiseImage(img_dical +'-image.fits', img_dical +'-residual.fits', )
+            i2 = makeNoiseImage(img_dical +'-residual-smooth.fits', img_dical +'-residual.fits',low=True, )
+            makeCombMask(i1, i2, clip1=5, clip2=7, outname=mask3)
 
     if 'ddcal' in steps:
 # Cluster
@@ -634,10 +671,10 @@ def main(msin, steps='all', outbase=None, cfgfile='imcal.yml', force=False):
 
         render(img_ddsub_1+'-image.fits', aomodel, out=img_ddcal_1+'-image.fits')
 
-        smoothImage(img_ddsub_1+'-image.fits')
-        makeNoiseImage(img_ddcal_1 +'-image.fits', img_ddsub_1 +'-residual.fits', )
-        makeNoiseImage(img_ddcal_1 +'-image-smooth.fits', img_ddsub_1 +'-residual.fits',low=True, )
-        makeCombMask(outname=mask4, clip1=3.5, clip2=5, )
+        smoothImage(img_ddcal_1+'-image.fits')
+        i1 = makeNoiseImage(img_ddcal_1 +'-image.fits', img_ddsub_1 +'-residual.fits', )
+        i2 = makeNoiseImage(img_ddcal_1 +'-image-smooth.fits', img_ddsub_1 +'-residual.fits',low=True, )
+        makeCombMask(i1, i2, clip1=3.5, clip2=5, outname=mask4,)
 
         if not force and os.path.exists(img_ddsub_2+'-image.fits'):
             pass
@@ -649,13 +686,14 @@ def main(msin, steps='all', outbase=None, cfgfile='imcal.yml', force=False):
         render(img_ddsub_2+'-image.fits', aomodel, out=img_ddcal_2+'-image.fits', )
 
 # test facet imaging:
-    # ds9_file = 'ddfacets.reg'
-    # ddvis = outbase + '_ddvis.MS'
-    # h5_ddvis = 'ddsols.h5'
-    # clustered_sdb = img_dical+'-clustered.sourcedb'
-    # # ddvis = ddecal(dical3, clustered_sdb, msout=ddvis, subtract=False, h5out=h5_ddvis, **cfg['ddcal'])
-    # # write_ds9(ds9_file, h5_ddvis, img_ddcal+'-image.fits')
-    # wsclean(ddvis, fitsmask=mask3, save_source_list=False, outname='img-facet', **cfg['clean6'],)
+    if 'facet' in steps:
+        ds9_file = 'ddfacets.reg'
+        ddvis = outbase + '_ddvis.MS'
+        h5_ddvis = 'ddsols.h5'
+        clustered_sdb = img_dical+'-clustered.sourcedb'
+        # ddvis = ddecal(dical3, clustered_sdb, msout=ddvis, subtract=False, h5out=h5_ddvis, **cfg['ddcal'])
+        # write_ds9(ds9_file, h5_ddvis, img_ddcal_2+'-image.fits')
+        wsclean(ddvis, fitsmask=mask3, save_source_list=False, outname='img-facet', **cfg['facet_clean'],)
 
 
     return 0
@@ -673,13 +711,13 @@ if __name__ == "__main__":
                         dest='configfile', help='Config file', type=str)
     parser.add_argument('-o', '--outbase', default=None, help='output prefix', type=str)
     parser.add_argument('-s', '--steps', default='all', help='steps to run. Example: "nvss,mask,dical,ddcal"', type=str)
-    parser.add_argument('-f', '--force', action='store_false', help='Overwrite the existing files')
+    parser.add_argument('-f', '--force', action='store_true', help='Overwrite the existing files')
 
     args = parser.parse_args()
     configfile = args.configfile or \
         os.path.join(os.path.dirname(os.path.realpath(__file__)), 'imcal.yml')
     # msin = args.msin
-    main(args.msin, outbase=args.outbase, steps=args.steps, cfgfile=configfile)
+    main(args.msin, outbase=args.outbase, steps=args.steps, cfgfile=configfile, force=args.force)
     # extime = Time.now() - t0
     # print("Execution time: {:.1f} min".format(extime.to("minute").value))
 

imcal.yml

@@ -23,7 +23,14 @@ clean0: # initial clean
     max_over_thresh: 250 # the threshold for initial CLEAN is set to image_max/max_over_thresh
 # TODO add fixed threshold of 100 uJy
 
-############################## SPLIT ##########################################
+############################## PREFLAG ##########################################
+
+preflag: # See DP3 steps.preflag
+    baseline: # e.g. '[[RTB,*]]'
+    reltime: # e.g. '10:20+-20m'
+    chan: # null
+    blmin: # in meters
+
 clean1: # wsclean setup
     imagesize: 3072
     pixelsize: 3
@@ -86,7 +93,7 @@ clean4:
 
 ####################### CLUSTERING #######################
 cluster:
-    nbright: 80 # number of brightest clean components (CC) to check for artefacts
+    nbright: 100 # number of brightest clean components (CC) to check for artefacts
     boxsize: 250 # the boxsize around CC in pixels where to check for artefacts
     nclusters: 10 # number of clusters ('auto' -- to set automatically)
 # the following is only for 'auto' and 'manual' mathods:
@@ -114,16 +121,16 @@ clean5:
     clip_model_level: null
     kwstring: '-use-wgridder -parallel-deconvolution 1400 -parallel-gridding 8 -deconvolution-channels 3 -weight briggs -1.5'
 
-
-# clean6:
-#     imagesize: 3072
-#     pixelsize: 3
-#     multifreq: 6
-#     automask: 3.5
-#     autothresh: 0.5
-#     multiscale: True
-#     clip_model_level: null
-#     kwstring: '-facet-regions ddfacets.reg -apply-facet-solutions ddsols.h5 amplitude000,phase000 -use-wgridder -parallel-deconvolution 1400 -parallel-gridding 8 -deconvolution-channels 3 -weight briggs -1.5'
+### test for facet clean
+facet_clean:
+    imagesize: 3072
+    pixelsize: 3
+    multifreq: 6
+    automask: 3.5
+    autothresh: 0.5
+    multiscale: True
+    clip_model_level: null
+    kwstring: '-facet-regions ddfacets.reg -apply-facet-solutions ddsols.h5 amplitude000,phase000 -use-wgridder -parallel-deconvolution 1400 -parallel-gridding 8 -deconvolution-channels 3 -weight briggs -1.5'
 
 
 ### END