DAP Python Tutorial

Back to MaNGA tutorials

Disclaimer: This tutorial teaches you how to look at the DAP output files using standard python packages. However, we highly recommend you consider using Marvin, a python package designed specifically for downloading, visualizing, and analyzing MaNGA data.

The goal of this introductory tutorial will be to show you the basics of loading and visualizing MaNGA DAP products. We assume you're using ipython or something similar. Let's first import the packages we'll need and turn on interactive plotting:

import numpy
from astropy.io import fits
from matplotlib import pyplot
pyplot.ion()

Loading a MAPS file

Now, read in a maps file. For this tutorial, we'll use the maps file for MaNGA object (plate-ifu) 7443-12703.

hdu = fits.open(dir+'manga-7443-12703-MAPS-HYB10-GAU-MILESHC.fits.gz')

where "dir" is a string specifying where this data is located on your computer. This file contains several extensions. To see them all, type:

hdu.info()

For a full description of these extensions, see the DAP documentation.

Making an emission line map and applying masks

Let's make a simple map of Hα flux. The Hα flux measurements are stored in the EMLINE_GFLUX extension. If we examine the size of the data in this extension:

hdu['EMLINE_GFLUX'].data.shape

we'll see that it has multiple layers. Each layer holds the flux measurements for a different emission line. We can see which index corresponds to which line by looking in the header:

hdu['EMLINE_GFLUX'].header

Hα is in channel 19, although since we're in python where indices start at 0, we want index 18:

flux_halpha = hdu['EMLINE_GFLUX'].data[18,:,:]

It may be more convenient to make a dictionary that maps the different line names to their corresponding index:

emline = {}
for k, v in hdu['EMLINE_GFLUX'].header.items():
    if k[0] == 'C':
        try:
            i = int(k[1:])-1
        except ValueError:
            continue
        emline[v] = i
print(emline) 

So now we can select the Hα flux map using:

flux_halpha = hdu['EMLINE_GFLUX'].data[emline['Ha-6564'],:,:]

Now we can plot the map with a colorbar:

pyplot.clf()
pyplot.imshow(hdu['EMLINE_GFLUX'].data[emline['Ha-6564'],:,:],cmap='inferno', origin='lower', interpolation='nearest')
pyplot.colorbar(label=r'H$\alpha$ flux ($1\times10^{-17}$ erg s$^{-1}$ spaxel$^{-1}$ cm$^{-2}$)')
pyplot.show()
 An example map of Hα emission.
An example map of Hα emission.

Typically the maps contain some spaxels with unreliable mesaurements. The MAPS files contains extensions which identify these spaxels. The header of the EMLINE_GFLUX extension tells you which extension holds its mask.

mask_extension = hdu['EMLINE_GFLUX'].header['QUALDATA']

We can use this extension to make a masked image

masked_image = numpy.ma.array(hdu['EMLINE_GFLUX'].data[emline['Ha-6564'],:,:],
                              mask=hdu[mask_extension].data[emline['Ha-6564'],:,:]>0)
pyplot.imshow(masked_image, origin='lower', cmap='inferno', interpolation='nearest')
pyplot.show()
 An example Hα emission map with a mask applied to remove unreliable data.
An example Hα emission map with a mask applied to remove unreliable data.

Plotting a velocity field and creating your own masks

Now, let's use the same basic procedure to plot the ionized gas velocity field.

mask_ext = hdu['EMLINE_GVEL'].header['QUALDATA']
gas_vfield = numpy.ma.MaskedArray(hdu['EMLINE_GVEL'].data[emline['Ha-6564'],:,:], mask=hdu[mask_ext].data[emline['Ha-6564'],:,:] > 0)
pyplot.clf()
pyplot.imshow(gas_vfield, origin='lower', interpolation='nearest', vmin=-125, vmax=125, cmap='RdBu_r')
pyplot.colorbar()
 An ionized gas velocity field
An ionized gas velocity field

Depending on your science, you may want to only examine regions with very well-detected emission lines. To do this we need to read in the flux uncertainties, which are given as inverse variances. We can find the proper extension with

