Mapping Galaxy Properties with MaNGA + HST

---

Learning Goals

By the end of this tutorial, you will:

• Understand how to search the MAST Archive and download SDSS MaNGA data using astroquery.mast

• Plot galaxy properties including H-\(\alpha\) emission line flux and stellar velocities with MaNGA

• Understand how to search for HST data complementing the MaNGA observations

• Create maps of galaxy properties by combining HST and MaNGA data

Table of Contents

• Introduction

• Imports

• Accessing MaNGA data at MAST

  • Querying All MaNGA data

  • Searching for a specific galaxy

  • Downloading MaNGA data products

  • Plotting velocity and flux maps from MaNGA

• Searching for HST data of this galaxy

  • Coordinate search using astroquery.mast

  • Colorizing HST images with astropy

• Combining MaNGA and HST data

  • Creating an H-alpha Emission Map

  • Plotting the stellar velocity field

• Additional Resources

  • How to Cite

  • About This Notebook

Introduction

The Mapping Nearby Galaxies at Apache Point Observatory (MaNGA) survey provides optical-wavelength integral field unit (IFU) spectroscopy for over 10,000 nearby galaxies. MaNGA collected data between 2014 - 2020 as part of the Sloan Digital Sky Survey (SDSS-IV) project. MaNGA data is now available at the Mikulski Archive for Space Telescopes (MAST) through the SDSS Legacy Archive at MAST.

In this notebook tutorial, we will demonstrate how to access MaNGA data at MAST using Python. One MaNGA observation, an interacting galaxy pair named MaNGA 7443-12703 will be used to demonstrate the basics of how to download and plot MaNGA data. We will then combine this MaNGA data with HST observations, also accessible from MAST, to study this galaxy pair in detail, exploring how its gas and stars are moving as the galaxies merge together.

Imports

The main packages we’re using for this notebook and their use-cases are:

• astropy.coordinates for handling astronomical coordinates

• astroquery.mast Observations for searching the MAST archive

• astropy.visualization for colorizing images

• astropy.wcs WCS for handling spatial footprints

• astropy.units for working with astronomical units

• astropy.io fits for accessing FITS files

• matplotlib.pyplot for plotting data

• numpy to handle array functions

• PIL for plotting and handling preview (png/jpg) images

• reproject for coordinate transformations and projections

%matplotlib inline

from astropy.coordinates import SkyCoord
from astroquery.mast import Observations
from astropy.visualization import make_lupton_rgb
from astropy.wcs import WCS

import astropy.units as u
import astropy.io.fits as fits
import matplotlib.pyplot as plt
import numpy as np
import PIL
import reproject

This cell updates some of the settings in matplotlib to use larger font sizes in the figures:

#Update Plotting Parameters
params = {'axes.labelsize': 12, 'xtick.labelsize': 12, 'ytick.labelsize': 12, 
          'text.usetex': False, 'lines.linewidth': 1,
          'axes.titlesize': 18, 'font.family': 'serif', 'font.size': 12}
plt.rcParams.update(params)

---

Accessing MaNGA data at MAST

The SDSS Legacy Archive at MAST hosts all of the science-ready data products from the SDSS-IV MaNGA Survey, which includes data for over 10,000 different nearby galaxies taken with the Apache Point Observatory SDSS-2.5m telescope. This notebook will demonstrate how to search and download MaNGA data using MAST.

Image Credit: Dana Berry / SkyWorks Digital Inc., David Law, and the SDSS collaboration.

Querying all MaNGA data

Searching for MaNGA data is straightforward with astroquery.mast. In this example, we use Observations.query_criteria and search for provenance_name = 'MaNGA'. This will return a table describing all of the MaNGA data hosted by the MAST archive.

Other useful search parameters for MaNGA data might include:

• obs_collection = 'SDSS': searches for all SDSS data

• target_classification = 'GALAXY': searches for only galaxies

• obs_id to search for specific targets

# Search for MaNGA data
manga_obs_list = Observations.query_criteria(provenance_name='MaNGA')

# Display First Ten Entries in Table
manga_obs_list[:10]

