Commit 85bc50da authored by EPINAT Benoit's avatar EPINAT Benoit
Browse files

Correction of variance weighting with smoothing

parent 845a8e78
......@@ -380,16 +380,20 @@ def remove_nanvar(var):
"""
cnan = np.isnan(var)
nanmap = cnan.prod(axis=0)
indnan = np.where(nanmap == 1)
indnum = np.where(nanmap == 0)
qvar = var[:, indnum[0], indnum[1]].reshape(var.shape[0], np.size(indnum[0]))
medianvar = np.median(qvar, axis=1)
cnan = np.isnan(var) # True where there are NaN
#nanmap = cnan.prod(axis=0) # 0 if there is at least one value
#indnan = np.where(nanmap == 1)
#indnum = np.where(nanmap == 0)
ncnan = np.logical_not(cnan) # True where there are values
nanmap = ncnan.prod(axis=0) # 0 if there is at least one NaN
indnan = np.where(nanmap == 0) # indicates the pixels where there is at least one NaN
indnum = np.where(nanmap == 1) # indicates the pixels where there are no NaN
qvar = var[:, indnum[0], indnum[1]].reshape(var.shape[0], np.size(indnum[0])) # spectra with no NaN
medianvar = np.median(qvar, axis=1) # median variance spectrum
if np.size(indnan[0]) == 0:
logger.info('Variance cube contains NaN and has been corrected')
for ind in range(np.size(indnan[0])):
var[:, indnan[0][ind], indnan[1][ind]] = medianvar
var[:, indnan[0][ind], indnan[1][ind]] = medianvar # If there is one NaN, the whole variance spectrum is replaced by the median variance spectrum
return var
......@@ -417,7 +421,7 @@ def clipcube(cube, hdr, xy=3, sclip=3, npass=10):
count = 0
ok = False
#while count < npass | ok:
while (count < npass) & not(ok):
while (count < npass) & (not(ok)):
count += 1
im1 = ndimage.median_filter(im, xy)
diff = im - im1
......@@ -454,6 +458,34 @@ def spatialsmooth(cube, hdr, fwhm):
return cube1, hdr
def var_spatialsmooth(cube, hdr, fwhm):
"""This function performs a Gaussian spatial smoothing on the cube
Parameters
----------
cube: array
input datacube
hdr: fits header
input data cube header
fwhm: float
full width half maximum of the 2D Gaussian kernel
"""
cube1 = np.zeros(cube.shape, dtype='>f4') # initialization of output cube
sigma = fwhm / (2. * np.sqrt(2. * np.log(2.)))
# For variance, we have to convolve with the square of the kernel.
# In 2D, it is equivalent to convolve by a kernel of width sigma2 = sigma / np.sqrt(2), hence FWHM2 = FWHM/np.sqrt
# and to renormalise the amplitude by multiplying by a factor fac = 1 / (2 * sigma) / np.sqrt(np.pi)
sigma2 = sigma / np.sqrt(2)
fac = 1 / (4 * sigma**2 * np.pi)
for z in range(cube.shape[0]): # loop on the z range
cube1[z, :] = ndimage.gaussian_filter(cube[z, :], sigma)
cube1 *= fac
hdr.append(('SSMOOTH', fwhm, 'FWHM in pixels of 2D Gaussian spatial smoothing'))
return cube1, hdr
def spectralsmooth(cube, hdr, fwhm):
"""This function performs a Gaussian spectral smoothing on the cube
......@@ -477,6 +509,36 @@ def spectralsmooth(cube, hdr, fwhm):
return cube1, hdr
def var_spectralsmooth(cube, hdr, fwhm):
"""This function performs a Gaussian spectral smoothing on the cube
Parameters
----------
cube: array
input datacube
hdr: fits header
input data cube header
fwhm: float
full width half maximum of the 1D Gaussian kernel
"""
cube1 = np.zeros(cube.shape, dtype='>f4') # initialization of output cube
sigma = fwhm / (2. * np.sqrt(2. * np.log(2.)))
# For variance, we have to convolve with the square of the kernel.
# In 1D, it is equivalent to convolve by a kernel of width sigma2 = sigma / np.sqrt(2), hence FWHM2 = FWHM/np.sqrt
# and to renormalise the amplitude by multiplying by a factor fac = 1 / (2 * sigma) / np.sqrt(np.pi)
sigma2 = sigma / np.sqrt(2)
fac = 1 / (2 * sigma * np.sqrt(np.pi))
for x in range(cube.shape[2]): # loop on the x range
for y in range(cube.shape[1]): # loop on the y range
cube1[:, y, x] = ndimage.gaussian_filter(cube[:, y, x], sigma2)
cube1 *= fac
hdr.append(('WSMOOTH', fwhm, 'FWHM in pixels of Gaussian spectral smoothing'))
# XXX add Hanning, Gauss, ...
return cube1, hdr
class line:
"""This class is the basic for defining a line. A line is defined by its name, its wavelength and the reference line to which it is attached if the ratio has to be constrained.
......@@ -1310,9 +1372,10 @@ def camel(filename, plot=False, debug=False, free_ratio=False, multi_comp=1, mul
hdr = writedata(cube, hdr, conf['OUTPUT'] + '_cube.fits')
if var is not None:
if var.ndim == 3:
var, varhdr = spectralsmooth(var, varhdr, conf['WSMOOTH'])
# variance is divided by the integral of the kernel to account for the variance lowering compared to the signal when smoothed
var /= ((conf['WSMOOTH'] / (2. * np.sqrt(2. * np.log(2.)))) * np.sqrt(2 * np.pi))
#var, varhdr = spectralsmooth(var, varhdr, conf['WSMOOTH'])
var, varhdr = var_spectralsmooth(var, varhdr, conf['WSMOOTH'])
## variance is divided by the integral of the kernel to account for the variance lowering compared to the signal when smoothed
#var /= ((conf['WSMOOTH'] / (2. * np.sqrt(2. * np.log(2.)))) * np.sqrt(2 * np.pi))
varhdr = writedata(var, varhdr, conf['OUTPUT'] + '_var.fits')
if (conf['SSMOOTH'] != 0) & (conf['SSMOOTH'] is not None):
......@@ -1322,9 +1385,10 @@ def camel(filename, plot=False, debug=False, free_ratio=False, multi_comp=1, mul
hdr = writedata(cube, hdr, conf['OUTPUT'] + '_cube.fits')
if var is not None:
if var.ndim == 3:
var, varhdr = spatialsmooth(var, varhdr, conf['SSMOOTH'])
# variance is divided by the integral of the kernel to account for the variance lowering compared to the signal when smoothed
var /= ((conf['SSMOOTH'] / (2. * np.sqrt(2. * np.log(2.)))) ** 2 * 2 * np.pi)
#var, varhdr = spatialsmooth(var, varhdr, conf['SSMOOTH'])
var, varhdr = var_spatialsmooth(var, varhdr, conf['SSMOOTH'])
## variance is divided by the integral of the kernel to account for the variance lowering compared to the signal when smoothed
#var /= ((conf['SSMOOTH'] / (2. * np.sqrt(2. * np.log(2.)))) ** 2 * 2 * np.pi)
varhdr = writedata(var, varhdr, conf['OUTPUT'] + '_var.fits')
wlh = muse_whitelight_image(cube, hdr, conf['OUTPUT'] + '_whitelight.fits')
if smooth:
......
......@@ -135,7 +135,8 @@ def create_config(path, cubefile, varfile, catfile, lines, suffixeout, commonw=T
hdu = fits.PrimaryHDU(data=wlim, header=hdr)
hdul1 = fits.HDUList(hdu)
hdr = hdul1[0].header
w = wcs.WCS(hdr, hdul1)
#w = wcs.WCS(hdr, hdul1)
w = wcs.WCS(hdr, hdul1, naxis=2) # In order that it also works for MXDF
#w = wcs.WCS(hdr[:28], hdul1)
logger.debug(' Naxis %d ' % (w.naxis))
......
......@@ -318,6 +318,36 @@ def spectralsmooth(data, hdr, fwhm):
return data1, hdr
def var_spectralsmooth(cube, hdr, fwhm):
"""This function performs a Gaussian spectral smoothing on the cube
Parameters
----------
cube: array
input datacube
hdr: fits header
input data cube header
fwhm: float
full width half maximum of the 1D Gaussian kernel
"""
cube1 = np.zeros(cube.shape, dtype='>f4') # initialization of output cube
sigma = fwhm / (2. * np.sqrt(2. * np.log(2.)))
# For variance, we have to convolve with the square of the kernel.
# In 1D, it is equivalent to convolve by a kernel of width sigma2 = sigma / np.sqrt(2), hence FWHM2 = FWHM/np.sqrt
# and to renormalise the amplitude by multiplying by a factor fac = 1 / (2 * sigma) / np.sqrt(np.pi)
sigma2 = sigma / np.sqrt(2)
fac = 1 / (2 * sigma) / np.sqrt(np.pi)
for x in range(cube.shape[2]): # loop on the x range
for y in range(cube.shape[1]): # loop on the y range
cube1[:, y, x] = ndimage.gaussian_filter(cube[:, y, x], sigma2)
cube1 *= fac
hdr.append(('WSMOOTH', fwhm, 'FWHM in pixels of Gaussian spectral smoothing'))
# XXX add Hanning, Gauss, ...
return cube1, hdr
class line:
"""This class is the basic for defining a line. A line is defined by its name, its wavelength and the reference line to which it is attached if the ratio has to be constrained.
......@@ -1050,9 +1080,10 @@ def samel(filename, plot=False, debug=False, free_ratio=False, multi_comp=1, mul
conf['OUTPUT'] = conf['OUTPUT'] + '_wsmooth'
hdr = writedata(data, hdr, conf['OUTPUT'] + '_spectrum.fits')
if var is not None:
var, varhdr = spectralsmooth(var, varhdr, conf['WSMOOTH'])
# variance is divided by the integral of the kernel to account for the variance lowering compared to the signal when smoothed
var /= ((conf['WSMOOTH'] / (2. * np.sqrt(2. * np.log(2.)))) * np.sqrt(2 * np.pi))
#var, varhdr = spectralsmooth(var, varhdr, conf['WSMOOTH'])
## variance is divided by the integral of the kernel to account for the variance lowering compared to the signal when smoothed
#var /= ((conf['WSMOOTH'] / (2. * np.sqrt(2. * np.log(2.)))) * np.sqrt(2 * np.pi))
var, varhdr = var_spectralsmooth(var, varhdr, conf['WSMOOTH'])
varhdr = writedata(var, varhdr, conf['OUTPUT'] + '_var.fits')
if smooth:
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment