diff --git a/cluster_to_reg.py b/cluster_to_reg.py new file mode 100644 index 0000000000000000000000000000000000000000..e37fcd2ce79c90cb9cece5d5e1f34c26cab553e3 --- /dev/null +++ b/cluster_to_reg.py @@ -0,0 +1,432 @@ +#!/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 diff --git a/imcal.py b/imcal.py index 2cba5d1325a755dc4b944428ad78ddf012aa2bdd..cb6fe2414dbfd10f5a102ce7a5f8ce2e4f1fd28f 100755 --- a/imcal.py +++ b/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,52 +275,63 @@ 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)) - 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]) - 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)]) + fig = plt.figure(figsize=[20, 15]) + 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 == 'amplitude000' : + ax.set_ylim(0,2) + else : + 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])]) - if i == 10: - ax.set_xlabel('Time (hrs)') - - 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') - 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') - except: - logging.error('No phase solutions found') + ax.legend(['XX','YY']) + if i == 10: + ax.set_xlabel('Time (hrs)') + return fig, ax + + if outname is not None: + try: + fig1, ax1 = plot_sols(h5param, 'amplitude000') + fig1.savefig(f'{outname}_amp.png') + except: + logging.error('No amplitude solutions found') + + try: + fig2, ax2 = plot_sols(h5param, 'phase000') + fig2.savefig(f'{outname}_phase.png') + except: + logging.error('No phase solutions found') def remove_model_components_below_level(model, level=0.0, out=None): @@ -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']): - ms_split = split_ms(msin, msout_path=ms_split, **cfg['split1']) - else: - ms_split = msin + 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']) + + 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') -# Clean + DIcal + 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 diff --git a/imcal.yml b/imcal.yml index 33c01140a6625d12f52b68aec4001981ad82db3a..24cce7a9b8e5103811f9a6756267053ddd4b8f5b 100644 --- a/imcal.yml +++ b/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' diff --git a/makemask/makeCombMaskFits.c b/makemask/makeCombMaskFits.c new file mode 100644 index 0000000000000000000000000000000000000000..f057703903d47639bd460bfbff59eac16101f7ce --- /dev/null +++ b/makemask/makeCombMaskFits.c @@ -0,0 +1,107 @@ +#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); + + + + +} diff --git a/makemask/makeMaskFits.c b/makemask/makeMaskFits.c index 3060cf9dd8a43da1b6678eb28b27b7808e722877..8f49affa28b031f652905bde857c2744f3a1ad79 100644 --- a/makemask/makeMaskFits.c +++ b/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]); diff --git a/makemask/makeNoiseMapFits.c b/makemask/makeNoiseMapFits.c index e420f18f0a4009fcb6ccec9e4b9de3746605f237..b8c74560e1123f9490f87d7a283bd9caa1106fdf 100644 --- a/makemask/makeNoiseMapFits.c +++ b/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; @@ -45,6 +45,11 @@ int main(int argc, char **argv) { status = 0; 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) ; diff --git a/makemask/makeNoiseMapFitsLow.c b/makemask/makeNoiseMapFitsLow.c new file mode 100644 index 0000000000000000000000000000000000000000..a040d7cbba16a995e4d2975ef3f92493e1540838 --- /dev/null +++ b/makemask/makeNoiseMapFitsLow.c @@ -0,0 +1,157 @@ +#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); + + + +}