Using SimSpin v2.0.0

In this walk-through, we are going to demonstrate how you can post-process SimSpin FITS files using observational software pPXF.

This example will show you how to:

Produce a FITS file with SimSpin. For a longer demonstration and specific details about producing these files, check out this link.
Generate kinematic fits for each spectra in the data cube.
Construct a line-of-sight (LOS) velocity and LOS velocity dispersion map of the observation.

For this walk-through, you will need an installation of Python3.x, R, SimSpin and an installation of pPXF. Assuming you have a suiatble installation of Python, you can simply install pPXF using pip install ppxf. Similarly, SimSpin can be installed using the command devtools::install_github('kateharborne/SimSpin') within R.

Producing a FITS file with SimSpin

From R, we can run the following script to generate and observation of the example galaxy included within the package.

library(SimSpin)

## Loading required package: foreach

example_galaxy = system.file("extdata",  # using the example particle file included in SimSpin 
                             "SimSpin_example_Gadget", 
                             package = "SimSpin")
example_spectra = make_simspin_file(filename      = example_galaxy, 
                                    template      = "EMILES", 
                                    write_to_file = FALSE)
SAMI = SimSpin::telescope(type="SAMI")            # defining the telescope parameters
strategy = SimSpin::observing_strategy(z = 0.05)  # defining the observation parameters
temp_loc = tempdir()                              # temporary location for saving FITS
output_FITS_path = paste0(temp_loc, "_example.FITS")

example_cube = SimSpin::build_datacube(simspin_file = example_spectra, # generating the data cube observation
                                       telescope = SAMI,
                                       observing_strategy = strategy,
                                       write_fits = TRUE,              # saving the output to a FITS file
                                       output_location = output_FITS_path) # at this PATH location

## FITS file written to: /var/folders/45/wfwdmf2d3fdc9_ffnwcjsgrhhqjvjp/T//RtmpBeKBXY_example.FITS

Below, we plot an example of the physical LOS velocities and velocity dispersion of the particles within the model. These maps are measured from the physical particle velocities within each spaxel region of the simulation, and are not the same as the fitted values that we expect to recover, as these do not include observational effects.

A single pixel is also highlighted, for which we shown the spectrum contained in that spaxel below.

par(options(scipen = 0), xpd=FALSE, ps=12, cex=1)  

magicaxis::magplot(example_cube$observation$wave_seq,
                   example_cube$spectral_cube[15,15,]/median(example_cube$spectral_cube[15,15,]),
                   type="l", lty=1, lwd=2, col ="blue",
                   ylab= "Normalised Flux", xlab = "Wavelength, Ang")

Having generated this FITS file, we can then analyse this data cube using pPXF.

Analysing a spectrum using pPXF

For the analysis of the stellar kinematics, we switch to using python in order to use the pPXF fitting routines. Let’s start taking that single spectrum from the cube above and fitting using the templates within pPXF.

Here, we have modified the “ppxf_example_kinematics_sauron.py” script included with the pPXF distribution for use with the SimSpin cube produced using the “SAMI” telescope. Comments in the chunk below are taken from the original Cappellari example, with additions relating to the analysis of SAMI/SimSpin data.


# We begin by importing all of the packages required for this analysis. 

import glob
from os import path
from astropy.io import fits
from scipy import ndimage
import numpy as np
from ppxf import ppxf
import ppxf.ppxf_util as util
import matplotlib.pyplot as plt

c = 299792.458 # the speed of light

# Next, we open the FITS file for analysis

hdu = fits.open(r.output_FITS_path) # opening the FITS file described in the "R" chunk above 
gal_lin = hdu[0].data # pulling in the spectral data cube
h1 = hdu[0].header    # pulling in the header information

lamRange1 = h1['CRVAL3'] + (np.array([-h1['CRPIX3'], h1['CRPIX3']])*h1['CDELT3']) # wavelength range
FWHM_gal = 2.65  # SAMI has an instrumental resolution FWHM of 2.65A.

# If the galaxy is at a significant redshift (z >= 0.03), one would need to apply
# a large velocity shift in pPXF to match the template to the galaxy spectrum.
# This would require a large initial value for the velocity (V >= 1e4 km/s)
# in the input parameter START = [V,sig]. 
#
# An alternative solution consists of bringing the galaxy spectrum roughly to the
# rest-frame wavelength, before calling pPXF. In practice there is no
# need to modify the spectrum before the usual LOG_REBIN, given that a
# redshift corresponds to a linear shift of the log-rebinned spectrum.
# One just needs to compute the wavelength range in the rest-frame
# and adjust the instrumental resolution of the galaxy observations.
#
# This is done with the following three lines:
#
z = 0.05                    # Initial estimate of the galaxy redshift
lamRange1 = lamRange1/(1+z) # Compute approximate restframe wavelength range
FWHM_gal = FWHM_gal/(1+z)   # Adjust resolution in Angstrom

