MaNGA Python Tutorial
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()
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()
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()