Searching for a specific galaxy

Let’s narrow down the search to one particular galaxy: MaNGA galaxy 7443-12703, also known as “Mrk 848B”. This is an interacting galaxy pair currently undergoing a merger event: the two galaxies are colliding and will eventually merge together into one large galaxy. Studying merger events like this one help us understand how galaxies grow and evolve over time.

We can search for this galaxy in particular using the obs_id keyword:

# Search for MaNGA galaxy 7443-12703
manga_obs_list = Observations.query_criteria(provenance_name='MaNGA',
                                             obs_id='sdss_manga_7443-12703')

# Display results
manga_obs_list

Downloading MaNGA data products

List all of the data products availble for this galaxy using Observations.get_product_list().

There are 13 total files available for this galaxy, which include include a preview image (12703.png), the 3D MaNGA spectral cube (manga-7443-12703-LOGCUBE.fits.gz), and the MaNGA MAP file (manga-7443-12703-MAPS-HYB10-MILESHC-MASTARSSP.fits.gz) which contains measurements and parameters from the MaNGA data analysis pipeline.

More information on all of these products can be seen in the search results table:

# Get list of products (files) associated with this observation 
manga_products_list = Observations.get_product_list(manga_obs_list)

# Show Products List
manga_products_list

Now we will download some of these files using Observations.download_products(). To start, let’s download the preview image to see what this galaxy pair looks like. The parameters used are:

• manga_products_list: the list of products from the previous cell

• productType='PREVIEW': limit the list to preview files only

• flat=True: download the data into the current directory

The download will print a status message when completed.

print(Observations.download_products(manga_products_list, productType='PREVIEW', flat=True))

The Right Ascension (RA) and declination (Dec) coordinates for this image are contained in the file as a World Coordinate System (WCS) axis. Using astropy.wcs.WCS(), we can retrieve this information from the file and use it in the plot.

def get_wcs_from_png(file):
    """ Retrieves WCS information from png file """
    with PIL.Image.open(file) as ii:
        ww = {}
        for k, v in ii.info.items():
            try:
                ww[k] = float(v)
            except ValueError:
                ww[k] = v
        ww.pop("WCSAXES")
        return WCS(ww)   

Now let’s plot this preview image to see what these galaxies look like!

# Open image file
manga_preview = np.asarray(PIL.Image.open('12703.png'))

# Set up plot using WCS projection from the file
plt.subplot(projection=get_wcs_from_png('12703.png'))

# Show the image
plt.imshow(manga_preview)

# Label title and axes
plt.title('MaNGA Preview Image')
plt.xlabel('Right Ascension')
plt.ylabel('Declination')

plt.show()

Next, let’s download the Minimum Recommended Products. These are the recommended files to download for this MaNGA galaxy. You can select the Minimum Recommended Products using mrp_only=True. This will download 4 files.

# Download Minimum Recommended Products
Observations.download_products(manga_products_list, mrp_only=True, flat=True)

Plotting velocity and flux maps from MaNGA

Now that the MaNGA data has been downloaded, we can plot some of the galaxy properties. The ManGA MAP file (manga-7443-12703-MAPS-HYB10-MILESHC-MASTARSSP.fits.gz) contains the output and results from the MaNGA Data Analysis Pipeline (DAP). This example is based off of the MaNGA DAP Python Tutorial, and plots the g-band flux, the H-alpha emission line flux, and the stellar velocity field from MaNGA.

The full description of contents of the the MaNGA MAP file is available here.

Let’s open the file and view some basic information:

# Open MaNGA MAP file
manga_map = fits.open('manga-7443-12703-MAPS-HYB10-MILESHC-MASTARSSP.fits.gz')

# Print information on file extensions
manga_map.info()

To make the plot in the next cell, we will use three different extensions of the MAP file. Here is the description of each of these extensions from the data model:

• SPX_MFLUX

  • Description: “g-band-weighted mean flux, not corrected for Galactic extinction or internal attenuation.”