ivar_extension = hdu['EMLINE_GFLUX'].header['ERRDATA']

Let's now calculate the S/N per spaxel and use it to reject spaxels where the line flux has S/N < 10

snr_map = hdu['EMLINE_GFLUX'].data[emline['Ha-6564'],:,:]*numpy.sqrt(hdu[ivar_extension].data[emline['Ha-6564'],:,:])
sncut = 10
gas_vfield_alt = numpy.ma.MaskedArray(hdu['EMLINE_GVEL'].data[emline['Ha-6564'],:,:],
                              mask=(hdu[mask_ext].data[emline['Ha-6564'],:,:]>0) | (snr_map[:,:] < sncut))
pyplot.clf()
pyplot.imshow(gas_vfield_alt, origin='lower', interpolation='nearest', vmin=-125, vmax=125, cmap='RdBu_r')
pyplot.colorbar()
 An ionized gas velocity field limited to regions where the Hα flux has S/N>10.
An ionized gas velocity field limited to regions where the Hα flux has S/N>10.

An important note regarding gas velocity fields: The velocities of the emission lines are tied together, so for example, the velocity of the [OIII]-5007 line is the same as the Hα line, as are the uncertainties. You cannot reduce the uncertainty on the measured velocity by averaging the velocities of several lines together.

Plotting stellar velocity dispersion

Next we'll make a map of the stellar velocity dispersion. First, we'll access the raw dispersion measurements:

disp_raw = hdu['STELLAR_SIGMA'].data

However, we need to correct the measured dispersion for the instrumental resolution, which is reported in a different extension":

disp_inst = hdu['STELLAR_SIGMACORR'].data

Now, let's apply the correction and plot the results (also removing masked values). The calculation below will ignore any points where the correction is larger than the measured dispersion:

disp_stars_corr = numpy.sqrt(numpy.square(disp_raw) - numpy.square(disp_inst))
mask_ext = hdu['STELLAR_SIGMA'].header['QUALDATA']
disp_stars_final = numpy.ma.MaskedArray(disp_stars_corr,mask = (hdu[mask_ext].data > 0)) 
pyplot.clf()
pyplot.imshow(disp_stars_final,origin='lower',interpolation='nearest',cmap='RdBu_r')
pyplot.colorbar()
 Stellar velocity dispersion after correction for instrumental spectral resolution.
Stellar velocity dispersion after correction for instrumental spectral resolution.

Plotting spectral-index measurements

Next we'll show how to plot spectral indices and correct them for velocity dispersion. The spectral indices are stored in the "SPECINDEX" extension of the maps file. Like the emission line measurements, We can create a dictionary which makes accessing different spectral indices easier. We'll also track their units, which will be relevant shortly:

spec_ind = {}
spec_unit = numpy.empty(numpy.shape(hdu['SPECINDEX'].data)[0],dtype=object)
for k, v in hdu['SPECINDEX'].header.items():
    if k[0] == 'C':
        try:
            i = int(k[1:])-1
        except ValueError:
            continue
        spec_ind[v] = i
    if k[0] == 'U':
        try:
            i = int(k[1:])-1
        except ValueError:
            continue
        spec_unit[i] = v.strip()

Let's make a masked array holding the Hβ spectral index

mask_ext = hdu['SPECINDEX'].header['QUALDATA']
hb_raw = numpy.ma.MaskedArray(hdu['SPECINDEX'].data[spec_ind['Hb'],:,:],mask=hdu[mask_ext].data[spec_ind['Hb'],:,:]>0)

The spectral index measurements need to be corrected for velocity dispersion. Keep in mind that the way in which the corrections are applied depends on whether the units are angstroms or magnitudes. Hβ is in Angstroms:

corr = hdu['SPECINDEX_CORR'].data[spec_ind['Hb'],:,:]
hb_corr = hb_raw*corr
pyplot.clf()
pyplot.imshow(hb_corr,origin='lower', cmap='inferno', interpolation='nearest')
pyplot.colorbar(label=spec_unit[spec_ind['Hb']])
 Hβ spectral index measurement after applying the correction for velocity dispersion.
