Commit eb0e5af9 authored by EPINAT Benoit's avatar EPINAT Benoit
Browse files

Unit handling in fits headers for output

parent 9d57e78e
......@@ -295,7 +295,7 @@ def muse_whitelight_image(cube, hdr, filename):
"""
data = np.sum(cube, axis=0)
# Suppressing header keywords about spectral axis
hdrw = copy.deepcopy(hdr)
if 'CTYPE3' in hdr.keys():
del hdrw['CTYPE3']
......@@ -315,8 +315,22 @@ def muse_whitelight_image(cube, hdr, filename):
del hdrw['CD2_3']
if 'CD1_3' in hdr.keys():
del hdrw['CD1_3']
wlh = writedata(data, hdrw, filename)
return wlh
if 'CDELT3' in hdr.keys():
del hdrw['CDELT3']
# Recovery of important information from header
if 'BUNIT' in hdr.keys():
bunit = hdr['BUNIT']
else:
bunit = ''
wave, lconvfac, sclp = waveindgen(hdr)
# We put a meaningful quantity and unit (integration of signal)
data = np.sum(cube, axis=0) * lconvfac * sclp
hdrw['BUNIT'] = '{} {}'.format(bunit, 'Angstrom')
wlih = writedata(data, hdrw, filename)
return wlih
def cutcube(cube, hdr, conf):
......@@ -623,6 +637,8 @@ def waveindgen(hdr):
lconvfac = 1e1
if cunit.strip().lower() in {'centimeters', 'cm'}:
lconvfac = 1e8
if cunit.strip().lower() in {'millimeters', 'mm'}:
lconvfac = 1e7
if 'CD3_3' in hdr.keys():
sclp = hdr['CD3_3']
elif 'CDELT3' in hdr.keys():
......@@ -828,7 +844,7 @@ def buildmaps(conf, cube, hdr, var=None, plot=False, debug=False, free_ratio=Fal
#els.append(line(conf['EXTRAL'][i][5], float(conf['EXTRAL'][i][0]), fit=True, ref=conf['EXTRAL'][i][1], th=float(conf['EXTRAL'][i][2]), low=float(conf['EXTRAL'][i][3]), up=float(conf['EXTRAL'][i][4])))
#Construction of wavelength index
wave0, lconvfac, sclp = waveindgen(hdr) # wavelength array in Angstrom and convertion factor
wave0, lconvfac, sclp = waveindgen(hdr) # wavelength array in Angstrom and conversion factor
#Condition on the wavelength to be within the redshift cuts
#XXX try zmin > zmax?
......@@ -1218,6 +1234,10 @@ def buildmaps(conf, cube, hdr, var=None, plot=False, debug=False, free_ratio=Fal
hh = writedata(econtcube, hdr, conf['OUTPUT'] + '_econtcube.fits')
#Images
if 'BUNIT' in hdr.keys():
bunit = hdr['BUNIT']
else:
bunit = ''
if 'CTYPE3' in hdr.keys():
del hdr['CTYPE3']
if 'CRVAL3' in hdr.keys():
......@@ -1247,39 +1267,49 @@ def buildmaps(conf, cube, hdr, var=None, plot=False, debug=False, free_ratio=Fal
#paramscube
#perrorcube
hdr['BUNIT'] = ''
hh = writedata(dofmap, hdr, conf['OUTPUT'] + '_dof' + suff + '.fits')
hh = writedata(fnormmap, hdr, conf['OUTPUT'] + '_fnorm' + suff + '.fits')
hh = writedata(statusmap, hdr, conf['OUTPUT'] + '_status' + suff + '.fits')
hh = writedata(chi2map, hdr, conf['OUTPUT'] + '_chi2' + suff + '.fits')
hh = writedata(dispcontmap, hdr, conf['OUTPUT'] + '_dispcont' + suff + '.fits')
hh = writedata(snrmap, hdr, conf['OUTPUT'] + '_snr' + suff + '.fits')
hdr['BUNIT'] = bunit
hh = writedata(dispcontmap, hdr, conf['OUTPUT'] + '_dispcont' + suff + '.fits')
for comp in range(multi_comp):
if multi_comp == 1:
suffixe = suff
else:
suffixe = suff + '_comp%i' % (comp + 1)
hdr['BUNIT'] = ''
hh = writedata(zmap[comp * els.index + comind, :, :], hdr, conf['OUTPUT'] + '_z' + suffixe + '.fits')
hh = writedata(ezmap[comp * els.index + comind, :, :], hdr, conf['OUTPUT'] + '_ez' + suffixe + '.fits')
hdr['BUNIT'] = 'KM/S'
hh = writedata(rvmap[comp * els.index + comind, :, :], hdr, conf['OUTPUT'] + '_vel' + suffixe + '.fits')
hh = writedata(ervmap[comp * els.index + comind, :, :], hdr, conf['OUTPUT'] + '_evel' + suffixe + '.fits')
for i in range(els.index):
if (not linefit[i]):
continue
hdr['BUNIT'] = 'Angstrom'
hh = writedata(wavemaps[comp * els.index + i, :, :], hdr, conf['OUTPUT'] + '_wave' + suffixe + '_' + linename[i] + '.fits')
hh = writedata(ewavemaps[comp * els.index + i, :, :], hdr, conf['OUTPUT'] + '_ewave' + suffixe + '_' + linename[i] + '.fits')
hh = writedata(wavedispmaps[comp * els.index + i, :, :], hdr, conf['OUTPUT'] + '_wavedisp' + suffixe + '_' + linename[i] + '.fits')
hh = writedata(ewavedispmaps[comp * els.index + i, :, :], hdr, conf['OUTPUT'] + '_ewavedisp' + suffixe + '_' + linename[i] + '.fits')
hdr['BUNIT'] = bunit
hh = writedata(intmaps[comp * els.index + i, :, :], hdr, conf['OUTPUT'] + '_int' + suffixe + '_' + linename[i] + '.fits')
hh = writedata(eintmaps[comp * els.index + i, :, :], hdr, conf['OUTPUT'] + '_eint' + suffixe + '_' + linename[i] + '.fits')
hdr['BUNIT'] = '{} {}'.format(bunit, 'Angstrom')
hh = writedata(fluxmaps[comp * els.index + i, :, :], hdr, conf['OUTPUT'] + '_flux' + suffixe + '_' + linename[i] + '.fits')
hh = writedata(efluxmaps[comp * els.index + i, :, :], hdr, conf['OUTPUT'] + '_eflux' + suffixe + '_' + linename[i] + '.fits')
hdr['BUNIT'] = ''
hh = writedata(snrmaps[comp * els.index + i, :, :], hdr, conf['OUTPUT'] + '_snr' + suffixe + '_' + linename[i] + '.fits')
hdr['BUNIT'] = 'KM/S'
if not(conf['COMMW']):
hh = writedata(dispmaps[comp * els.index + i, :, :], hdr, conf['OUTPUT'] + '_disp' + suffixe + '_' + linename[i] + '.fits')
hh = writedata(edispmaps[comp * els.index + i, :, :], hdr, conf['OUTPUT'] + '_edisp' + suffixe + '_' + linename[i] + '.fits')
if conf['COMMW']:
hdr['BUNIT'] = 'KM/S'
hh = writedata(dispmaps[comp * els.index + comind, :, :], hdr, conf['OUTPUT'] + '_disp' + suffixe + '.fits')
hh = writedata(edispmaps[comp * els.index + comind, :, :], hdr, conf['OUTPUT'] + '_edisp' + suffixe + '.fits')
......@@ -1351,7 +1381,7 @@ def camel(filename, plot=False, debug=False, free_ratio=False, multi_comp=1, mul
if var.ndim == 3:
var = remove_nanvar(var)
varhdr = writedata(var, varhdr, conf['OUTPUT'] + '_var_cut.fits')
wlh = muse_whitelight_image(cube, hdr, conf['OUTPUT'] + '_cut_whitelight.fits')
wlih = muse_whitelight_image(cube, hdr, conf['OUTPUT'] + '_cut_whitelight.fits')
if cut & (not(clip)) & (not(smooth)):
return
......@@ -1360,7 +1390,7 @@ def camel(filename, plot=False, debug=False, free_ratio=False, multi_comp=1, mul
logger.info('Performing %3.1f-sigma clipping using box of size %i. Number of pass: %i' % (conf['SCLIP'], conf['XYCLIP'], conf['NCLIP']))
cube, hdr = clipcube(cube, hdr, xy=conf['XYCLIP'], sclip=conf['SCLIP'], npass=conf['NCLIP'])
hdr = writedata(cube, hdr, conf['OUTPUT'] + '_cube_cut_clip.fits')
wlh = muse_whitelight_image(cube, hdr, conf['OUTPUT'] + '_cut_clip_whitelight.fits')
wlih = muse_whitelight_image(cube, hdr, conf['OUTPUT'] + '_cut_clip_whitelight.fits')
if clip & (not(smooth)):
return
......@@ -1390,7 +1420,7 @@ def camel(filename, plot=False, debug=False, free_ratio=False, multi_comp=1, mul
## 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')
wlih = muse_whitelight_image(cube, hdr, conf['OUTPUT'] + '_whitelight.fits')
if smooth:
return
......@@ -1496,7 +1526,7 @@ def main(argv):
parser.add_argument('--SCLIP', action="store", dest="SCLIP", default=None, type=float, help="sigma-clipping threshold to clean the cube")
parser.add_argument('--XYCLIP', action="store", dest="XYCLIP", default=None, type=int, help="width of the box to compute the median for the sigma-clipping")
parser.add_argument('--NCLIP', action="store", dest="NCLIP", default=None, type=int, help="number of pass for clipping")
parser.add_argument('--LSF', action="store", dest="SPSF", default=None, type=float, help="spectral PSF width in angstroms (dispersion of the gaussian)")
parser.add_argument('--LSF', action="store", dest="SPSF", default=None, type=float, help="spectral PSF width in angstroms (dispersion of the Gaussian)")
parser.add_argument('--WSMOOTH', action="store", dest="WSMOOTH", default=None, type=float, help="width for Gaussian spectral smoothing in pixel") # options? gaussian uniquement?
parser.add_argument('--SSMOOTH', action="store", dest="SSMOOTH", default=None, type=float, help="width for Gaussian spatial smoothing in pixel")
parser.add_argument('--XMIN', action="store", dest="XMIN", default=None, type=int, help="lower cut window in x direction of the cube in pixel")
......
File mode changed from 100644 to 100755
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