• The EMLINE_GFLUX

  • Description: “gaussian profile integrated flux from a combined continuum+emission-line fit. The flux ratio of the [NeIII], [OIII], [OI], [NII], and [S III] lines are fixed and cannot be treated as independent measurements. The emission-line fluxes account for Galactic reddening using the E(B-V) (copied to the DAP primary headers, see the EBVGAL header keyword) value provided by the DRP header and assuming an O’Donnell (1994, ApJ, 422, 158) reddening law; however, no attenuation correction is applied for dust internal to the galaxy. See here for more information.”

  • Specifically, we will be using the H-alpha emission line flux from this extension

• STELLAR_VEL

  • Description: “Line-of-sight stellar velocity, relative to the input guess redshift (given as cz by the keyword SCINPVEL in the header of the PRIMARY extension, and most often identical to the NSA redshift).”

#=========================================
# Set up Plot
#=========================================
plt.figure(figsize=(12, 6))
# Grab the WCS information from the header
manga_wcs = WCS(manga_map['SPX_MFLUX'].header)
ax1 = plt.subplot(131, projection=manga_wcs)
ax2 = plt.subplot(132, projection=manga_wcs)
ax3 = plt.subplot(133, projection=manga_wcs)

for ax in [ax1, ax2, ax3]:
    ax.set_xlabel('Right Ascension')
    ax.set_ylabel('Declination')


#=========================================
# Subplot 1: MaNGA Flux Map
#=========================================
# The 'SPX_MFLUX' ext contains the g-band-weighted mean flux
manga_flux = manga_map['SPX_MFLUX'].data
manga_flux[manga_flux == 0] = np.nan # mask for quality data

# Plot image
im = ax1.imshow(manga_flux, cmap='Greys_r')
ax1.set_title('MaNGA g-band Flux')
plt.colorbar(im, label=r'Flux [$10^{-17}$ erg/s/cm$^{2}$/Å/spaxel]',
             orientation='horizontal')
             

#=========================================
# Subplot 2: MaNGA H-Alpha Emission Map
#=========================================

# Define the emission line indexes for this ext
emline = {}
for k, v in manga_map['EMLINE_GFLUX'].header.items():
    if k[0] == 'C':
        try:
            i = int(k[1:])-1
        except ValueError:
            continue
        emline[v] = i
        
# The 'EMLINE_GFLUX' ext contains the emission line measurements
h_alpha_flux = np.copy(manga_map['EMLINE_GFLUX'].data[emline['Ha-6564']])
h_alpha_flux[h_alpha_flux == 0] = np.nan # mask for quality data

# Plot image
im = ax2.imshow(h_alpha_flux, cmap='magma', vmin=0, vmax=200)
plt.colorbar(im, label=r'H-$\alpha$ Emission [$10^{-17}$ erg/s/cm$^{2}$/spaxel]',
             orientation='horizontal')

ax2.set_title(r'H-$\alpha$ Emission')

#=========================================
# Subplot 3: MaNGA Stellar Velocity Field
#=========================================
# The 'STELLAR_VEL' ext contains the stellar velocity measurements
qual_mask = manga_map['STELLAR_VEL'].header['QUALDATA'] # mask for quality data
velocity_map = np.ma.MaskedArray(manga_map['STELLAR_VEL'].data,
                                 mask=manga_map[qual_mask].data > 0)

# plot Image
im = ax3.imshow(velocity_map, interpolation='nearest', 
                vmin=-125, vmax=125, cmap='RdBu_r')

ax3.set_title(r'Stellar Velocity')
plt.colorbar(im, label=r'Stellar Velocity [km/s]',
             orientation='horizontal')

#plt.subplots_adjust(hspace=0)

plt.tight_layout()
plt.show()

---

Searching for HST observations of this galaxy

Coordinate search using astroquery.mast

Now let’s search for HST obersvations of this same galaxy. Similar to before, we will be using Observations.query_criteria() to search the MAST archive, but this time, we will search for HST observations (obs_collection='HST') near the coordinates of this MaNGA galaxy pair.

