def spatial_smooth(self, npix, output="smoothed.fits", test=False, **kwargs):
"""Applies Gaussian filter of std=npix in both spatial directions
and writes it to disk as a new MUSE Cube.
Notes: the STAT cube is not touched.
Parameters
----------
npix : int
Std of Gaussian kernel in spaxel units.
output : str, optional
Name of the output file
test : bool, optional
Whether to check for flux being conserved
**kwargs are passed down to scipy.ndimage.gaussian_filter()
Return
------
Writes a new file to disk.
"""
if not isinstance(npix, int):
raise ValueError("npix must be integer.")
cube_new = copy.deepcopy(self.cube)
ntot = len(self.cube)
for wv_ii in range(ntot):
print('{}/{}'.format(wv_ii + 1, ntot))
image_aux = self.cube[wv_ii, :, :]
smooth_ii = ma.MaskedArray(ndimage.gaussian_filter(image_aux, sigma=npix, **kwargs))
smooth_ii.mask = image_aux.mask | np.isnan(smooth_ii)
# test the fluxes are conserved
if test:
gd_pix = ~smooth_ii.mask
try:
med_1 = np.nansum(smooth_ii[gd_pix])
med_2 = np.nansum(image_aux[gd_pix])
print(med_1, med_2, (med_1 - med_2) / med_1)
np.testing.assert_allclose(med_1, med_2, decimal=4)
except AssertionError:
import pdb
pdb.set_trace()
cube_new[wv_ii, :, :] = smooth_ii
# import pdb; pdb.set_trace()
hdulist = fits.open(self.filename)
hdulist[1].data = cube_new.data
prihdr = hdulist[0].header
comment = 'Spatially smoothed with a Gaussian kernel of sigma={} spaxels (by MuseCube)'.format(npix)
print(comment)
prihdr['history'] = comment
hdulist.writeto(output, clobber=True)
print("MuseCube: new smoothed cube written to {}".format(output))
评论列表
文章目录
