<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Spatially-Resolved Mass-Metallicity Relation
We're going to construct the spatially-resolved mass-metallicity relation (MZR) for a MaNGA galaxy, where mass refers to stellar mass and metallicity refers to gas-phase oxygen abundance.
Roadmap
Compute metallicity.
Select spaxels that are
star-forming, 
not flagged as "bad data," and
above a signal-to-noise ratio threshold.
Compute stellar mass surface density.
Plot metallicity as a function of stellar mass surface density.
Step1:  Load Maps for Galaxy
Import the Marvin Maps class from marvin.tools.maps and initialize a Maps object for the galaxy 8077-6104.
Step2:  Measure Metallicity
Pettini & Pagel (2004) N2 metallicity calibration
We are going to use the N2 metallicity calibration (their Equation 1) from Pettini & Pagel (2004)
Step3:  Calculate the necessary line ratio.
Marvin can do map arithmetic, which propagates the inverse variances and masks, so you can just do +, -, *, /, and ** operations as normal.  (Note
Step4:  Finally, calculate the metallicity.
Step5:  Select Spaxels
Using the BPT Diagram to select star-forming spaxels
Metallicity indicators only work for star-forming spaxels, so we need a way to select only these spaxels.
The classic diagnostic diagram for classify the emission from galaxies (or galactic sub-regions) as star-forming or non-star-forming (i.e., from active galactic nuclei (AGN) or evolved stars) was originally proposed in Baldwin, Phillips, & Terlevich (1981) and is known as the BPT diagram.
The BPT diagram uses ratios of emission lines to separate thermal and non-thermal emission.
The classic BPT diagram uses [OIII]5007 / Hbeta vs. [NII]6583 / Halpha, but there are several versions of the BPT diagram that use different lines ratios.
BPT Diagrams with Marvin
Let's use Marvin's maps.get_bpt() method to make BPT diagrams for this galaxy.
red line
Step6:  The BPT masks are dictionaries of dictionaries of a boolean (True/False) arrays.  We are interested in the spaxels that are classified as star-forming in all three BPT diagrams are designated as True, which is designated with the global key.  Print this mask.
Step7:  Masks
MaNGA (and SDSS generally) use bitmasks to communicate data quality.
Marvin has built-in methods to convert from the bitmasks integer values to individual bits or labels and to create new masks by specifying a set of labels.
Show the mask schema with n2.pixmask.schema.
Step8:  Select non-star-forming spaxels (from the BPT mask) and set their mask value to the DAP's DONOTUSE value with the n2.pixmask.labels_to_value() method.  Note that we are selecting spaxels that we want from the BPT mask (i.e., True is a spaxel to keep), whereas we are using the pixmask to select spaxels that we want to exclude (i.e., True is a spaxel to ignore).
Step9:  Select spaxels classified by the DAP as bad data according to the masks for spaxels with no IFU coverage, with unreliable measurements, or otherwise unfit for science.  Use the n2.pixmask.get_mask method.
Step10:  Select spaxels with signal-to-noise ratios (SNRs) > 3 on both [NII] 6585 and Halpha.
ha.ivar = inverse variance = $\frac{1}{\sigma^2}$, where $\sigma$ is the error.
Step11:  Do a bitwise (binary) OR to create a master mask of spaxels to ignore.
Step12:  Plot the Metallicity Map
Plot the map of metallicity using the plot() method from your Marvin Map metallicity object.  Also, mask undesirable spaxels and label the colorbar.
Note
Step13:  Compute Stellar Mass Surface Density
Read in spaxel stellar mass measurements from the Firefly spectral fitting catalog (Goddard et al. 2017).
The Firefly stellar population fitting results file is large (1.8 GB), so we have extracted only the measurements needed for our purposes.
Summary of the Firefly Value Added Catalog
Datamodel of the Firefly Value Added Catalog
Convert spaxel angular size to a physical scale in pc.
Divide stellar mass by area to get stellar surface mass density.
Read in stellar masses
Use pandas to read in the csv file with stellar masses.
Step14:  Plot stellar mass map using ax.imshow().  MaNGA maps are oriented such that you want to specify origin='lower'.  Also include a labelled colorbar.
Step15:  Calculate physical size of a spaxel
MaNGA's maps (and data cubes) have a spaxel size of 0.5 arcsec. Let's convert that into a physical scale for our galaxy.
Step16:  Get the redshift of the galaxy from the maps.nsa attribute.
Step17:  We'll use the small angle approximation to estimate the physical scale
Step18:  Rearrange the small angle formula to solve for the scale ($\frac{d}{\theta}$) in pc / arcsec.
Step19:  Now convert the spaxel size from arcsec to parsecs and calculate the area of a spaxel.
Step20:  Finally, we simply divide the stellar mass by the area to get the stellar mass surface density $\Sigma_\star$ in units of $\frac{M_\odot}{pc^2}$.
Step21:  Let's plot metallicity as a function of $\Sigma_\star$!  Remember to apply the mask.  Also set the axis range to be [0, 4, 8, 8.8].
Step22:  MaNGA Spatially-Resolved Mass-Metallicity Relation
We have constructed the spatially-resolved MZR for one galaxy, but we are interested in understanding the evolution of galaxies in general, so we want to repeat this exercise for many galaxies.  In Barrera-Ballesteros et al. (2016), Jorge Barrera-Ballesteros (who gave a talk at Pitt in November 2017) did just this, and here is the analogous figure for 653 disk galaxies.
<img src="images/barrera-ballesteros_local_mzr.png" style="width
Step23:  Remake the spatially-resolved MZR plot for our galaxy showing the he best fit line from Barrera-Ballesteros et al. (2016).
<ASSISTANT_TASK:>
Python Code:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
import os
from os.path import join
path_notebooks = os.path.abspath('.')
path_data = join(path_notebooks, 'data')
Explanation: Spatially-Resolved Mass-Metallicity Relation
We're going to construct the spatially-resolved mass-metallicity relation (MZR) for a MaNGA galaxy, where mass refers to stellar mass and metallicity refers to gas-phase oxygen abundance.
Roadmap
Compute metallicity.
Select spaxels that are
star-forming, 
not flagged as "bad data," and
above a signal-to-noise ratio threshold.
Compute stellar mass surface density.
Plot metallicity as a function of stellar mass surface density.
End of explanation
from marvin.tools.maps import Maps
# REMOVE FROM NOTEBOOK
filename = '/Users/andrews/hacks/galaxies-mzr/data/manga-8077-6104-MAPS-SPX-GAU-MILESHC.fits.gz'
maps = Maps(filename=filename)
# maps = Maps('8077-6104')
Explanation: Load Maps for Galaxy
Import the Marvin Maps class from marvin.tools.maps and initialize a Maps object for the galaxy 8077-6104.
End of explanation
nii = maps.emline_gflux_nii_6585
ha = maps.emline_gflux_ha_6564
Explanation: Measure Metallicity
Pettini & Pagel (2004) N2 metallicity calibration
We are going to use the N2 metallicity calibration (their Equation 1) from Pettini & Pagel (2004):
12 + log(O/H) = 8.90 + 0.57 $\times$ log( $\frac{F([NII])}{F(H\alpha)}$ ).
One of the benefits of this calibration is that the required lines are very close in wavelength, so the reddening correction is negligible.
Get [NII] 6585 and Halpha flux maps from the Marvin Maps object.  Note: MaNGA (and Marvin) use the wavelengths of lines in vaccuum, whereas they are usually reported in air, hence the slight offsets.
End of explanation
n2 = nii / ha
logn2 = np.log10(n2)
Explanation: Calculate the necessary line ratio.
Marvin can do map arithmetic, which propagates the inverse variances and masks, so you can just do +, -, *, /, and ** operations as normal.  (Note: taking the log of a Marvin Map will work for the values but the inverse variance propagation does not correctly propagate the inverse variance yet.)
End of explanation
oh = 8.90 + 0.57 * logn2
Explanation: Finally, calculate the metallicity.
End of explanation
masks_bpt, __, __ = maps.get_bpt()
Explanation: Select Spaxels
Using the BPT Diagram to select star-forming spaxels
Metallicity indicators only work for star-forming spaxels, so we need a way to select only these spaxels.
The classic diagnostic diagram for classify the emission from galaxies (or galactic sub-regions) as star-forming or non-star-forming (i.e., from active galactic nuclei (AGN) or evolved stars) was originally proposed in Baldwin, Phillips, & Terlevich (1981) and is known as the BPT diagram.
The BPT diagram uses ratios of emission lines to separate thermal and non-thermal emission.
The classic BPT diagram uses [OIII]5007 / Hbeta vs. [NII]6583 / Halpha, but there are several versions of the BPT diagram that use different lines ratios.
BPT Diagrams with Marvin
Let's use Marvin's maps.get_bpt() method to make BPT diagrams for this galaxy.
red line: maximal starbust (Kewley et al 2001) -- everything to the right is non-star-forming.
dashed black line: conservative star-forming cut (Kauffmann et al. 2003) -- everything to the left is star-forming.
Line ratios that fall in between these two lines are designated "Composite" with contributions from both star-forming and non-star-forming emission.
blue line: separates non-star-forming spaxels into Seyferts and LINERs.
Seyferts are a type of AGNs.
LINERs (Low Ionization Nuclear Emission Regions) are not always nuclear (LIER is a better acronym) and not always AGN (oftern hot evolved stars).
Sometimes these diagnostic diagrams disagree with each other, hence the "Ambiguous" designation.
Try using maps.get_bpt? to read the documentation on how to use this function.
End of explanation
masks_bpt['sf']['global']
Explanation: The BPT masks are dictionaries of dictionaries of a boolean (True/False) arrays.  We are interested in the spaxels that are classified as star-forming in all three BPT diagrams are designated as True, which is designated with the global key.  Print this mask.
End of explanation
n2.pixmask.schema
Explanation: Masks
MaNGA (and SDSS generally) use bitmasks to communicate data quality.
Marvin has built-in methods to convert from the bitmasks integer values to individual bits or labels and to create new masks by specifying a set of labels.
Show the mask schema with n2.pixmask.schema.
End of explanation
mask_non_sf = ~masks_bpt['sf']['global'] * n2.pixmask.labels_to_value('DONOTUSE')
Explanation: Select non-star-forming spaxels (from the BPT mask) and set their mask value to the DAP's DONOTUSE value with the n2.pixmask.labels_to_value() method.  Note that we are selecting spaxels that we want from the BPT mask (i.e., True is a spaxel to keep), whereas we are using the pixmask to select spaxels that we want to exclude (i.e., True is a spaxel to ignore).
End of explanation
mask_bad_data = n2.pixmask.get_mask(['NOCOV', 'UNRELIABLE', 'DONOTUSE'])
Explanation: Select spaxels classified by the DAP as bad data according to the masks for spaxels with no IFU coverage, with unreliable measurements, or otherwise unfit for science.  Use the n2.pixmask.get_mask method.
End of explanation
min_snr = 3.
mask_nii_low_snr = (np.abs(nii.value * np.sqrt(nii.ivar)) < min_snr)
mask_ha_low_snr = (np.abs(ha.value * np.sqrt(ha.ivar)) < min_snr)
Explanation: Select spaxels with signal-to-noise ratios (SNRs) > 3 on both [NII] 6585 and Halpha.
ha.ivar = inverse variance = $\frac{1}{\sigma^2}$, where $\sigma$ is the error.
End of explanation
mask = mask_non_sf | mask_bad_data | mask_nii_low_snr | mask_ha_low_snr
Explanation: Do a bitwise (binary) OR to create a master mask of spaxels to ignore.
End of explanation
fig, ax = oh.plot(mask=mask, cblabel='12+log(O/H)')
Explanation: Plot the Metallicity Map
Plot the map of metallicity using the plot() method from your Marvin Map metallicity object.  Also, mask undesirable spaxels and label the colorbar.
Note: solar metallicity is about 8.7.
End of explanation
import pandas as pd
mstar = pd.read_csv(join(path_data, 'manga-{}_mstar.csv'.format(maps.plateifu)))
Explanation: Compute Stellar Mass Surface Density
Read in spaxel stellar mass measurements from the Firefly spectral fitting catalog (Goddard et al. 2017).
The Firefly stellar population fitting results file is large (1.8 GB), so we have extracted only the measurements needed for our purposes.
Summary of the Firefly Value Added Catalog
Datamodel of the Firefly Value Added Catalog
Convert spaxel angular size to a physical scale in pc.
Divide stellar mass by area to get stellar surface mass density.
Read in stellar masses
Use pandas to read in the csv file with stellar masses.
End of explanation
fig, ax = plt.subplots()
p = ax.imshow(mstar, origin='lower')
ax.set_xlabel('spaxel')
ax.set_ylabel('spaxel')
cb = fig.colorbar(p)
cb.set_label('log(Mstar) [M$_\odot$]')
Explanation: Plot stellar mass map using ax.imshow().  MaNGA maps are oriented such that you want to specify origin='lower'.  Also include a labelled colorbar.
End of explanation
spaxel_size = 0.5  # [arcsec]
# or programmatically:
# spaxel_size = float(maps.getCube().header['CD2_2']) * 3600
Explanation: Calculate physical size of a spaxel
MaNGA's maps (and data cubes) have a spaxel size of 0.5 arcsec. Let's convert that into a physical scale for our galaxy.
End of explanation
redshift = maps.nsa['z']
Explanation: Get the redshift of the galaxy from the maps.nsa attribute.
End of explanation
c = 299792  # speed of light [km/s]
H0 = 70  # [km s^-1 Mpc^-1]
D = c * redshift / H0  # approx. distance to galaxy [Mpc]
Explanation: We'll use the small angle approximation to estimate the physical scale:
$\theta = \mathrm{tan}^{-1}(\frac{d}{D}) \approx \frac{206,265 \, \mathrm{arcsec}}{1 \, \mathrm{radian}} \frac{d}{D}$,
where
$\theta$ is the angular size of the object (in our case spaxel) in arcsec,
$d$ is the diameter of the object (spaxel), and
$D$ is the angular diameter distance.
The distance (via the Hubble Law --- which is fairly accurate for low redshift objects) is
$D \approx \frac{cz}{H_0}$,
where
$c$ is the speed of light in km/s,
$z$ is the redshift, and
$H_0$ is the Hubble constant in km/s/Mpc.
Calculate $D$.
End of explanation
scale = 1 / 206265 * D * 1e6  # 1 radian = 206265 arcsec [pc / arcsec]
Explanation: Rearrange the small angle formula to solve for the scale ($\frac{d}{\theta}$) in pc / arcsec.
End of explanation
spaxel_area = (scale * spaxel_size)**2  # [pc^2]
Explanation: Now convert the spaxel size from arcsec to parsecs and calculate the area of a spaxel.
End of explanation
sigma_star = np.log10(10**mstar / spaxel_area)  # [Msun / pc^2]
Explanation: Finally, we simply divide the stellar mass by the area to get the stellar mass surface density $\Sigma_\star$ in units of $\frac{M_\odot}{pc^2}$.
End of explanation
fig, ax = plt.subplots(figsize=(6, 6))
ax.scatter(sigma_star.values[mask == 0], oh.value[mask == 0], alpha=0.15)
ax.set_xlabel('log(Mstar) [M$_\odot$]')
ax.set_ylabel('12+log(O/H)')
ax.axis([0, 4, 8.0, 8.8])
Explanation: Let's plot metallicity as a function of $\Sigma_\star$!  Remember to apply the mask.  Also set the axis range to be [0, 4, 8, 8.8].
End of explanation
# fitting formula
aa = 8.55
bb = 0.014
cc = 3.14
xx = np.linspace(1, 3, 1000)
yy = aa + bb * (xx - cc) * np.exp(-(xx - cc))
Explanation: MaNGA Spatially-Resolved Mass-Metallicity Relation
We have constructed the spatially-resolved MZR for one galaxy, but we are interested in understanding the evolution of galaxies in general, so we want to repeat this exercise for many galaxies.  In Barrera-Ballesteros et al. (2016), Jorge Barrera-Ballesteros (who gave a talk at Pitt in November 2017) did just this, and here is the analogous figure for 653 disk galaxies.
<img src="images/barrera-ballesteros_local_mzr.png" style="width: 400px;"/>
The best fit line from Barrera-Ballesteros et al. (2016) is given in the next cell.
End of explanation
fig, ax = plt.subplots(figsize=(6, 6))
ax.scatter(sigma_star.values[mask == 0], oh.value[mask == 0], alpha=0.15)
ax.plot(xx, yy)
ax.set_xlabel('log(Mstar) [M$_\odot$]')
ax.set_ylabel('12+log(O/H)')
ax.axis([0, 4, 8.0, 8.8])
Explanation: Remake the spatially-resolved MZR plot for our galaxy showing the he best fit line from Barrera-Ballesteros et al. (2016).
End of explanation
<END_TASK>