Commit bb630310 authored by bepinat's avatar bepinat
Browse files

Variance corrections, whitelight image creation, new keywords to only create...

Variance corrections, whitelight image creation, new keywords to only create cut/clip/smoothed cubes, changes in create_config.py to use command line
parent a6eab647
......@@ -232,6 +232,36 @@ def writedata(data, hdr, filename):
hdulist.writeto(filename, clobber=True)
return hdulist[0].header
def muse_whitelight_image(cube, hdr, filename):
"""This function creates and writes a white light image
Parameters
----------
data: ndarray
cube to be collapsed
hdr: fits header
fits header of the cube
filename: string
name of the output
"""
data = np.sum(cube, axis=0)
hdrw = copy.deepcopy(hdr)
if 'CTYPE3' in hdr.ascard.keys(): del hdrw['CTYPE3']
if 'CRVAL3' in hdr.ascard.keys(): del hdrw['CRVAL3']
if 'CRPIX3' in hdr.ascard.keys(): del hdrw['CRPIX3']
if 'CUNIT3' in hdr.ascard.keys(): del hdrw['CUNIT3']
if 'CD3_3' in hdr.ascard.keys(): del hdrw['CD3_3']
if 'CD3_2' in hdr.ascard.keys(): del hdrw['CD3_2']
if 'CD3_1' in hdr.ascard.keys(): del hdrw['CD3_1']
if 'CD2_3' in hdr.ascard.keys(): del hdrw['CD2_3']
if 'CD1_3' in hdr.ascard.keys(): del hdrw['CD1_3']
wlh = writedata(data, hdrw, filename)
return wlh
def cutcube(cube, hdr, conf):
"""This function cuts the cube according to the limits requested
......@@ -562,7 +592,11 @@ def buildmaps(conf, cube, hdr, var=None, plot=False, debug=False, free_ratio=Fal
if var == None: var = np.ones(cube.shape)
if var.ndim == 1: var = np.tile(var.reshape(var.size, 1, 1), (1, cube.shape[1], cube.shape[2]))
#if var.ndim == 3: weight = ((1. / var[indok,:,:]) / np.sum(1. / var[indok,:,:], axis=0)) * indok.size
if var.ndim == 3: weight = ((1. / var) / np.sum(1. / var, axis=0)) * var.shape[0]
# weight inverse variance
#if var.ndim == 3: weight = ((1. / var) / np.sum(1. / var, axis=0)) * var.shape[0]
# weight inverse standard deviation
#var += 0
if var.ndim == 3: weight = ((1. / np.sqrt(var)) / np.sum(1. / np.sqrt(var), axis=0)) * var.shape[0]
#Redshift step
zstep = (conf['DFIT'] * 1e3 / ct.c + zcosmo) / (1 - conf['DFIT'] * 1e3 / ct.c) - zcosmo
......@@ -800,7 +834,7 @@ def buildmaps(conf, cube, hdr, var=None, plot=False, debug=False, free_ratio=Fal
intmaps = np.zeros(np.append(els.index, cube.shape[1:]))
eintmaps = np.zeros(np.append(els.index, cube.shape[1:]))
fluxmaps = np.zeros(np.append(els.index, cube.shape[1:]))
efluxmaps = np.zeros(np.append(els.index, cube.shape[1:]))
efluxmaps = np.zeros(np.append(els.index, cube.shape[1:]))
snrmaps = np.zeros(np.append(els.index, cube.shape[1:]))
# dzmap and edzmap are intermediate products
......@@ -839,7 +873,21 @@ def buildmaps(conf, cube, hdr, var=None, plot=False, debug=False, free_ratio=Fal
efluxmaps[i,:,:] = np.sqrt(np.pi * 2) * np.sqrt((eintmaps[i,:,:] * np.sqrt((linewave[i] * dzmap) ** 2 + conf['SPSF'] ** 2)) ** 2 + (linewave[i] * dzmap / np.sqrt((linewave[i] * dzmap) ** 2 + conf['SPSF'] ** 2) * intmaps[i,:,:] * edzmap) ** 2)
snrmap = np.sum(weight[indok,:,:] * (modcube[indok,:,:] - contcube[indok,:,:]), axis = 0) * sclp * lconvfac / (np.sqrt(np.pi * 2) * np.sqrt((dzmap * np.median(wave) / (1 + zmap)) ** 2 + conf['SPSF'] ** 2) * np.sqrt( np.sum(weight[indok,:,:] * (cube[indok,:,:] - modcube[indok,:,:]) ** 2, axis=0) / np.sum(weight[indok,:,:], axis=0)))
# XXX Pour les raies individuelles, ne devrait-on pas calculer le flux dans la raie uniquement?
# XXX ajouter un cas pour lequel on ajoute un niveau seuil a la variance (e.g. sky spectrum pour lequel on veut rajouter un bruit de detecteur, ou si on estime que le niveau bas est trop bas pour que les weights soient realistes)
# case SNR weighted by inverse of variance / standard deviation (normalised)
#snrmap = np.sum(weight[indok,:,:] * (modcube[indok,:,:] - contcube[indok,:,:]), axis = 0) * sclp * lconvfac / (np.sqrt(np.pi * 2) * np.sqrt((dzmap * np.median(wave) / (1 + zmap)) ** 2 + conf['SPSF'] ** 2) * np.sqrt( np.sum(weight[indok,:,:] * (cube[indok,:,:] - modcube[indok,:,:]) ** 2, axis=0) / np.sum(weight[indok,:,:], axis=0)))
# case SNR weighted by inverse of variance / standard deviation (normalised) + normalised to number of elements
snrmap = np.sum(weight[indok,:,:] * (modcube[indok,:,:] - contcube[indok,:,:]), axis = 0) * sclp * lconvfac * np.size(indok) / (np.sum(weight[indok,:,:], axis=0)) /(np.sqrt(np.pi * 2) * np.sqrt((dzmap * np.median(wave) / (1 + zmap)) ** 2 + conf['SPSF'] ** 2) * np.sqrt( np.sum(weight[indok,:,:] * (cube[indok,:,:] - modcube[indok,:,:]) ** 2, axis=0) / np.sum(weight[indok,:,:], axis=0)))
# case SNR using standard deviation as estimate of noise (might be incorrect after smoothing?)
#snrmap = np.sum((modcube[indok,:,:] - contcube[indok,:,:]) / np.sqrt(var[indok,:,:]), axis = 0) * sclp * lconvfac / (np.sqrt(np.pi * 2) * np.sqrt((dzmap * np.median(wave) / (1 + zmap)) ** 2 + conf['SPSF'] ** 2)) # * np.sqrt( np.sum((cube[indok,:,:] - modcube[indok,:,:]) ** 2 / np.sqrt(var[indok,:,:]), axis=0) / np.sum(var[indok,:,:], axis=0)))
# case SNR without variance. No weight, noise estimated as the dispersion of the residuals
#snrmap = np.sum((modcube[indok,:,:] - contcube[indok,:,:]), axis = 0) * sclp * lconvfac / (np.sqrt(np.pi * 2) * np.sqrt((dzmap * np.median(wave) / (1 + zmap)) ** 2 + conf['SPSF'] ** 2) * np.sqrt( np.sum((cube[indok,:,:] - modcube[indok,:,:]) ** 2, axis=0) / np.size(indok)))
#Output writing
#Cubes
......@@ -895,7 +943,8 @@ def buildmaps(conf, cube, hdr, var=None, plot=False, debug=False, free_ratio=Fal
return rvmap
def camel(filename, plot=False, debug=False, free_ratio=False, multi_ext=False, config=None, ftol=1e-7, xtol=1e-7, gtol=1e-7, factor=10, maxiter=200):
def camel(filename, plot=False, debug=False, free_ratio=False, multi_ext=False, config=None, cut=False, clip=False, smooth=False, ftol=1e-7, xtol=1e-7, gtol=1e-7, factor=10, maxiter=200):
"""CAMEL (Cube Analysis: Moment maps of Emission Lines) computes moment maps
Parameters
......@@ -956,12 +1005,16 @@ def camel(filename, plot=False, debug=False, free_ratio=False, multi_ext=False,
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')
if cut & (not(clip)) & (not(smooth)): return
#Clipping the cube
if (conf['SCLIP'] != 0) & (conf['SCLIP'] != None):
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')
if clip & (not(smooth)): return
#Smoothing the cubes
if (conf['WSMOOTH'] != 0) & (conf['WSMOOTH'] != None):
......@@ -972,6 +1025,8 @@ def camel(filename, plot=False, debug=False, free_ratio=False, multi_ext=False,
if var != 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['SSMOOTH'] / (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'] != None):
......@@ -982,7 +1037,11 @@ def camel(filename, plot=False, debug=False, free_ratio=False, multi_ext=False,
if var != None:
if var.ndim == 3:
var, varhdr = spatialsmooth(var, varhdr, conf['SSMOOTH'])
varhdr = writedata(var, varhdr, conf['OUTPUT'] + '_var.fits')
# 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'] + '_ssmooth_whitelight.fits')
if smooth: return
#Creating the maps
rv = buildmaps(conf, cube, hdr, var=var, plot=plot, debug=debug, free_ratio=free_ratio, ftol=ftol, xtol=xtol, gtol=gtol, factor=factor, maxiter=maxiter)
......@@ -1089,6 +1148,9 @@ def main(argv):
parser.add_argument('--XMAX', action="store", dest="XMAX", default=None, type=int, help="upper cut window in y direction of the cube in pixel")
parser.add_argument('--YMAX', action="store", dest="YMAX", default=None, type=int, help="lower cut window in z direction of the cube in pixel")
parser.add_argument('--ZMAX', action="store", dest="ZMAX", default=None, type=int, help="upper cut window in z direction of the cube in pixel")
parser.add_argument('--cut', action="store_true", dest="cut", default=False, help="keword to stop after cutting the cube")
parser.add_argument('--clip', action="store_true", dest="clip", default=False, help="keword to stop after clipping the cube")
parser.add_argument('--smooth', action="store_true", dest="smooth", default=False, help="keword to stop after smoothing the cube")
#parser.add_argument('--MFIT', action="store", dest="MFIT", default=None)
#parser.add_argument('--THRES', action="store", dest="THRES", default=None, type=float)
......@@ -1105,7 +1167,7 @@ def main(argv):
opts = vars(options) # this is a dictionnary
opts['OII3726'] = opts['OII3729']
camel(options.filename, plot=options.plot, debug=options.debug, free_ratio=options.free_ratio, multi_ext=options.multi_ext, ftol=options.ftol, xtol=options.xtol, gtol=options.gtol, maxiter=options.maxiter, factor=options.factor, config=opts)
camel(options.filename, plot=options.plot, debug=options.debug, free_ratio=options.free_ratio, multi_ext=options.multi_ext, ftol=options.ftol, xtol=options.xtol, gtol=options.gtol, maxiter=options.maxiter, factor=options.factor, config=opts, cut=options.cut, clip=options.clip, smooth=options.smooth)
return
if __name__ == "__main__":
......
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
import os, sys
import argparse
import ipdb
import numpy as np
import pyfits as pf
import astropy.wcs as wcs
......@@ -34,7 +37,7 @@ def create_config(path, cubefile, varfile, catfile, lines, suffixeout, commonw=T
specifies if a common width has to be used for the various lines
dz: float
interval of redshift around the systemic redshift (from the catalogue) to search for lines and for cutting the cube spectrally
dxyz: int
dxy: int
spatial extent around the object for cutting the cube (in pixels)
initw: float
initial guess for the line width (dispersion) in km/s
......@@ -58,7 +61,8 @@ def create_config(path, cubefile, varfile, catfile, lines, suffixeout, commonw=T
gaussian spatial smooth FWHM in pixels
'''
cat = np.genfromtxt(path + catfile, names=True)
#ipdb.set_trace()
cat = np.genfromtxt(path + catfile, names=True, case_sensitive=False)
try:
hdul = pf.open(path + cubefile)
......@@ -93,6 +97,7 @@ def create_config(path, cubefile, varfile, catfile, lines, suffixeout, commonw=T
if 'CRVAL3' in hdr.ascard.keys(): del hdr['CRVAL3']
if 'CRPIX3' in hdr.ascard.keys(): del hdr['CRPIX3']
if 'CUNIT3' in hdr.ascard.keys(): del hdr['CUNIT3']
if 'CDELT3' in hdr.ascard.keys(): del hdr['CDELT3']
if 'CD3_3' in hdr.ascard.keys(): del hdr['CD3_3']
if 'CD3_2' in hdr.ascard.keys(): del hdr['CD3_2']
if 'CD3_1' in hdr.ascard.keys(): del hdr['CD3_1']
......@@ -117,13 +122,14 @@ def create_config(path, cubefile, varfile, catfile, lines, suffixeout, commonw=T
obj = int(line['ID'])
suffixe = suffixeout + '_' + str(obj) + '_' + lines
fileout = suffixe + '.config'
fileout = suffixe + '/' + suffixe + '.config'
logger.info('Writing config ' + fileout)
print(path + fileout)
fo = open(path + fileout, 'w')
fo.write('FITSFILE= ' + cubefile + ' / path of the fits file containing the cube\n')
fo.write('OUTPUT = ' + suffixe + '/' + suffixe + ' / where to write the output (default: ./out.fits)\n')
fo.write('SKYFILE = ' + varfile + ' / path of the file containing the sky spectrum or cube\n')
fo.write('FITSFILE= ' + '../' + cubefile + ' / path of the fits file containing the cube\n')
fo.write('OUTPUT = ' + suffixe + ' / where to write the output (default: ./out.fits)\n')
fo.write('SKYFILE = ' + '../' + varfile + ' / path of the file containing the sky spectrum or cube\n')
lbda = []
if 'ha' in lines:
......@@ -161,10 +167,10 @@ def create_config(path, cubefile, varfile, catfile, lines, suffixeout, commonw=T
lbda.append(4101.73)
fo.write('COMMW = ' + str(commonw) + ' / do we use a common width (default: FALSE)\n')
fo.write('REDSHIFT= ' + str(line['z']) + ' / initial redshift of the galaxy\n')
fo.write('REDSHIFT= ' + str(line['Z']) + ' / initial redshift of the galaxy\n')
if dz is not None:
fo.write('REDMIN = ' + str(line['z'] - dz) + ' / minimum redshift\n')
fo.write('REDMAX = ' + str(line['z'] + dz) + ' / maximum redshift\n')
fo.write('REDMIN = ' + str(line['Z'] - dz) + ' / minimum redshift\n')
fo.write('REDMAX = ' + str(line['Z'] + dz) + ' / maximum redshift\n')
else:
logger.error('dz must be not None and non-zero')
......@@ -179,9 +185,9 @@ def create_config(path, cubefile, varfile, catfile, lines, suffixeout, commonw=T
fo.write('WSMOOTH = ' + str(wsmooth) + ' / spectral smoothing: 1 for Hanning, 2 for gaussian (FWHM=3 pixels), set 0 for no smoothing (default:0)\n')
fo.write('SSMOOTH = ' + str(ssmooth) + ' / gaussian spatial smooth FWHM in pixel, set 0 for no smoothing (default:0)\n')
ra = coord.Angle(line['ra'], unit=u.degree)
dec = coord.Angle(line['dec'], unit=u.degree)
logger.debug(' converting ra,dec %f,%f ' % (line['ra'],line['dec']) )
ra = coord.Angle(line['RA'], unit=u.degree)
dec = coord.Angle(line['DEC'], unit=u.degree)
logger.debug(' converting ra,dec %f,%f ' % (line['RA'],line['DEC']) )
xc, yc = np.round(w.wcs_world2pix(ra, dec, 0))
#xc, yc, _ = np.round(w.wcs_world2pix(ra, dec, 0, 0))
......@@ -198,8 +204,8 @@ def create_config(path, cubefile, varfile, catfile, lines, suffixeout, commonw=T
fo.write('YMIN = None / Y minimum value for spatial cut of the input cube (pixel)\n')
fo.write('YMAX = None / Y maximum value for spatial cut of the input cube (pixel)\n')
if dz is not None:
minlbda = np.min(np.array(lbda)) * (1 + line['z'] - 1.5 * dz)
maxlbda = np.max(np.array(lbda)) * (1 + line['z'] + 1.5 * dz)
minlbda = np.min(np.array(lbda)) * (1 + line['Z'] - 1.5 * dz)
maxlbda = np.max(np.array(lbda)) * (1 + line['Z'] + 1.5 * dz)
zmin = np.round((minlbda - crval) / cdelt + crpix - 1)
zmax = np.round((maxlbda - crval) / cdelt + crpix - 1)
fo.write('ZMIN = ' + str(zmin) + ' / Z minimum value for spectral cut of the input cube (spectral channel, integer)\n')
......@@ -211,21 +217,47 @@ def create_config(path, cubefile, varfile, catfile, lines, suffixeout, commonw=T
fo.close()
def main():
def main_byhand():
gr = '116'
#gr = '83'
#gr = '28'
#gr = '32'
path = '/home/bepinat/Instruments/MUSE/groups/kinematics/gr' + gr + '/'
cubefile = 'DATACUBE_FINAL_ZAP_COSMOS-GR' + gr + '_newWCS.fits'
varfile = cubefile
catfile = 'catalog_gr' + gr + '.txt'
suffixeout = 'gr' + gr
lines = 'o3hb'
lines = 'o2'
suffixeout = 'gr' + gr
path = '/home/bepinat/Instruments/MUSE/groups/kinematics/gr' + gr + '/' + lines + '/'
cubefile = '../DATACUBE_FINAL_ZAP_COSMOS-GR' + gr + '_newWCS.fits'
varfile = cubefile
catfile = '../catalog_gr' + gr + '.txt'
create_config(path, cubefile, varfile, catfile, lines, suffixeout)
def main(argv):
parser = argparse.ArgumentParser()
parser.add_argument('--lines', '-l', action="store", dest="lines", default="", help="name of the lines")
parser.add_argument('--suffixe', '-s', action="store", dest="suffixeout", default="", help="suffixe")
parser.add_argument('--path', '-p', action="store", dest="pathb", default="", help="initial path")
parser.add_argument('--cat', '-c', action="store", dest="catfile", default="", help="path of the catalogue with respect to initial path")
parser.add_argument('--data', '-d', action="store", dest="cubefile", default="", help="path of the cube with respect to initial path")
parser.add_argument('--var', '-v', action="store", dest="varfile", default=None, help="path of the variance with respect to initial path")
parser.add_argument('--dxy', action="store", dest="dxy", default='15', help="size in pixels around the center")
try:
options = parser.parse_args()
except:
sys.exit(2)
opts = vars(options) # this is a dictionnary
#print(opts)
if opts['varfile'] is None:
opts['varfile'] = opts['cubefile']
create_config(opts['pathb'], opts['cubefile'], opts['varfile'], opts['catfile'], opts['lines'], opts['suffixeout'], dxy=int(opts['dxy']))
if __name__ == "__main__":
main()
main(sys.argv[1:])
#main_byhand()
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