Download this as a Jupyter notebook!

GZ3D VAC usage within Marvin

This Notebook will demonstrate how to access and use GZ3D masks.

For more information on GZ3D see the VAC documentation website or Masters et al. (2021), which is the paper describing the project.

We will start by importing both `marvin <https://sdss-marvin.readthedocs.io/en/latest/>`__ and some `matplotlib <https://matplotlib.org/>`__ functions

import marvin
marvin.config.download = True
marvin.config.release = 'DR17'

%matplotlib inline
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec

print('Marvin and Matplotlib imported. Warnings about missing paths/files are typical, so no need to be concerned')

[INFO]: No release version set. Setting default to MPL-11
[WARNING]: redis-py works best with hiredis. Please consider installing (UserWarning)

Marvin and Matplotlib imported. Warnings about missing paths/files are typical, so no need to be concerned

The cell below just increases matplotlib font size, tick size, and line thickness to make for better plots.

style = {
    'lines.linewidth': 1.5,
    'font.size': 16.0,
    'axes.linewidth': 1.5,
    'axes.titlesize': 'large',
    'axes.labelsize': 'large',
    'xtick.major.size': 6,
    'xtick.minor.size': 4,
    'xtick.major.width': 1.5,
    'xtick.minor.width': 1.5,
    'xtick.major.pad': 6,
    'xtick.minor.pad': 6,
    'xtick.labelsize': 'medium',
    'ytick.major.size': 6,
    'ytick.minor.size': 4,
    'ytick.major.width': 1.5,
    'ytick.minor.width': 1.5,
    'ytick.major.pad': 6,
    'ytick.minor.pad': 6,
    'ytick.labelsize': 'medium',
    'legend.fancybox': True,
    'legend.fontsize': 'small',
    'legend.scatterpoints': 5,
    'legend.loc': 'best',
    'figure.titlesize': 'large',
}

plt.style.use(style)

Accessing the Files in the VAC

The GZ3D data can be accessed either through the `marvin.tools.Cube <https://sdss-marvin.readthedocs.io/en/latest/tools/cube.html#marvin-cube>`__ object or the `marvin.tools.Maps <https://sdss-marvin.readthedocs.io/en/latest/tools/maps.html>`__ object using the .vacs.gz3d attribute. This attribute will return a special GZ3DTarget class instance that will aid in processing and using the GZ3D VAC fits files. Calling print on this instance will show summary information about the GZ3D classifications.

For the example we include below, manga_id = 1-604761 this should return the below, which summarizes some data about the galaxy, how many GZ3D volunteers classified on each of the three tasks (center/star identification; drawing around the bar; drawing around the spirals), as well as how many clusters are returned for galaxy centers and foreground stars (if present). In this example all three tasks return masks. Some MaNGA galaxies do not have bars and/or spiral arms and will only have center/star counts.

Subject info:
    subject id: 14742283
    manga id: 1-604761
    ra: 113.47227478027344
    dec: 37.02590560913086
    ifu size: 127
Classification counts:
    15 center/star, 15 non_blank
    15 spiral, 14 non_blank
    15 bar, 15 non_blank
Cluster counts:
    1 center(s)
    3 star(s)

manga_id = '1-604761'

cube = marvin.tools.Cube(manga_id)
data = cube.vacs.gz3d
print(data)

[WARNING]: mismatch between file version=MPL-11 and object release=DR17. Setting object release to MPL-11 (MarvinUserWarning)

Subject info:
    subject id: 14742283
    manga id: 1-604761
    ra: 113.47227478027344
    dec: 37.02590560913086
    ifu size: 127
Classification counts:
    15 center/star, 15 non_blank
    15 spiral, 14 non_blank
    15 bar, 15 non_blank
Cluster counts:
    1 center(s)
    3 star(s)

Plotting the image and feature masks

Each GZ3D fits file includes the original (SDSS) image shown to the volunteers. A hexagon was overlayed to indicate roughly the region where spectra would be taken by the MaNGA IFU. Unfortunately the hexagon that was shown is slightly smaller than the true IFU size due to a bug in the code that made the cutout images (magenta), the correct IFU will be over layed (grey) by default when displaying the image. For those interested, the source of this bug was not accounting for the space between IFU fibers when making the image for GZ3D (it assumed each fiber was 2 arcsec on the sky when in reality they are 2.5 arcsec on the sky).

