casa6
Public
Commits
| 151 151 | nbin = 30 |
| 152 152 | |
| 153 153 | npixx, npixy = uvgrid.shape[0], uvgrid.shape[1] |
| 154 154 | cx, cy = npixx//2, npixy//2 |
| 155 155 | |
| 156 156 | uvlen_m_grid = uvlen_m.reshape([npixx, npixy]) |
| 157 157 | |
| 158 158 | # Generate grid of radii |
| 159 159 | x = np.arange(npixx) - cx |
| 160 160 | y = np.arange(npixy) - cy |
| 161 - | rad = np.sqrt(x**2 + y**2).astype(int) |
| 161 + | #rad = np.sqrt(x**2 + y**2).astype(int) |
| 162 162 | |
| 163 163 | yy, xx = np.meshgrid(x, y) |
| 164 - | radgrid = np.sqrt(xx**2 + yy**2).astype(int) |
| 164 + | #radgrid = np.sqrt(xx**2 + yy**2).astype(int) |
| 165 165 | |
| 166 166 | # Create log-spaced annuli to account for reducing UV coverage with radius |
| 167 167 | # Minimum annulus is 5px |
| 168 168 | annuli = np.logspace(0, np.log10(uvlen_m_grid.max()), nbin) |
| 169 169 | annuli = np.round(annuli).astype(int) |
| 170 170 | |
| 171 171 | radial_mean = np.zeros(nbin) |
| 172 172 | radial_mad = np.zeros(nbin) |
| 173 173 | |
| 174 174 | ann_min = 0 |
| 259 259 | radamp[rval]+=absval*wgtgrid[k,j] |
| 260 260 | radamp2[rval]+=absval*absval*wgtgrid[k,j] |
| 261 261 | radwght[rval]+=wgtgrid[k,j] |
| 262 262 | if(np.max(radwght)==0.0): |
| 263 263 | #empty channel |
| 264 264 | return |
| 265 265 | radamp[radwght!=0]=radamp[radwght!=0]/radwght[radwght!=0] |
| 266 266 | radamp2[radwght!=0]=radamp2[radwght!=0]/radwght[radwght!=0] |
| 267 267 | maxsenspos=np.argmax(radwght) |
| 268 268 | #normalize radweight |
| 269 - | normrdwght=radwght/np.max(radwght) |
| 269 + | #normrdwght=radwght/np.max(radwght) |
| 270 270 | sig=np.sqrt(np.abs(radamp2-radamp*radamp)) |
| 271 271 | #nescale relative sigmas by number of weights att the point |
| 272 - | medsig=np.median(sig[sig !=0]) |
| 272 + | #medsig=np.median(sig[sig !=0]) |
| 273 273 | sigtouse=sig[maxsenspos] |
| 274 274 | #sig[normrdwght!=0]=sig[normrdwght!=0]/normrdwght[normrdwght!=0] |
| 275 275 | xvalnz=xval[(sig !=0.0) & (radamp !=0)] |
| 276 276 | radampnz=radamp[(sig !=0) & (radamp !=0)] |
| 277 277 | fitnz=curve_fit(self.resid_cinco, xvalnz, radampnz) |
| 278 278 | ### |
| 279 279 | #print('min max of sig and max sens one', np.min(sig), np.max(sig), sigtouse) |
| 280 280 | signz=sig[sig !=0.0] |
| 281 281 | sig=np.interp(xval, xvalnz, signz) |
| 282 282 | #print( 'corvar ', fitnz[1]) |
| 364 364 | uvgrid_cont = np.zeros((npoints, npoints), dtype=np.complex128) |
| 365 365 | wgtgrid_cont = np.zeros((npoints, npoints), dtype=np.float64) |
| 366 366 | |
| 367 367 | for pol in range(npol): |
| 368 368 | for chan in range(nchan): |
| 369 369 | dat = tb.getcolslice('DATA', [pol, chan], [pol, chan], [1, 1]) |
| 370 370 | flg = tb.getcolslice('FLAG', [pol, chan], [pol, chan], [1, 1]) |
| 371 371 | wgt = tb.getcolslice('WEIGHT_SPECTRUM', [pol, chan], [pol, chan], [1, 1]) |
| 372 372 | |
| 373 373 | if dat.size == 0 or flg.size == 0 or wgt.size == 0: |
| 374 - | casalog.post(f"Zero size array read. Skipping.", |
| 374 + | casalog.post("Zero size array read. Skipping.", |
| 375 375 | 'WARN', |
| 376 376 | 'task_msuvbinflag') |
| 377 377 | continue |
| 378 378 | |
| 379 379 | dat_grid = dat[0, 0, :].reshape([npoints, npoints]) |
| 380 380 | flg_grid = flg[0, 0, :].reshape([npoints, npoints]) |
| 381 381 | wgt_grid = wgt[0, 0, :].reshape([npoints, npoints]) |
| 382 382 | |
| 383 383 | # Flag the data as necessary |
| 384 384 | dat_grid = dat_grid * ~flg_grid |
| 450 450 | #plt.savefig('uvgrid.jpg', bbox_inches='tight') |
| 451 451 | plt.show() |
| 452 452 | |
| 453 453 | for pol in range(npol): |
| 454 454 | for chan in range(nchan): |
| 455 455 | dat = np.asarray(tb.getcolslice('DATA', [pol, chan], [pol, chan], [1, 1])) |
| 456 456 | flg = np.asarray(tb.getcolslice('FLAG', [pol, chan], [pol, chan], [1, 1])) |
| 457 457 | wgt = np.asarray(tb.getcolslice('WEIGHT_SPECTRUM', [pol, chan], [pol, chan], [1, 1])) |
| 458 458 | |
| 459 459 | if dat.size == 0 or flg.size == 0 or wgt.size == 0: |
| 460 - | casalog.post(f"Zero size array read. Skipping.", |
| 460 + | casalog.post("Zero size array read. Skipping.", |
| 461 461 | 'WARN', |
| 462 462 | 'task_msuvbinflag') |
| 463 463 | continue |
| 464 464 | |
| 465 465 | # Do the flagging and write back |
| 466 466 | flg_new = self.apply_flags(dat, flg, uvlen_m, radial_mean, radial_mad, annuli, nsigma=sigma) |
| 467 467 | |
| 468 468 | tb.putcolslice('DATA', dat, [pol, chan], [pol, chan]) |
| 469 469 | tb.putcolslice('FLAG', flg_new, [pol, chan], [pol, chan]) |
| 470 470 | tb.putcolslice('WEIGHT_SPECTRUM', wgt, [pol, chan], [pol, chan]) |
| 503 503 | a = dat[k, 0, :].reshape([nx, nx]) |
| 504 504 | f = flg[k, 0, :].reshape([nx, nx]) |
| 505 505 | w = wgt[k, 0, :].reshape([nx, nx]) |
| 506 506 | self.calc_radial_profile_and_fit(a,w,f,sigma) |
| 507 507 | if(self.doplot): |
| 508 508 | #pl.show() |
| 509 509 | pl.savefig(f'rad_{self.binnedvis}_c{c}_p{k}.jpg') |
| 510 510 | #input("Press Enter to continue...") |
| 511 511 | ######### |
| 512 512 | of = np.sum(flg[:,0,:]) |
| 513 - | casa.logpost('AFTER chan %d number of unflagged points: %d max: %f' % (c, nx*nx*npol-of, np.max(np.abs(dat[:,0,:])),), |
| 513 + | casalog.post('AFTER chan %d number of unflagged points: %d max: %f' % (c, nx*nx*npol-of, np.max(np.abs(dat[:,0,:])),), |
| 514 514 | 'DEBUG', |
| 515 515 | 'task_msuvbinflag') |
| 516 516 | #print (f'AFTER chan {c} number of unflagged points: {nx*nx*npol-of} max: {np.max(np.abs(dat[:,0,:]))}') |
| 517 517 | #print ('=======================================================================') |
| 518 518 | ######## |
| 519 519 | if(not self.dryrun): |
| 520 520 | tb.putcolslice('FLAG', flg, [0,c], [npol-1,c], [1,1]) |
| 521 521 | tb.putcolslice('DATA', dat, [0,c], [npol-1,c], [1,1]) |
| 522 522 | tb.putcolslice('WEIGHT_SPECTRUM', wgt, [0, c], [npol-1, c], [1, 1]) |
| 523 523 | |
| 524 524 | tb.done() |
| 525 525 | ######################################################################################### |
| 526 526 | def flag_gradient(self) -> None: |
| 527 + | sigma=3 |
| 527 528 | tb.open(self.binnedvis, nomodify=False) |
| 528 529 | msuvbinkey = tb.getkeyword('MSUVBIN') |
| 529 530 | nx = msuvbinkey['nx'] |
| 530 531 | ny = msuvbinkey['ny'] |
| 531 532 | if(nx != ny): |
| 532 533 | raise Exception('Do not deal with non square gridded vis') |
| 533 534 | npol = msuvbinkey['npol'] |
| 534 535 | nchan = msuvbinkey['nchan'] |
| 535 536 | |
| 536 537 | for c in range(nchan): |
| 567 568 | tb.done() |
| 568 569 | |
| 569 570 | |
| 570 571 | |
| 571 572 | |
| 572 573 | |
| 573 574 | def locate_von(grid, radius=1, scale=0.3): |
| 574 575 | flagpoints = [] |
| 575 576 | gridpoints = [] |
| 576 577 | npoints = len(grid) |
| 577 - | gradients = gradient(10, 10, grid.data) |
| 578 + | gradients = np.array(np.gradient(grid.data)) |
| 578 579 | du = gradients[0] |
| 579 580 | dv = gradients[1] |
| 580 581 | th = np.sqrt(np.max(du)*np.max(du) + np.max(dv)*np.max(dv)) |
| 581 582 | print ('Max grad', th, np.max(du), np.max(dv)) |
| 582 583 | scale = th * scale |
| 583 584 | |
| 584 585 | for i in range(int(radius), int(npoints - radius)): |
| 585 586 | for j in range(int(radius), int(npoints - radius)): |
| 586 587 | if(grid.mask[i,j]==False): |
| 587 588 | du_up = du[i + radius][j] |