Some new Tom's code

cluster_to_reg.py

+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+from __future__ import division, absolute_import, print_function
+
+
+
+
+import numpy as np
+from scipy.spatial import Voronoi
+
+import polygon as Polygon
+
+def test():
+
+    RadiusTot=1.
+    Poly=np.array([[-RadiusTot,-RadiusTot],
+                   [+RadiusTot,-RadiusTot],
+                   [+RadiusTot,+RadiusTot],
+                   [-RadiusTot,+RadiusTot]])*1
+
+    Line=np.array([[0.,0.],
+                   [5.,5.]])
+
+    print((CutLineInside(Poly,Line)))
+
+def GiveABLin(P0,P1):
+    x0,y0=P0
+    x1,y1=P1
+    a=(y1-y0)/(x1-x0)
+    b=y1-a*x1
+    return b,a
+
+def GiveB(P0,P1):
+    x0,y0=P0
+    x1,y1=P1
+
+    B=np.array([x0-x1,y0-y1])
+    B/=np.sqrt(np.sum(B**2))
+    if B[0]<0: B=-B
+    return B
+
+def CutLineInside(Poly,Line):
+    P=Polygon.Polygon(Poly)
+
+    dx=1e-4
+    PLine=np.array(Line.tolist()+Line.tolist()[::-1]).reshape((4,2))
+    #PLine[2,0]+=dx
+    #PLine[3,0]+=2*dx
+    PLine[2:,:]+=np.random.randn(2,2)*1e-6
+    P0=Polygon.Polygon(PLine)
+    PP=np.array(P0&P)[0].tolist()
+    PP.append(PP[0])
+
+    B0=GiveB(Line[0],Line[1])
+    B=[GiveB(PP[i],PP[i+1]) for i in range(len(PP)-1)]
+
+    PLine=[]
+    for iB in range(len(B)):
+        d=np.sum((B[iB]-B0)**2)
+        print((d,PP[iB],PP[iB+1]))
+        if d==0:
+            PLine.append([PP[iB],PP[iB+1]])
+
+
+    LOut=np.array(PLine[0])
+
+    return LOut
+
+
+def voronoi_finite_polygons_2d(vor, radius=None):
+    """
+    Reconstruct infinite voronoi regions in a 2D diagram to finite
+    regions.
+    Parameters
+    ----------
+    vor : Voronoi
+        Input diagram
+    radius : float, optional
+        Distance to 'points at infinity'.
+    Returns
+    -------
+    regions : list of tuples
+        Indices of vertices in each revised Voronoi regions.
+    vertices : list of tuples
+        Coordinates for revised Voronoi vertices. Same as coordinates
+        of input vertices, with 'points at infinity' appended to the
+        end.
+    """
+
+    if vor.points.shape[1] != 2:
+        raise ValueError("Requires 2D input")
+
+    new_regions = []
+    new_vertices = vor.vertices.tolist()
+
+    center = vor.points.mean(axis=0)
+    if radius is None:
+        radius = vor.points.ptp().max()
+
+    # Construct a map containing all ridges for a given point
+    all_ridges = {}
+    for (p1, p2), (v1, v2) in zip(vor.ridge_points, vor.ridge_vertices):
+        all_ridges.setdefault(p1, []).append((p2, v1, v2))
+        all_ridges.setdefault(p2, []).append((p1, v1, v2))
+
+    # Reconstruct infinite regions
+    for p1, region in enumerate(vor.point_region):
+        vertices = vor.regions[region]
+
+        if all(v >= 0 for v in vertices):
+            # finite region
+            new_regions.append(vertices)
+            continue
+
+
+        # reconstruct a non-finite region
+        if p1 not in list(all_ridges.keys()):
+            new_regions.append(vertices)
+            continue
+        ridges = all_ridges[p1]
+        new_region = [v for v in vertices if v >= 0]
+
+        for p2, v1, v2 in ridges:
+            if v2 < 0:
+                v1, v2 = v2, v1
+            if v1 >= 0:
+                # finite ridge: already in the region
+                continue
+
+            # Compute the missing endpoint of an infinite ridge
+
+            t = vor.points[p2] - vor.points[p1] # tangent
+            t /= np.linalg.norm(t)
+            n = np.array([-t[1], t[0]])  # normal
+
+            midpoint = vor.points[[p1, p2]].mean(axis=0)
+            direction = np.sign(np.dot(midpoint - center, n)) * n
+
+            V=vor.vertices[v2]
+            R=np.sqrt(np.sum(V**2))
+
+            # while R>0.1:
+            #     V*=0.9
+            #     R=np.sqrt(np.sum(V**2))
+            #     #print R
+
+            vor.vertices[v2][:]=V[:]
+
+            ThisRad=radius
+            far_point = vor.vertices[v2] + direction * radius
+            R=np.sqrt(np.sum(far_point**2))
+            ThisRad=R
+
+            # while R>1:
+            #     ThisRad*=0.9
+            #     far_point = vor.vertices[v2] + direction * ThisRad
+            #     R=np.sqrt(np.sum(far_point**2))
+            #     print "=============="
+            #     print R,np.sqrt(np.sum(vor.vertices[v2]**2))
+            #     print vor.vertices[v2]
+            #     print direction
+            #     print ThisRad
+            #     #if R>1000:
+            #     #    stop
+
+            # RadiusTot=.3
+            # Poly=np.array([[-RadiusTot,-RadiusTot],
+            #                [+RadiusTot,-RadiusTot],
+            #                [+RadiusTot,+RadiusTot],
+            #                [-RadiusTot,+RadiusTot]])*1
+            # stop
+            # PT.CutLineInside(Poly,Line)
+
+
+            new_region.append(len(new_vertices))
+            new_vertices.append(far_point.tolist())
+
+        # sort region counterclockwise
+        vs = np.asarray([new_vertices[v] for v in new_region])
+        c = vs.mean(axis=0)
+        angles = np.arctan2(vs[:,1] - c[1], vs[:,0] - c[0])
+        new_region = np.array(new_region)[np.argsort(angles)]
+
+        # finish
+        new_regions.append(new_region.tolist())
+
+    regions, vertices=new_regions, np.asarray(new_vertices)
+
+    return regions, vertices
+
+
+HEADER = '\033[95m'
+OKBLUE = '\033[94m'
+OKGREEN = '\033[92m'
+WARNING = '\033[93m'
+FAIL = '\033[91m'
+ENDC = '\033[0m'
+bold='\033[1m'
+nobold='\033[0m'
+Separator="================================%s=================================="
+silent=0
+
+def Str(strin0,col="red",Bold=True):
+    if silent==1: return strin0
+    strin=str(strin0)
+    if col=="red":
+        ss=FAIL
+    if col=="green":
+        ss=OKGREEN
+    elif col=="yellow":
+        ss=WARNING
+    elif col=="blue":
+        ss=OKBLUE
+    elif col=="green":
+        ss=OKGREEN
+    elif col=="white":
+        ss=""
+    ss="%s%s%s"%(ss,strin,ENDC)
+    if Bold: ss="%s%s%s"%(bold,ss,nobold)
+    return ss
+
+def Sep(strin=None,D=1):
+    if D!=1:
+        return Str(Separator%("="*len(strin)))
+    else:
+        return Str(Separator%(strin))
+
+def Title(strin,Big=False):
+    print()
+    print()
+    if Big: print((Sep(strin,D=0)))
+    print((Sep(strin)))
+    if Big: print((Sep(strin,D=0)))
+    print()
+
+def disable():
+    HEADER = ''
+    OKBLUE = ''
+    OKGREEN = ''
+    WARNING = ''
+    FAIL = ''
+    ENDC = ''
+
+
+class ClassCoordConv():
+    def __init__(self,rac,decc):
+
+        rarad=rac
+        decrad=decc
+        self.rarad=rarad
+        self.decrad=decrad
+        cos=np.cos
+        sin=np.sin
+        mrot=np.array([[cos(rarad)*cos(decrad), sin(rarad)*cos(decrad),sin(decrad)],[-sin(rarad),cos(rarad),0.],[-cos(rarad)*sin(decrad),-sin(rarad)*sin(decrad),cos(decrad)]]).T
+        vm=np.array([[0.,0.,1.]]).T
+        vl=np.array([[0.,1., 0.]]).T
+        vn=np.array([[1., 0, 0.]]).T
+        self.vl2=np.dot(mrot,vl)
+        self.vm2=np.dot(mrot,vm)
+        self.vn2=np.dot(mrot,vn)
+        self.R=np.array([[cos(decrad)*cos(rarad),cos(decrad)*sin(rarad),sin(decrad)]]).T
+
+
+    def lm2radec(self,l_list,m_list):
+
+        ra_list=np.zeros(l_list.shape,dtype=np.float)
+        dec_list=np.zeros(l_list.shape,dtype=np.float)
+
+        for i in range(l_list.shape[0]):
+            l=l_list[i]
+            m=m_list[i]
+            if (l_list[i]==0.)&(m_list[i]==0.):
+                ra_list[i]=self.rarad
+                dec_list[i]=self.decrad
+                continue
+            Rp=self.R+self.vl2*l+self.vm2*m-(1.-np.sqrt(1.-l**2-m**2))*self.vn2
+            dec_list[i]=np.arcsin(Rp[2])
+            ra_list[i]=np.arctan(Rp[1]/Rp[0])
+            if Rp[0]<0.: ra_list[i]+=np.pi
+
+        return ra_list,dec_list
+
+    def radec2lm(self,ra,dec):
+        l = np.cos(dec) * np.sin(ra - self.rarad)
+        m = np.sin(dec) * np.cos(self.decrad) - np.cos(dec) * np.sin(self.decrad) * np.cos(ra - self.rarad)
+        return l,m
+
+
+
+
+
+class VoronoiToReg():
+    def __init__(self,rac,decc):
+        self.rac=rac
+        self.decc=decc
+        self.CoordMachine=ClassCoordConv(rac,decc)
+
+    def ToReg(self,regFile,xc,yc,radius=0.1,Col="red"):
+        print("Writing voronoi in: %s"%Str(regFile,col="blue"))
+        f=open(regFile,"w")
+        f.write("# Region file format: DS9 version 4.1\n")
+        ss0='global color=green dashlist=8 3 width=1 font="helvetica 10 normal roman" select=1 highlite=1 dash=0'
+        ss1=' fixed=0 edit=1 move=1 delete=1 include=1 source=1\n'
+
+        f.write(ss0+ss1)
+        f.write("fk5\n")
+
+        CoordMachine=self.CoordMachine
+
+        xy=np.zeros((xc.size,2),np.float32)
+        xy[:,0]=xc
+        xy[:,1]=yc
+        vor = Voronoi(xy)
+        regions, vertices = voronoi_finite_polygons_2d(vor,radius=radius)
+
+
+        for region in regions:
+            polygon0 = vertices[region]
+            P=polygon0.tolist()
+            polygon=np.array(P+[P[0]])
+            for iline in range(polygon.shape[0]-1):
+
+                x0,y0=CoordMachine.lm2radec(np.array([polygon[iline][0]]),np.array([polygon[iline][1]]))
+                x1,y1=CoordMachine.lm2radec(np.array([polygon[iline+1][0]]),np.array([polygon[iline+1][1]]))
+
+                x0*=180./np.pi
+                y0*=180./np.pi
+                x1*=180./np.pi
+                y1*=180./np.pi
+
+
+                f.write("line(%f,%f,%f,%f) # line=0 0 color=%s dash=1\n"%(x0,y0,x1,y1,Col))
+                #f.write("line(%f,%f,%f,%f) # line=0 0 color=red dash=1\n"%(x1,y0,x0,y1))
+
+        f.close()
+
+
+    def VorToReg(self,regFile,vor,radius=0.1,Col="red"):
+        print("Writing voronoi in: %s"%Str(regFile,col="blue"))
+
+        f=open(regFile,"w")
+        f.write("# Region file format: DS9 version 4.1\n")
+        ss0='global color=green dashlist=8 3 width=1 font="helvetica 10 normal roman" select=1 highlite=1 dash=0'
+        ss1=' fixed=0 edit=1 move=1 delete=1 include=1 source=1\n'
+
+        f.write(ss0+ss1)
+        f.write("fk5\n")
+
+        CoordMachine=self.CoordMachine
+
+        regions, vertices = vor.regions,vor.vertices
+
+
+        for region in regions:
+            if len(region)==0: continue
+            polygon0 = vertices[region]
+            P=polygon0.tolist()
+            polygon=np.array(P+[P[0]])
+            for iline in range(polygon.shape[0]-1):
+
+                x0,y0=CoordMachine.lm2radec(np.array([polygon[iline][0]]),np.array([polygon[iline][1]]))
+                x1,y1=CoordMachine.lm2radec(np.array([polygon[iline+1][0]]),np.array([polygon[iline+1][1]]))
+
+                x0*=180./np.pi
+                y0*=180./np.pi
+                x1*=180./np.pi
+                y1*=180./np.pi
+
+                f.write("line(%f,%f,%f,%f) # line=0 0 color=%s dash=1\n"%(x0,y0,x1,y1,Col))
+                #f.write("line(%f,%f,%f,%f) # line=0 0 color=red dash=1\n"%(x1,y0,x0,y1))
+
+        f.close()
+
+    def PolygonToReg(self,regFile,LPolygon,radius=0.1,Col="red",labels=None):
+        print("Writing voronoi in: %s"%Str(regFile,col="blue"))
+
+        f=open(regFile,"w")
+        f.write("# Region file format: DS9 version 4.1\n")
+        ss0='global color=green dashlist=8 3 width=1 font="helvetica 10 normal roman" select=1 highlite=1 dash=0'
+        ss1=' fixed=0 edit=1 move=1 delete=1 include=1 source=1\n'
+
+        f.write(ss0+ss1)
+        f.write("fk5\n")
+
+        CoordMachine=self.CoordMachine
+
+
+        for iFacet,polygon0 in zip(list(range(len(LPolygon))),LPolygon):
+            #polygon0 = vertices[region]
+            P=polygon0.tolist()
+            if len(polygon0)==0: continue
+            polygon=np.array(P+[P[0]])
+            ThisText=""
+            if labels is not None:
+                lmean0=np.mean(polygon[:,0])
+                mmean0=np.mean(polygon[:,1])
+
+                lmean,mmean,ThisText=labels[iFacet]
+                # print "!!!!======"
+                # print lmean0,mmean0
+                # print lmean,mmean
+
+                xm,ym=CoordMachine.lm2radec(np.array([lmean]),np.array([mmean]))
+                xm*=180./np.pi
+                ym*=180./np.pi
+                f.write("point(%f,%f) # text={%s} point=circle 5 color=red width=2\n"%(xm,ym,ThisText))
+
+            for iline in range(polygon.shape[0]-1):
+
+
+                L0,M0=np.array([polygon[iline][0]]),np.array([polygon[iline][1]])
+                x0,y0=CoordMachine.lm2radec(L0,M0)
+                L1,M1=np.array([polygon[iline+1][0]]),np.array([polygon[iline+1][1]])
+                x1,y1=CoordMachine.lm2radec(L1,M1)
+
+                x0*=180./np.pi
+                y0*=180./np.pi
+                x1*=180./np.pi
+                y1*=180./np.pi
+
+                # print "===================="
+                # print "[%3.3i] %f %f %f %f"%(iline,x0,y0,x1,y1)
+                # print "       %s"%str(L0)
+                # print "       %s"%str(L1)
+                # print "       %s"%str(M0)
+                # print "       %s"%str(M1)
+                f.write("line(%f,%f,%f,%f) # line=0 0 color=%s dash=1 \n"%(x0,y0,x1,y1,Col))
+
+                #f.write("line(%f,%f,%f,%f) # line=0 0 color=red dash=1\n"%(x1,y0,x0,y1))
+
+        f.close()
\ No newline at end of file

