Commit ec6af203 authored by bepinat's avatar bepinat

bug correction + variance creation from data in case no variance cube exist

parent bd8f5722
...@@ -181,10 +181,25 @@ def readcubes(conf, multi_ext=False): ...@@ -181,10 +181,25 @@ def readcubes(conf, multi_ext=False):
var = hdul[0].data var = hdul[0].data
varhdr = hdul[0].header varhdr = hdul[0].header
else: else:
print('toto')
hdul = pf.open(conf['FITSFILE']) hdul = pf.open(conf['FITSFILE'])
cube = hdul[0].data cube = hdul[0].data
hdr = hdul[0].header hdr = hdul[0].header
if (not conf['SKYFILE']): if conf['SKYFILE'] == conf['FITSFILE']:
var = np.zeros(cube.shape)
if conf['XMIN'] == None:
conf['XMIN'] = 0
if conf['YMIN'] == None:
conf['YMIN'] = 0
if conf['XMAX'] == None:
conf['XMAX'] = cube.shape[2] - 1
if conf['YMAX'] == None:
conf['YMAX'] = cube.shape[1] - 1
for i in np.arange(cube.shape[0]):
var[i,:,:] = np.var(cube[i,conf['YMIN']:conf['YMAX'],conf['XMIN']:conf['XMAX']])
varhdr = hdul[0].header
var = writedata(var, varhdr, conf['OUTPUT'] + '_variance.fits')
elif (not conf['SKYFILE']):
var = None var = None
varhdr = '' varhdr = ''
else: else:
...@@ -229,7 +244,7 @@ def writedata(data, hdr, filename): ...@@ -229,7 +244,7 @@ def writedata(data, hdr, filename):
""" """
hdu = pf.PrimaryHDU(data, hdr) hdu = pf.PrimaryHDU(data, hdr)
hdulist = pf.HDUList(hdu) hdulist = pf.HDUList(hdu)
hdulist.writeto(filename, clobber=True) hdulist.writeto(filename, clobber=True, output_verify='fix')
return hdulist[0].header return hdulist[0].header
...@@ -476,7 +491,7 @@ def waveindgen(hdr): ...@@ -476,7 +491,7 @@ def waveindgen(hdr):
''' '''
cunit = hdr['CUNIT3'] cunit = hdr['CUNIT3']
if cunit.strip().lower() in {'microns', 'micrometers', 'mum'}: lconvfac=1e4 if cunit.strip().lower() in {'microns', 'micron', 'micrometers', 'mum'}: lconvfac=1e4
if cunit.strip().lower() in {'angstroms', 'angstrom'}: lconvfac=1 if cunit.strip().lower() in {'angstroms', 'angstrom'}: lconvfac=1
if cunit.strip().lower() in {'nanometers', 'nm'}: lconvfac=1e3 if cunit.strip().lower() in {'nanometers', 'nm'}: lconvfac=1e3
if 'CDELT3' in hdr.ascard.keys(): if 'CDELT3' in hdr.ascard.keys():
...@@ -750,6 +765,7 @@ def buildmaps(conf, cube, hdr, var=None, plot=False, debug=False, free_ratio=Fal ...@@ -750,6 +765,7 @@ def buildmaps(conf, cube, hdr, var=None, plot=False, debug=False, free_ratio=Fal
fa['spectrum'] = cube[indok,y,x] # Spectrum fa['spectrum'] = cube[indok,y,x] # Spectrum
fa['err'] = np.sqrt(var[indok,y,x]) # Error fa['err'] = np.sqrt(var[indok,y,x]) # Error
#fa['weight'] = np.sqrt(weight[indok,y,x]) # Weight #fa['weight'] = np.sqrt(weight[indok,y,x]) # Weight
if (np.min(fa['spectrum']) == np.max(fa['spectrum'])) & (np.min(fa['spectrum']) == 0): continue # No need for fitting if (np.min(fa['spectrum']) == np.max(fa['spectrum'])) & (np.min(fa['spectrum']) == 0): continue # No need for fitting
minfnorm = np.infty minfnorm = np.infty
...@@ -906,6 +922,7 @@ def buildmaps(conf, cube, hdr, var=None, plot=False, debug=False, free_ratio=Fal ...@@ -906,6 +922,7 @@ def buildmaps(conf, cube, hdr, var=None, plot=False, debug=False, free_ratio=Fal
if 'CD3_1' in hdr.ascard.keys(): del hdr['CD3_1'] if 'CD3_1' in hdr.ascard.keys(): del hdr['CD3_1']
if 'CD2_3' in hdr.ascard.keys(): del hdr['CD2_3'] if 'CD2_3' in hdr.ascard.keys(): del hdr['CD2_3']
if 'CD1_3' in hdr.ascard.keys(): del hdr['CD1_3'] if 'CD1_3' in hdr.ascard.keys(): del hdr['CD1_3']
if 'CDELT3' in hdr.ascard.keys(): del hdr['CDELT3']
suffixe = '' suffixe = ''
if conf['COMMW']: if conf['COMMW']:
...@@ -942,8 +959,8 @@ def buildmaps(conf, cube, hdr, var=None, plot=False, debug=False, free_ratio=Fal ...@@ -942,8 +959,8 @@ def buildmaps(conf, cube, hdr, var=None, plot=False, debug=False, free_ratio=Fal
hh = writedata(edispmaps[i,:,:], hdr, conf['OUTPUT'] + '_edisp' + suffixe + '_' + linename[i] + '.fits') hh = writedata(edispmaps[i,:,:], hdr, conf['OUTPUT'] + '_edisp' + suffixe + '_' + linename[i] + '.fits')
if conf['COMMW']: if conf['COMMW']:
cond = np.where(linefit) cond = np.where(linefit)
hh = writedata(dispmaps[cond[0]:,:], hdr, conf['OUTPUT'] + '_disp' + suffixe + '.fits') hh = writedata(dispmaps[cond[0][0],:,:], hdr, conf['OUTPUT'] + '_disp' + suffixe + '.fits')
hh = writedata(edispmaps[cond[0],:,:], hdr, conf['OUTPUT'] + '_edisp' + suffixe + '.fits') hh = writedata(edispmaps[cond[0][0],:,:], hdr, conf['OUTPUT'] + '_edisp' + suffixe + '.fits')
return rvmap return rvmap
...@@ -999,6 +1016,8 @@ def camel(filename, plot=False, debug=False, free_ratio=False, multi_ext=False, ...@@ -999,6 +1016,8 @@ def camel(filename, plot=False, debug=False, free_ratio=False, multi_ext=False,
logger.info("Variance size problem") logger.info("Variance size problem")
return return
#XXX when no variance, I should add an option to compute the variance from the cube itself (as done in the IDL code)
#Cutting the cubes #Cutting the cubes
logger.info("Cutting data cube") logger.info("Cutting data cube")
cube, hdr = cutcube(cube, hdr, conf) cube, hdr = cutcube(cube, hdr, conf)
......
...@@ -103,6 +103,7 @@ def create_config(path, cubefile, varfile, catfile, lines, suffixeout, commonw=T ...@@ -103,6 +103,7 @@ def create_config(path, cubefile, varfile, catfile, lines, suffixeout, commonw=T
if 'CD3_1' in hdr.ascard.keys(): del hdr['CD3_1'] if 'CD3_1' in hdr.ascard.keys(): del hdr['CD3_1']
if 'CD2_3' in hdr.ascard.keys(): del hdr['CD2_3'] if 'CD2_3' in hdr.ascard.keys(): del hdr['CD2_3']
if 'CD1_3' in hdr.ascard.keys(): del hdr['CD1_3'] if 'CD1_3' in hdr.ascard.keys(): del hdr['CD1_3']
if '' in hdr.ascard.keys(): del hdr['']
if 'COMMENT' in hdr.ascard.keys(): del hdr['COMMENT'] if 'COMMENT' in hdr.ascard.keys(): del hdr['COMMENT']
if 'HISTORY' in hdr.ascard.keys(): del hdr['HISTORY'] if 'HISTORY' in hdr.ascard.keys(): del hdr['HISTORY']
...@@ -116,6 +117,7 @@ def create_config(path, cubefile, varfile, catfile, lines, suffixeout, commonw=T ...@@ -116,6 +117,7 @@ def create_config(path, cubefile, varfile, catfile, lines, suffixeout, commonw=T
logger.debug(' Naxis %d ' % (w.naxis)) logger.debug(' Naxis %d ' % (w.naxis))
lc = np.size(cat) lc = np.size(cat)
if lc == 1: cat = [cat]
logger.debug(" catatalog length %s " % (str(lc)) ) logger.debug(" catatalog length %s " % (str(lc)) )
for line in cat: for line in cat:
......
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