ddf-pipeline
ddf-pipeline copied to clipboard
random pixels in mosaic
Hi, the mosaic image from mosaic.py contains a small number of pixels outside of the beam factor used for blanking which are not set to NaN. While this is not an issue for LoTSS where the mosaic is cut out again, this can be annoying when working with a general mosaic image and e.g. a source finder.
I suspect the issue here might be that the image data contains very small numbers. The beam factor is computed as the ratio of apparent to intrinsic, and if both the intrinsic and apparent are ~0, then this is a div by zero. However, this is just a guess, I did not explicitly check this.
The relevant lines seem to be here.
The solution might be to filter out local outliers and replace them by the local mean - this was kinda slow for me with basic scipy.ndimage, a faster (~1 min) solution was to explicitly fit a 2d Gaussian as primary beam to the app/int ratio, sample this Gaussian and then cut on the samples.
app[i].data=np.divide(app[i].data,hdus[i].data)
from astropy.modeling import functional_models, fitting
G2d = functional_models.Gaussian2D(1., len(app[i].data)/2, len(app[i].data[0])/2,len(app[i].data)/4,
len(app[i].data[0])/4, fixed={'amplitude': True})
fitter = fitting.LevMarLSQFitter()
yp, xp = np.indices(app[i].data.shape)
fit_G2d = fitter(G2d, xp, yp, app[i].data)
beam_smooth = fit_G2d(xp, yp)
app[i].data[beam_smooth<threshold]=0 # could also use the beam_smooth values instead of app/int ratio
Not sure if this should be added to mosaic.py or if this is the best way to do this, I just wanted to leave this here as reference.
Scipy's binary erosion function might be a faster way to remove them at the expense of shrinking the perimeter of the mosaic by a pixel.
Thanks, the thing I'd be worried about with this solution is the assumption that the beam is well characterized by a Gaussian, but I expect you're right about the nature of the problem. I used to have some code that did this reasonably efficiently (awimager used to have the same problem) so let me see if I can dig it out. Or we could use Ian's solution...
Hey, thanks for the quick reply. Yeah, I think my solution has a few problems, Ian's is indeed fast and stable, this might be the way to go.