update on code: remove pyfits for astropy.io.fits

camel.py

-#!/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")

create_config.py

-#!/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)