# Now we initialise the velocity scale and rebin the wavelength space into a 
# logarithmic scale. 
galaxy, logLam1, velscale = util.log_rebin(lamRange1, gal_lin[:, 15, 15]) 
galaxy = galaxy/np.median(galaxy)  # Normalize spectrum to avoid numerical issues
noise  = np.full_like(galaxy, 0.01) # Assume constant noise per pixel 

# We also need to pull in the template spectra with which we will fit the observed cube.
# pPXF is distributed included some template files, specifically a subset of the MILES templates. 
# For more information, check out: https://www-astro.physics.ox.ac.uk/~mxc/software/#ppxf

file_dir = '/Users/00094926/.pyenv/versions/3.9.1/envs/SimSpin/lib/python3.9/site-packages/ppxf'
vazdekis = glob.glob(file_dir + '/miles_models/Mun1.30Z*.fits') 
# where the variable "file_dir" describes the location of your pPXF installation
FWHM_tem = 2.51  # Vazdekis+10 spectra have a constant resolution FWHM of 2.51A.
velscale_ratio = 2  # adopts 2x higher spectral sampling for templates than for galaxy

# Extract the wavelength range and logarithmically rebin one spectrum to determine
# the size needed for the array which will contain the template spectra.
hdu = fits.open(vazdekis[0])
ssp = hdu[0].data
h2 = hdu[0].header
lamRange2 = h2['CRVAL1'] + np.array([0., h2['CDELT1']*(h2['NAXIS1'] - 1)])
sspNew, logLam2, velscale_temp = util.log_rebin(lamRange2, ssp, velscale=velscale/velscale_ratio)
templates = np.empty((sspNew.size, len(vazdekis)))

# Convolve the whole Vazdekis library of spectral templates
# with the quadratic difference between the instrument and the
# Vazdekis instrumental resolution. Logarithmically rebin
# and store each template as a column in the array TEMPLATES.

# Quadratic sigma difference in pixels Vazdekis --> SAMI
# The formula below is rigorously valid if the shapes of the
# instrumental spectral profiles are well approximated by Gaussians.
#
FWHM_dif = np.sqrt(FWHM_gal**2 - FWHM_tem**2)
sigma = FWHM_dif/2.355/h2['CDELT1']  # Sigma difference in pixels

for j, file in enumerate(vazdekis):
    hdu = fits.open(file)
    ssp = hdu[0].data
    ssp = ndimage.gaussian_filter1d(ssp, sigma)
    sspNew, logLam2, velscale_temp = util.log_rebin(lamRange2, ssp, velscale=velscale/velscale_ratio)
    if np.median(sspNew) != 0:
        templates[:, j] = sspNew/np.median(sspNew)  # Normalizes templates
    else:
        templates[:, j] = 0.

# Now fitting the spectra at the highlighted spaxel:
dv = (np.mean(logLam2[:velscale_ratio]) - logLam1[0])*c

start = [100., 200.]  # (km/s), starting guess for [V, sigma]

pp = ppxf.ppxf(templates, galaxy, noise, velscale, start,
               plot=True, moments=2, quiet=False,
               degree=4, vsyst=dv, velscale_ratio = velscale_ratio)

##  Best Fit:       Vel     sigma
##  comp.  0:       -80       144
## chi2/DOF: 16.91; degree = 4; mdegree = 0
## method = capfit; Jac calls: 5; Func calls: 17; Status: 2
## linear_method = lsq_box; Nonzero Templates (>0.1%):  4 / 150

plt.show()

Here, we can see that pPXF can successfully recover the kinematics of the spectra. We now need to run this process for every spaxel in the cube in order to recover the LOS velocity and LOS velocity dispersion at each spaxel and construct a similar velocity and dispersion map to those above.

Analysing the SimSpin data cube using pPXF

To do this, we have written two functions that perform this analysis in a looped fashion. These functions are shown below, but are also publically available as a python script on GitHub here.