# Retrieve RA and DEC of MaNGA observations
ra = manga_obs_list['s_ra'][0]
dec = manga_obs_list['s_dec'][0]
# make a SkyCoord object from these coordinates
coord = SkyCoord(ra=ra*u.deg, dec=dec*u.deg)
print(coord)

# Search for HST observations based on coordinates
hst_obs = Observations.query_criteria(# Search by coordinates
                                      coordinates=coord, 
                                      # Search for HST observations
                                      obs_collection='HST',
                                      # Select only Science observations (not calibration files)
                                      intentType='science',
                                      # Select calibrated reduced observations
                                      provenance_name='CAL*') 

# Display Results
hst_obs

The table of results shown above tells us that there are three observations of this galaxy from HST, in three different filters: F160W, F435W, and F814W. Let’s see what data products are available associated with these observations:

hst_products = Observations.get_product_list(hst_obs)

# Filter product list 
hst_products = Observations.filter_products(hst_products,
                                            # Select science files
                                            productType="SCIENCE",
                                            # Recommended Products
                                            productGroupDescription="Minimum Recommended Products", 
                                            # Select DRZ files - the calibrated combined images
                                            productSubGroupDescription='DRZ')
hst_products

With our products selected, we can proceed to download. We’ll use flat=True to put them all into the same directory.

# Download products
Observations.download_products(hst_products, flat=True)

Plot the HST images

Let’s plot the three HST images.

# Open files
f160w = fits.open('ia1e42010_drz.fits')
f435w = fits.open('j9cv55010_drz.fits')
f814w = fits.open('j9cv55020_drz.fits')


f160w_wcs = WCS(f160w[1].header)
f435w_wcs = WCS(f435w[1].header)
f814w_wcs = WCS(f814w[1].header)

#=========================================
# Set up Plot
#=========================================
fig = plt.figure(figsize=(15, 5))
ax1 = plt.subplot(131, projection=f160w_wcs)
ax2 = plt.subplot(132, projection=f435w_wcs)
ax3 = plt.subplot(133, projection=f814w_wcs)

for ax in [ax1, ax2, ax3]:
    ax.set_xlabel('Right Ascension')
    ax.set_ylabel('Declination')

#=========================================
# Subplot 1: F160W filter image
#=========================================
# Note - this image is upside-down compared to the other two
fits_file = f160w
ax1.imshow(fits_file[1].data, cmap='Greys_r', vmin=0, vmax=2.5, origin='lower')
ax1.set_title(f"{fits_file[0].header['INSTRUME']}/{fits_file[0].header['FILTER']}")

#=========================================
# Subplot 2: F435W filter image
#=========================================
fits_file = f435w
wcs = WCS(fits_file[1].header)
ax2.imshow(fits_file[1].data, cmap='Greys_r', vmin=0, vmax=0.5, origin='lower')
ax2.set_title(f"{fits_file[0].header['INSTRUME']}/{fits_file[0].header['FILTER2']}")

#=========================================
# Subplot 3: F814W filter image
#=========================================
fits_file = f814w
wcs = WCS(fits_file[1].header)
ax3.imshow(fits_file[1].data, cmap='Greys_r', vmin=0, vmax=0.5, origin='lower')
ax3.set_title(f"{fits_file[0].header['INSTRUME']}/{fits_file[0].header['FILTER2']}")

plt.tight_layout()
plt.show()

Colorizing HST images with astropy

We can combine all three HST images here and colorize them using the make_lupton_rgb() function from astropy. This function takes image data from three filter images, and combines them as an RGB image. In this example, we will order by wavelength and map the IR filter (F160W) to the red channel (‘r’), the F814W filter to the green channel (‘g’), and the F435W filter to the blue channel (‘b’).

Before we colorize the images, however, we need to reproject them onto the same coordinate system. We will resample the F160W and F814W images to the same projection as the F435W image, because the F435W image is the largest.

# Print image shapes:
print(f"F160W shape: {f160w[1].data.shape}")
print(f"F814W shape: {f814w[1].data.shape}")
print(f"F435W shape: {f435w[1].data.shape}")

