Source code for brutus_ppxf

# -*- coding: utf-8 -*-
'''
 brutus: a set of Python modules to process datacubes from integral field spectrographs.\n
 Copyright (C) 2016,  F.P.A. Vogt
 
 -----------------------------------------------------------------------------------------
 
 This file contains functions to call ppxf to fit the continuum in a given spectra.
 It is based on the examples provided by M. Cappellari within ppxf istelf
 (e.g. ppxf_kinematics_example_sauron.py & ppxf_population_example_sdss.py), but has been
 modified to account for a wavelength dependant spectral resolution of the data. 

 Basically, this is designed with MUSE in mind.

 Created February 2016, F.P.A. Vogt - frederic.vogt@alumni.anu.edu.au
'''
# ----------------------------------------------------------------------------------------

#from __future__ import print_function

from astropy.io import fits as pyfits
from scipy import ndimage
import numpy as np
import glob
import matplotlib.pyplot as plt
from time import clock
import datetime
import warnings

import os
import scipy.signal as signal
import sys

from ppxf import ppxf
import ppxf_util as util

import brutus_tools
from brutus_metadata import *
        
# ----------------------------------------------------------------------------------------
[docs]def sl_resample(fn, z=0, inst = 'MUSE', sl_name='MILES_ppxf_default', do_cap=True): ''' This function resamples a given spectrum (i.e. from a stellar library) to match theresolution of the instrument that took the data. :param: fn: filename of stellar template :param: z: redshift of the target observations :param: inst: which instrument are we using :param: sl_name: which templates do we use ? :do_cap: do the convolution using M.Cappellari in-built function Returns the convolved spectrum ''' if not(os.path.isfile(fn)): sys.exit('"%s": not a valid file.' % fn) # Load the ssp spectra hdu = pyfits.open(fn) ssp_raw = hdu[0].data ssp_header = hdu[0].header ssp_out = np.zeros_like(ssp_raw) hdu.close() # The wavelength step and wavelength array dlam = ssp_header['CDELT1'] lams = np.arange(0,ssp_header['NAXIS1'],1)*ssp_header['CDELT1'] + \ ssp_header['CRVAL1'] # What is the target resolution I want to achieve ? if inst: R_target = brutus_tools.inst_resolution(inst=inst, get_ff=False, show_plot=False) fwhm_target = lams*(z+1)/R_target(lams*(z+1)) else: # Assuming MUSE - probably a bad idea that I will forget about ... Issue a warning # to be on the safe side. warnings.warn('No instrument specified. Assuming MUSE for the ppxf input spectra') R_target = brutus_tools.inst_resolution(inst='MUSE', get_ff=False, show_plot=False) fwhm_target = lams*(z+1)/R_target(lams*(z+1)) # WARNING: what if the resolution of the models is NOT sufficient ? if np.any((fwhm_target - sl_models[sl_name]['fwhm']) < 0): sys.exit('WARNING: stellar templates have too high fwhm given the redshift '+ 'of the target ! Ignoring z for now ... you really should use '+ 'other models !') # Use the ppxf-in-built gaussian_fliter1d home made by M. Cappellari (10%faster !) if do_cap: # Compute the correction in FWHM to be applied. fwhm_diff = np.sqrt(fwhm_target**2 - sl_models[sl_name]['fwhm']**2) # Transform this from Angstroem to resolution elements fwhm_diff = fwhm_diff/dlam # In resolution elements # Compute the corresponding sigma for the gaussian profile sigma = (2*np.sqrt(2*np.log(2)))**-1 * fwhm_diff # Feed this to M. Cappellari elegant routine. N.A.: I checked that function, and # I can reproduce its results (see below) doing things my way. ssp_out = util.gaussian_filter1d(ssp_raw,sigma) else: # do things my way - this is perfectly equivalent, but slower ... # What is the max FWHM in resolution elements ? fwhm_pix_max = np.max(fwhm_target/dlam) # Very well, let's use this to define my box size. Make it 8 times # the max FWHM to correctly include the wings of the gaussian kernel. window_size = 8*np.ceil(fwhm_pix_max) # Make sure the window size is odd (does it matter ?) if window_size % 2 ==0: window_size +=1 window_size = np.int(window_size) # Ok, start the convolution. I do it semi-manually, for cases where it # varies as a function of wavelength (i.e. MUSE) for i in range(len(ssp_raw)): this_lam = lams[i] # The wavelength of the i-th element this_fwhm_target = fwhm_target[i] # Compute the correction in FWHM to be applied. fwhm_diff = np.sqrt(this_fwhm_target**2 - sl_models[sl_name]['fwhm']**2) # Transform this from Angstroem to resolution elements fwhm_diff = fwhm_diff/dlam # Compute the corresponding sigma for the gaussian profile sigma = (2*np.sqrt(2*np.log(2)))**-1 * fwhm_diff # In res. elements # Here's my gaussian kernel w = signal.gaussian(window_size,sigma,sym = True) # If the window lands outside the spectrum, do nothing ... if (i < window_size/2.) or (i > len(ssp_raw) - window_size/2.): ssp_out[i] = 0.0 else: # Here, I am using correlate, because it is exactly what # np.convolve does. It's a way to get the same numbers out, # (down to the machine precision limit). ssp_out[i] = np.correlate(w/w.sum(), ssp_raw[i-(window_size-1)/2: i+(window_size-1)/2+1], mode='valid') ''' # Some comparison with similar methods, to see if I # really understand what I do ... out_spec2 = np.convolve(w/w.sum(), raw_spec[:,1][i-(window_size-1)/2: i+(window_size-1)/2+1], mode='valid') # -> identical. Makes sense, since np.convolve wraps around # np.correlate. # Next, check the "official" gaussian_filter1d out_spec3 = ndimage.gaussian_filter1d(raw_spec[:,1][i-(window_size-1)/2: i+(window_size-1)/2+1], sigma) # -> different at 10**-6 level. Results are different # because I have a larger window size compared to what # gaussian_filter1d does by default. # But checking closer with the function inside # gaussian_filter1d, the results are actually identical, # when I use the same weights ! out_spec4 = ndimage.filters.correlate1d(raw_spec[:,1][i-(window_size-1)/2: i+(window_size-1)/2+1], w/w.sum()) import pdb pdb.set_trace() # Bottom line I'm doing things correctly. Slightly more # accurately than in the example of M. Cappellari. ''' ''' plt.close(1) plt.figure(1) plt.plot(lams,ssp_raw,'k-') plt.plot(lams,ssp_out,'r-') plt.show() ''' return ssp_out
[docs]def setup_spectral_library(velscale, inst, sl_name, fit_lims, z): ''' This prepares the set of stellar templates to be fed to ppxf. :param: velscale - the velocity scale for the log-rebinning :param: inst - the instrument that took the data, e.g. 'MUSE' :param: sl_name - the name of the stellar library to use :param: fit_lims - the spectral range for the fit (data & templates overlap) :param: z - redshift of target - needed for correct reshaping of fwhm_target. :return: ...lots of stuf... ''' # Fetch the stellar template files fns = glob.glob(os.path.join(sl_models[sl_name]['sl_loc'],'*')) fns.sort() # Extract the wavelength range and logarithmically rebin one spectrum # to the same velocity scale of the data spectra, to determine # the size needed for the array which will contain the template spectra. hdu = pyfits.open(fns[0]) ssp = hdu[0].data h2 = hdu[0].header lams_temp = np.arange(0,h2['NAXIS1'],1)*h2['CDELT1'] + h2['CRVAL1'] mask = (lams_temp >= fit_lims[0]) & (lams_temp <= fit_lims[1]) ssp = ssp[mask] lams_temp = lams_temp[mask] lam_range_temp = np.array([lams_temp[0],lams_temp[-1]]) log_lam, ssp_new, velscale = brutus_tools.log_rebin(lam_range_temp, ssp, velscale=velscale) # Create a three dimensional array to store the two dimensional grid of model spectra nAges = sl_models[sl_name]['nAges'] nMetal = sl_models[sl_name]['nMetal'] templates = np.empty((ssp_new.size, nAges, nMetal)) # These are the array where we want to store the characteristics of each SSP model logAge_grid = np.empty((nAges, nMetal)) metal_grid = np.empty((nAges, nMetal)) # These are the characteristics of the adopted rectangular grid of SSP models logAge = sl_models[sl_name]['logAge'] metal = sl_models[sl_name]['metal'] # Here we make sure the spectra are sorted in both [M/H] # and Age along the two axes of the rectangular grid of templates. metal_str = sl_models[sl_name]['metal_str'] for k, mh in enumerate(metal_str): files = [s for s in fns if 'Z'+mh in s] for j, filename in enumerate(files): # Resample the spectra to match the instrumental resolution ssp = sl_resample(filename, z=z, inst=inst, sl_name=sl_name, do_cap=True) ssp = ssp[mask] # Log-rebin the spectrum log_lam, ssp_new, velscale = brutus_tools.log_rebin(lam_range_temp, ssp, velscale=velscale) templates[:, j, k] = ssp_new # Templates are *not* normalized here logAge_grid[j, k] = logAge[j] metal_grid[j, k] = metal[k] return templates, lam_range_temp, logAge_grid, metal_grid, log_lam
#------------------------------------------------------------------------------
[docs]def ppxf_population(specerr, templates=None, velscale=None, start = None, goodpixels = None, plot=False, moments=4, degree=-1, vsyst=None, clean=False, mdegree=10, regul=None): ''' Function that calls ppxf, designed to be called from multiprocessing (i.e. only require 1 argument). Everything else is taken care of outside of this. ''' galaxy = specerr[0] noise = specerr[1] # Only run ppxf if I have a non-nan's spectrum. Even one nan's breaks ppxf ... sigh. if not(np.any(np.isnan(galaxy))): pp = ppxf(templates, galaxy, noise, velscale, start, goodpixels=goodpixels, plot=False, moments=moments, degree=degree, vsyst=vsyst, clean=clean, mdegree=mdegree, regul=regul, quiet=True) # In a perfect world, I would return pp. But that is a very massive variable. # 1 row = 300 spaxels = 5GB => 1 MUSE cube = 1.5TB !!! # Instead, just store what I need/trust. The best fit, the galaxy spectra (for # sanity checks) and the sol array (with the kinematics information, that I # sort-of trust. return [pp.bestfit, pp.galaxy, pp.sol] else: return None
#------------------------------------------------------------------------------