Hβ spectral index measurement after applying the correction for velocity dispersion.

Identifying unique bins

The spaxels are binned in different ways depending on the measurement being made (the DAP documentation provides more information). This binning means that two spaxels can belong the same bin, and therefore a derived quantity at those locations will be identical. The BINID extension provides information about which spaxels are in which bins. There are 5 channels providing the IDs of spaxels associated with

  • 0: each binned spectrum. Any spaxel with BINID=-1 as not included in any bin.
  • 1: any binned spectrum with an attempted stellar kinematics fit.
  • 2: any binned spectrum with emission-line moment measurements.
  • 3: any binned spectrum with an attempted emission-line fit.
  • 4: any binned spectrum with spectral-index measurements.

For any analysis, you'll want to extract the unique spectra and/or maps values. For instance, to find the indices of the unique bins where stellar kinematics were fit:

bin_indx = hdu['binid'].data[1,:,:]
unique_bins, unique_indices = tuple(map(lambda x : x[1:], numpy.unique(bin_indx.ravel(), return_index=True)))

Let's now use this information to plot the position of each unique bin and color-code it by the measured stellar velocity:

pyplot.clf()

# Get the x and y coordinates and the stellar velocities
x = hdu['BIN_LWSKYCOO'].data[0,:,:].ravel()[unique_indices]
y = hdu['BIN_LWSKYCOO'].data[1,:,:].ravel()[unique_indices]
v = numpy.ma.MaskedArray(hdu['STELLAR_VEL'].data.ravel()[unique_indices],
                         mask=hdu['STELLAR_VEL_MASK'].data.ravel()[unique_indices] > 0)

pyplot.scatter(x, y, c=v, vmin=-150, vmax=150, cmap='RdBu', marker='.', s=30, lw=0, zorder=3)
pyplot.colorbar()

Positions of each unique stellar velocity measurement from the binned spectra, color-coded by value.
Positions of each unique stellar velocity measurement from the binned spectra, color-coded by value.

Extract a binned spectrum and its model fit

The model cube files provide detailed information about the output the binned spaxels and the model fitting. We can use these files to compare an individual bin's measured spectrum, model fit, model residuals, and so on. However, one needs to be careful when comparing binned spectra with the model fits. Specifically, there are two types of files with different binning schemes:

  • VOR10: The spectra are voronoi binned to S/N~10. Stellar and emission line parameters are estimated from those bins.
  • HYB10: The spectra are again voronoi binned to S/N~10, and the stellar parameters are calculated using these voronoi binned spectra. However, emission line parameters are measured using the individual 0.5"x0.5" spaxels.

If you are using the HYB10 files and want to compare the best fitting model (including stellar continuum and emission lines) to the data, you need to compare the models to the individual spectra (measured in 0.5"x0.5" spaxels) from the DRP LOGCUBE files, not the binned spectra in the HYB10 files. If you are comparing the best-fitting stellar continuum models to the data, you should use the binned spectra within the HYB10 files.

Below are a few examples using both VOR10 and HYB10 files to demonstrate the proper way to compare the models and data using these two files types.

VOR10 Files

We'll start with the VOR10 files where comparing the models and data is simpler. First, let's read in the BIN_SNR extension from the VOR10 maps file and find the bin with the highest S/N.

hdu = fits.open(path+'manga-7443-12703-MAPS-VOR10-GAU-MILESHC.fits.gz')
snr = numpy.ma.MaskedArray(hdu['BIN_SNR'].data, mask=hdu['BINID'].data[0,:,:] < 0)
j,i = numpy.unravel_index(snr.argmax(), snr.shape)

Next we'll load the modelcube file for this galaxy, pull out the binned spectrum at its location, and then plot the binned spectrum, the full best fit model, the model stellar continuum, the model emission lines, and the residuals.