To plot the original image you can use the plot_image method and to plot the various masks you can use the plot_masks method. Below are some examples using these methods that shows how to pass in existing axis objects for embedding these in subplots. Note: plot_masks uses matplotlib’s `GridSpec <https://matplotlib.org/stable/tutorials/intermediate/gridspec.html>`__ internally when colorbars for each of the masks, so the example below will use this method for creating subplots but normal axis objects can be passed in on subplot_spec and it will work.

If you pass in a your own axis object make sure the projection is set to the correct WCS for the image (see example below).

#Setting up the plot
plt.figure(figsize=(23, 5))
gs = gridspec.GridSpec(1, 3, wspace=0.5)
ax0 = plt.subplot(gs[0], projection=data.wcs)

#Plot the orignal image shown to GZ3D volunteers
data.plot_image(color_grid='grey', ax=ax0)

#Plot a heat map showing vote counts for the four different tasks in GZ3D. We call these feature masks.
data.plot_masks(color_grid='grey', subplot_spec=gs[1]);

#As above, but overlay on the image
data.plot_masks(show_image=True, color_grid='grey', subplot_spec=gs[2]);

Image vs Spaxel masks

The GZ3D feature masks are created in the pixel coordinates of the image shown to the volunteers (which have ~0.1 arcsecond pixels), but to be useful for MaNGA maps they are projected into spaxel space (a spaxel is a spatial pixel in a MaNGA data cube, which are gridded on 0.5arcsecond spaxels). The image masks are found on the attributes bar_mask, spiral_mask, star_mask, and center_mask while the spaxel masks are found on bar_mask_spaxel, spiral_mask_spaxel, star_mask_spaxel, and center_mask_spaxel. Below are some plots showing how this process is done.

The spaxel masks are all the same size and orientation as the manga.tools.maps.Maps arrays and can be used directly as masks for those arrays.

Note: The WCS for MaNGA maps (and the spaxel masks) has the y-axis flipped compared to the image (and original mask) WCS. In the plot below the ylim of the axis is inverted to keep the orientation the same. The Dec values are correct, the origin is just flipped.

fig = plt.figure(figsize=(14, 5))
gs = gridspec.GridSpec(1, 2, wspace=0.4)

ax0 = plt.subplot(gs[0])
ax0.imshow(data.image)
# add the correct IFU hexagon to the image plot
ax0.add_patch(data.get_hexagon(edgecolor='C7'))

# plot the lines for the spaxel grid
v_grid, h_grid = data.get_spaxel_grid()
ax0.vlines(v_grid[0], ymin=v_grid[1][0], ymax=v_grid[1][1], color='white', alpha=0.2)
ax0.hlines(h_grid[1], xmin=h_grid[0][0], xmax=h_grid[0][1], color='white', alpha=0.2)

# zoom in on the image to show the grid clearly
ax0.set_xlim(70, 455)
ax0.set_ylim(70, 455)

# plot the spaxel masks
data.plot_masks(color_grid='grey', spaxel_masks=True, subplot_spec=gs[1])
ax1 = fig.axes[1]

# flip the y-axis to match the orientation of other plots
ax1.set_ylim(*ax1.get_ylim()[::-1]);

Placing GZ3D mask contours on marvin maps

You may wish to overlay contours at some appropriate vote threshold onto a MaNGA map.

To get an image in the same orientation as the maps you can use `marvin.tools.Image <https://sdss-marvin.readthedocs.io/en/latest/tools/image.html>`__, plot it along side the spaxel masks, and then use the masks as contours on top of any marvin map plot.

For these plots we demonstrate how to show one contour level, which in this example corresponds to the contour within which four or more volunteers marked a spaxel as belonging to each of the masks.

#Get the Image from Marvin
image = marvin.tools.Image(manga_id)

# Set up the Multiplot.
fig = plt.figure(figsize=(15, 12))
gs = gridspec.GridSpec(2, 2, wspace=0.4, hspace=0.6)

# Plot the original image as stored in Marvin
ax0 = plt.subplot(gs[0, 0], projection=image.wcs)
data._set_up_axes(ax0, color_grid='grey')
ax0.imshow(image.data)

# The masks plotted on the spaxel (0.5" spatial pixel) grid, which matches the MaNGA cube gridding
ax1 = plt.subplot(gs[0, 1])
data.plot_masks(color_grid='grey', spaxel_masks=True, subplot_spec=gs[1])

