From c7a66e8ccabb7a57c674ef6b78a0f336b71de643 Mon Sep 17 00:00:00 2001
From: AK <kutkin@gmail.com>
Date: Mon, 9 Aug 2021 13:35:42 +0000
Subject: [PATCH] fixed wrong mask file usage

---
 cluster.py | 20 +++++++++++---------
 imcal.py   | 27 ++++++++++++++-------------
 imcal.yml  | 14 +++++---------
 3 files changed, 30 insertions(+), 31 deletions(-)

diff --git a/cluster.py b/cluster.py
index 67aa4a1..a2c2cbc 100755
--- a/cluster.py
+++ b/cluster.py
@@ -378,12 +378,6 @@ def auto_clustering(fig, ax, df, wcs, resid_data, pix_arcmin_scale, nbright,
 
         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)
-        # print ellipse_mean, ecc, amaj, numpix
-
-
-
-        # if abs(ring_mean/resid_mean) > 1.7e5:
-        # if abs(ring_mean/small_resid_mean) > 1e4:
 
         if nclusters=='auto':
             if abs(ellipse_mean/rms) > 1.4:
@@ -451,12 +445,20 @@ def auto_clustering(fig, ax, df, wcs, resid_data, pix_arcmin_scale, nbright,
 
         return final_clusters
 
+# TODO
+def voronoy_clustering():
+    """
+    Use Voronoy clustering instead of fixed radius around sources
+    """
+    pass
+
+
 
 def main(img, resid, model, auto=True, add_manual=False, nclusters=5, boxsize=250,
          nbright=80, cluster_radius=5, cluster_overlap=1.6):
 
-    path = os.path.split(os.path.abspath(img))[0]
-    output = os.path.splitext(img)[0]+'-clustered.txt'
+    imgbase = os.path.splitext(img)[0]
+    output = imgbase + '-clustered.txt'
 
     df = pd.read_csv(model, skipinitialspace=True)
     df['ra'] = df.Ra.apply(ra2deg)
@@ -491,7 +493,7 @@ def main(img, resid, model, auto=True, add_manual=False, nclusters=5, boxsize=25
     if clusters:
         write_df(df, clusters, output=output)
     fig.tight_layout()
-    fig.savefig(path+'/clustering.png')
+    fig.savefig(imgbase+'-clustering.png')
     return output
 
 
diff --git a/imcal.py b/imcal.py
index f7e5236..cd6a1e4 100755
--- a/imcal.py
+++ b/imcal.py
@@ -326,13 +326,13 @@ def main(msin, outbase=None, cfgfile='imcal.yml'):
     img1 = outbase + '_1'
     img2 = outbase + '_2'
     img3 = outbase + '_3'
-    img_final = outbase
+    img_final = outbase + '-dical'
     img_ddsub = outbase + '-ddsub'
     img_ddcal = outbase + '-ddcal'
 
-    mask0 = 'mask0.fits'
-    mask1 = 'mask1.fits'
-    mask2 = 'mask2.fits'
+    mask0 = outbase + '-mask0.fits'
+    mask1 = outbase + '-mask1.fits'
+    mask2 = outbase + '-mask2.fits'
 
     model1 = outbase + '_model1.sourcedb'
     model2 = outbase + '_model2.sourcedb'
@@ -364,45 +364,46 @@ def main(msin, outbase=None, cfgfile='imcal.yml'):
         threshold = img_max/cfg['clean0']['max_over_thresh']
         wsclean(ms_split, outname=img0, automask=None, save_source_list=False, multifreq=False, mgain=None,
             kwstring=f'-threshold {threshold}')
-    if not os.path.exists('mask0.fits')    :
-        create_mask(img0 +'-image.fits', img0 +'-residual.fits', clipval=20, outname='mask0.fits')
+        create_mask(img0 +'-image.fits', img0 +'-residual.fits', clipval=20, outname=mask0)
 
 # clean1
     if not os.path.exists(img1 +'-image.fits') and (not os.path.exists(img1 +'-MFS-image.fits')):
-        wsclean(ms_split, outname=img1, **cfg['clean1']) # fast shallow clean
+        wsclean(ms_split, fitsmask=mask0, outname=img1, **cfg['clean1']) # fast shallow clean
 
     if not os.path.exists(model1):
         makesourcedb(img1+'-sources.txt', out=model1)
-
+# dical1
     if not os.path.exists(dical1):
         dical1 = dical(ms_split, model1, msout=dical1, h5out=h5_1, **cfg['dical1'])
         view_sols(h5_1, outname='dical1')
 # clean2
     if (not os.path.exists(img2 +'-image.fits')) and (not os.path.exists(img2 +'-MFS-image.fits')):
-        wsclean(dical1, outname=img2, **cfg['clean2'])
-        create_mask(img2 +'-image.fits', img2 +'-residual.fits', clipval=7, outname='mask1.fits')
+        wsclean(dical1, fitsmask=mask0, outname=img2, **cfg['clean2'])
+        create_mask(img2 +'-image.fits', img2 +'-residual.fits', clipval=7, outname=mask1)
 
     if not os.path.exists(model2):
         makesourcedb(img2+'-sources.txt', out=model2)
 
+# dical2
     if not os.path.exists(dical2):
         dical2 = dical(dical1, model2, msout=dical2, h5out=h5_2, **cfg['dical2'])
         view_sols(h5_2, outname='dical2')
 # clean3
     if (not os.path.exists(img3 +'-image.fits')) and (not os.path.exists(img3 +'-MFS-image.fits')):
-        wsclean(dical2, outname=img3, **cfg['clean3'])
-        create_mask(img3 +'-image.fits', img3 +'-residual.fits', clipval=5, outname='mask2.fits')
+        wsclean(dical2, fitsmask=mask1, outname=img3, **cfg['clean3'])
+        create_mask(img3 +'-image.fits', img3 +'-residual.fits', clipval=5, outname=mask2)
 
     if not os.path.exists(model3):
         makesourcedb(img3+'-sources.txt', out=model3)
 
+# dical3
     if not os.path.exists(dical3):
         dical3 = dical(dical2, model3, msout=dical3, h5out=h5_3, **cfg['dical3'])
         view_sols(h5_3, outname='dical3')
 
 # clean4
     if (not os.path.exists(img_final +'-image.fits')) and (not os.path.exists(img_final +'-MFS-image.fits')):
-        wsclean(dical3, outname=img_final, **cfg['clean4'])
+        wsclean(dical3, fitsmask=mask2, outname=img_final, **cfg['clean4'])
 
 
 # Cluster
diff --git a/imcal.yml b/imcal.yml
index 6fc2213..997fc12 100644
--- a/imcal.yml
+++ b/imcal.yml
@@ -24,7 +24,6 @@ clean1: # wsclean setup
     mgain: 0
     automask: 20
     autothresh: 5
-    fitsmask: 'mask0.fits'
     multiscale: False
     kwstring: '-use-wgridder -parallel-deconvolution 1400' # use this for additional wsclean options, e.g. '-weight uniform -use-idg' 
 
@@ -32,7 +31,7 @@ dical1: # DPPP setup for direction independent calibration
     solint: 20
     mode: 'phaseonly'
     uvlambdamin: 500 # Ignore baselines / channels with UV < uvlambdamin wavelengths.
-#    cal_nchan: 31 # number of chans with the same solutions
+    cal_nchan: 30 # number of chans with the same solutions
 
 clean2:
     imagesize: 3072
@@ -40,7 +39,6 @@ clean2:
     multifreq: 8  
     automask: 10
     autothresh: 5
-    fitsmask: 'mask0.fits'
     multiscale: True
     kwstring: '-use-wgridder -parallel-deconvolution 1400 -parallel-gridding 8 -deconvolution-channels 3'
     
@@ -48,7 +46,7 @@ dical2:
     solint: 1
     mode: 'phaseonly'
     uvlambdamin: 500 # Ignore baselines / channels with UV < uvlambdamin wavelengths.
-    cal_nchan: 31 # number of chans with the same solutions
+    cal_nchan: 30 # number of chans with the same solutions
 
 clean3:
     imagesize: 3072
@@ -56,7 +54,6 @@ clean3:
     multifreq: 8  
     automask: 7
     autothresh: 3.5
-    fitsmask: 'mask1.fits'
     multiscale: True
     kwstring: '-use-wgridder -parallel-deconvolution 1400 -parallel-gridding 8 -deconvolution-channels 3'
     
@@ -64,7 +61,7 @@ dical3:
     solint: 800
     mode: 'diagonal'
     uvlambdamin: 500 # Ignore baselines / channels with UV < uvlambdamin wavelengths.
-    cal_nchan: 31 # number of chans with the same solutions
+    cal_nchan: 30 # number of chans with the same solutions.
     
 clean4:
     imagesize: 3072
@@ -73,7 +70,6 @@ clean4:
     automask: 4.5
     autothresh: 0.5
     multiscale: True
-    fitsmask: 'mask2.fits'
     kwstring: '-use-wgridder -parallel-deconvolution 1400 -parallel-gridding 8 -deconvolution-channels 3 -weight briggs 0.0'
 
 ####################### CLUSTERING #######################
@@ -88,9 +84,9 @@ cluster:
 
 ######################## DD CALIBRATION #######################
 ddcal: # see DPPP/DDECal documentation
-    solint: 120 # Solution interval in timesteps (1 ~ 30sec for Apertif).
+    solint: 120 # Solution interval in timesteps (1 corresponds to 30 seconds for Apertif).
     mode: 'diagonal' # Type of constraint to apply. 
-    nfreq: 15 # Number of channels in each channel block, for which the solution is assumed to be constant.
+    nfreq: 30 # Number of channels in each channel block, for which the solution is assumed to be constant.
     startchan: 0
     nchan: 0
     uvlambdamin: 500 # Ignore baselines / channels with UV < uvlambdamin wavelengths.
-- 
GitLab