diff --git a/camel.py b/camel.py index 3caf5aa7124c61160652b08c23d4f17cf4185219..d4affc3731e22b998fd1524d4b9d6efd50ccb0e9 100755 --- a/camel.py +++ b/camel.py @@ -1,4 +1,4 @@ -#!/usr/bin/env python3 +#!/usr/bin/env python # -*- coding: utf-8 -*- """ @@ -23,7 +23,8 @@ import time import os, sys import copy import numpy as np -import pyfits as pf +#import pyfits as pf +import astropy.io.fits as fits from cap_mpfit import mpfit # python 2.7 & 3! #from mpfit import mpfit # python 2.7... import ipdb @@ -50,7 +51,7 @@ def compute_rv(z, zcosmo=None): """ - if zcosmo == None: + if zcosmo is None: v = ((z + 1) ** 2 - 1) / ((z + 1) ** 2 + 1) * ct.c * 1e-3 else: v = (z - zcosmo) / (1 + z) * ct.c * 1e-3 @@ -70,7 +71,7 @@ def compute_fwhm(dz, z=None, zcosmo=0): """ - if z == None: z = zcosmo + if z is None: z = zcosmo dv = dz * (1 + zcosmo) / (1 + z) ** 2 * ct.c *1e-3 return dv @@ -87,7 +88,7 @@ def readconfig(filename, config): """ conf = {} - if filename != None: + if filename is not None: try: data = open(filename) except: @@ -122,7 +123,7 @@ def readconfig(filename, config): value = int(value) except: value = value - if (key == 'EXTRAL') & (value != None): + if (key == 'EXTRAL') & (value is not None): if count_extral == 0: conf[key] = np.array([value]) else: @@ -132,9 +133,9 @@ def readconfig(filename, config): conf[key] = value # add input options to configuration dictionnary + check needed keywords - if config != None: + if config is not None: for key in config: - if config[key] != None: + if config[key] is not None: conf[key] = config[key] elif not(key in conf.keys()): conf[key] = None @@ -148,8 +149,8 @@ def readconfig(filename, config): if 'OII' in conf.keys(): conf['OII3729'] = conf['OII3726'] = conf.pop('OII') - if conf['SPSF'] == None: conf['SPSF'] = 0. - if conf['COMMW'] == None: conf['COMMW'] = False + if conf['SPSF'] is None: conf['SPSF'] = 0. + if conf['COMMW'] is None: conf['COMMW'] = False # XXX or give an error message but in that case, must test that the keywords have correct values return conf @@ -166,7 +167,7 @@ def readcubes(conf, multi_ext=False): """ if multi_ext: - hdul = pf.open(conf['FITSFILE']) + hdul = fits.open(conf['FITSFILE']) cube = hdul[1].data hdr = hdul[1].header if conf['SKYFILE'] == conf['FITSFILE']: @@ -176,22 +177,22 @@ def readcubes(conf, multi_ext=False): var = None varhdr = '' else: - hdul = pf.open(conf['SKYFILE']) + hdul = fits.open(conf['SKYFILE']) var = hdul[0].data varhdr = hdul[0].header else: - hdul = pf.open(conf['FITSFILE']) + hdul = fits.open(conf['FITSFILE']) cube = hdul[0].data hdr = hdul[0].header if conf['SKYFILE'] == conf['FITSFILE']: var = np.zeros(cube.shape, dtype='>f4') - if conf['XMIN'] == None: + if conf['XMIN'] is None: conf['XMIN'] = 0 - if conf['YMIN'] == None: + if conf['YMIN'] is None: conf['YMIN'] = 0 - if conf['XMAX'] == None: + if conf['XMAX'] is None: conf['XMAX'] = cube.shape[2] - 1 - if conf['YMAX'] == None: + if conf['YMAX'] is 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']]) @@ -201,7 +202,7 @@ def readcubes(conf, multi_ext=False): var = None varhdr = '' else: - hdul = pf.open(conf['SKYFILE']) + hdul = fits.open(conf['SKYFILE']) var = hdul[0].data varhdr = hdul[0].header return cube, hdr, var, varhdr @@ -217,7 +218,7 @@ def testcubes(cube, var): input variance data """ - if var == None: + if var is None: return if var.ndim == 3: if var.shape != cube.shape: @@ -240,8 +241,8 @@ def writedata(data, hdr, filename): name of the output """ - hdu = pf.PrimaryHDU(data, hdr) - hdulist = pf.HDUList(hdu) + hdu = fits.PrimaryHDU(data, hdr) + hdulist = fits.HDUList(hdu) hdulist.writeto(filename, clobber=True, output_verify='fix') return hdulist[0].header @@ -288,17 +289,17 @@ def cutcube(cube, hdr, conf): """ - if conf['XMIN'] == None: + if conf['XMIN'] is None: conf['XMIN'] = 0 - if conf['YMIN'] == None: + if conf['YMIN'] is None: conf['YMIN'] = 0 - if conf['ZMIN'] == None: + if conf['ZMIN'] is None: conf['ZMIN'] = 0 - if conf['XMAX'] == None: + if conf['XMAX'] is None: conf['XMAX'] = cube.shape[2] - 1 - if conf['YMAX'] == None: + if conf['YMAX'] is None: conf['YMAX'] = cube.shape[1] - 1 - if conf['ZMAX'] == None: + if conf['ZMAX'] is None: conf['ZMAX'] = cube.shape[0] - 1 conf['ZMIN'], conf['YMIN'], conf['XMIN'] = np.max([[conf['ZMIN'], conf['YMIN'], conf['XMIN']], [0, 0, 0]],axis=0) @@ -429,7 +430,7 @@ class line: self.name = name self.wave = wave self.fit = fit - if ref == None: ref = name + if ref is None: ref = name self.ref = ref if ref == name: self.low = 1 @@ -465,8 +466,8 @@ class lines: self.append(line('SII6716', 6716.44, ref='SII6731', low=0.45, up=1.45, th=1.)) self.append(line('OIII5007', 5006.843, ref='OIII4959', low=2.7, up=3.3, th=3.)) self.append(line('OIII4959', 4958.911)) - self.append(line('OIII4363', 4363., ref='OIII4959')) # XXX low=, up=, th= - self.append(line('OII3729', 3728.80, ref='OII3726', low=0.3, up=1.5, th=1.)) + self.append(line('OIII4363', 4363.21, ref='OIII4959')) # XXX low=, up=, th= + self.append(line('OII3729', 3728.80, ref='OII3726', low=0.35, up=1.5, th=1.)) self.append(line('OII3726', 3726.04)) self.append(line('HEI4471', 4471.)) self.append(line('HEI5876', 5876., ref='HEI4471', low=2.5, th=2.5)) @@ -619,12 +620,12 @@ def build_parinfo_lines(els, linename, conf, free_ratio, comind, multi_comp=1): parinfo[2 * els.index + i]['tied'] = 'p[%i]'%(2 * els.index + j + multi_comp * len(parinfo)) # line intensity is tied to reference #Ratio limits parinfo[3 * els.index + i]['limited'][0] = 1 - if els.lines[name].low != None: + if els.lines[name].low is not None: parinfo[3 * els.index + i]['limits'][0] = els.lines[name].low - if els.lines[name].up != None: + if els.lines[name].up is not None: parinfo[3 * els.index + i]['limited'][1] = 1 parinfo[3 * els.index + i]['limits'][1] = els.lines[name].up - if els.lines[name].th != None: + if els.lines[name].th is not None: parinfo[3 * els.index + i]['value'] = els.lines[name].th #params[3 * els.index + i] = els.lines[name].th else: @@ -648,12 +649,12 @@ def build_parinfo_lines(els, linename, conf, free_ratio, comind, multi_comp=1): #if (not free_ratio) & (i!=j): #parinfo[2 * els.index + i]['tied'] = 'p[%i]'%(2 * els.index + j + multi_comp * len(parinfo)) #parinfo[3 * els.index + i]['limited'][0] = 1 - #if els.lines[name].low != None: + #if els.lines[name].low is not None: #parinfo[3 * els.index + i]['limits'][0] = els.lines[name].low - #if els.lines[name].up != None: + #if els.lines[name].up is not None: #parinfo[3 * els.index + i]['limited'][1] = 1 #parinfo[3 * els.index + i]['limits'][1] = els.lines[name].up - #if els.lines[name].th != None: + #if els.lines[name].th is not None: #parinfo[3 * els.index + i]['value'] = params[3 * els.index + i] = els.lines[name].th #else: #params[3 * els.index + i] = 1 @@ -678,7 +679,7 @@ def buildmaps(conf, cube, hdr, var=None, plot=False, debug=False, free_ratio=Fal if lline in conf.keys(): els.lines[lline].fit = conf[lline] - if conf['EXTRAL'] != None: + if conf['EXTRAL'] is not None: for i in np.arange(np.shape(conf['EXTRAL'])[0]): els.append(line('EXTRAL%i'%(i+1), conf['EXTRAL'][i], fit=True)) @@ -698,10 +699,10 @@ def buildmaps(conf, cube, hdr, var=None, plot=False, debug=False, free_ratio=Fal linewave = np.array([els.lines[i].wave for i in els.lines])[argsorted] linename = np.array([els.lines[i].name for i in els.lines])[argsorted] indok = np.zeros(0, dtype = int) - firstind = np.zeros(linefit.size, dtype='>f4') - lastind = np.zeros(linefit.size, dtype='>f4') - firstind0 = np.zeros(linefit.size, dtype='>f4') - lastind0 = np.zeros(linefit.size, dtype='>f4') + firstind = np.zeros(linefit.size, dtype=int) + lastind = np.zeros(linefit.size, dtype=int) + firstind0 = np.zeros(linefit.size, dtype=int) + lastind0 = np.zeros(linefit.size, dtype=int) for i in range(linefit.size): if (not linefit[i]): continue linf = np.max([(zmin - zrange / 2. + 1) * linewave[i] , np.min(wave0)]) @@ -747,8 +748,8 @@ def buildmaps(conf, cube, hdr, var=None, plot=False, debug=False, free_ratio=Fal #lastind[i] = condo[0][np.size(condo) - 1] #Weights - #if var == None: var = np.ones([indok.size, cube.shape[1], cube.shape[2]]) - if var == None: var = np.ones(cube.shape, dtype='>f4') + #if var is None: var = np.ones([indok.size, cube.shape[1], cube.shape[2]]) + if var is None: var = np.ones(cube.shape, dtype='>f4') 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 # weight inverse variance @@ -791,6 +792,12 @@ def buildmaps(conf, cube, hdr, var=None, plot=False, debug=False, free_ratio=Fal #Continuum parametres parbase = {'value':0., 'fixed':0, 'limited':[0,0], 'limits':[0.,0.], 'tied':'', 'mpmaxstep':0, 'mpminstep':0, 'parname':''} + if conf['DGCTNUM'] == -1: # case without continuum + conf['DGCTNUM'] = 0 + parbase['fixed'] = 1 + contfixed = True + else: + contfixed = False contind = [i for i in range(len(parinfo), len(parinfo) + conf['DGCTNUM'] + 1)] for i in range(conf['DGCTNUM'] + 1): parinfo.append(copy.deepcopy(parbase)) @@ -849,7 +856,8 @@ def buildmaps(conf, cube, hdr, var=None, plot=False, debug=False, free_ratio=Fal #if parinfo[2 * els.index + i]['fixed'] == 0: p[2 * els.index + i] = 2 * std #Inititalisation of continuum (degree 0) - p[contind[0]] = np.median(fa['spectrum']) + if not contfixed: + p[contind[0]] = np.median(fa['spectrum']) #p[els.index * 4] = np.median(fa['spectrum']) for z in np.arange(zzmin, zzmax, zstep): #p[comind] = z @@ -1159,7 +1167,7 @@ def camel(filename, plot=False, debug=False, free_ratio=False, multi_comp=1, mul logger.info("Cutting data cube") cube, hdr = cutcube(cube, hdr, conf) hdr = writedata(cube, hdr, conf['OUTPUT'] + '_cube_cut.fits') - if var != None: + if var is not None: logger.info("Cutting variance") var, varhdr = cutcube(var, varhdr, conf) if var.ndim == 3: @@ -1169,7 +1177,7 @@ def camel(filename, plot=False, debug=False, free_ratio=False, multi_comp=1, mul if cut & (not(clip)) & (not(smooth)): return #Clipping the cube - if (conf['SCLIP'] != 0) & (conf['SCLIP'] != None): + if (conf['SCLIP'] != 0) & (conf['SCLIP'] is not 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') @@ -1177,24 +1185,24 @@ def camel(filename, plot=False, debug=False, free_ratio=False, multi_comp=1, mul if clip & (not(smooth)): return #Smoothing the cubes - if (conf['WSMOOTH'] != 0) & (conf['WSMOOTH'] != None): + if (conf['WSMOOTH'] != 0) & (conf['WSMOOTH'] is not None): 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') - if var != None: + 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['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): + if (conf['SSMOOTH'] != 0) & (conf['SSMOOTH'] is not None): 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') - if var != None: + 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 @@ -1255,7 +1263,7 @@ def main(argv): """ """ parser = argparse.ArgumentParser(description=usage()) - parser.add_argument('--file', '-f', action="store", dest="filename", default=None, help="name of the configuration file") + parser.add_argument('--file', '-f', action="store", dest="filename", default='', help="name of the configuration file") parser.add_argument('--plot', '-p', action="store_true", dest="plot", default=False, help="keyword to plot fits during the process (not implemented yet)") parser.add_argument('--debug', '-d', action="store_true", dest="debug", default=False, help="keyword to show debug information (not implemented yet)") parser.add_argument('--free_line_ratio', '-r', action="store_true", dest="free_ratio", default=False, help="keyword to allow unconstrained line ratio") diff --git a/create_config.py b/create_config.py index 782774485c1491deae545c98112d9a47539eb735..59ba271ae3e163084c683d77c81f888aa6ddfed8 100755 --- a/create_config.py +++ b/create_config.py @@ -1,11 +1,12 @@ -#!/usr/bin/env python3 +#!/usr/bin/env python # -*- coding: utf-8 -*- import os, sys import argparse import ipdb import numpy as np -import pyfits as pf +#import pyfits as pf +import astropy.io.fits as fits import astropy.wcs as wcs import astropy.coordinates as coord import astropy.units as u @@ -65,7 +66,7 @@ def create_config(path, cubefile, varfile, catfile, lines, suffixeout, commonw=T cat = np.genfromtxt(path + catfile, names=True, case_sensitive=False) try: - hdul = pf.open(path + cubefile) + hdul = fits.open(path + cubefile) logger.info('Using cube %s ' % (path + cubefile) ) except: logger.info('Not able to read cube %s !' % (path + cubefile) ) @@ -92,24 +93,23 @@ def create_config(path, cubefile, varfile, catfile, lines, suffixeout, commonw=T except: cdelt = hdr['CDELT3'] logger.debug(" wcs check: %s %d %.4f %.6f " % (cunit,crpix,crval,cdelt)) - - if 'CTYPE3' in hdr.ascard.keys(): del hdr['CTYPE3'] - 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'] - if 'CD2_3' in hdr.ascard.keys(): del hdr['CD2_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 'HISTORY' in hdr.ascard.keys(): del hdr['HISTORY'] + if 'CTYPE3' in hdr.keys(): del hdr['CTYPE3'] + if 'CRVAL3' in hdr.keys(): del hdr['CRVAL3'] + if 'CRPIX3' in hdr.keys(): del hdr['CRPIX3'] + if 'CUNIT3' in hdr.keys(): del hdr['CUNIT3'] + if 'CDELT3' in hdr.keys(): del hdr['CDELT3'] + if 'CD3_3' in hdr.keys(): del hdr['CD3_3'] + if 'CD3_2' in hdr.keys(): del hdr['CD3_2'] + if 'CD3_1' in hdr.keys(): del hdr['CD3_1'] + if 'CD2_3' in hdr.keys(): del hdr['CD2_3'] + if 'CD1_3' in hdr.keys(): del hdr['CD1_3'] + if '' in hdr.keys(): del hdr[''] + if 'COMMENT' in hdr.keys(): del hdr['COMMENT'] + if 'HISTORY' in hdr.keys(): del hdr['HISTORY'] wlim = np.zeros(np.shape(cube)[1:]) - hdu = pf.PrimaryHDU(data=wlim, header=hdr) - hdul1 = pf.HDUList(hdu) + hdu = fits.PrimaryHDU(data=wlim, header=hdr) + hdul1 = fits.HDUList(hdu) hdr = hdul1[0].header w = wcs.WCS(hdr, hdul1) #w = wcs.WCS(hdr[:28], hdul1)