# Overlaying the masks as contours on the H-alpha map
ax2 = plt.subplot(gs[1, 0])
data.maps.emline_gflux_ha_6564.plot(ax=ax2, fig=fig)
ax2.contour(data.spiral_mask_spaxel, levels=[4], zorder=20, colors='C0')
ax2.contour(data.bar_mask_spaxel, levels=[4], zorder=20, colors='C1')
ax2.contour(data.center_mask_spaxel, levels=[4], zorder=20, colors='C2')
ax2.contour(data.star_mask_spaxel, levels=[4], zorder=20, colors='C4')

# Overlaying the masks as countours on the stellar velocity map
ax3 = plt.subplot(gs[1, 1])
data.maps.stellar_vel.plot(ax=ax3, fig=fig)
ax3.contour(data.spiral_mask_spaxel, levels=[4], zorder=20, colors='C0')
ax3.contour(data.bar_mask_spaxel, levels=[4], zorder=20, colors='C1')
ax3.contour(data.center_mask_spaxel, levels=[4], zorder=20, colors='C2')
ax3.contour(data.star_mask_spaxel, levels=[4], zorder=20, colors='C4');

[WARNING]: Warning: converting a masked element to nan. (UserWarning)
[WARNING]: Warning: converting a masked element to nan. (UserWarning)
[WARNING]: Warning: converting a masked element to nan. (UserWarning)

Plotting BPT diagrams

The GZ3D class has a special method for plotting BPT diagrams (a popular emission line ratio diagnostics plot for spectra from galaxies). This class method is called plot_bpt, and it produced BPT diagrams that are color coded by the feature mask each spaxel belongs to.

There are three kinds of BPT diagrams that can be made by passing in bpt_kind as either "log_nii_ha", "log_sii_ha", or "log_oi_ha" (where you can see the names of the emission lines used in those names). Other keywords are passed to matplotlib’s scatter function (in this example we pass the keyword, s, that controls the scatter point size).

Note: When making a legend sometimes there is transparency on the legend symbols, see the code below for an easy way to remove the alpha effect.

#Set up the plot grid
plt.figure(figsize=(18, 5))
gs = gridspec.GridSpec(1, 3, wspace=0.4)

#Plot the 'log_nii_ha' type BPT diagram in the first plot and include a legend with the feature identification
ax0 = plt.subplot(gs[0])
data.plot_bpt(ax=ax0, bpt_kind='log_nii_ha', s=3)
legend = ax0.legend(loc=3)

# remove transparency from legend
for lh in legend.legendHandles:
    lh.set_alpha(1)

#Plot another BPT diagram in the second plot
ax1 = plt.subplot(gs[1])
data.plot_bpt(ax=ax1, bpt_kind='log_sii_ha', s=3)

#Plot another BPT diagram in the third plot
ax2 = plt.subplot(gs[2])
data.plot_bpt(ax=ax2, bpt_kind='log_oi_ha', s=3);

Obtaining the mean spectra from spaxels in a feature mask

The get_mean_spectra method will calculate the weighted mean spectra inside each mask (weighted by classification counts). After running this method new attributes will become available for mean_spiral, mean_bar, and mean_center that will return marvin.tools.quantities.spectrum objects with the mean spectra. The error for these mean spectra is fully calculated and accounts for spacial covariance as outlined in Westfall et al. (2019) (the paper describing the MaNGA data analysis pipeline).

By using the inv=True keyword, get_mean_spectra will also calculate the mean spectra of the inverted masks (e.g. mean of all not spiral spaxels, etc..). This can be used for calculating difference spectrum to look for excess inside various masks over the background. This calculation can take a few minutes to finish for galaxies with larger IFUs.

data.get_mean_spectra(inv=True)

[WARNING]: `np.float` is a deprecated alias for the builtin `float`. To silence this warning, use `float` by itself. Doing this will not modify any behavior and is safe. If you specifically wanted the numpy scalar type, use `np.float64` here. (DeprecationWarning)
[WARNING]: `np.int` is a deprecated alias for the builtin `int`. To silence this warning, use `int` by itself. Doing this will not modify any behavior and is safe. When replacing `np.int`, you may wish to use e.g. `np.int64` or `np.int32` to specify the precision. If you wish to review your current use, check the release note link for additional information. (DeprecationWarning)

plt.figure(figsize=(18, 10))
gs = gridspec.GridSpec(3, 1, hspace=0.4)

# Mean spectra for Spiral, Bar, and Center masks
ax0 = plt.subplot(gs[0])

ax0.plot(data.mean_spiral.wavelength, data.mean_spiral, color='C0', label='Spiral')
ax0.plot(data.mean_bar.wavelength, data.mean_bar, color='C1', label='Bar')
ax0.plot(data.mean_center.wavelength, data.mean_center, color='C2', label='Center')

