Commit 9ce864e2 authored by bepinat's avatar bepinat

message

parent daf55af2
......@@ -87,7 +87,11 @@ def readconfig(filename, config):
conf = {}
if filename != None:
data = open(filename)
try:
data = open(filename)
except:
logger.info('Camel is not able to open ' + filename)
data = open(filename)
count_extral = 0
for raw in data:
keyvalcom = raw.split('= ')
......@@ -138,7 +142,7 @@ def readconfig(filename, config):
# XXX or give an error message but in that case, must test that the keywords have correct values
## XXX more keywords? add default values?
needed_keys = ['FITSFILE', 'OUTPUT', 'SKYFILE', 'HALPHA', 'NII6548', 'NII6583', 'SII6716', 'SII6731', 'OIII4959', 'OIII5007', 'HBETA', 'OII3729', 'OII3726', 'COMMW', 'REDSHIFT', 'REDMIN', 'REDMAX', 'INITW', 'WMIN', 'WMAX', 'DFIT', 'DGCTNUM', 'MFIT', 'SCLIP', 'XYCLIP', 'NCLIP', 'SPSF', 'WSMOOTH', 'SSMOOTH', 'THRES', 'MEDIAN', 'FITSPSF', 'XMIN', 'YMIN', 'ZMIN', 'XMAX', 'YMAX', 'ZMAX']
needed_keys = ['FITSFILE', 'OUTPUT', 'SKYFILE', 'HALPHA', 'NII6548', 'NII6583', 'SII6716', 'SII6731', 'OIII4959', 'OIII5007', 'HBETA', 'OII3729', 'OII3726', 'EXTRAL', 'COMMW', 'REDSHIFT', 'REDMIN', 'REDMAX', 'INITW', 'WMIN', 'WMAX', 'DFIT', 'DGCTNUM', 'MFIT', 'SCLIP', 'XYCLIP', 'NCLIP', 'SPSF', 'WSMOOTH', 'SSMOOTH', 'THRES', 'MEDIAN', 'FITSPSF', 'XMIN', 'YMIN', 'ZMIN', 'XMAX', 'YMAX', 'ZMAX']
keys = conf.keys()
for key in needed_keys:
#if keys.isdisjoint([key]): # c'est du python3
......@@ -289,7 +293,7 @@ def remove_nanvar(var):
qvar = var[:, indnum[0], indnum[1]].reshape(var.shape[0], np.size(indnum[0]))
medianvar = np.median(qvar, axis=1)
if np.size(indnan[0]) == 0:
print('Variance cube contains NaN and has been corrected')
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
return var
......@@ -706,7 +710,7 @@ def buildmaps(conf, cube, hdr, var=None, plot=False, debug=False, free_ratio=Fal
fac_prev = np.int((counter / number_of_pixels) * 100)
sys.stdout.write("Progress => {:3d}%\r".format(fac_prev))
sys.stdout.flush()
#if (counter / 10.) == np.int(counter / 10.): print('pixel %i / %i'%(counter, number_of_pixels))
if (counter / 10.) == np.int(counter / 10.): logger.debug('pixel %i / %i'%(counter, number_of_pixels))
counter += 1
fa['spectrum'] = cube[indok,y,x] # Spectrum
......@@ -742,7 +746,7 @@ def buildmaps(conf, cube, hdr, var=None, plot=False, debug=False, free_ratio=Fal
fit = copy.copy(fit0)
#print(z,tf - ti)
if minfnorm == np.infty: continue
#print(fit.params[2 * els.index + 1 + 12], parinfo[2 * els.index + 1 + 12]['limits'], fit.params[2 * els.index + 1 + 13])
#logger.debug(fit.params[2 * els.index + 1 + 12], parinfo[2 * els.index + 1 + 12]['limits'], fit.params[2 * els.index + 1 + 13])
paramscube[:, y, x] = fit.params[:]
perrorcube[:, y, x] = fit.perror[:]
fnormmap[y, x] = fit.fnorm
......@@ -751,7 +755,7 @@ def buildmaps(conf, cube, hdr, var=None, plot=False, debug=False, free_ratio=Fal
modcube[indok, y, x] = elspectrum(fit.params, wave=wave, lines=linewave, psf=conf['SPSF'], degc=conf['DGCTNUM'])
tff = time.time()
#print(tff - tii)
logger.debug('Execution time: %3.1f' %(tff - tii))
#print(tt)
#Creating outputs
......@@ -917,51 +921,51 @@ def camel(filename, plot=False, debug=False, free_ratio=False, multi_ext=False,
t0 = time.time()
#Reading of configuration file
print("Reading configuration file " + filename)
logger.info("Reading configuration file " + filename)
try:
conf = readconfig(filename, config)
except:
print("Configuration file problem")
logger.info("Configuration file problem")
return
#Reading data
try:
print("Reading data cube " + conf['FITSFILE'] + " and variance " + conf['SKYFILE'])
logger.info("Reading data cube " + conf['FITSFILE'] + " and variance " + conf['SKYFILE'])
except:
print("Reading data cube " + conf['FITSFILE'] + ". No variance data.")
logger.info("Reading data cube " + conf['FITSFILE'] + ". No variance data.")
try:
cube, hdr, var, varhdr = readcubes(conf, multi_ext=multi_ext)
except:
print("Input cube problem")
logger.info("Input cube problem")
return
try:
testcubes(cube, var)
except:
print("Variance shape: %i; Data cube shape: %i"%(np.shape(var), cube.shape))
print("Variance size problem")
logger.info("Variance shape: %i; Data cube shape: %i"%(np.shape(var), cube.shape))
logger.info("Variance size problem")
return
#Cutting the cubes
print("Cutting data cube")
logger.info("Cutting data cube")
cube, hdr = cutcube(cube, hdr, conf)
hdr = writedata(cube, hdr, conf['OUTPUT']+'_cube_cut.fits')
hdr = writedata(cube, hdr, conf['OUTPUT'] + '_cube_cut.fits')
if var != None:
print("Cutting variance")
logger.info("Cutting variance")
var, varhdr = cutcube(var, varhdr, conf)
if var.ndim == 3:
var = remove_nanvar(var)
varhdr = writedata(var, varhdr, conf['OUTPUT']+'_var_cut.fits')
varhdr = writedata(var, varhdr, conf['OUTPUT'] + '_var_cut.fits')
#Clipping the cube
if (conf['SCLIP'] != 0) & (conf['SCLIP'] != None):
print('Performing %3.1f-sigma clipping using box of size %i. Number of pass: %i'%(conf['SCLIP'], conf['XYCLIP'], conf['NCLIP']))
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')
#Smoothing the cubes
if (conf['WSMOOTH'] != 0) & (conf['WSMOOTH'] != None):
print('Performing %3.1f pixels 1D spectral smoothing'%conf['WSMOOTH'])
logger.info('Performing %3.1f pixels 1D spectral smoothing'%conf['WSMOOTH'])
cube, hdr = spectralsmooth(cube, hdr, conf['WSMOOTH'])
conf['OUTPUT'] = conf['OUTPUT'] + '_wsmooth'
hdr = writedata(cube, hdr, conf['OUTPUT'] + '_cube.fits')
......@@ -971,7 +975,7 @@ def camel(filename, plot=False, debug=False, free_ratio=False, multi_ext=False,
varhdr = writedata(var, varhdr, conf['OUTPUT'] + '_var.fits')
if (conf['SSMOOTH'] != 0) & (conf['SSMOOTH'] != None):
print('Performing %3.1f pixels 2D gaussian spatial smoothing'%conf['SSMOOTH'] )
logger.info('Performing %3.1f pixels 2D gaussian spatial smoothing'%conf['SSMOOTH'] )
cube, hdr = spatialsmooth(cube, hdr, conf['SSMOOTH'])
conf['OUTPUT'] = conf['OUTPUT'] + '_ssmooth'
hdr = writedata(cube, hdr, conf['OUTPUT'] + '_cube.fits')
......
......@@ -60,11 +60,25 @@ def create_config(path, cubefile, varfile, catfile, lines, suffixeout, commonw=T
cat = np.genfromtxt(path + catfile, names=True)
hdul = pf.open(path + cubefile)
logger.info(' using cube %s ' % (path + cubefile) )
try:
hdul = pf.open(path + cubefile)
logger.info('Using cube %s ' % (path + cubefile) )
except:
logger.info('Not able to read cube %s !' % (path + cubefile) )
return
try:
cube = hdul[1].data
hdr = hdul[1].header
logger.info('Using extension 1 of cube %s ' % (path + cubefile))
except:
try:
cube = hdul[0].data
hdr = hdul[0].header
logger.info('Using extension 0 of cube %s ' % (path + cubefile))
except:
logger.info('Problem in data stored in extensions 0 and 1 of cube %s ' % (path + cubefile))
cube = hdul[1].data
hdr = hdul[1].header
cunit = hdr['CUNIT3']
crpix = hdr['CRPIX3']
crval = hdr['CRVAL3']
......@@ -84,12 +98,15 @@ def create_config(path, cubefile, varfile, catfile, lines, suffixeout, commonw=T
if 'CD3_1' in hdr.ascard.keys(): del hdr['CD3_1']
if 'CD2_3' in hdr.ascard.keys(): del hdr['CD2_3']
if 'CD1_3' in hdr.ascard.keys(): del hdr['CD1_3']
if 'COMMENT' in hdr.ascard.keys(): del hdr['COMMENT']
if 'HISTORY' in hdr.ascard.keys(): del hdr['HISTORY']
wlim = np.zeros(np.shape(cube)[1:])
hdu = pf.PrimaryHDU(data=wlim, header=hdr)
hdul1 = pf.HDUList(hdu)
hdr = hdul1[0].header
w = wcs.WCS(hdr[:28], hdul1)
w = wcs.WCS(hdr, hdul1)
#w = wcs.WCS(hdr[:28], hdul1)
logger.debug(' Naxis %d ' % (w.naxis))
......@@ -146,10 +163,10 @@ def create_config(path, cubefile, varfile, catfile, lines, suffixeout, commonw=T
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')
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')
logger.error('dz must be not None and non-zero')
fo.write('INITW = ' + str(initw) + ' / initial linewidth in km/s during the line fitting process \n')
fo.write('WMIN = ' + str(wmin) + ' / minimum linewidth in km/s\n')
......@@ -171,25 +188,25 @@ def create_config(path, cubefile, varfile, catfile, lines, suffixeout, commonw=T
logger.debug(' converting x,y %f %f ' % (xc,yc) )
if dxy is not None:
fo.write('XMIN = ' + str(xc - dxy) + ' / X minimum value for spatial cut of the input cube (pixel)\n')
fo.write('XMAX = ' + str(xc + dxy) + ' / X maximum value for spatial cut of the input cube (pixel)\n')
fo.write('YMIN = ' + str(yc - dxy) + ' / Y minimum value for spatial cut of the input cube (pixel)\n')
fo.write('YMAX = ' + str(yc + dxy) + ' / Y maximum value for spatial cut of the input cube (pixel)\n')
fo.write('XMIN = ' + str(xc - dxy) + ' / X minimum value for spatial cut of the input cube (pixel)\n')
fo.write('XMAX = ' + str(xc + dxy) + ' / X maximum value for spatial cut of the input cube (pixel)\n')
fo.write('YMIN = ' + str(yc - dxy) + ' / Y minimum value for spatial cut of the input cube (pixel)\n')
fo.write('YMAX = ' + str(yc + dxy) + ' / Y maximum value for spatial cut of the input cube (pixel)\n')
else:
fo.write('XMIN = None / X minimum value for spatial cut of the input cube (pixel)\n')
fo.write('XMAX = None / X maximum value for spatial cut of the input cube (pixel)\n')
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')
fo.write('XMIN = None / X minimum value for spatial cut of the input cube (pixel)\n')
fo.write('XMAX = None / X maximum value for spatial cut of the input cube (pixel)\n')
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)
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')
fo.write('ZMAX = ' + str(zmax) + ' / Z maximum value for spectral cut of the input cube (spectral channel, integer)')
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')
fo.write('ZMAX = ' + str(zmax) + ' / Z maximum value for spectral cut of the input cube (spectral channel, integer)')
else:
fo.write('ZMIN = None / Z minimum value for spectral cut of the input cube (spectral channel, integer)\n')
fo.write('ZMAX = None / Z maximum value for spectral cut of the input cube (spectral channel, integer)')
fo.write('ZMIN = None / Z minimum value for spectral cut of the input cube (spectral channel, integer)\n')
fo.write('ZMAX = None / Z maximum value for spectral cut of the input cube (spectral channel, integer)')
fo.close()
......
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