# Reproject all three HST images into the same frame (using F435W image as base)
# This may take about a minute.
r, _ = reproject.reproject_interp(f160w[1], f435w[1].header)
g, _ = reproject.reproject_interp(f814w[1], f435w[1].header)
b, _ = reproject.reproject_interp(f435w[1], f435w[1].header)

# Colorize image using the three filters
hst_image = make_lupton_rgb(r*0.1, g, b*2.5, Q=4, stretch=0.75)

Our reprojection and colorization is complete! Let’s now plot the image.

# Make plot
plt.figure(figsize=(5, 10))
ax = plt.subplot(projection=f435w_wcs)

# Plot image
ax.imshow(hst_image)

# Zoom in
ax.set_xlim(1800, 2600)
ax.set_ylim(2000, 2900)

# Label Plot
ax.set_xlabel('Right Ascension')
ax.set_ylabel('Declination')
ax.set_title('HST Image (colorized)')
plt.show()

This image looks great! There are certainly a few places where bright spots in one image are affecting the final result; this can be seen as the green speckles across the lower third of the image.

---

Combining MaNGA and HST data

Now that we have downloaded both MaNGA and HST data, let’s combine them into one plot to map some of properties of this merging pair of galaxies!

The key to combining the HST and MaNGA data is the World Coordinate System (WCS) transformations done by astropy. Using information from the file headers, astropy.wcs() will calculate the RA and DEC coordinates corresponding to each pixel in the image, allowing us to project both datasets on the same figure.

# Store WCS for easy coordinate transformations
hst_wcs = WCS(f435w[1].header)
manga_wcs = WCS(manga_map[3].header)

#=========================================
# Set up Plot
#=========================================
plt.figure(figsize=(15, 5))

ax1 = plt.subplot(131, projection=hst_wcs)
ax2 = plt.subplot(132, projection=hst_wcs)
ax3 = plt.subplot(133, projection=hst_wcs)

for ax in [ax1, ax2, ax3]:
    # Zoom in
    ax.set_xlim(1850, 2550)
    ax.set_ylim(2300, 3000)
    # Label Axes
    ax.set_xlabel('RA')
    ax.set_ylabel('Dec')
    # Plot Grid
    ax.grid(color='white', ls='dotted', alpha=0.5)

    
#=========================================
# Subplot 1: HST Image
#=========================================
ax1.imshow(hst_image)
ax1.set_title('HST Image (colorized)')

#=========================================
# Subplot 2: MaNGA H-alpha Emission
#=========================================
ax2.set_title(r'MaNGA H-$\alpha$ Flux' + '\n' + r'[$10^{-17}$ erg/s/cm$^{2}$/spaxel]')

h_alpha_flux = np.copy(manga_map['EMLINE_GFLUX'].data[emline['Ha-6564']])
h_alpha_flux[h_alpha_flux == 0] = np.nan # mask for quality data


# Plot H-alpha flux
ax2.imshow(np.log10(h_alpha_flux),
           transform=ax2.get_transform(manga_wcs),
           cmap='magma', vmin=0, vmax=2)

# Plot MaNGA contours
contour_levels = [1, 10, 20, 50, 100]
contour_labels = [str(c) for c in contour_levels]
fmt = {}
for label_level, label_string in zip(contour_levels, contour_labels):
    fmt[label_level] = label_string

contours = ax2.contour(h_alpha_flux,
                       transform=ax2.get_transform(manga_wcs),
                       levels=contour_levels, colors='white')
ax2.clabel(contours, contours.levels, inline=True, fmt=fmt, fontsize=12)


#=========================================
# Subplot 3: HST + MaNGA
#=========================================
ax3.set_title(r'HST + MaNGA H-$\alpha$ Contours')

# Plot background HST image
ax3.imshow(hst_image)
    
# Plot MaNGA contours
contours = ax3.contour(h_alpha_flux,
                       transform=ax3.get_transform(manga_wcs),
                       levels=contour_levels, colors='white')
ax3.clabel(contours, contours.levels, inline=True, fmt=fmt, fontsize=8)


plt.tight_layout()
plt.show()