# Function for generating the list of templates at the same resolution as the observation.
def degrade_templates(template_dir, FWHM_templates, FWHM_observation, velscale_obs, vel_ratio=1):
    """
    Function that combines pPXF "log_rebin" and "gaussian_filter1d" convolution 
    to return a list of template spectra that match the resolution and gridding of
    the observed data. 

    :param template_dir: A character string describing the path to the directory
        containing the raw template files. 
    :param FWHM_templates: A float describing the resolution of the templates.
    :param FWHM_observation: A float describing the resolution of the observation.
        This can be different from the instrumental resolution if the galaxy is at
        a signifcant redshift. 
    :param vescale_obs: A float describing the velocity scale in km/s per pixels of
        the observed spectra.
    :param vel_ratio: A float describing if you would like the velocity scale of 
        the templates to be sampled at a higher rate. 

    :return: An array of templates that have been degraded and rebinned to match
        the observations.

    """

    vazdekis = glob.glob(template_dir) # returns list of all template files contained within directory
    FWHM_tem = FWHM_templates
    FWHM_gal = FWHM_observation
    velscale = velscale_obs

    hdu = fits.open(vazdekis[0]) # opening the first template file in the list
    ssp = hdu[0].data            #  and measuring the wavelength range and SSP 
    h2 = hdu[0].header           #  size for initiallising the template array.
    lamRange2 = h2['CRVAL1'] + np.array([0., h2['CDELT1']*(h2['NAXIS1'] - 1)]) # Gives the range of the spectra wavelengths
    sspNew, logLam2, velscale_temp = util.log_rebin(lamRange2, ssp, velscale=velscale/vel_ratio)
    templates = np.empty((sspNew.size, len(vazdekis)))
    FWHM_dif = np.sqrt(FWHM_gal**2 - FWHM_tem**2)
    sigma = FWHM_dif/2.355/h2['CDELT1']  # Sigma difference in resolutions

    for j, file in enumerate(vazdekis):
        hdu = fits.open(file)
        ssp = hdu[0].data
        ssp = ndimage.gaussian_filter1d(ssp, sigma)
        sspNew, logLam2, velscale_temp = util.log_rebin(lamRange2, ssp, velscale=velscale/vel_ratio)
        if np.median(sspNew) != 0:
            templates[:, j] = sspNew/np.median(sspNew)  # Normalizes templates
        else:
            templates[:, j] = 0.

    return templates, logLam2

def apply_pPXF(spectra_file, fwhm_instrument, redshift, noise, vel_ratio=1, template_dir=None, fwhm_templates=None, templates=None, logLam2=None):
    """
    Function to fit galaxy spectra using pPXF.

    If provided with the template directory and resolution, this function will run the degrade_templates()
    function in order to match the resolution of the observation and the templates. Else, if the degraded
    templates are already provided, this function will just use the templates provided to fit the observed 
    galaxy spectra. 

    :param spectra_file: A string describing the path to the SimSpin spectra file.
    :param fwhm_instrument: A float describing the resolution of the observed spectra.
    :param redshift: A float describing the redshift, z, of the observed galaxy.
    :param noise: A float describing the level of noise within the observed spectra.
    :param vel_ratio: A float describing the sampling rate of the templates relative
        to the observed spectra. Default value is 1.
    
    If wishing to degrade templates to the appropriate resolution, 

    :param template_dir: A string describing the path to the directory in which the template files are 
        located. Default is None.
    :param fwhm_templates: A float describing the resolution of the template spectra. Default 
        is None.

    If you have already degraded the templates to the appropriate resolution for a previous 
    fit, you can provide these variables directly to the function to avoid recalculating the
    comupationally expensive degradation and rebinning:

    :param templates: A matrix containing the rebinned and degraded templates. Default is None.
    :param logLam2: An array describing the wavelength labels of the templates. Default is None.

    """
    hdu = fits.open(spectra_file)
    dim = hdu[0].data.shape
    mid = np.array([round(dim[1]/2), round(dim[2]/2)])
    
    gal_lin = hdu[0].data # pulling in the spectra for each pixel
    h1 = hdu[0].header
    lamRange1 = h1['CRVAL3'] + (np.array([-h1['CRPIX3'], h1['CRPIX3']])*h1['CDELT3']) # wavelength range
    FWHM_gal = fwhm_instrument  # SAMI has an instrumental resolution FWHM of 2.65A.
    
    z = redshift      # Initial estimate of the galaxy redshift
    lamRange1 = lamRange1/(1+z) # Compute approximate restframe wavelength range
    FWHM_gal = FWHM_gal/(1+z)   # Adjust resolution in Angstrom

    galaxy, logLam1, velscale = util.log_rebin(lamRange1, gal_lin[:, mid[0], mid[1]]) 

    if not templates or not logLam2:
        assert isinstance(template_dir, str) and isinstance(fwhm_templates, float), 'Please provide the path to the template directory and template resolution' 
        templates, logLam2 = degrade_templates(template_dir, fwhm_templates, FWHM_observation=FWHM_gal, velscale_obs=velscale, vel_ratio=vel_ratio)

    velocity = np.empty([dim[1],dim[2]])
    dispersion = np.empty([dim[1],dim[2]])

    for x in range(0,dim[1]):
            for y in range(0,dim[2]):
                
                if all(np.isnan(gal_lin[:,x,y])):
                    velocity[x,y] = None
                    dispersion[x,y] = None
                    
                else:
                    galaxy, logLam1, velscale = util.log_rebin(lamRange1, gal_lin[:,x,y])
                    galaxy = galaxy/np.median(galaxy)    # Normalize spectrum to avoid numerical issues
                    noise  = np.full_like(galaxy, noise) # Assume constant noise per pixel here
                    
                    if velscale_ratio > 1:
                        dv = (np.mean(logLam2[:velscale_ratio]) - logLam1[0])*c  # km/s
                    else:
                        dv = (logLam2[0] - logLam1[0])*c  # km/s

                    start = [100, 200.]  # (km/s), starting guess for [V, sigma]

                    pp = ppxf.ppxf(templates, galaxy, noise, velscale, start,
                                plot=False, moments=2, quiet=True,
                                degree=4, vsyst=dv, velscale_ratio = vel_ratio)
                                
                    velocity[x,y] = pp.sol[0]
                    dispersion[x,y] = pp.sol[1]

    return velocity, dispersion

