Source
1
-
import os
2
-
import numpy as np
3
-
import numba as nb
1
+
4
2
import casatools
5
3
from casatasks.private.imagerhelpers.input_parameters import saveparams2last
6
4
from casatasks.private.imagerhelpers.msuvbinflag_algorithms import UVGridFlag
7
5
import time
8
-
from scipy.optimize import curve_fit
9
-
import pdb
10
-
import matplotlib.pyplot as pl
11
-
ms = casatools.ms()
12
-
tb = casatools.table()
13
-
me = casatools.measures()
14
-
qa = casatools.quanta()
15
-
ia = casatools.image()
16
-
im = casatools.imager()
17
-
#####debug plots
18
-
doplot=True
19
-
dryrun=False
20
-
############
21
-
pl.ion()
22
-
@nb.njit(cache=True)
23
-
def populate_grid(uvw, stokesI, uvgrid, uvgrid_npt, deltau, deltav, npix):
24
-
for ii in range(len(uvw[0])):
25
-
uidx = int(np.round(uvw[0][ii]//deltau + npix//2))
26
-
vidx = int(np.round(uvw[1][ii]//deltav + npix//2))
27
-
28
-
uvgrid[uidx, vidx] += stokesI[ii]
29
-
uvgrid_npt[uidx, vidx] += 1
30
-
31
-
return uvgrid, uvgrid_npt
32
-
33
-
34
-
def mad(inpdat):
35
-
"""
36
-
Calculate the STD via MAD for the input data
37
-
38
-
Inputs:
39
-
inpdat Input numpy array
40
-
41
-
Returns:
42
-
std Calculate the std via mad
43
-
"""
44
-
45
-
med = np.median(inpdat)
46
-
mad = np.median(np.abs(inpdat - med))
47
-
48
-
return 1.4826 * mad
49
-
50
-
51
-
def msuvbin_to_uvgrid(ms, npix, deltau, deltav):
52
-
tb.open(ms)
53
-
uvw = tb.getcol('UVW')
54
-
data = tb.getcol('DATA')
55
-
tb.close()
56
-
57
-
umin, umax = np.min(uvw[0]), np.max(uvw[0])
58
-
vmin, vmax = np.min(uvw[1]), np.max(uvw[1])
59
-
wmin, wmax = np.min(uvw[2]), np.max(uvw[2])
60
-
61
-
uvals = np.linspace(umin, umax, npix)
62
-
vvals = np.linspace(vmin, vmax, npix)
63
-
64
-
uvgrid = np.zeros((npix, npix), dtype=np.complex128)
65
-
uvgrid_npt = np.zeros((npix, npix), dtype=int)
66
-
67
-
stokesI = 0.5*(data[0] + data[1])
68
-
stokesI = np.squeeze(stokesI)
69
-
70
-
uvgrid, uvgrid_npt = populate_grid(uvw, stokesI, uvgrid, uvgrid_npt, deltau, deltav, npix)
71
-
# Average per uv cell
72
-
idx = np.where(uvgrid_npt != 0)
73
-
uvgrid[idx] = uvgrid[idx]/uvgrid_npt[idx]
74
-
75
-
return uvgrid, uvgrid_npt
76
-
77
-
78
-
def resid_cube(x, a, b, c, d):
79
-
return a*x**3 + b*x**2 + c*x + d
80
-
def resid_cinco(x, a, b, c, d, e, f):
81
-
return a*x**5 + b*x**4 + c*x**3 + d*x**2+ e*x +f
82
-
83
-
def fit_radial_profile(xvals, yvals, ystd, deg=3, clip_sigma=3):
84
-
"""
85
-
Fit the radial profile with a polynomial
86
-
"""
87
-
88
-
print(f"deg {deg}, clip_sigma {clip_sigma}")
89
-
print(f"xvals shape {xvals.shape}, yvals shape {yvals.shape}")
90
-
print(f"non-zero x {np.count_nonzero(xvals)}, non-zero y {np.count_nonzero(yvals)}")
91
-
92
-
idx = np.where(yvals != 0)
93
-
xvals = xvals[idx]
94
-
yvals = yvals[idx]
95
-
ystd = ystd[idx]
96
-
#print(f"xvals {xvals}, yvals {yvals}")
97
-
98
-
# Fit the radial profile with a polynomial
99
-
#pfit, pcov = curve_fit(resid_function, xvals, yvals, sigma=ystd, p0=[1, 1, 1, 1])
100
-
101
-
pfit = np.polyfit(xvals, yvals, deg=2)
102
-
yfit = np.polyval(pfit, np.linspace(xvals.min(), xvals.max(), 100))
103
-
104
-
has_converged = False
105
-
while not has_converged:
106
-
# Outlier rejection while fitting
107
-
resid = yvals - np.polyval(pfit, xvals)
108
-
resid_mad = mad(resid)
109
-
resid_med = np.median(resid)
110
-
idx = np.where((resid > resid_med - clip_sigma*resid_mad) & (resid < resid_med + clip_sigma*resid_mad))
111
-
112
-
if len(idx[0]) == len(xvals):
113
-
has_converged = True
114
-
continue
115
-
116
-
xvals = xvals[idx]
117
-
yvals = yvals[idx]
118
-
ystd = ystd[idx]
119
-
120
-
#pfit = np.polyfit(xvals, yvals, deg=deg)
121
-
122
-
pfit, pcov = curve_fit(resid_cube, xvals, yvals, sigma=ystd, p0=[1, 1, 1, 1])
123
-
yfit = np.polyval(pfit, np.linspace(xvals.min(), xvals.max(), 100))
124
-
125
-
#import matplotlib.pyplot as plt
126
-
127
-
#fig, ax = plt.subplots()
128
-
#ax.plot(xvals, yvals, '-o', label='Data')
129
-
#ax.plot(np.linspace(xvals.min(), xvals.max(), 100), yfit, '--', label='Fit')
130
-
#ax.set_xlabel('Radius')
131
-
#ax.set_ylabel('Intensity')
132
-
#ax.legend()
133
-
#plt.tight_layout()
134
-
#plt.show()
135
-
136
-
return pfit
137
-
138
-
139
-
def calc_radial_profile_ann(uvgrid, deltau, deltav):
140
-
"""
141
-
Calculate the annular average of the uvgrid for every radius,
142
-
and fit the 1D radial profile with a polynomial.
143
-
"""
144
-
145
-
nbin = 30
146
-
147
-
npixx, npixy = uvgrid.shape[0], uvgrid.shape[1]
148
-
cx, cy = npixx//2, npixy//2
149
-
150
-
print(f"npixx {npixx}, npixy {npixy}, centerx {cx} centery {cy}")
151
-
152
-
# Generate grid of radii
153
-
x = np.arange(npixx) - cx
154
-
y = np.arange(npixy) - cy
155
-
rad = np.sqrt(x**2 + y**2).astype(int)
156
-
157
-
yy, xx = np.meshgrid(x, y)
158
-
radgrid = np.sqrt(xx**2 + yy**2).astype(int)
159
-
160
-
print(f"x shape {x.shape}, y shape {y.shape}, rad shape {rad.shape}")
161
-
print(f"x min {x.min()}, x max {x.max()}")
162
-
print(f"y min {y.min()}, y max {y.max()}")
163
-
print(f"rad min {rad.min()}, rad max {rad.max()}")
164
-
165
-
# Create log-spaced annuli to account for reducing UV coverage with radius
166
-
# Minimum annulus is 5px
167
-
annuli = np.logspace(np.log10(5), np.log10(rad.max()), nbin)
168
-
annuli = np.round(annuli).astype(int)
169
-
170
-
radial_mean = np.zeros(nbin)
171
-
radial_mad = np.zeros(nbin)
172
-
173
-
ann_min = 0
174
-
for idx, ann in enumerate(annuli):
175
-
ridx = np.where((radgrid >= ann_min) & (radgrid < ann))
176
-
uvgrid_sel = uvgrid[ridx]
177
-
uvgrid_sel = uvgrid_sel[np.abs(uvgrid_sel) != 0]
178
-
179
-
if len(uvgrid_sel) == 0:
180
-
radial_mean[idx] = 0
181
-
radial_mad[idx] = 0.
182
-
else:
183
-
radial_mean[idx] = np.mean(uvgrid_sel)
184
-
radial_mad[idx] = mad(uvgrid_sel)
185
-
186
-
print(f" ann_left {ann_min} ann_right {ann} radial_mean {radial_mean[idx]} radial_mad {radial_mad[idx]}")
187
-
print('-------------------------------')
188
-
189
-
ann_min = ann
190
-
191
-
return uvgrid, radial_mean, radial_mad, annuli
192
-
193
-
194
-
def calc_radial_profile_pix(uvgrid, deltau, deltav):
195
-
"""
196
-
Calculate the radial per-pixel average of the uvgrid for every radius,
197
-
and fit the 1D radial profile with a polynomial.
198
-
"""
199
-
200
-
npixx, npixy = uvgrid.shape[0], uvgrid.shape[1]
201
-
cx, cy = npixx//2, npixy//2
202
-
203
-
print(f"npixx {npixx}, npixy {npixy}, centerx {cx} centery {cy}")
204
-
205
-
# Generate radial values from 0 to max
206
-
x = np.arange(cx)
207
-
y = np.arange(cy)
208
-
rad = np.sqrt(x**2 + y**2).astype(int)
209
-
210
-
print(f"x shape {x.shape}, y shape {y.shape}, rad shape {rad.shape}")
211
-
print(f"rad [0] {rad[0]}, rad [-1] {rad[-1]}")
212
-
print(f"rad {rad}")
213
-
214
-
# Generate grid of radii
215
-
x = np.arange(npixx) - cx
216
-
y = np.arange(npixy) - cy
217
-
yy, xx = np.meshgrid(x, y)
218
-
radgrid = np.sqrt(xx**2 + yy**2).astype(int)
219
-
220
-
radial_mean = np.zeros(np.max([cx, cy]))
221
-
radial_mad = np.zeros(np.max([cx, cy]))
222
-
223
-
for idx, rr in enumerate(rad):
224
-
if idx == len(rad) - 1:
225
-
ridx = np.where((radgrid >= rad[idx]))
226
-
else:
227
-
ridx = np.where((radgrid > rad[idx]) & (radgrid <= rad[idx+1]))
228
-
229
-
uvgrid_sel = uvgrid[ridx]
230
-
uvgrid_sel = uvgrid_sel[np.abs(uvgrid_sel) != 0]
231
-
232
-
if len(uvgrid_sel) == 0:
233
-
radial_mean[idx] = 0
234
-
radial_mad[idx] = 0.
235
-
else:
236
-
radial_mean[idx] = np.mean(uvgrid_sel)
237
-
radial_mad[idx] = mad(uvgrid_sel)
238
-
239
-
return uvgrid, radial_mean, radial_mad, rad*deltau
240
-
def calc_radial_profile_and_fit(uvgrid, wgtgrid, flggrid, nsigma):
241
-
npixx, npixy = uvgrid.shape[0], uvgrid.shape[1]
242
-
cx, cy = npixx//2, npixy//2
243
-
#print(f"npixx {npixx}, npixy {npixy}, centerx {cx} centery {cy}")
244
-
245
-
# Generate radial values from 0 to max
246
-
x = np.arange(cx)
247
-
y = np.arange(cy)
248
-
rad = np.sqrt(x**2 + y**2).astype(int)
249
-
npoints=int(np.max(rad))+1
250
-
radamp=np.zeros(npoints)
251
-
radamp2=np.zeros(npoints)
252
-
radwght=np.zeros(npoints)
253
-
xval=np.array(range(npoints), dtype='float')
254
-
for j in range(1,npixy):
255
-
yval2=(j-cy)*(j-cy)
256
-
for k in range(1,npixx):
257
-
rval=int(np.sqrt((k-cx)*(k-cx)+yval2))
258
-
if(wgtgrid[k,j] > 0.0):
259
-
absval=np.abs(uvgrid[k,j])
260
-
radamp[rval]+=absval*wgtgrid[k,j]
261
-
radamp2[rval]+=absval*absval*wgtgrid[k,j]
262
-
radwght[rval]+=wgtgrid[k,j]
263
-
if(np.max(radwght)==0.0):
264
-
#empty channel
265
-
return
266
-
radamp[radwght!=0]=radamp[radwght!=0]/radwght[radwght!=0]
267
-
radamp2[radwght!=0]=radamp2[radwght!=0]/radwght[radwght!=0]
268
-
maxsenspos=np.argmax(radwght)
269
-
#normalize radweight
270
-
normrdwght=radwght/np.max(radwght)
271
-
sig=np.sqrt(np.abs(radamp2-radamp*radamp))
272
-
#nescale relative sigmas by number of weights att the point
273
-
medsig=np.median(sig[sig !=0])
274
-
sigtouse=sig[maxsenspos]
275
-
#sig[normrdwght!=0]=sig[normrdwght!=0]/normrdwght[normrdwght!=0]
276
-
xvalnz=xval[(sig !=0.0) & (radamp !=0)]
277
-
radampnz=radamp[(sig !=0) & (radamp !=0)]
278
-
fitnz=curve_fit(resid_cinco, xvalnz, radampnz)
279
-
###
280
-
print('min max of sig and max sens one', np.min(sig), np.max(sig), sigtouse)
281
-
signz=sig[sig !=0.0]
282
-
sig=np.interp(xval, xvalnz, signz)
283
-
#print( 'corvar ', fitnz[1])
284
-
radfit=resid_cinco(xval, *fitnz[0])
285
-
#radamp=np.ma.array(radamp, mask=(radwght == 0))
286
-
#radfit=np.ma.array(radfit, mask=(radwght == 0))
287
-
max_rad_idx=np.where(xval==np.max(xvalnz))[0][0]
288
-
if(doplot):
289
-
#pl.figure()
290
-
ax1=pl.subplot(211)
291
-
292
-
#pl.plot(xval, radfit+sig,'+')
293
-
ax1.errorbar(xval[0:max_rad_idx], radfit[0:max_rad_idx], yerr=sig[0:max_rad_idx], ecolor='lime', fmt='none', label='sigma')
294
-
ax1.plot(xval[0:max_rad_idx],radamp[0:max_rad_idx],'o',color='magenta', label='mean radial value')
295
-
ax1.plot(xval[0:max_rad_idx], radfit[0:max_rad_idx], 'k', label='fitted radial value')
296
-
ax1.set_ylabel('Amplitude')
297
-
ax1.set_xlabel('uvdist in pix')
298
-
ax1.legend()
299
-
if(doplot):
300
-
ax2=pl.subplot(212)
301
-
for j in range(npixy):
302
-
for k in range(npixx):
303
-
#sweep over all points
304
-
# if points are not already flagged
305
-
yval2=(j-cy)*(j-cy)
306
-
if(not flggrid[k,j]):
307
-
r=int(np.sqrt(yval2 +(k-cx)*(k-cx)))
308
-
#if(r < npoints and np.abs(uvgrid[k,j]) > (radfit[r]+nsigma*max(medsig, sig[r]))):
309
-
if(r < npoints and (np.abs(uvgrid[k,j]) > (radfit[r]+nsigma*sigtouse))):
310
-
if(doplot):
311
-
ax2.plot(r, np.abs(uvgrid[k,j]), 'go')
312
-
uvgrid[k,j]=0
313
-
wgtgrid[k,j]=0
314
-
flggrid[k,j]=True
315
-
else:
316
-
if(doplot):
317
-
ax2.plot(r, np.abs(uvgrid[k,j]), 'b+')
318
-
319
-
320
-
321
-
322
-
@nb.njit(cache=True)
323
-
def apply_flags(dat, flg, uvw, radial_mean, radial_mad, annuli, nsigma=5.):
324
-
"""
325
-
Apply flags based on the radial profile to the input data column
326
-
"""
327
-
328
-
print(f"dat.shape {dat.shape}, flg.shape {flg.shape}, uvw.shape {uvw.shape}")
329
-
330
-
nrow = uvw.shape[1]
331
-
332
-
for rr in range(nrow):
333
-
uvlen = np.sqrt(uvw[0][rr]**2 + uvw[1][rr]**2)
334
-
annidx = np.searchsorted(annuli, uvlen)
335
-
336
-
if np.abs(dat[..., rr]) > radial_mean[annidx] + nsigma*radial_mad[annidx]:
337
-
flg[..., rr] = True
338
-
339
-
return flg
340
-
341
-
342
-
343
-
def accumulate_continuum_grid(tb, npol, nchan, npoints, deltau, deltav):
344
-
"""
345
-
If the input msuvbin has multiple channels, loop over them to
346
-
accumulate on a single grid. This allows for a better estimate of the
347
-
radial profile from a "fuller" UV grid before flagging outliers
348
-
per-plane.
349
-
350
-
Inputs:
351
-
tb Table object - must be open
352
-
npol Number of polarizations'
353
-
nchan Number of channels
354
-
npoints Number of points in the grid
355
-
deltau U spacing in lambda
356
-
deltav V spacing in lambda
357
-
358
-
Returns:
359
-
uvgrid Accumulated UV grid
360
-
uvnpt Number of points per UV cell
361
-
"""
362
-
363
-
uvgrid_cont = np.zeros((npoints, npoints), dtype=np.complex128)
364
-
wgtgrid_cont = np.zeros((npoints, npoints), dtype=np.float64)
365
-
366
-
for pol in range(npol):
367
-
for chan in range(nchan):
368
-
dat = tb.getcolslice('DATA', [pol, chan], [pol, chan], [1, 1])
369
-
flg = tb.getcolslice('FLAG', [pol, chan], [pol, chan], [1, 1])
370
-
wgt = tb.getcolslice('WEIGHT_SPECTRUM', [pol, chan], [pol, chan], [1, 1])
371
-
uvw = tb.getcol('UVW')
372
-
373
-
if dat.size == 0 or flg.size == 0 or wgt.size == 0:
374
-
print(f"Zero size array read. Skipping.")
375
-
continue
376
-
377
-
dat_grid = dat[0, 0, :].reshape([npoints, npoints])
378
-
flg_grid = flg[0, 0, :].reshape([npoints, npoints])
379
-
wgt_grid = wgt[0, 0, :].reshape([npoints, npoints])
380
-
381
-
# Flag the data as necessary
382
-
#dat_grid = dat_grid * ~flg_grid
383
-
384
-
uvgrid_cont += dat_grid
385
-
wgtgrid_cont += wgt_grid
386
-
387
-
return uvgrid_cont, wgtgrid_cont
388
-
389
-
def flagViaBin_radial(binnedvis, sigma=5):
390
-
pdb.set_trace()
391
-
tb.open(binnedvis, nomodify=False)
392
-
msuvbinkey = tb.getkeyword('MSUVBIN')
393
-
394
-
# msuvbinkey.keys()
395
-
# Out[5]: dict_keys(['csys', 'nchan', 'npol', 'numvis', 'nx', 'ny', 'sumweight'])
396
-
397
-
# in radian
398
-
dra = msuvbinkey['csys']['direction0']['cdelt'][0]
399
-
ddec = msuvbinkey['csys']['direction0']['cdelt'][1]
400
-
401
-
nx = msuvbinkey['nx']
402
-
ny = msuvbinkey['ny']
403
-
404
-
# in radian
405
-
ra_extent = dra*nx
406
-
dec_extent = ddec*ny
407
-
408
-
# in Lambda
409
-
deltau = 1./ra_extent
410
-
deltav = 1./dec_extent
411
-
412
-
npol = msuvbinkey['npol']
413
-
nchan = msuvbinkey['nchan']
414
-
415
-
npoints = min(nx, ny)
416
-
417
-
dat = tb.getcol('UVW')
418
-
419
-
# Accumulate all channels in a single grid
420
-
uvgrid_cont, wgtgrid_cont = accumulate_continuum_grid(tb, npol, nchan, npoints, deltau, deltav)
421
-
# Calculate the radial profile
422
-
uvgrid, radial_mean, radial_mad, annuli = calc_radial_profile_ann(uvgrid_cont, deltau, deltav)
423
-
#radial_fit = fit_radial_profile(annuli, np.abs(radial_mean), radial_mad, deg=2)
424
-
425
-
for mean, mad, ann in zip(radial_mean, radial_mad, annuli):
426
-
print(f"mean {mean}, mad {mad}, ann {ann}")
427
-
428
-
# XXX : For debugging only. Remove later
429
-
if doplot:
430
-
import matplotlib.pyplot as plt
431
-
from matplotlib.colors import LogNorm
432
-
fig, ax = plt.subplots()
433
-
ax.plot(annuli[2:], np.abs(radial_mean)[2:], '-o', label='data')
434
-
#ax.plot(np.linspace(annuli.min(), annuli.max(), 100), np.polyval(radial_fit, np.linspace(annuli.min(), annuli.max(), 100)), label='fit')
435
-
ax.fill_between(annuli, np.abs(radial_mean)-radial_mad, np.abs(radial_mean)+radial_mad, alpha=0.5)
436
-
ax.set_xlabel('Radius')
437
-
ax.set_ylabel('Intensity (Jy)')
438
-
ax.set_title('Radial Mean')
439
-
#ax.set_yscale('symlog', linthresh=1e-9)
440
-
ax.legend()
441
-
plt.tight_layout()
442
-
443
-
plt.savefig('radprof.jpg', bbox_inches='tight')
444
-
445
-
fig, ax = plt.subplots(1, 1)
446
-
uvgrid_shape = uvgrid_cont.shape
447
-
ax.imshow(np.abs(uvgrid_cont), origin='lower', norm=LogNorm(vmin=1e-12, vmax=1), extent=[-uvgrid_shape[0]//2, uvgrid_shape[0]//2, -uvgrid_shape[1]//2, uvgrid_shape[1]//2])
448
-
ax.set_title('UV grid')
449
-
plt.tight_layout()
450
-
#plt.savefig('uvgrid.jpg', bbox_inches='tight')
451
-
plt.show()
452
-
453
-
for pol in range(npol):
454
-
for chan in range(nchan):
455
-
dat = np.asarray(tb.getcolslice('DATA', [pol, chan], [pol, chan], [1, 1]))
456
-
flg = np.asarray(tb.getcolslice('FLAG', [pol, chan], [pol, chan], [1, 1]))
457
-
wgt = np.asarray(tb.getcolslice('WEIGHT_SPECTRUM', [pol, chan], [pol, chan], [1, 1]))
458
-
uvw = np.asarray(tb.getcol('UVW'))
459
-
460
-
print(f"pol {pol}, chan {chan}")
461
-
print(f"data shape {dat.shape}, flag shape {flg.shape}, weight shape {wgt.shape} uvw shape {uvw.shape}")
462
-
463
-
if dat.size == 0 or flg.size == 0 or wgt.size == 0:
464
-
print(f"Zero size array read. Skipping.")
465
-
continue
466
-
467
-
dat_grid = dat[0, 0, :].reshape([npoints, npoints])
468
-
flg_grid = flg[0, 0, :].reshape([npoints, npoints])
469
-
wgt_grid = wgt[0, 0, :].reshape([npoints, npoints])
470
-
471
-
# Flag the data as necessary
472
-
dat_grid = dat_grid * ~flg_grid
473
-
474
-
print("Old flagged points: ", np.count_nonzero(flg))
475
-
476
-
# Do the flagging and write back
477
-
flg_new = apply_flags(dat, flg, uvw, radial_mean, radial_mad, annuli, nsigma=sigma)
478
-
479
-
print("Flagged points: ", np.count_nonzero(flg_new))
480
-
481
-
tb.putcolslice('DATA', dat, [pol, chan], [pol, chan])
482
-
tb.putcolslice('FLAG', flg_new, [pol, chan], [pol, chan])
483
-
tb.putcolslice('WEIGHT_SPECTRUM', wgt, [pol, chan], [pol, chan])
484
6
485
-
tb.clearlocks()
486
-
tb.done()
487
-
def flag_radial_per_plane(binnedvis, sigma=5):
488
-
tb.open(binnedvis, nomodify=False)
489
-
msuvbinkey = tb.getkeyword('MSUVBIN')
490
-
nx = msuvbinkey['nx']
491
-
ny = msuvbinkey['ny']
492
-
if(nx != ny):
493
-
raise Exception('Do not deal with non square gridded vis')
494
-
npol = msuvbinkey['npol']
495
-
nchan = msuvbinkey['nchan']
496
-
497
-
for c in range(nchan):
498
-
dat = tb.getcolslice('DATA', [0,c], [npol-1,c], [1,1])
499
-
flg = tb.getcolslice('FLAG', [0, c], [npol-1, c], [1, 1])
500
-
wgt = tb.getcolslice('WEIGHT_SPECTRUM', [0, c], [npol-1, c], [1, 1])
501
-
#########
502
-
of = np.sum(flg[:,0,:])
503
-
print (f'BEFORE chan {c} number of unflagged points: {nx*nx*npol-of} max: {np.max(np.abs(dat[:,0,:]))}')
504
-
########
505
-
for k in range(npol):
506
-
if(doplot):
507
-
pl.clf()
508
-
ax1=pl.subplot(211)
509
-
ax1.set_title(f'radial mean Amp and fit for chan {c} and pol {k} ')
510
-
a = dat[k, 0, :].reshape([nx, nx])
511
-
f = flg[k, 0, :].reshape([nx, nx])
512
-
w = wgt[k, 0, :].reshape([nx, nx])
513
-
calc_radial_profile_and_fit(a,w,f,sigma)
514
-
if(doplot):
515
-
#pl.show()
516
-
pl.savefig(f'rad_{binnedvis}_c{c}_p{k}.jpg')
517
-
#input("Press Enter to continue...")
518
-
#########
519
-
of = np.sum(flg[:,0,:])
520
-
print (f'AFTER chan {c} number of unflagged points: {nx*nx*npol-of} max: {np.max(np.abs(dat[:,0,:]))}')
521
-
print ('=======================================================================')
522
-
########
523
-
if(not dryrun):
524
-
tb.putcolslice('FLAG', flg, [0,c], [npol-1,c], [1,1])
525
-
tb.putcolslice('DATA', dat, [0,c], [npol-1,c], [1,1])
526
-
tb.putcolslice('WEIGHT_SPECTRUM', wgt, [0, c], [npol-1, c], [1, 1])
527
-
528
-
tb.done()
529
7
530
8
@saveparams2last(multibackup=True)
531
9
def msuvbinflag(
532
10
######## binned ms
533
11
binnedvis, #= 'thegridded_from_msuvbin.ms',
534
12
######## flagging algorithm option
535
13
method, #= 'radial','RegionalMean','Gradient'
536
14
sigma, #=5
537
15
sizeRegion, #=20
538
16
radius, #=1.0
539
-
ignorPoint, #= True
540
17
doplot=False
541
18
) -> None:
542
19
flagger=UVGridFlag(binnedvis, doplot)
543
-
if method == 'radial_mean':
20
+
if method == 'radial_mean_annular':
544
21
tic = time.time()
545
22
flagger.flagViaBin_radial(sigma=sigma)
546
23
toc = time.time()
547
24
print("MSuvbinflag radial method Running time :", toc - tic)
548
25
elif method == 'radial_per_plane':
549
-
pdb.set_trace()
550
26
tic = time.time()
551
27
flagger.flag_radial_per_plane(sigma=sigma)
552
28
toc = time.time()
553
29
print("MSuvbinflag radial method Running time :", toc - tic)
554
30
elif method == 'regionalMean':
555
31
print("No bueno")
556
32
#tic = time.time()
557
33
#flagViaBin_regionalMean(binnedvis=binnedvis, sizeRegion=20, sigma=5,ignorPoint=True)
558
34
#toc = time.time()
559
35
#print("MSuvbinflag other method Running time:", toc - tic)