This plot encapsulates everything we have learned so far, showing the colorized HST image (left), the H-alpha flux from MaNGA (middle), and combining them both into a single image (right).

Creating an H-alpha Emission Map

We can isolate the last panel and plot it by itself, too.

plt.figure(figsize=(10, 10))
ax = plt.subplot(projection=hst_wcs)
ax.set_xlim(1950, 2550)
ax.set_ylim(2350, 2900)
plt.xlabel(r'Right Ascension')
plt.ylabel(r'Declination')
plt.title(r'HST Image + H-$\alpha$ Emission from MaNGA')

# Plot the HST image
ax.imshow(hst_image)


# Plot MaNGA contours
contour_levels = [1, 10, 20, 50, 100]
contour_labels = [str(c) for c in contour_levels]
contour_labels[0] = r'1 $\times10^{-17}$ erg/s/cm$^{2}$/spaxel'
fmt = {}
for label_level, label_string in zip(contour_levels, contour_labels):
    fmt[label_level] = label_string

contours = ax.contour(h_alpha_flux,
                      transform=ax.get_transform(manga_wcs),
                      levels=contour_levels, colors='white')
ax.clabel(contours, contours.levels, inline=True, fmt=fmt, fontsize=14)

plt.savefig('manga_halpha_map.png')
plt.show()

In this figure, we can see that the H-alpha flux from MaNGA follows the spiral arms of the top galaxy, and the emission is strongest near the center of both galaxies.

Plotting the stellar velocity field

Last but not least, let’s do the same thing with the stellar velocity! The stellar velocity field from MaNGA will show which parts of the galaxies are redshifted (moving away from us, which we will plot in red), and which parts are blueshifted (moving toward us, which we will plot in blue).

plt.figure(figsize=(14, 10))
ax = plt.subplot(projection=WCS(f435w[1].header))

ax.imshow(f435w[1].data, vmin=0, vmax=0.5, cmap='Greys_r',
          transform=ax.get_transform(WCS(f435w[1].header)))

ax.set_xlim(1950, 2550)
ax.set_ylim(2350, 2900)

bin_indx = manga_map['BINID'].data[1]
unique_bins, unique_indices = tuple(map(lambda x: x[1:],
                                        np.unique(bin_indx.ravel(), return_index=True)))

x_pix = np.array([x for y in range(74) for x in range(74)])[unique_indices]
y_pix = np.array([y for y in range(74) for x in range(74)])[unique_indices]
v_map = velocity_map.ravel()[unique_indices]

# Get the luminosity-weighted x and y coordinates of the unique bins
im = ax.scatter(x_pix, y_pix, c=v_map,
                marker='.', s=100, lw=0,
                vmin=-150, vmax=150, cmap='seismic', 
                transform=ax.get_transform(manga_wcs))

plt.colorbar(im, label='Stellar Velocity [km/s]')


plt.xlabel(r'Right Ascension')
plt.ylabel(r'Declination')
plt.title(r'Stellar Velocity from MaNGA')

plt.savefig('manga_velocity_map.png')
plt.show()

This plot shows that this galaxy pair is rotating vertically (from the perspective of this image). The arm in the upper galaxy is blueshifted, while most of the lower galaxy, and the end of the arm is redshifted.

---

Congratulations! You have reached the end of this tutorial notebook. You have learned how to access and download MaNGA data from MAST, and combine it with HST images to map different properties for this merging pair of galaxies.

Additional Resources

Additional resources are linked below:

• SDSS Legacy Archive at MAST

• SDSS-IV MaNGA Survey

• astroquery.mast User Manual

• MAST API

Citations

If you use MaNGA data for published research, see the following links for information on which citations to include in your paper:

• Citing SDSS

• Citing MaNGA

• Citing MAST

• Citing astropy

About this Notebook

For questions or issues with this notebook, you can open a github issue or email archive@stsci.edu.

Authors: Julie Imig, Brian Cherinka
Keywords: Tutorial, SDSS, MaNGA, HST, galaxies
First published: October 2024
Last updated: October 2024

---

Top of Page