MaNGA Python Tutorial

Back to MaNGA tutorials

This tutorial will use the MaNGA galaxy 12-193481 (PLATE-IFU = 7443-12703). You can download the RSS file and datacube for this galaxy yourself with the following rsync commands in your terminal (replace LOCAL_PATH):

RSS

rsync -avz rsync://data.sdss.org/dr16/manga/spectro/redux/v2_4_3/7443/stack/manga-7443-12703-LOGRSS.fits.gz LOCAL_PATH

Cube

rsync -avz rsync://data.sdss.org/dr16/manga/spectro/redux/v2_4_3/7443/stack/manga-7443-12703-LOGCUBE.fits.gz LOCAL_PATH

This tutorial will first describe how to access RSS files with Python, and then continue with a description on how to access datacubes. Examples of plotting a spectrum and an Hα map are included.

The differences between the RSS and datacube files are here.

RSS

Import python modules and load data. Note that the extensions for the MaNGA RSS files are listed in its datamodel.

import numpy as np
import matplotlib.pyplot as plt
from astropy.io import fits
from matplotlib import cm

rss = fits.open('manga-7443-12703-LOGRSS.fits.gz')
flux_rss = rss['FLUX'].data
ivar_rss = rss['IVAR'].data
mask_rss = rss['MASK'].data
wave_rss = rss['WAVE'].data
xpos = rss['XPOS'].data
ypos = rss['YPOS'].data

Plot spectrum of fiber closest to (x=0, y=0) at H$\alpha$:

# Create masked array to skip plotting of bad pixels
bad_bits = (mask_rss != 0)
flux_rss_m = np.ma.array(flux_rss, mask=bad_bits)

# index element of central fiber for this particular data cube is 1313
ind_center = 1313    
plt.plot(wave_rss, flux_rss_m[ind_center])
plt.xlabel('$\lambda \, [\AA]$')
plt.ylabel(rss['FLUX'].header['BUNIT'])
plt.show()

rss_spec

Cube

Import python modules and load data, an overview of all extensions in the MaNGA data cube format is given in its datamodel.

import os
import numpy as np
import matplotlib.pyplot as plt
from astropy import wcs
from astropy.io import fits

cube = fits.open('manga-7443-12703-LOGCUBE.fits.gz')

# Re-order FLUX, IVAR, and MASK arrays from (wavelength, DEC, RA) to (RA, DEC, wavelength).
flux = np.transpose(cube['FLUX'].data, axes=(2, 1, 0))
ivar = np.transpose(cube['IVAR'].data, axes=(2, 1, 0))
mask = np.transpose(cube['MASK'].data, axes=(2, 1, 0))

wave = cube['WAVE'].data
flux_header = cube['FLUX'].header

Note: For convenience, we have reordered the axes of the data arrays to the intended ordering of (x,y,λ); see the discussion of array indexing on the Caveats page. In the flux, ivar, and mask arrays, (x=0, y=0) corresponds to the upper left if North is up and East is left.

Plot the central spectrum:

x_center = np.int(flux_header['CRPIX1']) - 1
y_center = np.int(flux_header['CRPIX2']) - 1

plt.plot(wave, flux[x_center, y_center])
plt.xlabel('$\lambda \, [\AA]$')
plt.ylabel(flux_header['BUNIT'])
plt.show()

cube_spec

Find array indices of a particular RA and DEC:

cubeWCS = wcs.WCS(flux_header)
ra = 229.525580000
dec = 42.7458420000
x_cube_coord, y_cube_coord, __ = cubeWCS.wcs_world2pix([[ra, dec, 1.]], 1)[0]
x_spaxel = np.int(np.round(x_cube_coord)) - 1
y_spaxel = np.int(np.round(y_cube_coord)) - 1

Plot an H$\alpha$ narrow band image from the datacube

Apply bitmasks and select region around H$\alpha$:

do_not_use = (mask & 2**10) != 0
flux_m = np.ma.array(flux, mask=do_not_use)

redshift = 0.0402719
ind_wave = np.where((wave / (1 + redshift) > 6550) & (wave / (1 + redshift) < 6680))[0]
halpha = flux_m[:, :, ind_wave].sum(axis=2)
im = halpha.T

# Convert from array indices to arcsec relative to IFU center
dx = flux_header['CD1_1'] * 3600.  # deg to arcsec
dy = flux_header['CD2_2'] * 3600.  # deg to arcsec
x_extent = (np.array([0., im.shape[0]]) - (im.shape[0] - x_center)) * dx * (-1)
y_extent = (np.array([0., im.shape[1]]) - (im.shape[1] - y_center)) * dy
extent = [x_extent[0], x_extent[1], y_extent[0], y_extent[1]]

Note: How to find the redshift of a galaxy from the drpall file

Generate plot:

plt.imshow(im, extent=extent, cmap=cm.YlGnBu_r, vmin=0.1, vmax=100, origin='lower', interpolation='none')
plt.colorbar(label=flux_header['BUNIT'])
plt.xlabel('arcsec')
plt.ylabel('arcsec')
plt.show()

Halpha_map