imcal.py

@@ -86,11 +86,21 @@ def wsclean(msin, datacolumn='DATA', outname=None, pixelsize=3, imagesize=3072,
     return 0
 
 
+def smoothImage(imgfits) :
+    """
+    Smoothe an image
+    """
+    cmd = f'smoFits.py {imgfits}'
+    logging.debug("Running command: %s", cmd)
+    subprocess.call(cmd, shell=True)
+    return
+
+
 def create_mask(imgfits, residfits, clipval, outname='mask.fits'):
     """
     Create mask using Tom's code (e-mail on 1 Jul 2021)
     """
-    cmd = f'makeNoiseMapFits {imgfits} {residfits} noise.fits noiseMap.fits'
+    cmd = f'makeNoiseMapFitsLow {imgfits} {residfits} noise.fits noiseMap.fits'
     logging.debug("Running command: %s", cmd)
     subprocess.call(cmd, shell=True)
     cmd = f'makeMaskFits noiseMap.fits {outname} {clipval}'
@@ -99,6 +109,30 @@ def create_mask(imgfits, residfits, clipval, outname='mask.fits'):
     return outname
 
 
+def makeNoiseImage(imgfits, residfits, low=False) :
+    """
+    Create mask using Tom's code (e-mail on 1 Jul 2021)
+    """
+    if low:
+      cmd = f'makeNoiseMapFitsLow {imgfits} {residfits} noiseLow.fits noiseMapLow.fits'
+    else :
+      cmd = f'makeNoiseMapFits {imgfits} {residfits} noise.fits noiseMap.fits'
+    logging.debug("Running command: %s", cmd)
+    subprocess.call(cmd, shell=True)
+    return
+
+
+def makeCombMask(ima1='noiseMap.fits', ima2='noiseMapLow.fits',
+		 clip1=5, clip2=7, outname='mask.fits') :
+    """
+    Create mask using Tom's code (e-mail on 1 Jul 2021)
+    """
+    cmd = f'makeCombMaskFits {ima1} {ima2} {outname} {clip1} {clip2}'
+    logging.debug("Running command: %s", cmd)
+    subprocess.call(cmd, shell=True)
+    return
+
+
 def get_image_max(msin):
     """
     Determine maximum image value for msin
@@ -187,8 +221,8 @@ def dical(msin, srcdb, msout=None, h5out=None, solint=1, startchan=0, split_ncha
            f'cal.solint={solint}',
            f'cal.parmdb={h5out}',
            f'cal.nchan={cal_nchan}',
-           f'cal.uvlambdamin={uvlambdamin}',
            'cal.applysolution=True',
+           'cal.blrange=[100,1000000]',
            'cal.type=gaincal',
            'steps=[cal]']
     if startchan or split_nchan:
@@ -199,9 +233,8 @@ def dical(msin, srcdb, msout=None, h5out=None, solint=1, startchan=0, split_ncha
     check_return_code(return_code)
     return msout
 
-
 def ddecal(msin, srcdb, msout=None, h5out=None, solint=120, nfreq=30,
-           startchan=0, nchan=0, minvisratio=0.6, mode='diagonal', uvlambdamin=500, subtract=True):
+           startchan=0, nchan=0,  mode='diagonal', uvlambdamin=500, subtract=True):
     """ Perform direction dependent calibration with DPPP """
     h5out = h5out or os.path.split(msin)[0] + '/ddcal.h5'
     msbase = os.path.basename(msin).split('.')[0]
@@ -218,13 +251,12 @@ def ddecal(msin, srcdb, msout=None, h5out=None, solint=120, nfreq=30,
           cal.subtract={subtract} \
           cal.propagatesolutions=true \
           cal.propagateconvergedonly=true \
-          cal.minvisratio={minvisratio} \
           cal.nchan={nfreq} \
           cal.uvlambdamin={uvlambdamin} \
           steps=[cal] \
           '.format(msin=msin, msout=msout, startchan=startchan, nchan=nchan, mode=mode,
             srcdb=srcdb, solint=solint, h5out=h5out, subtract=subtract, nfreq=nfreq,
-            minvisratio=minvisratio, uvlambdamin=uvlambdamin)
+            uvlambdamin=uvlambdamin)
     cmd = " ".join(cmd.split())
     logging.debug("Running command: %s", cmd)
     subprocess.call(cmd, shell=True)
@@ -243,50 +275,61 @@ def phase_shift(msin, new_center, msout=None):
 def view_sols(h5param, outname=None):
     """ read and plot the gains """
     def plot_sols(h5param, key):
-        print('AAA')
-        figs = []
-        axs = []
         with h5py.File(h5param, 'r') as f:
             grp = f['sol000/{}'.format(key)]
             data = grp['val'][()]
             time = grp['time'][()]
             ants = ['RT2','RT3','RT4','RT5','RT6','RT7','RT8','RT9','RTA','RTB','RTC','RTD']
-            freq_avg_gains = np.nanmean(data, axis=1)  # average by frequency
-            print(dict(f['sol000/amplitude000/val'].attrs.items())) # h5 attributes
-            for ipol, pol in enumerate(['XX', 'YY']):
             fig = plt.figure(figsize=[20, 15])
-                fig.suptitle('Freq. averaged {} gain solutions ({})'.format(key.rstrip('000'), pol))
+            fig.suptitle('Freq. averaged {} gain solutions'.format(key.rstrip('000')))
             for i, ant in enumerate(ants):
                 ax = fig.add_subplot(4, 3, i+1)
                 ax.set_title(ant)
-                    if key.startswith('phase'):
-                        ax.plot((time-time[0])/3600.0, freq_avg_gains[:, i,...,ipol]*180.0/np.pi, alpha=0.7)
-                        ax.set_ylim([-180,180])
+                if key == 'amplitude000' :
+                   ax.set_ylim(0,2)
                 else :
-                        ax.plot((time-time[0])/3600.0, freq_avg_gains[:, i,...,ipol], alpha=0.7)
-                    if key.startswith('amplitude'):
-                        ax.set_ylim([-0.1, np.max(freq_avg_gains)])
+                   ax.set_ylim(-180,180)
+                if len(data.shape) == 5: # several directions
+                    # a = ax.imshow(data[:,:,i,1,0].T, aspect='auto')
+                    # plt.colorbar(a)
+                    gavg = np.nanmean(data, axis=1)
+                    if key == 'amplitude000' :
+                      ax.plot((time-time[0])/60.0, gavg[:, i, :, 0], alpha=0.7)
+                      ax.plot((time-time[0])/60.0, gavg[:, i, :, 1], alpha=0.7)
+                    else :
+                      ax.plot((time-time[0])/60.0, 360.0/np.pi*gavg[:, i, :, 0], alpha=0.7)
+                      ax.plot((time-time[0])/60.0, 360.0/np.pi*gavg[:, i, :, 1], alpha=0.7)
+
+                elif len(data.shape) == 4: # a single direction
+                    if key == 'amplitude000' :
+                      gavg = np.nanmean(data,axis=1)
+#                      ax.plot((time-time[0])/3600.0, data[:, 0, i, 0], alpha=0.7)
+                      ax.plot((time-time[0])/3600.0, gavg[:,  i, 0], alpha=0.7,label='XX')
+                      ax.plot((time-time[0])/3600.0, gavg[:,  i, 0], alpha=0.7,label='YY')
+                    else :
+                      gavg = np.nanmean(data,axis=1)
+#                      ax.plot((time-time[0])/3600.0, 360.0/np.pi*data[:, 0, i, 0], alpha=0.7)
+#                      ax.plot((time-time[0])/3600.0, 360.0/np.pi*data[:,3 , i, 0], alpha=0.7)
+                      ax.plot((time-time[0])/3600.0, 360.0/np.pi*gavg[:,  i, 0], alpha=0.7,label='XX')
+                      ax.plot((time-time[0])/3600.0, 360.0/np.pi*gavg[:,  i, 1], alpha=0.7,label='YY')
+
+
                     if i == 0:
-                        ax.legend(['c{}'.format(_) for _ in range(data.shape[-2])])
+                      ax.legend(['XX','YY'])
                 if i == 10:
                     ax.set_xlabel('Time (hrs)')
+        return fig, ax
 
-                figs.append(fig)
-                axs.append(ax)
-        return figs, axs
-
-    try:
-        (fig1, fig2), (ax1, ax2) = plot_sols(h5param, 'amplitude000')
     if outname is not None:
-            fig1.savefig(f'{outname}_amp_XX.png')
-            fig2.savefig(f'{outname}_amp_YY.png')
+        try:
+            fig1, ax1 = plot_sols(h5param, 'amplitude000')
+            fig1.savefig(f'{outname}_amp.png')
         except:
             logging.error('No amplitude solutions found')
+
         try:
-        (fig1, fig2), (ax1, ax2) = plot_sols(h5param, 'phase000')
-        if outname is not None:
-            fig1.savefig(f'{outname}_phase_XX.png')
-            fig2.savefig(f'{outname}_phase_YY.png')
+            fig2, ax2 = plot_sols(h5param, 'phase000')
+            fig2.savefig(f'{outname}_phase.png')
         except:
             logging.error('No phase solutions found')
 
@@ -318,28 +361,6 @@ def remove_model_components_below_level(model, level=0.0, out=None):
     return out
 
 
-def apply_ddcal_sols(msin, parmdb, msout='.', msout_datacolumn='CORRECTED_DATA', amp_interp = 'nearest', ph_interp = 'nearest'):
-    """ apply calibration solution from h5file """
-
-    with h5py.File(parmdb, 'r') as f:
-        directions = list(dict(f['sol000/source']).keys())
-
-    dir_names = [i.decode().strip('[]') for i in directions]
-    dp3args = [f'msin={msin}', 'msin.datacolumn=DATA', 'msout=.', 'msout.datacolumn=CORRECTED_DATA',
-               f'steps={dir_names}']
-    for d in dir_names:
-        dp3args += [f'{d}.type=applycal',
-                    f'{d}.parmdb={parmdb}',
-                    f'{d}.steps=[ph,amp]',
-                    f'{d}.amp.correction=amplitude000',
-                    f'{d}.ph.correction=phase000',
-                    f'{d}.amp.interpolation={amp_interp}',
-                    f'{d}.ph.interpolation={ph_interp}',
-                    f'{d}.direction={[d]}',
-                    ]
-    execute_dppp(dp3args)
-
-
 def main(msin, outbase=None, cfgfile='imcal.yml'):
     msin = msin.rstrip('/')
     logging.info('Processing {}'.format(msin))
@@ -367,54 +388,69 @@ def main(msin, outbase=None, cfgfile='imcal.yml'):
     mask0 = outbase + '-mask0.fits'
     mask1 = outbase + '-mask1.fits'
     mask2 = outbase + '-mask2.fits'
+    mask3 = outbase + '-mask3.fits'
+    mask4 = outbase + '-mask4.fits'
 
+    nvssMod = outbase + '_nvss.sourcedb'
     model1 = outbase + '_model1.sourcedb'
     model2 = outbase + '_model2.sourcedb'
     model3 = outbase + '_model3.sourcedb'
 
+    dical0 = outbase + '_dical0.MS'
     dical1 = outbase + '_dical1.MS'
     dical2 = outbase + '_dical2.MS'
     dical3 = outbase + '_dical3.MS'
     ddsub = outbase + '_ddsub.MS'
 
+    h5_0 = outbase + '_dical0.h5'
     h5_1 = outbase + '_dical1.h5'
     h5_2 = outbase + '_dical2.h5'
     h5_3 = outbase + '_dical3.h5'
-    h5_dd = outbase + '_dd.h5'
+    h5_dd = outbase + '_ddcal.h5'
 
 
     if os.path.exists(img_ddcal+'-image.fits'):
-        logging.warning('The final image exists!..')
-        # return 0
+        logging.info('The final image exists. Exiting...')
+        return 0
 
-    if (not os.path.exists(ms_split)) and (cfg['split1']['startchan'] or cfg['split1']['nchan']):
+    if cfg['nvss']['doit'] :
+#    if (not os.path.exists(ms_split)) and (cfg['split1']['startchan'] or cfg['split1']['nchan']):
        ms_split = split_ms(msin, msout_path=ms_split, **cfg['split1'])
-    else:
-        ms_split = msin
 
-# Clean + DIcal
+       print('------------------makesource')
+       makesourcedb('nvss-model.txt', out=nvssMod)
+       print('-------------nvssCal')
+       dical(ms_split, nvssMod, msout=dical0, h5out=h5_0, **cfg['nvssCal'])
+       view_sols(h5_0, outname=msbase+'_sols_dical0')
+
+    else :
+       ms_split = split_ms(msin, msout_path=dical0, **cfg['split1'])
 
     if not os.path.exists(img0 +'-image.fits') and (not os.path.exists(img0 +'-MFS-image.fits')):
-        img_max = get_image_max(ms_split)
+        img_max = get_image_max(dical0)
         threshold = img_max/cfg['clean0']['max_over_thresh']
-        wsclean(ms_split, outname=img0, automask=None, save_source_list=False, multifreq=False, mgain=None,
+        threshold = max(threshold,0.001)
+        wsclean(dical0, outname=img0, automask=None, save_source_list=False, multifreq=False, mgain=None,
             kwstring=f'-threshold {threshold}')
-        create_mask(img0 +'-image.fits', img0 +'-residual.fits', clipval=20, outname=mask0)
+        create_mask(img0 +'-image.fits', img0 +'-residual.fits', clipval=10, outname=mask0)
 
 # clean1
     if not os.path.exists(img1 +'-image.fits') and (not os.path.exists(img1 +'-MFS-image.fits')):
-        wsclean(ms_split, fitsmask=mask0, outname=img1, **cfg['clean1']) # fast shallow clean
+        wsclean(dical0, 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'])
+        dical1 = dical(dical0, model1, msout=dical1, h5out=h5_1, **cfg['dical1'])
         view_sols(h5_1, outname=msbase+'_sols_dical1')
 # clean2
     if (not os.path.exists(img2 +'-image.fits')) and (not os.path.exists(img2 +'-MFS-image.fits')):
         wsclean(dical1, fitsmask=mask0, outname=img2, **cfg['clean2'])
-        create_mask(img2 +'-image.fits', img2 +'-residual.fits', clipval=7, outname=mask1)
+        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)
 
     if not os.path.exists(model2):
         makesourcedb(img2+'-sources.txt', out=model2)
@@ -426,7 +462,11 @@ def main(msin, outbase=None, cfgfile='imcal.yml'):
 # clean3
     if (not os.path.exists(img3 +'-image.fits')) and (not os.path.exists(img3 +'-MFS-image.fits')):
         wsclean(dical2, fitsmask=mask1, outname=img3, **cfg['clean3'])
-        create_mask(img3 +'-image.fits', img3 +'-residual.fits', clipval=5, outname=mask2)
+        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)
+
 
     if not os.path.exists(model3):
         makesourcedb(img3+'-sources.txt', out=model3)
@@ -439,6 +479,10 @@ def main(msin, outbase=None, cfgfile='imcal.yml'):
 # clean4
     if (not os.path.exists(img_final +'-image.fits')) and (not os.path.exists(img_final +'-MFS-image.fits')):
         wsclean(dical3, fitsmask=mask2, outname=img_final, **cfg['clean4'])
+        smoothImage(img_final+'-residual.fits')
+        makeNoiseImage(img_final +'-image.fits', img_final +'-residual.fits')
+        makeNoiseImage(img_final +'-residual-smooth.fits', img_final +'-residual.fits',low=True)
+        makeCombMask(outname=mask3,clip1=5,clip2=7)
 
 
 # Cluster
@@ -451,24 +495,29 @@ def main(msin, outbase=None, cfgfile='imcal.yml'):
 # DDE calibration + peeling everything
     if (not os.path.exists(ddsub)):
         ddsub, h5out = ddecal(dical3, clustered_sdb, msout=ddsub, h5out=h5_dd, **cfg['ddcal'])
+
 # view the solutions and save figure
         view_sols(h5_dd, outname=msbase+'_sols_ddcal')
 
-
-# Apply DDEcal soulutions to dical3
-
-    # apply_ddcal_sols(dical3, h5_dd) # does not work -- image is bad. Can you apply all solution to all directions with DP3?
-    # wsclean(dical3, datacolumn='CORRECTED_DATA', outname='test_apply_ddcal', **cfg['clean4'])
-    # sys.exit('planned')
-
     if (not os.path.exists(img_ddsub+'-image.fits')):
-        wsclean(ddsub, outname=img_ddsub, **cfg['clean4'])
+        wsclean(ddsub, fitsmask=mask3,outname=img_ddsub, **cfg['clean5'])
+#TAO        wsclean(ddsub,outname=img_ddsub, **cfg['clean5'])
 
     aomodel = bbs2model(img_final+'-sources.txt', img_final+'-model.ao')
 
     render(img_ddsub+'-image.fits', aomodel, out=img_ddcal+'-image.fits')
 
+    smoothImage(img_ddsub+'-image.fits')
+    makeNoiseImage(img_ddcal +'-image.fits', img_ddsub +'-residual.fits')
+    makeNoiseImage(img_ddsub +'-image-smooth.fits', img_ddsub +'-residual.fits',low=True)
+    makeCombMask(outname=mask4,clip1=3.5,clip2=5)
+
+    if (not os.path.exists(img_ddsub+'-2-image.fits')):
+        wsclean(ddsub, fitsmask=mask4,outname=img_ddsub+'-2', **cfg['clean5'])
+#TAO        wsclean(ddsub,outname=img_ddsub, **cfg['clean5'])
 
+    aomodel = bbs2model(img_final+'-sources.txt', img_final+'-model.ao')
+    render(img_ddsub+'-2-image.fits', aomodel, out=img_ddcal+'-2-image.fits')
 
     return 0
 

imcal.yml

@@ -11,9 +11,17 @@
 ####################### IMAGING #######################
 
 split1:
-    startchan: 40 # start channel to split from
+    startchan: 0 # start channel to split from
     nchan: 0 # 0 means till the end
 
+nvssCal: # DPPP setup for direction independent calibration
+    solint: 60
+    mode: 'phaseonly'
+    uvlambdamin: 500 # Ignore baselines / channels with UV < uvlambdamin wavele$
+
+nvss:
+    doit : False
+
 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
@@ -23,11 +31,11 @@ clean1: # wsclean setup
     pixelsize: 3
     multifreq: 0
     mgain: 0
-    automask: 20
+    automask: 10
     autothresh: 5
     multiscale: False
     clip_model_level: null # use a float number to clip the model to above this level. null is None.
-    kwstring: '-use-wgridder -parallel-deconvolution 1400 -weight briggs 0.0' # use this for additional wsclean options, e.g. '-weight uniform -use-idg'
+    kwstring: '-use-wgridder -parallel-deconvolution 1400 -weight briggs -1.5' # use this for additional wsclean options, e.g. '-weight uniform -use-idg'
 
 dical1: # DPPP setup for direction independent calibration
     solint: 20
@@ -38,13 +46,13 @@ dical1: # DPPP setup for direction independent calibration
 clean2:
     imagesize: 3072
     pixelsize: 3
-    mgain: 0
+    mgain: 0.8
     multifreq: 6
     automask: 10
     autothresh: 5
-    multiscale: True
+    multiscale: False
     clip_model_level: null
-    kwstring: '-use-wgridder -parallel-deconvolution 1400 -parallel-gridding 8 -deconvolution-channels 3 -weight briggs 0.0'
+    kwstring: '-use-wgridder -parallel-deconvolution 1400 -parallel-gridding 8 -deconvolution-channels 3 -weight briggs -1.5'
 
 dical2:
     solint: 1
@@ -58,12 +66,12 @@ clean3:
     multifreq: 6
     automask: 7
     autothresh: 3.5
-    multiscale: True
+    multiscale: False
     clip_model_level: null
-    kwstring: '-use-wgridder -parallel-deconvolution 1400 -parallel-gridding 8 -deconvolution-channels 3 -weight briggs 0.0'
+    kwstring: '-use-wgridder -parallel-deconvolution 1400 -parallel-gridding 8 -deconvolution-channels 3 -weight briggs -1.5'
 
 dical3:
-    solint: 240
+    solint: 120
     mode: 'diagonal'
     uvlambdamin: 500 # Ignore baselines / channels with UV < uvlambdamin wavelengths.
     cal_nchan: 72 # number of chans with the same solutions.
@@ -82,31 +90,29 @@ clean4:
 cluster:
     nbright: 80 # number of brightest clean components (CC) to check for artefacts
     boxsize: 250 # the boxsize around CC in pixels where to check for artefacts
-    nclusters: 6 # number of clusters ('auto' -- to set automatically)
-    clustering_method: 'Voronoi'
+    nclusters: 10 # number of clusters ('auto' -- to set automatically)
 # the following is only for 'auto' and 'manual' mathods:
     cluster_radius: 5 # arcmin
-    cluster_overlap: 2 # if lower than 2 clusters can intersect
+    cluster_overlap: 1.6 # if lower than 2 clusters can intersect
+    clustering_method :  voronoi # or auto
     add_manual: False
 
 ######################## DD CALIBRATION #######################
 ddcal: # see DPPP/DDECal documentation
-    solint: 60 # Solution interval in timesteps (1 corresponds to 30 seconds for Apertif).
+    solint: 120 # Solution interval in timesteps (1 corresponds to 30 seconds for Apertif).
     mode: 'diagonal' # Type of constraint to apply.
     nfreq: 72 # Number of channels in each channel block, for which the solution is assumed to be constant.
     startchan: 0
-    minvisratio: 0.6
     nchan: 0
     uvlambdamin: 500 # Ignore baselines / channels with UV < uvlambdamin wavelengths.
 
 clean5:
     imagesize: 3072
     pixelsize: 3
-    multifreq: 8
-    automask: 4
-    autothresh: 0.5
+    multifreq: 6
+    automask: 3.5
+    autothresh: 1.0
     multiscale: True
-    fitsmask: null
     clip_model_level: null
     kwstring: '-use-wgridder -parallel-deconvolution 1400 -parallel-gridding 8 -deconvolution-channels 3 -weight briggs -1.5'
 

makemask/makeCombMaskFits.c

+#include <string.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+
+/*
+  protgram to compute the 
+*/
+#include "fitsio.h"
+
+
+
+int main(int argc, char **argv) {
+
+  fitsfile *inImage1;      /* pointer to input image */
+  fitsfile *inImage2;      /* pointer to input image */
+
+  fitsfile *outMask;      /* pointer to output mask */
+
+  int status, morekeys, hdutype;
+  int  i, j, ind;
+  int naxis;
+  int xlen, ylen;
+    
+  int maxdim;
+  long  nelements, index;
+  long naxes[10];
+  int bitpix = SHORT_IMG ; /* output image is ints; BITPIX = -16*/
+    
+  float *imArray1;
+  float *imArray2;
+  int   *maskArray;
+     
+  float clipVal1;
+  float clipVal2;
+  int nullval ;
+  int anynull;
+   
+  status = 0;
+  maxdim = 10;
+
+  if (argc != 6) {
+    puts("makeCombMaskFits image1 image2 outmask clip1 clip2 ");
+    exit(1);
+  }
+  clipVal1 = atof(argv[4]);
+  clipVal2 = atof(argv[5]);
+ 
+          
+  /* open the existing image */
+  fits_open_file(&inImage1, argv[1], READONLY, &status) ;
+  fits_open_file(&inImage2, argv[2], READONLY, &status) ;
+
+    
+  /* get number of dimensions */
+  fits_get_img_dim(inImage1, &naxis, &status);
+  /* get size of image. assume x and y are axes 0 and 1 */  
+  fits_get_img_size(inImage1, maxdim, naxes, &status);
+
+  /* number of pixels in image */
+  nelements = naxes[0]*naxes[1];
+  xlen = naxes[0];
+  ylen = naxes[1];
+
+  /* allocare memeory for the arrays */
+  imArray1 = (float *)malloc(nelements*sizeof(float));
+  imArray2 = (float *)malloc(nelements*sizeof(float));
+  maskArray = (int *)malloc(nelements*sizeof(int));
+
+  /* don't check for null values */
+  nullval = 0;
+
+  /* read the images */
+  fits_read_img(inImage1, TFLOAT, 1, nelements, &nullval,
+                imArray1, &anynull, &status) ;
+  fits_read_img(inImage2, TFLOAT, 1, nelements, &nullval,
+                imArray2, &anynull, &status) ;
+
+  /* loop over image */
+  for (i = 0; i < xlen; i++) {
+    for (j = 0; j < ylen;  j++) {
+      index = i + xlen*(j) ;
+      if (imArray1[index] > clipVal1 || imArray2[index] > clipVal2) {
+         maskArray[index]= 1;
+      } else {
+        maskArray[index]= 0;
+      }
+    }
+  }
+
+  /* delete output mask file */
+  remove(argv[3]);
+  /* create output mask file, copy header and write maskfits
+     file */
+  fits_create_file(&outMask, argv[3], &status);
+  fits_copy_hdu(inImage1, outMask, 0, &status);
+  fits_write_img(outMask,TINT,1,nelements,maskArray,&status);
+
+  /* close all files */
+  fits_close_file(inImage1,&status);
+  fits_close_file(inImage2,&status);
+  fits_close_file(outMask,&status);
+
+
+
+
+}

makemask/makeMaskFits.c

@@ -6,7 +6,7 @@
 /*
   protgram to compute the 
 */
-#include <fitsio.h>
+#include "fitsio.h"
 
 
 
@@ -35,6 +35,12 @@ int main(int argc, char **argv) {
    
   status = 0;
   maxdim = 10;
+  
+  if (argc != 4) {
+    puts("makeMaskFits inimage outmask clip");
+    exit(1);
+  }
+
   clipVal = atof(argv[3]);
  
           

makemask/makeNoiseMapFits.c

@@ -6,7 +6,7 @@
 /*
   protgram to compute the 
 */
-#include <fitsio.h>
+#include "fitsio.h"
 
 int comp(const void *a,const void *b) {
 float *x = (float *) a;
@@ -46,6 +46,11 @@ int main(int argc, char **argv) {
   maxdim = 10;
   boxsize = 13;
 
+  if (argc != 5) {
+    puts("makeNoiseMapFits inImage inResidual outNoise ourNoiseMap");
+    exit(1);
+  }
+        
   /* open the existing image */
   fits_open_file(&inImage, argv[1], READONLY, &status) ;
   /* open the residual image */

makemask/makeNoiseMapFitsLow.c

+#include <string.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+
+/*
+  protgram to compute the 
+*/
+#include "fitsio.h"
+
+int comp(const void *a,const void *b) {
+float *x = (float *) a;
+float *y = (float *) b;
+// return *x - *y; // this is WRONG...
+if (*x < *y) return -1;
+else if (*x > *y) return 1; return 0;
+}
+
+int main(int argc, char **argv) {
+
+  fitsfile *inImage;      /* pointer to input image */
+  fitsfile *inRes;        /* pointer to input residual */
+
+  fitsfile *outNoise;     /* pointer to output noise image */
+  fitsfile *outSN;        /* pointer to output S/N image */
+
+  int status, morekeys, hdutype;
+  int ii, jj, i, j, n, m, ind;
+  int naxis;
+  int xlen, ylen;
+    
+  int maxdim;
+  long  nelements, index, numP;;
+  long naxes[10];
+  int bitpix = FLOAT_IMG ; /* output image is floats; BITPIX = -32*/
+    
+  float *imArray, *resArray, *noiseArray, *snArray;
+  float tmpAr[4096];
+    
+  int boxsize;
+  int nullval ;
+  int anynull;
+  float sig;
+  
+  status = 0;
+  maxdim = 10;
+  boxsize = 50;
+
+  if (argc != 5) {
+    puts("makeNoiseMapFitsLow.c inImage inResidual outNoise outNoiseMap");
+    exit(1);
+  }
+
+  /* open the existing image */
+  fits_open_file(&inImage, argv[1], READONLY, &status) ;
+  /* open the residual image */
+  fits_open_file(&inRes, argv[2], READONLY, &status) ;
+    
+  /* get number of dimensions */
+  fits_get_img_dim(inImage, &naxis, &status);
+  /* get size of image. assume x and y are axes 0 and 1 */  
+  fits_get_img_size(inImage, maxdim, naxes, &status);
+
+  /* number of pixels in image */
+  nelements = naxes[0]*naxes[1];
+  xlen = naxes[0];
+  ylen = naxes[1];
+
+
+  /* allocare memeory for the arrays */
+  imArray = (float *)malloc(nelements*sizeof(float));
+  resArray = (float *)malloc(nelements*sizeof(float));
+  noiseArray = (float *)malloc(nelements*sizeof(float));
+  snArray = (float *)malloc(nelements*sizeof(float));
+
+
+  /* don't check for null values */
+  nullval = 0;
+
+  /* read the image */
+  fits_read_img(inImage, TFLOAT, 1, nelements, &nullval,
+                imArray, &anynull, &status) ;
+  /* read the residual */
+  fits_read_img(inRes, TFLOAT, 1, nelements, &nullval,
+               resArray, &anynull, &status) ;
+
+
+  /* set outpit noise array to zero */
+  for (ii = 0; ii < nelements; ii++) {
+    noiseArray[ii] = 0.0;
+  }
+
+
+  /* number of pixels in the box for which we compute the median */
+  numP = (2*boxsize+1)*(2*boxsize+1);
+  /* compute index that will point to the median */
+  numP = numP/2;
+
+  /* loop over image in steps of 5 pixels */
+  for (i = boxsize; i < xlen-boxsize-1; i=i+9) {
+    for (j = boxsize; j < ylen-boxsize-1; j= j+9) {
+      /* get the median for the current box */
+
+      /* reset index */
+      ind = 0;
+      /* loop over box, putting absolute value of the image value in temp array */
+      for (m = -boxsize; m <= boxsize; m=m+5) {
+        for (n = -boxsize; n <= boxsize; n=n+5) {
+          index = i+n + xlen*(j+m);
+          tmpAr[ind] = fabs(resArray[index]);
+          ind++;
+        }
+      }
+
+      /* sort the temp array and take the central value */
+      qsort(tmpAr,ind,sizeof(float),comp);
+      /* multiply by 1.4826 to turn the median absolute value into sigma */
+      sig = tmpAr[ind/2]*1.4826;
+      /* and put this value in the output noise image */
+      for (m = -4; m < 5; m++) {
+        for (n = -4; n < 5; n++) {
+          index = i+m + xlen*(j+n) ;
+          noiseArray[index]= sig;
+          /* if not zero, compute the S/N, i.e. the image value / noise value */
+          if (sig > 0.0) {
+            snArray[index]= imArray[index]/sig;
+          }
+
+        }
+      }
+    }
+  }
+
+  /* delete output noise file */
+  remove(argv[3]);
+  /* create output noise file, copy header and write noise values into fits
+     file */
+  fits_create_file(&outNoise, argv[3], &status);
+  fits_copy_hdu(inRes, outNoise, 0, &status);
+  fits_write_img(outNoise,TFLOAT,1,nelements,noiseArray,&status);
+
+  /* same as above, but now for S/N image */
+  remove(argv[4]);
+  
+  fits_create_file(&outSN, argv[4], &status);
+  fits_copy_hdu(inRes, outSN, 0, &status);
+  fits_write_img(outSN,TFLOAT,1,nelements,snArray,&status);
+
+  /* close all files */
+  fits_close_file(inImage,&status);
+  fits_close_file(inRes,&status);
+  fits_close_file(outNoise,&status);
+  fits_close_file(outSN,&status);
+
+
+
+}