from __future__ import print_function, division, absolute_import
from marvin.tools import maps
from .base import VACMixIn, VACTarget
import numpy as np
import matplotlib.patches as patches
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
import matplotlib as mpl
import scipy.linalg as sl
import json
import marvin.utils.dap.bpt as bpt
import marvin
from astropy.io import fits
from astropy.wcs import WCS
from astropy.table import Table
from scipy.interpolate import RectBivariateSpline
from scipy.spatial import distance_matrix
from marvin.tools.quantities.spectrum import Spectrum
from marvin import log
LUT = {7: 3, 19: 5, 37: 7, 61: 9, 91: 11, 127: 13}
spaxel_grid = {7: 24, 19: 34, 37: 44, 61: 54, 91: 64, 127: 74}
[docs]def convert_json(table, column_name):
'''Unpacks the JSON column of a table
Paramters:
table (astropy.table.Table):
An astropy table
column_name (str):
The name of the column made up of JSON strings
The input table is updated in place by appending `_string` to
the end of the JSON column name and adding a new column with
`_list` on the end with the list representation of the same
column.
'''
# this unpacks the json column of a table
new_col = [json.loads(i) for i in table[column_name]]
table.rename_column(column_name, '{0}_string'.format(column_name))
table['{0}_list'.format(column_name)] = new_col
[docs]def non_blank(table, *column_name):
'''Count how many non-blank classifications are in the given columns
of the input table.
Paramters:
table (astropy.table.Table):
An astropy table with Zooniverse classifications
column_name(s) (str):
One or multiple column names
Returns:
non_blank (int):
The total number of non-blank classifications across *all*
input column names (combined with "logical or" in a single
row).
'''
for cdx, c in enumerate(column_name):
if cdx == 0:
non_blank = np.array([len(i) > 0 for i in table[c]])
else:
non_blank = non_blank | np.array([len(i) > 0 for i in table[c]])
return non_blank.sum()
[docs]def cov_to_ellipse(cov, pos, nstd=1, **kwargs):
'''Create a covariance ellipse given an covariance matrix and postion
Paramters:
cov (numpy.array):
2x2 covariance matrix
pos (numpy.array):
1x2 center position of the ellipse
Keywords:
nstd (int):
Number of standard deviations to make the output ellipse (Default=1)
kwargs:
All other keywords are passed to matplotlib.patches.Ellipse
Returns:
ellipse (matplotlib.patches.Ellipse):
matplotlib ellipse patch object
'''
eigvec, eigval, V = sl.svd(cov, full_matrices=False)
# the angle the first eigenvector makes with the x-axis
theta = np.degrees(np.arctan2(eigvec[1, 0], eigvec[0, 0]))
# full width and height of ellipse, not radius
# the eigenvalues are the variance along the eigenvectors
width, height = 2 * nstd * np.sqrt(eigval)
return patches.Ellipse(xy=pos, width=width, height=height, angle=theta, **kwargs)
[docs]def alpha_overlay(C_a, a_a, C_b, a_b=None):
'''Take a base color (C_a), an alpha map (a_a), background image (C_b), and optional
background alpha map (a_b) and overlay them.
Paramters:
C_a (numpy.array):
1x3 RGB array for the base color to be overlayed
a_a (numpy.array):
NxM array of alpha values for each postion on an image
C_b (numpy.array):
1x3 RGB array for the background color or NxMx3 RGB array
for a background image
a_b (numpy.array):
NxM array of alpha values for the background color/image
Returns:
c_out (numpy.array):
NxMx3 RGB array containing the alpha overlayed image.
'''
if a_b is None:
a_b = np.ones(a_a.shape)
c_a = np.array([a_a.T] * 3).T * C_a
c_b = np.array([a_b.T] * 3).T * C_b
c_out = c_a + ((1 - a_a.T) * c_b.T).T
return c_out
[docs]def alpha_maps(maps, colors=None, vmin=0, vmax=15, background_image=None):
'''Take a list of color masks and base color values
and make an alpha-mask overlay image.
Parameters:
maps (list):
List of masks to use as alpha maps
Keywords:
colors (list):
What matplotlib color to use for each of the input maps
(defaults to standard MPL color cycle)
vmin (int):
Value in the maps at or below this value will be 100% transparent
vmax (int):
Value in the maps at or above this value will be 100% opaque
background_image (numpy.array):
RGB array to use as the background image (default solid white)
Returns:
overlay (numpy.array):
RGB array with each map overlayed on each other with alpha
transparency.
'''
norm = mpl.colors.Normalize(vmin=vmin, vmax=vmax, clip=True)
iter_cycle = iter(mpl.rcParams['axes.prop_cycle'])
for mdx, m in enumerate(maps):
if colors is None:
c = next(iter_cycle)['color']
else:
c = colors[mdx]
base_color = np.array(mpl.colors.to_rgb(c))
norm_map = norm(m)
if mdx == 0:
if background_image is None:
background_color = np.ones(3)
else:
background_color = background_image
background_color = alpha_overlay(base_color, norm_map, background_color)
return background_color
[docs]def make_alpha_bar(color, vmin=-1, vmax=15):
'''Make a matplotlib color bar for a alpha mask of a single color
Parameters:
color (string):
A matplotlib color (any format matplotlib accepts)
Keywords:
vmin (int):
The minimum value for the colorbar. Default value is -1
to ensure the labels show up correctly when used with
plot_alpha_bar.
vmax (int):
The maximum value for the colorbar. Default value is 15.
Returns:
colormap (mpl.colors.ListedColormap):
The colormap for the colorbar
norm (mpl.colors.Normalize):
The normalization for the color bar
'''
norm = mpl.colors.Normalize(vmin=vmin, vmax=vmax, clip=True)
a_a = norm(range(vmin, vmax))
C_a = np.array(mpl.colors.to_rgb(color))
new_cm = alpha_overlay(C_a, a_a, np.ones(3))
return mpl.colors.ListedColormap(new_cm), norm
[docs]def make_alpha_color(count, color, vmin=1, vmax=15):
'''Give a matplotlib color and alpha channel proportional to
the input count value.
Parameters:
count (int):
The count value used to select an alpha value
color (string):
A matplotlib color (any format matplotlib accepts)
Keywords:
vmin (int):
The count value to be associated with transparent.
Default is 1.
vmax (int):
The count value to be associated with opaque. Default
is 15.
Returns:
alpha_color (tuple):
An rgba tuple for the new alpha color
'''
norm = mpl.colors.Normalize(vmin=vmin, vmax=vmax, clip=True)
return mpl.colors.to_rgb(color) + (norm(count), )
[docs]def plot_alpha_bar(color, grid, ticks=[]):
'''Display and alpha colorbar on a plot grid.
Parameters:
color (string):
A matplotlib color (any format matplotlib accepts)
grid (matplotlib.gridspec.SubplotSpec):
A gridspec subplot specification to place the color bar in
Keywords:
ticks (list):
A list of tick value for the colorbar
Returns:
ax_bar (matplotlib.axes.Axes):
Matplotlib axis object for the colorbar
colorbar (mpl.colorbar.ColorbarBase):
Matplotlib colorbar object
'''
bar, norm = make_alpha_bar(color)
ax_bar = plt.subplot(grid)
cb = mpl.colorbar.ColorbarBase(ax_bar, cmap=bar, norm=norm, orientation='vertical', ticks=ticks)
cb.outline.set_linewidth(0)
return ax_bar, cb
[docs]def plot_alpha_scatter(x, y, mask, color, ax, snr=None, sf_mask=None, value=True, **kwargs):
'''Make a scatter plot where each x-y point has and alpha transparency
set by the values in a count mask array.
Parmeters:
x (numpy.array):
1-D numpy array with x-values to be plotted
y (numpy.array or spectral line object from a Marvin Maps cube):
1-D numpy array with y-values to be plotted
mask (numpy.array):
1-D numpy array with mask array containing the "count" value for each
(x,y) data point
color (string):
A matplotlib color (any format matplotlib accepts) used for the
base color of the data points
ax (matplotlib.axes.Axes):
The maplotlib axes to use for the plot
Keywords:
snr (float):
Minimum signal to noise ratio to use as a cutoff for the y-values.
Defaults to `None`. Only used if `value=True`
sf_mask (numpy.array):
1-D numpy array with A star formation region mask that is 1 when there
is star formation in a spaxel and 0 otherwise. If passed in only spxels
where this mask is 1 will be plotted.
value (bool):
If True y is a spectral line object from a Marvin Maps cube, otherwise
y is assumed to be regular np.array object.
**kwargs:
All other keywords are passed forward to matplotlib's scatter plot
function.
Returns:
scatter (matplotlib.collections.PathCollection):
A maplotlib scatter plot object
'''
idx = mask > 0
if value:
idx = idx & (y.value > 0)
if (value) and (snr is not None):
idx = idx & (y.snr > snr)
if sf_mask is not None:
idx = idx & sf_mask
c = mpl.colors.to_rgb(color)
c_a = np.array([c + (i, ) for i in mask[idx] / 15])
c_a[c_a > 1] = 1
if idx.sum() > 0:
if value:
return ax.scatter(x[idx], y.value[idx], c=c_a, edgecolor=c_a, **kwargs)
else:
return ax.scatter(x[idx], y[idx], c=c_a, edgecolor=c_a, **kwargs)
return None
[docs]class GZ3DVAC(VACMixIn):
'''Provides access to the Galaxy Zoo 3D spaxel masks.
VAC name: Galaxy Zoo: 3D
URL: https://www.sdss.org/dr17/data_access/value-added-catalogs/?vac_id=galaxy-zoo-3d
Description: Galaxy Zoo: 3D (GZ: 3D) made use of a project on the Zooniverse platform to
crowdsource spaxel masks locating galaxy centers, foreground stars, bars and spirals
in the SDSS images of MaNGA target galaxies. These masks (available for use within Marvin)
can be used to pick out spectra, or map quantities associated with the different
structures. See Masters et al. 2021 for more information, advice on useage and examples.
Authors: Coleman Krawczyk, Karen Masters and the rest of the Galaxy Zoo 3D Team.
'''
name = 'gz3d'
description = 'Return object for working with Galaxy Zoo 3D data masks'
version = {'DR17': 'v4_0_0', 'MPL-11': 'v4_0_0'}
display_name = 'Galaxy Zoo 3D'
url = 'https://www.sdss.org/dr17/data_access/value-added-catalogs/?vac_id=galaxy-zoo-3d'
include = (marvin.tools.cube.Cube, marvin.tools.maps.Maps)
[docs] def set_summary_file(self, release):
''' Sets the path to the GalaxyZoo3D summary file '''
self.path_params = {'ver': self.version[release]}
self.summary_file = self.get_path('mangagz3dmetadata', path_params=self.path_params)
self.center_summary_file = self.get_path('mangagz3dcenters', path_params=self.path_params)
self.stars_summary_file = self.get_path('mangagz3dstars', path_params=self.path_params)
[docs] def get_target(self, parent_object):
'''Find the GZ3D data based on the manga ID'''
mangaid = parent_object.mangaid
if parent_object.__class__ == marvin.tools.cube.Cube:
cube = parent_object
maps = parent_object.getMaps()
else:
cube = parent_object.getCube()
maps = parent_object
if not self.file_exists(self.summary_file):
self.summary_file = self.download_vac('mangagz3dmetadata', path_params=self.path_params)
self.center_summary_file = self.download_vac('mangagz3dcenters', path_params=self.path_params)
self.stars_summary_file = self.download_vac('mangagz3dstars', path_params=self.path_params)
summary_table = Table.read(self.summary_file, hdu=1)
# Table adds extra spaces to short strings, these need to be stripped off
gz3d_mangaids = np.array([mid.strip() for mid in summary_table['MANGAID']])
idx = gz3d_mangaids == mangaid
if idx.sum() > 0:
file_name = summary_table[idx]['file_name'][0].strip()
self.path_params.update(file_name=file_name)
self.gz3d_filename = self.get_path('mangagz3d', path_params=self.path_params)
if not self.file_exists(self.gz3d_filename):
self.gz3d_filename = self.download_vac('mangagz3d', path_params=self.path_params)
return GZ3DTarget(self.gz3d_filename, cube, maps)
log.info('There is no GZ3D data for this mangaid: {0}'.format(mangaid))
return None
[docs]class GZ3DTarget(object):
'''A customized class to open and display GZ3D spaxel masks
Parameters:
filename (str):
Path to the GZ3D fits file
cube (marvin.tools.cube.Cube):
Marvin Cube object
maps (marvin.tools.maps.Maps):
Mavin Maps object
Attributes:
hdulist (list):
List containing the 11 HDUs present in the GZ3D fits file (see <url> for full data model)
wcs (astropy.wcs):
WCS object for the GZ3D masks (e.g. HDU[1] to HDU[4])
image (numpy.array):
The galaxy image shown to GZ3D volunteers
center_mask (numpy.array):
Pixel mask (same shape as image) of the clustering results for the galaxy center(s). Each identified
center is represented by a 2 sigma ellipse of clustered points with the value of the pixels inside
the ellipse equal to the number of points belonging to that cluster.
star_mask (numpy.array):
Pixel mask (same shape as image) of the clustering results for forground star(s). Each identified
star is represented by a 2 sigma ellipse of clustered points with the value of the pixels inside
the ellipse equal to the number of points belonging to that cluster.
spiral_mask (numpy.array):
Pixel mask (same shape as image) of the spiral arm location(s). The value for the pixels is the number
of overlapping polygons at that location.
bar_mask (numpy.array):
Pixel mask (same shape as image) of the bar location. The value for the pixels is the number of
overlapping polygons at that location.
metadata (astropy.Table):
Table containing metadata about the galaxy.
ifu_size (int):
Size of IFU
center_clusters (astropy.Table):
Position for identified galaxy center(s) in both image pixels and (RA, DEC)
num_centers (int):
Number of galaxy centers identified
star_clusters (astropy.Table):
Position for identified forground star(s) in both image pixels and (RA, DEC)
num_stars (int):
Number of forground stars identified
center_star_classifications (astropy.Table):
Raw GZ3D classifications for center(s) and star(s)
num_center_star_classifications (int):
Total number of classifications made for either center(s) or star(s)
num_center_star_classifications_non_blank (int):
Total number of non-blank classifications made for either center(s) or star(s)
spiral_classifications (astropy.Table):
Raw GZ3D classifications for spiral arms
num_spiral_classifications (int):
Total number of spiral arm classifications made
num_spiral_classifications_non_blank (int):
Total number of non-blank spiral arm classifications made
bar_classifications (astropy.Table):
Raw GZ3D classifications for bars
num_bar_classifications (int):
Total number of bar classifications made
num_bar_classifications_non_blank (int):
Total number of non-blank bar classifications made
cube (marvin.tools.cube.Cube):
Marvin Cube object
maps (marvin.tools.maps.Maps)
Marvin Maps object
center_mask_spaxel (numpy.array):
The center_mask projected into spaxel space
star_mask_spaxel (numpy.array):
The star_mask projected into spaxel space
spiral_mask_spaxel (numpy.array):
The spiral_mask projected into spaxel space
bar_mask_spaxel (numpy.array):
The bar_mask projected into spaxel space
other_mask_spaxel (numpy.array):
A mask for spaxel not contained in any of the other spaxel masks
'''
def __init__(self, filename, cube, maps):
'''Set useful paramters and process the GZ3D masks'''
# get the subject id from the filename
self.subject_id = filename.split('/')[-1].split('_')[-1].split('.')[0]
self.cube = cube
self.maps = maps
self.mean_bar = None
self.mean_spiral = None
self.mean_center = None
self.mean_not_bar = None
self.mean_not_spiral = None
self.mean_not_center = None
self.log_oiii_hb = None
self.log_nii_ha = None
self.log_sii_ha = None
self.log_oi_ha = None
self.dis = None
# read in the fits file
with fits.open(filename) as hdulist:
self.hdulist = hdulist
# grab the wcs
self.wcs = WCS(self.hdulist[1].header)
self._process_images()
# read in metadata
self.metadata = Table(self.hdulist[5].data)
self.ifu_size = int(self.metadata['IFUDESIGNSIZE'][0])
self._process_clusters()
self._process_clusters_classifications()
self._process_spiral_classifications()
self._process_bar_classifications()
self._process_all_spaxel_masks()
self._get_bpt()
def _process_images(self):
'''Extract the data from the fits file and give it useful names'''
# read in images
self.image = self.hdulist[0].data
self.center_mask = self.hdulist[1].data
self.star_mask = self.hdulist[2].data
self.spiral_mask = self.hdulist[3].data
self.bar_mask = self.hdulist[4].data
def _process_clusters(self):
'''Format the cluster tables form the fits file as Astropy tables and count the number of rows'''
# read in center clusters
self.center_clusters = Table(self.hdulist[6].data)
self.num_centers = len(self.center_clusters)
# read in star clusters
self.star_clusters = Table(self.hdulist[7].data)
self.num_stars = len(self.star_clusters)
def _process_clusters_classifications(self):
'''Format the Zooniverse point classifications as Astropy tables and count the number of unique non-blank classifications'''
# read in center and star classifications
self.center_star_classifications = Table(self.hdulist[8].data)
self.num_center_star_classifications = len(self.center_star_classifications)
convert_json(self.center_star_classifications, 'center_points')
convert_json(self.center_star_classifications, 'star_points')
self.num_center_star_classifications_non_blank = non_blank(self.center_star_classifications, 'center_points_list', 'star_points_list')
def _process_spiral_classifications(self):
'''Format the Zooniverse spiral arm classifications as Astropy tables and count the number of unique non-blank classifications'''
# read in spiral classifications
self.spiral_classifications = Table(self.hdulist[9].data)
self.num_spiral_classifications = len(self.spiral_classifications)
convert_json(self.spiral_classifications, 'spiral_paths')
self.num_spiral_classifications_non_blank = non_blank(self.spiral_classifications, 'spiral_paths_list')
def _process_bar_classifications(self):
'''Format the Zooniverse bar classifications as Astropy tables and count the number of unique non-blank classifications'''
# read in bar classifications
self.bar_classifications = Table(self.hdulist[10].data)
self.num_bar_classifications = len(self.bar_classifications)
convert_json(self.bar_classifications, 'bar_paths')
self.num_bar_classifications_non_blank = non_blank(self.bar_classifications, 'bar_paths_list')
[docs] def center_in_pix(self):
'''Return the center of the IFU in image coordinates'''
return self.wcs.wcs_world2pix(np.array([[self.metadata['ra'][0], self.metadata['dec'][0]]]), 1)[0]
[docs] def get_hexagon(self, correct_hex=True, edgecolor='magenta'):
'''Get the IFU hexagon in image as a matplotlib polygon for plotting
Paramters:
correct_hex (bool, default=True):
If True it returns the correct IFU hexagon, if False it returns the hexagon shown
to the GZ3D volunteers (this was slightly too small due to a bug when producing the
original images for the project).
edgecolor (matplotlib color):
What color to make the hexagon.
Returns:
hexagon (matplotlib.patches.RegularPolygon):
A matplotlib patch object of the IFU hexagon returned in image coordinates.
'''
# the spacing should be ~0.5 arcsec not 0, and it should not be rotated by np.sqrt(3) / 2
if correct_hex:
# each hex has a total diameter of 2.5 arcsec on the sky (only 2 of it is a fiber)
diameter = 2.5 / 0.099
# the radius for mpl is from the center to each vertex, not center to side
r = LUT[self.ifu_size] * diameter / 2
else:
# this was me being wrong about the hexagon params forgetting about the space between fibers
diameter = 2.0 / 0.099
# The factor of 1.1 was to try (and fail) to account for the space between fibers :(
r = 1.1 * LUT[self.ifu_size] * diameter / 2
c = self.center_in_pix()
return patches.RegularPolygon(c, 6, r, fill=False, orientation=np.deg2rad(30), edgecolor=edgecolor, linewidth=0.8)
def _get_ellipse_list(self, table):
'''Convert table of x, y, var_x, var_y, var_x_y into matplotlib ellipse objects (one for each row)'''
ellip_list = []
for idx in range(len(table)):
pos = np.array([table['x'][idx], table['y'][idx]])
cov = np.array([[table['var_x'][idx], table['var_x_y'][idx]], [table['var_x_y'][idx], table['var_y'][idx]]])
ellip_list.append(cov_to_ellipse(cov, pos, nstd=2, edgecolor='k', facecolor='none', lw=1))
return ellip_list
[docs] def get_center_ellipse_list(self):
'''Return matplotlib ellipse objects for identified galaxy center(s)'''
return self._get_ellipse_list(self.center_clusters)
[docs] def get_star_ellipse_list(self):
'''Return matplotlib ellipse objects for identified forground star(s)'''
return self._get_ellipse_list(self.star_clusters)
def _get_spaxel_grid_xy(self, include_edges=False, grid_size=None):
'''Find the spaxel grid (in image coordinates) for the images in the fits file'''
if grid_size is None:
grid_size = self.cube.data['FLUX'].data.shape[1:]
one_grid = 0.5 / 0.099
c = self.center_in_pix()
grid_y = np.arange(grid_size[0] + include_edges) * one_grid
grid_x = np.arange(grid_size[1] + include_edges) * one_grid
grid_y = grid_y - np.median(grid_y) + c[0]
grid_x = grid_x - np.median(grid_x) + c[1]
return grid_x, grid_y
[docs] def get_spaxel_grid(self, grid_size=None):
'''Return the data needed to plot the spaxel grid over the GZ3D image'''
grid_x, grid_y = self._get_spaxel_grid_xy(include_edges=True, grid_size=grid_size)
v_line_x = np.vstack([grid_x, grid_x])
v_line_y = np.array([[grid_y[0]], [grid_y[-1]]])
h_line_x = np.array([[grid_x[0]], [grid_x[-1]]])
h_line_y = np.vstack([grid_y, grid_y])
return [(v_line_x, v_line_y), (h_line_x, h_line_y)]
def _get_spaxel_mask(self, mask, grid_size=None):
'''Resample GZ3D masks onto the spaxel grid using a Bivariate spline resampling'''
# assumes a 0.5 arcsec grid centered on the ifu's ra and dec
# use a Bivariate spline approximation to resample mask to the spaxel grid
xx = np.arange(mask.shape[1])
yy = np.arange(mask.shape[0])
s = RectBivariateSpline(xx, yy, mask)
grid_x, grid_y = self._get_spaxel_grid_xy(grid_size=grid_size)
# flip the output mask so the origin is the lower left of the image
s_mask = np.flipud(s(grid_x, grid_y))
# zero out small values
s_mask[s_mask < 0.5] = 0
return s_mask
def _process_all_spaxel_masks(self, grid_size=None):
'''Process all GZ3D masks onto the MaNGA spaxel grid and give them useful names'''
self.center_mask_spaxel = self._get_spaxel_mask(self.center_mask, grid_size=grid_size)
self.star_mask_spaxel = self._get_spaxel_mask(self.star_mask, grid_size=grid_size)
self.spiral_mask_spaxel = self._get_spaxel_mask(self.spiral_mask, grid_size=grid_size)
self.bar_mask_spaxel = self._get_spaxel_mask(self.bar_mask, grid_size=grid_size)
self.other_mask_spaxel = (self.spiral_mask_spaxel == 0) & (self.bar_mask_spaxel == 0) & (self.center_mask_spaxel == 0)
def _stack_spectra(self, mask_name, inv=False):
'''Stack multiple spectra withing a spaxel mask following Westfall et al. 2019 for covariance correction factors'''
mask = getattr(self, mask_name)
if inv:
mask = mask.max() - mask
mdx = np.where(mask > 0)
if len(mdx[0] > 0):
weights = mask[mdx]
spaxel_index = np.array(mdx).T
spectra = [s.flux for s in self.cube[mdx]]
# only keep spectra inside the IFU
in_ifu = np.array([not any(2**0 & s.mask) for s in spectra])
if in_ifu.sum() == 0:
return None
spectra = [spectra[i] for i in in_ifu.nonzero()[0]]
spaxel_index = spaxel_index[in_ifu]
weights = weights[in_ifu]
weights_total = weights.sum()
if len(spectra) == 1:
return spectra[0]
flux = np.array([sp.value for sp in spectra])
# the weighted mean
mean = (flux * weights[:, None]).sum(axis=0) / weights_total
# we need to handle covariance between spaxels when calculating
# uncertainties. We follow Westfall et al. 2019's method based in
# distance between spaxels
d = distance_matrix(spaxel_index, spaxel_index) / 1.92
roh = np.exp(-0.5 * d**2)
# only work with value where roh > 0.003
xx, yy = np.where(roh > 0.003)
sigma = np.array([sp.error.value for sp in spectra])
running_sum = np.zeros_like(sigma[0])
for idx, jdx in zip(xx, yy):
running_sum += roh[idx, jdx] * weights[idx] * weights[jdx] * sigma[idx] * sigma[jdx]
ivar = (weights_total**2) / running_sum
return Spectrum(
mean,
unit=spectra[0].unit,
wavelength=spectra[0].wavelength,
wavelength_unit=spectra[0].wavelength.unit,
pixmask_flag=spectra[0].pixmask_flag,
ivar=ivar
)
return None
[docs] def get_mean_spectra(self, inv=False):
'''Calculate the mean spectra inside each of the spaxel masks accounting
for covariance following Westfall et al. 2019's method based in distance
between spaxels.
Parameters:
inv (bool, default=False):
If true this function will also calculate the mean spectra
for each inverted mask. Useful if you want to make difference
spectra (e.g. <spiral> - <not spiral>).
Attributes:
mean_bar (marvin.tools.quantities.spectrum):
average spectra inside the bar mask
mean_spiral (marvin.tools.quantities.spectrum):
average spectra inside the spiral mask
mean_center (marvin.tools.quantities.spectrum):
average spectra inside the center mask
mean_not_bar (marvin.tools.quantities.spectrum):
average spectra outside the bar mask
mean_not_spiral (marvin.tools.quantities.spectrum):
average spectra outside the spiral mask
mean_not_center (marvin.tools.quantities.spectrum):
average spectra outside the center mask
'''
if self.mean_bar is None:
self.mean_bar = self._stack_spectra('bar_mask_spaxel')
if self.mean_spiral is None:
self.mean_spiral = self._stack_spectra('spiral_mask_spaxel')
if self.mean_center is None:
self.mean_center = self._stack_spectra('center_mask_spaxel')
if inv:
if self.mean_not_bar is None:
self.mean_not_bar = self._stack_spectra('bar_mask_spaxel', inv=True)
if self.mean_not_spiral is None:
self.mean_not_spiral = self._stack_spectra('spiral_mask_spaxel', inv=True)
if self.mean_not_center is None:
self.mean_not_center = self._stack_spectra('center_mask_spaxel', inv=True)
def _get_bpt(self, snr_min=3, oi_sf=False):
'''Use the `bpt` module to grab the information needed to make BPT plots color coded by the GZ3D masks'''
# Gets the necessary emission line maps
oiii = bpt.get_masked(self.maps, 'oiii_5008', snr=bpt.get_snr(snr_min, 'oiii'))
nii = bpt.get_masked(self.maps, 'nii_6585', snr=bpt.get_snr(snr_min, 'nii'))
ha = bpt.get_masked(self.maps, 'ha_6564', snr=bpt.get_snr(snr_min, 'ha'))
hb = bpt.get_masked(self.maps, 'hb_4862', snr=bpt.get_snr(snr_min, 'hb'))
sii = bpt.get_masked(self.maps, 'sii_6718', snr=bpt.get_snr(snr_min, 'sii'))
oi = bpt.get_masked(self.maps, 'oi_6302', snr=bpt.get_snr(snr_min, 'oi'))
self.log_oiii_hb = np.ma.log10(oiii / hb)
self.log_nii_ha = np.ma.log10(nii / ha)
self.log_sii_ha = np.ma.log10(sii / ha)
self.log_oi_ha = np.ma.log10(oi / ha)
sf_mask_nii = ((self.log_oiii_hb < bpt.kewley_sf_nii(self.log_nii_ha)) & (self.log_nii_ha < 0.05)).filled(False)
sf_mask_sii = ((self.log_oiii_hb < bpt.kewley_sf_sii(self.log_sii_ha)) & (self.log_sii_ha < 0.32)).filled(False)
sf_mask_oi = ((self.log_oiii_hb < bpt.kewley_sf_oi(self.log_oi_ha)) & (self.log_oi_ha < -0.59)).filled(False)
if oi_sf:
self.sf_mask = sf_mask_nii & sf_mask_sii & sf_mask_oi
else:
self.sf_mask = sf_mask_nii & sf_mask_sii
[docs] def get_distance(self):
'''Find the radial distance between each spaxel and the center of the galaxy'''
if self.dis is None:
cdx = np.unravel_index(self.center_mask_spaxel.argmax(), self.center_mask_spaxel.shape)
self.dis = np.zeros_like(self.center_mask_spaxel)
for yy in range(self.dis.shape[0]):
for xx in range(self.dis.shape[1]):
self.dis[yy, xx] = np.linalg.norm([yy - cdx[0], xx - cdx[1]])
def _set_up_axes(self, ax, color_grid=None):
'''Helper function to set RA and DEC ticks on plots'''
ra = ax.coords['ra']
dec = ax.coords['dec']
# add axis labels
ra.set_axislabel('RA')
dec.set_axislabel('Dec')
ra.set_major_formatter('d.ddd')
ra.ticklabels.set_rotation(90)
ra.ticklabels.set_rotation_mode('anchor')
ra.ticklabels.set_pad(15)
dec.set_major_formatter('d.ddd')
ra.display_minor_ticks(True)
dec.display_minor_ticks(True)
# add a coordinate grid to the image
if color_grid is not None:
ax.coords.grid(color=color_grid, alpha=0.5, linestyle='solid', lw=1.5)
[docs] def plot_image(self, ax=None, color_grid=None, correct_hex=True, hex_color='C7'):
'''Plot original GZ3D image that was shown to volunteers.
Keywords:
ax (matplotlib.axes.Axes):
Matplotlib axis object. This axis must have a WCS projection set e.g.
`ax = fig.add_subplot(111, projection=data.wcs)`. If not provided a new
figure and axis will be created with the correct projection.
color_grid (string):
A matplotlib color to use for the RA-DEC grid lines. Default `None`.
correct_hex (bool):
If set to true the correct MaNGA hexagon will be plotted on top of the
galaxy cutout (the hexagon in the image shown to the volunteers was slightly
too small due to a bug when producing the original images for the project).
hex_color (string):
A matplotlib color to use for the correct MaNGA hexagon if `correct_hex` is
True. Default is `'C7'`.
Returns:
ax (matplotlib.axes.Axes):
The matplotlib axis object for the resulting plot.
'''
if (ax is None):
fig = plt.figure()
ax = fig.add_subplot(111, projection=self.wcs)
try:
self._set_up_axes(ax, color_grid=color_grid)
except AttributeError:
raise TypeError('ax must have a WCS project, e.g. `ax = fig.add_subplot(111, projection=data.wcs)`')
if correct_hex:
ax.add_patch(self.get_hexagon(correct_hex=True, edgecolor=hex_color))
ax.imshow(self.image)
return ax
[docs] def plot_masks(
self,
colors=['C1', 'C0', 'C4', 'C2'],
color_grid=None,
hex=True,
hex_color='C7',
show_image=False,
subplot_spec=None,
spaxel_masks=False
):
'''Plot GZ3D masks
Keywords:
colors (list):
A list of matplotlib colors to use for each of the masks. The order of the list is:
[Bar, Spiral, Forground Stars, Galaxy Center(s)]. Default value is
`['C1', 'C0', 'C4', 'C2']`.
color_grid (string):
A matplotlib color to use for the RA-DEC grid lines. Default `None`.
hex (bool):
If `True` plot the MaNGA hexagon. Default value is `True`.
hex_color (string):
A matplotlib color to use for the correct MaNGA hexagon if `correct_hex` is
True. Default is `'C7'`.
show_image (bool):
If `True` plot the original galaxy image behind the masks. Default is `False`.
subplot_spec (matplotlib.gridspec.SubplotSpec):
A gridspec subplot specification for this plot. If `None` is provided a new
figure will be created.
spaxel_masks (bool):
If `True` use the masks projected on to the MaNGA spaxel grid, other wise
plot them on the pixel grid of the GZ3D image shown to the volunteers. Default
value is `False`.
'''
if subplot_spec is None:
fig = plt.figure()
# image axis
gs = gridspec.GridSpec(1, 2, width_ratios=[0.9, 0.1], wspace=0.01)
else:
gs = gridspec.GridSpecFromSubplotSpec(1, 2, width_ratios=[0.9, 0.1], wspace=0.01, subplot_spec=subplot_spec)
# color bar axis
gs_color_bars = gridspec.GridSpecFromSubplotSpec(1, 4, wspace=0, subplot_spec=gs[1])
# alpha overlay all masks with correct colors
if spaxel_masks:
mask_list = [self.bar_mask_spaxel, self.spiral_mask_spaxel, self.star_mask_spaxel, self.center_mask_spaxel]
ax1 = plt.subplot(gs[0], projection=self.maps.wcs)
else:
mask_list = [self.bar_mask, self.spiral_mask, self.star_mask, self.center_mask]
ax1 = plt.subplot(gs[0], projection=self.wcs)
if (show_image) and (not spaxel_masks):
all_masks = alpha_maps(mask_list, colors, background_image=self.image / 255)
else:
all_masks = alpha_maps(mask_list, colors)
self._set_up_axes(ax1, color_grid=color_grid)
ax1.imshow(all_masks)
# overlay IFU hexagon
if hex:
ax1.add_patch(self.get_hexagon(correct_hex=True, edgecolor=hex_color))
# plot center and star ellipses to better define ellipse shape
if (not spaxel_masks):
center_ellipse = self.get_center_ellipse_list()
for e, count in zip(center_ellipse, self.center_clusters['count']):
e.set_edgecolor(make_alpha_color(count, colors[3]))
ax1.add_artist(e)
star_ellipse = self.get_star_ellipse_list()
for e, count in zip(star_ellipse, self.star_clusters['count']):
e.set_edgecolor(make_alpha_color(count, colors[2]))
ax1.add_artist(e)
# make the legend
bar_patch = mpl.patches.Patch(color=colors[0], label='Bar')
spiral_patch = mpl.patches.Patch(color=colors[1], label='Spiral')
star_patch = mpl.patches.Patch(color=colors[2], label='Star')
center_patch = mpl.patches.Patch(color=colors[3], label='Center')
plt.legend(
handles=[bar_patch, spiral_patch, star_patch, center_patch],
ncol=2,
loc='lower center',
mode='expand'
)
# make the colorbars
ax_bar, cb_bar = plot_alpha_bar(colors[0], gs_color_bars[0])
ax_spiral, cb_spiral = plot_alpha_bar(colors[1], gs_color_bars[1])
ax_star, cb_star = plot_alpha_bar(colors[2], gs_color_bars[2])
ax_center, cb_center = plot_alpha_bar(colors[3], gs_color_bars[3])
ax_center.tick_params(axis=u'both', which=u'both', length=0)
tick_labels = np.arange(0, 16)
tick_locs = tick_labels - 0.5
cb_center.set_ticks(tick_locs)
cb_center.set_ticklabels(tick_labels)
cb_center.set_label('Count')
return gs
def _plot_bpt_boundary(self, ax, bpt_kind):
'''Plot the BPT boundary lines'''
if bpt_kind == 'log_nii_ha':
xx_sf_nii = np.linspace(-1.281, 0.045, int(1e4))
xx_comp_nii = np.linspace(-2, 0.4, int(1e4))
ax.plot(xx_sf_nii, bpt.kewley_sf_nii(xx_sf_nii), 'k--', zorder=90)
ax.plot(xx_comp_nii, bpt.kewley_comp_nii(xx_comp_nii), 'k-', zorder=90)
ax.set_xlim(-2, 0.5)
ax.set_ylim(-1.5, 1.3)
ax.set_xlabel(r'log([NII]/H$\alpha$)')
ax.set_ylabel(r'log([OIII]/H$\beta$)')
elif bpt_kind == 'log_sii_ha':
xx_sf_sii = np.linspace(-2, 0.315, int(1e4))
xx_agn_sii = np.array([-0.308, 1.0])
ax.plot(xx_sf_sii, bpt.kewley_sf_sii(xx_sf_sii), 'k-', zorder=90)
ax.plot(xx_agn_sii, bpt.kewley_agn_sii(xx_agn_sii), 'k-', zorder=90)
ax.set_xlim(-1.5, 0.5)
ax.set_ylim(-1.5, 1.3)
ax.set_xlabel(r'log([SII]/H$\alpha$)')
ax.set_ylabel(r'log([OIII]/H$\beta$)')
elif bpt_kind == 'log_oi_ha':
xx_sf_oi = np.linspace(-2.5, -0.7, int(1e4))
xx_agn_oi = np.array([-1.12, 0.5])
ax.plot(xx_sf_oi, bpt.kewley_sf_oi(xx_sf_oi), 'k-', zorder=90)
ax.plot(xx_agn_oi, bpt.kewley_agn_oi(xx_agn_oi), 'k-', zorder=90)
ax.set_xlim(-2.5, 0.0)
ax.set_ylim(-1.5, 1.3)
ax.set_xlabel(r'log([OI]/H$\alpha$)')
ax.set_ylabel(r'log([OIII]/H$\beta$)')
else:
raise AttributeError('bpt_kind must be one of "log_nii_ha", "log_sii_ha", or "log_oi_ha", {0} was given'.format(bpt_kind))
[docs] def plot_bpt(
self,
ax=None,
colors=['C1', 'C0', 'C4', 'C2'],
bpt_kind='log_nii_ha',
**kwargs
):
'''Plot a BPT diagram for a galaxy that colors the data points based on the GZ3D masks
Keywords:
ax (matplotlib.axes.Axes):
The matplotlib axis object to use for the plot. If `None` is provided a new
figure and axis is created for the plot.
colors (list):
A list of matplotlib colors to use for each of the masks. The order of the list is:
[Bar, Spiral, Forground Stars, Galaxy Center(s)]. Default value is
`['C1', 'C0', 'C4', 'C2']`.
bpt_kind (string):
The kind of BPT plot to make. This can be one of three values `'log_nii_ha'` (default),
`'log_sii_ha'`, or `'log_oi_ha'`.
kwargs:
All other keywords are pass forward to matplotlib's scatter plot function.
Returns:
ax (matplotlib.axes.Axes):
The matplotlib axis object for the resulting plot.
'''
if bpt_kind not in ["log_nii_ha", "log_sii_ha", "log_oi_ha"]:
raise AttributeError('bpt_kind must be one of "log_nii_ha", "log_sii_ha", or "log_oi_ha", {0} was given'.format(bpt_kind))
y = self.log_oiii_hb
x = getattr(self, bpt_kind)
mdx = ~(y.mask | x.mask)
if ax is None:
fig = plt.figure()
ax = fig.add_subplot(111)
s = kwargs.pop('s', 8)
odx = mdx & (self.other_mask_spaxel > 0)
ax.scatter(x[odx], y[odx], c='#c5c5c5', edgecolor='#c5c5c5', s=s, label='Other', **kwargs)
plot_alpha_scatter(x, y, self.spiral_mask_spaxel, colors[1], ax, s=s, sf_mask=mdx, snr=None, value=False, label='Spiral', **kwargs)
plot_alpha_scatter(x, y, self.bar_mask_spaxel, colors[0], ax, s=s, sf_mask=mdx, snr=None, value=False, label='Bar', **kwargs)
plot_alpha_scatter(x, y, self.star_mask_spaxel, colors[2], ax, s=s, sf_mask=mdx, snr=None, value=False, label='Star', **kwargs)
plot_alpha_scatter(x, y, self.center_mask_spaxel, colors[3], ax, s=s, sf_mask=mdx, snr=None, value=False, label='Center', **kwargs)
self._plot_bpt_boundary(ax, bpt_kind)
return ax
[docs] def polar_plot(
self,
x_unit='theta',
ax=None,
colors=['C1', 'C0', 'C4', 'C2'],
key='specindex_dn4000',
ylabel=r'D_{n}4000',
snr=3,
sf_only=False,
**kwargs
):
'''Make a plot of a MaNGA Map value vs. R or theta with the points color coded by
what GZ3D mask they belong to.
x_unit (string):
What x-value to plot against. Either `'theta'` (default) or `'radius'`.
ax (matplotlib.axes.Axes):
The matplotlib axis object to use for the plot. If `None` is provided a new
figure and axis is created for the plot.
colors (list):
A list of matplotlib colors to use for each of the masks. The order of the list is:
[Bar, Spiral, Forground Stars, Galaxy Center(s)]. Default value is
`['C1', 'C0', 'C4', 'C2']`.
key (string):
Name of the MaNGA Map attribute to plot. The default value is `'specindex_dn4000'`.
ylabel (string):
The `ylabel` to use for the plot (units will automatically be added to the label
based on the map being used).
snr (float):
The minimum signal to noise cutoff to use for the plot. The default value is `3`.
sf_only (bool):
If `True` only plot spaxes that are star forming. The default value is `False`.
kwargs:
All other keywords are pass forward to matplotlib's scatter plot function.
Returns:
ax (matplotlib.axes.Axes):
The matplotlib axis object for the resulting plot.
'''
title = []
s = kwargs.pop('s', 8)
if ax is None:
fig = plt.figure()
ax = fig.add_subplot(111)
if x_unit.lower() == 'theta':
x = self.maps.spx_ellcoo_elliptical_azimuth.value
ax.set_xticks([0, 90, 180, 270, 360])
ax.set_xlabel(r'$\theta$')
elif x_unit.lower() == 'radius':
r = self.maps.spx_ellcoo_elliptical_radius.value
r_50 = self.maps.nsa['elpetro_th50_r']
x = r / r_50
ax.set_xlabel(r'R / R$_{50}$')
else:
raise AttributeError('x_unit must be either `theta` or `radius`, {0} was given'.format(x_unit))
line = self.maps[key]
# other spaxel masks
odx = (self.other_mask_spaxel > 0) & (line.value > 0)
if snr is not None:
title.append('S/N > {0}'.format(snr))
odx = odx & (line.snr > snr)
sf_mask = None
if sf_only:
# star forming only
sf_mask = self.sf_mask
title.append('star forming only')
odx = odx & sf_mask
# plot scatter points
ax.scatter(x[odx], line.value[odx], c='#c5c5c5', edgecolor='#c5c5c5', s=s, **kwargs, label='Other')
plot_alpha_scatter(x, line, self.spiral_mask_spaxel, colors[1], ax, s=s, snr=snr, sf_mask=sf_mask, label='Spiral')
plot_alpha_scatter(x, line, self.bar_mask_spaxel, colors[0], ax, s=s, snr=snr, sf_mask=sf_mask, label='Bar')
plot_alpha_scatter(x, line, self.star_mask_spaxel, colors[2], ax, s=s, snr=snr, sf_mask=sf_mask, label='Star')
plot_alpha_scatter(x, line, self.center_mask_spaxel, colors[3], ax, s=s, snr=snr, sf_mask=sf_mask, label='Center')
if len(title) > 0:
ax.set_title(','.join(title))
if line.unit.to_string() != '':
ax.set_ylabel('$\\mathrm{{{}}}\\,[${}$]$'.format(ylabel, line.unit.to_string('latex')))
else:
ax.set_ylabel('$\\mathrm{{{}}}$'.format(ylabel))
return ax
def __str__(self):
'''A useful summary of the data in the GZ3D fits file'''
return '\n'.join([
'Subject info:',
' subject id: {0}'.format(self.subject_id),
' manga id: {0}'.format(self.metadata['MANGAID'][0]),
' ra: {0}'.format(self.metadata['ra'][0]),
' dec: {0}'.format(self.metadata['dec'][0]),
' ifu size: {0}'.format(self.ifu_size),
'Classification counts:',
' {0} center/star, {1} non_blank'.format(self.num_center_star_classifications, self.num_center_star_classifications_non_blank),
' {0} spiral, {1} non_blank'.format(self.num_spiral_classifications, self.num_spiral_classifications_non_blank),
' {0} bar, {1} non_blank'.format(self.num_bar_classifications, self.num_bar_classifications_non_blank),
'Cluster counts:',
' {0} center(s)'.format(self.num_centers),
' {0} star(s)'.format(self.num_stars)
])