ax0.set_xlabel('Wavelength')
ax0.set_ylabel('Flux')
ax0.legend()

# Mean spectra for the inverse masks
ax1 = plt.subplot(gs[1])

ax1.plot(data.mean_not_spiral.wavelength, data.mean_not_spiral, color='C0', label='Not Spiral')
ax1.plot(data.mean_not_bar.wavelength, data.mean_not_bar, color='C1', label='Not Bar')
ax1.plot(data.mean_not_center.wavelength, data.mean_not_center, color='C2', label='Not Center')

ax1.set_xlabel('Wavelength')
ax1.set_ylabel('Flux')
ax1.legend()

# Difference spectra for each mask
ax2 = plt.subplot(gs[2])

ax2.plot(data.mean_not_spiral.wavelength, data.mean_spiral - data.mean_not_spiral, color='C0', label='$\\Delta\\mathrm{{Spiral}}$')
ax2.plot(data.mean_not_bar.wavelength, data.mean_bar - data.mean_not_bar, color='C1', label='$\\Delta\\mathrm{{Bar}}$')
ax2.plot(data.mean_not_center.wavelength, data.mean_center - data.mean_not_center, color='C2', label='$\\Delta\\mathrm{{Center}}$')

ax2.set_xlabel('Wavelength')
ax2.set_ylabel('$\\Delta\\mathrm{{Flux}}$')
ax2.legend();

Polar plots (Radial and/or Azimuthal Averages)

The final plot types supported by the GZ3D class are polar plots. These can either be by azimuthal angle, theta or normalized radius, R/R_50, where R_50 is the effective radius, or the radius containing 50% of the light from the SDSS image - specifically it is NSA_ELPETRO_TH50_R in the DAP file (see data model). Specify which you want by changing the x_unit keyword from theta to radius. Any individual MaNGA map layer can be passed in under the key keyword and you will need to provide the correct ylabel (units will automatically be added to the ylabel based on the map chosen).

The default map is D_n 4000, a spectral index measurement which correlates with stellar population ages (larger values indicate older stars).

plt.figure(figsize=(15, 10))
gs = gridspec.GridSpec(2, 2, hspace=0.4, wspace=0.4)

# D_n 4000 vs radius
ax0 = plt.subplot(gs[0, 0])
data.polar_plot(x_unit='radius', ax=ax0)
legend0 = ax0.legend(loc=3)

# remove transparency from legend
for lh in legend0.legendHandles:
    lh.set_alpha(1)

# D_n 4000 vs theta
ax1 = plt.subplot(gs[0, 1])
data.polar_plot(x_unit='theta', ax=ax1)

# EW(H-alpha) vs radius
ax2 = plt.subplot(gs[1, 0])
data.polar_plot(x_unit='radius', key='emline_sew_ha_6564', ylabel='EW(H\\alpha)', ax=ax2)
legend2 = ax2.legend(loc=2)

# remove transparency from legend
for lh in legend2.legendHandles:
    lh.set_alpha(1)

# EW(H-alpha) vs theta
ax3 = plt.subplot(gs[1, 1])
data.polar_plot(x_unit='theta', key='emline_sew_ha_6564', ylabel='EW(H\\alpha)', ax=ax3);

[WARNING]: np.asscalar(a) is deprecated since NumPy v1.16, use a.item() instead (DeprecationWarning)

Note that the definition of the zero-point for azimuthal angle, theta, is different for each galaxy, as they are set by the orientation of elliptical fits to the photometry (for more information see the explanation of the elliptical polar coordinates of each spaxel from the galaxy center in the DAP Datamodel). An easy way to check what direction the zero-point in azimuth points in RA, Dec co-ordinates is by plotting the theta map along side the spaxel masks as below:

fig = plt.figure(figsize=(13, 4))
gs = gridspec.GridSpec(1, 2, wspace=0.4)

# theta map
ax0 = plt.subplot(gs[0], projection=data.maps.wcs)
im = ax0.imshow(data.maps.spx_ellcoo_elliptical_azimuth.value)
# change the default units and labels on the axes
data._set_up_axes(ax0, color_grid='white')
# incluce a colorbar for the theta map
fig.colorbar(im, ax=ax0)

# plot the spaxel masks
data.plot_masks(color_grid='grey', spaxel_masks=True, subplot_spec=gs[1])
ax1 = fig.axes[1]