hdu_cube = fits.open(dir+'manga-7443-12703-LOGCUBE-VOR10-GAU-MILESHC.fits.gz')
wave=hdu_cube['wave'].data
flux = numpy.ma.MaskedArray(hdu_cube['FLUX'].data[:,j,i],mask=hdu_cube['MASK'].data[:,j,i] > 0)
model = numpy.ma.MaskedArray(hdu_cube['MODEL'].data[:,j,i],mask=hdu_cube['MASK'].data[:,j,i]>0)
stellarcontinuum = numpy.ma.MaskedArray(hdu_cube['MODEL'].data[:,j,i] - hdu_cube['EMLINE'].data[:,j,i] - hdu_cube['EMLINE_BASE'].data[:,j,i], mask=hdu_cube['MASK'].data[:,j,i] > 0)
emlines = numpy.ma.MaskedArray(hdu_cube['EMLINE'].data[:,j,i],mask=hdu_cube['EMLINE_MASK'].data[:,j,i] > 0)
resid = flux-model-0.5
pyplot.clf()
pyplot.step(wave, flux, where='mid', color='k', lw=0.5,label='flux')
pyplot.plot(wave, model, color='r', lw=1,label='model')
pyplot.plot(wave, stellarcontinuum, color='g', lw=1,label='stellar cont.')
pyplot.plot(wave, emlines, color='b', lw=1,label='Emission lines')
pyplot.step(wave, resid, where='mid', color='0.5', lw=0.5,label='residuals')
pyplot.legend()
 Example of one binned spectrum in a MaNGA data cube.  Lines show the binned flux, full model fit, model stellar continuum, model emission lines, and model fit residuals
Example of one binned spectrum in a MaNGA data cube. Lines show the binned flux, full model fit, model stellar continuum, model emission lines, and model fit residuals

HYB10 Files

Identifying the bin with the highest S/N is done the same way as for the VOR10 files.

hdu = fits.open(dir+'manga-7443-12703-MAPS-HYB10-GAU-MILESHC.fits.gz')
snr = numpy.ma.MaskedArray(hdu['BIN_SNR'].data, mask=hdu['BINID'].data[0,:,:] < 0)
j,i = numpy.unravel_index(snr.argmax(), snr.shape)

Recall that although stellar parameters are measured using the voronoi bins, the emission line parameters are estimated using the individual spaxels. Therefore, if we want to compare the data to the full best-fitting model which includes stellar continuum and emission lines, we need to use the spectra from the DRP LOGCUBE file, not the DAP HYB10 LOGCUBE file. Let's compare the best-fitting model to the data at the position (j,i) found above:

hdu_cube = fits.open(dir+'manga-7443-12703-LOGCUBE-HYB10-GAU-MILESHC.fits.gz')  #DAP MODELCUBE file
hdu_drpcube = fits.open(dir+'manga-7443-12703-LOGCUBE.fits.gz') #DRP LOGCUBE file
flux = numpy.ma.MaskedArray(hdu_drpcube['FLUX'].data[:,j,i],mask=hdu_drpcube['MASK'].data[:,j,i] > 0)
model = numpy.ma.MaskedArray(hdu_cube['MODEL'].data[:,j,i],mask=hdu_cube['MASK'].data[:,j,i]>0)
pyplot.clf()
pyplot.step(wave, flux, where='mid', color='k', lw=0.5,label='flux')
pyplot.plot(wave, model, color='r', lw=1,label='model')
pyplot.legend()

If we just want to compare the best fitting stellar continuum to the data, we should use the binned spectra within the DAP HYB10 LOGCUBE file:

flux = numpy.ma.MaskedArray(hdu_cube['FLUX'].data[:,j,i],mask=hdu_cube['MASK'].data[:,j,i] > 0)
stellarcontinuum = numpy.ma.MaskedArray(hdu_cube['MODEL'].data[:,j,i] - hdu_cube['EMLINE'].data[:,j,i] - hdu_cube['EMLINE_BASE'].data[:,j,i], mask=hdu_cube['MASK'].data[:,j,i] > 0)
pyplot.clf()
pyplot.step(wave, flux, where='mid', color='k', lw=0.5,label='flux')
pyplot.plot(wave, stellarcontinuum, color='g', lw=1,label='stellar cont.')
pyplot.ylim(1.2,2.2)
pyplot.legend()

Now go use Marvin!