Having defined these functions, we can then use them to analyse the spectral data cube above.


kinematics = apply_pPXF(spectra_file=r.output_FITS_path, fwhm_instrument=2.65, redshift=0.05, noise=0.01, vel_ratio=2, template_dir='/Users/00094926/.pyenv/versions/3.9.1/envs/SimSpin/lib/python3.9/site-packages/ppxf/miles_models/Mun1.30Z*.fits', fwhm_templates=2.51)

Examining this output in the same manner that we plotted the particle properties at the start:

# Pulling the velocity and dispersion images from the output 
velocity_fit = py$kinematics[[1]]
dispersion_fit = py$kinematics[[2]]

# finding the range of values in map for colour map plotting
vel_val = max(c(abs(min(velocity_fit, na.rm = T)), abs(max(velocity_fit, na.rm = T)))) 
disp_val = c(floor(min(dispersion_fit, na.rm=T)/2), max(dispersion_fit, na.rm=T))

par(pty="s", fig=c(0,0.5,0,1), options(scipen = 999), xpd=FALSE, ps=12, cex=1)  
image_nan(z = t(velocity_fit), zlim = c(-vel_val,vel_val), col = velo_map_cols, na.color = "white",
          xaxt="n", yaxt="n", ann=FALSE, magmap=FALSE, family="mono", font=1)
fields::image.plot(legend.only = TRUE, zlim = c(-vel_val,vel_val), col = velo_map_cols,
                   horizontal = TRUE, family="serif", font=1,
                   legend.lab = expression("velocity"[LOS] * ", km s"^{-1}))

par(pty="s", fig=c(0.5,1,0,1), options(scipen = 999), xpd=FALSE, ps=12, cex=1, new=TRUE)  
image_nan(z = t(dispersion_fit), zlim = disp_val, col = disp_map_cols, na.color = "white", xaxt="n",
          yaxt="n", ann=FALSE, magmap=FALSE, family="mono", font=1)
fields::image.plot(legend.only = TRUE, zlim = c(disp_val), col = disp_map_cols,
                   horizontal = TRUE, family="serif", font=1,
                   legend.lab = expression("dispersion"[LOS] * ", km s"^{-1}))

The resulting images look very similar to those produced using the raw particle data. This is positive as it means we can be confident adding further observational constraints, such as seeing conditions, and trusting the output data as a result.

In this walk-through, we have demonstrated the ease with which SimSpin observations can be processed using observational software pPXF. This process has been packaged up and hosted on GitHub as SimSpin-post-process. We have shown you how to:

Generate an observation of a particle simulation of a galaxy.
Use observational pipeline, pPXF, to analyse the output SimSpin data cube.
Produce LOS velocity and LOS velocity dispersion maps using the pPXF fits.

If you are interested in learning more about the SimSpin mock-IFS software, please see the references below or check out our GitHub page. For questions, queries or complaints, please log an issue there.

References

SimSpin Harborne, Power & Robotham (2020) PASA, 37, e016; Astrophysics Source Code Library, record ascl:1903.006
pPXF Cappellari & Emsellem (2004) PASP, 116, 138C; Cappellari (2017) MNRAS, 466, 798, Astrophysics Source Code Library, record ascl:1210.002
MILES templates - Vazdekis et al., 2016, MNRAS, 463, 3409

SimSpin - Kinematic fitting using pPXF

Katherine Harborne

Using SimSpin v2.0.0

Producing a FITS file with SimSpin

Analysing a spectrum using pPXF

Analysing the SimSpin data cube using pPXF

References