Source code for marvin.utils.dap.bpt

#!/usr/bin/env python
# encoding: utf-8
#
# bpt.py
#
# Created by José Sánchez-Gallego on 19 Jan 2017.


from __future__ import division
from __future__ import print_function
from __future__ import absolute_import

from packaging.version import parse
import warnings

import matplotlib
import matplotlib.pyplot as plt
import numpy as np

from mpl_toolkits.axes_grid1 import ImageGrid

from marvin.core.exceptions import MarvinDeprecationWarning, MarvinError
from marvin.utils.plot import bind_to_figure


__ALL__ = ('get_snr', 'kewley_sf_nii', 'kewley_sf_sii', 'kewley_sf_oi',
           'kewley_comp_nii', 'kewley_agn_sii', 'kewley_agn_oi',
           'bpt_kewley06')


def get_snr(snr_min, emission_line, default=3):
    """Convenience function to get the minimum SNR for a certain emission line.

    If ``snr_min`` is a dictionary and ``emision_line`` is one of the keys,
    returns that value. If the emission line is not included in the dictionary,
    returns ``default``. If ``snr_min`` is a float, returns that value
    regardless of the ``emission_line``.

    """

    if not isinstance(snr_min, dict):
        return snr_min

    if emission_line in snr_min:
        return snr_min[emission_line]
    else:
        return default


def get_masked(maps, emline, snr=1):
    """Convenience function to get masked arrays without negative values."""

    gflux = maps['emline_gflux_' + emline]
    gflux_masked = gflux.masked

    # Masks spaxels with flux <= 0
    gflux_masked.mask |= (gflux_masked.data <= 0)

    # Masks all spaxels that don't reach the cutoff SNR
    gflux_masked.mask |= gflux.snr < snr
    gflux_masked.mask |= gflux.ivar == 0

    return gflux_masked


def _get_kewley06_axes(use_oi=True):
    """Creates custom axes for displaying Kewley06 plots."""

    fig = plt.figure(None, (8.5, 10))
    fig.clf()

    plt.subplots_adjust(top=0.99, bottom=0.08, hspace=0.01)

    # The axes for the three classification plots
    imgrid_kwargs = {'add_all': True} if parse(matplotlib.__version__) < parse('3.5.0') else {}
    grid_bpt = ImageGrid(fig, 211,
                         nrows_ncols=(1, 3) if use_oi else (1, 2),
                         direction='row',
                         axes_pad=0.1,
                         label_mode='L',
                         share_all=False, **imgrid_kwargs)

    # The axes for the galaxy display
    gal_bpt = ImageGrid(fig, 212, nrows_ncols=(1, 1))

    # Plots the classification boundary lines
    xx_sf_nii = np.linspace(-1.281, 0.045, int(1e4))
    xx_sf_sii = np.linspace(-2, 0.315, int(1e4))
    xx_sf_oi = np.linspace(-2.5, -0.7, int(1e4))

    xx_comp_nii = np.linspace(-2, 0.4, int(1e4))

    xx_agn_sii = np.array([-0.308, 1.0])
    xx_agn_oi = np.array([-1.12, 0.5])

    grid_bpt[0].plot(xx_sf_nii, kewley_sf_nii(xx_sf_nii), 'k--', zorder=90)
    grid_bpt[1].plot(xx_sf_sii, kewley_sf_sii(xx_sf_sii), 'r-', zorder=90)
    if use_oi:
        grid_bpt[2].plot(xx_sf_oi, kewley_sf_oi(xx_sf_oi), 'r-', zorder=90)

    grid_bpt[0].plot(xx_comp_nii, kewley_comp_nii(xx_comp_nii), 'r-', zorder=90)

    grid_bpt[1].plot(xx_agn_sii, kewley_agn_sii(xx_agn_sii), 'b-', zorder=80)
    if use_oi:
        grid_bpt[2].plot(xx_agn_oi, kewley_agn_oi(xx_agn_oi), 'b-', zorder=80)

    # Adds captions
    grid_bpt[0].text(-1, -0.5, 'SF', ha='center', fontsize=12, zorder=100, color='c')
    grid_bpt[0].text(0.5, 0.5, 'AGN', ha='left', fontsize=12, zorder=100)
    grid_bpt[0].text(-0.08, -1.2, 'Comp', ha='left', fontsize=12, zorder=100, color='g')

    grid_bpt[1].text(-1.2, -0.5, 'SF', ha='center', fontsize=12, zorder=100)
    grid_bpt[1].text(-1, 1.2, 'Seyfert', ha='left', fontsize=12, zorder=100, color='r')
    grid_bpt[1].text(0.3, -1, 'LINER', ha='left', fontsize=12, zorder=100, color='m')

    if use_oi:
        grid_bpt[2].text(-2, -0.5, 'SF', ha='center', fontsize=12, zorder=100)
        grid_bpt[2].text(-1.5, 1, 'Seyfert', ha='left', fontsize=12, zorder=100)
        grid_bpt[2].text(-0.1, -1, 'LINER', ha='right', fontsize=12, zorder=100)

    # Sets the ticks, ticklabels, and other details
    xtick_limits = ((-2, 1), (-1.5, 1), (-2.5, 0.5))
    axes = [0, 1, 2] if use_oi else [0, 1]

    for ii in axes:

        grid_bpt[ii].get_xaxis().set_tick_params(direction='in')
        grid_bpt[ii].get_yaxis().set_tick_params(direction='in')

        grid_bpt[ii].set_xticks(np.arange(xtick_limits[ii][0], xtick_limits[ii][1] + 0.5, 0.5))
        grid_bpt[ii].set_xticks(np.arange(xtick_limits[ii][0],
                                          xtick_limits[ii][1] + 0.1, 0.1), minor=True)
        grid_bpt[ii].set_yticks(np.arange(-1.5, 2.0, 0.5))
        grid_bpt[ii].set_yticks(np.arange(-1.5, 1.6, 0.1), minor=True)

        grid_bpt[ii].grid(which='minor', alpha=0.2)
        grid_bpt[ii].grid(which='major', alpha=0.5)

        grid_bpt[ii].set_xlim(xtick_limits[ii][0], xtick_limits[ii][1])

        grid_bpt[ii].set_ylim(-1.5, 1.6)
        if use_oi:
            grid_bpt[ii].set_ylim(-1.5, 1.8)

        grid_bpt[ii].spines['top'].set_visible(True)

        if ii in [0, 1]:
            if not use_oi and ii == 1:
                continue
            grid_bpt[ii].get_xticklabels()[-1].set_visible(False)

    grid_bpt[0].set_ylabel(r'log([OIII]/H$\beta$)')

    grid_bpt[0].set_xlabel(r'log([NII]/H$\alpha$)')
    grid_bpt[1].set_xlabel(r'log([SII]/H$\alpha$)')
    if use_oi:
        grid_bpt[2].set_xlabel(r'log([OI]/H$\alpha$)')

    gal_bpt[0].grid(False)

    return fig, grid_bpt, gal_bpt[0]


[docs]def kewley_sf_nii(log_nii_ha): """Star forming classification line for log([NII]/Ha).""" return 0.61 / (log_nii_ha - 0.05) + 1.3
[docs]def kewley_sf_sii(log_sii_ha): """Star forming classification line for log([SII]/Ha).""" return 0.72 / (log_sii_ha - 0.32) + 1.3
[docs]def kewley_sf_oi(log_oi_ha): """Star forming classification line for log([OI]/Ha).""" return 0.73 / (log_oi_ha + 0.59) + 1.33
[docs]def kewley_comp_nii(log_nii_ha): """Composite classification line for log([NII]/Ha).""" return 0.61 / (log_nii_ha - 0.47) + 1.19
[docs]def kewley_agn_sii(log_sii_ha): """Seyfert/LINER classification line for log([SII]/Ha).""" return 1.89 * log_sii_ha + 0.76
[docs]def kewley_agn_oi(log_oi_ha): """Seyfert/LINER classification line for log([OI]/Ha).""" return 1.18 * log_oi_ha + 1.30
[docs]def bpt_kewley06(maps, snr_min=3, return_figure=True, use_oi=True, **kwargs): """Returns a classification of ionisation regions, as defined in Kewley+06. Makes use of the classification system defined by `Kewley et al. (2006) <https://ui.adsabs.harvard.edu/#abs/2006MNRAS.372..961K/abstract>`_ to return classification masks for different ionisation mechanisms. If ``return_figure=True``, produces and returns a matplotlib figure with the classification plots (based on Kewley+06 Fig. 4) and the 2D spatial distribution of classified spaxels (i.e., a map of the galaxy in which each spaxel is colour-coded based on its emission mechanism). While it is possible to call this function directly, its normal use will be via the :func:`~marvin.tools.maps.Maps.get_bpt` method. Parameters: maps (a Marvin :class:`~marvin.tools.maps.Maps` object) The Marvin Maps object that contains the emission line maps to be used to determine the BPT classification. snr_min (float or dict): The signal-to-noise cutoff value for the emission lines used to generate the BPT diagram. If ``snr_min`` is a single value, that signal-to-noise will be used for all the lines. Alternatively, a dictionary of signal-to-noise values, with the emission line channels as keys, can be used. E.g., ``snr_min={'ha': 5, 'nii': 3, 'oi': 1}``. If some values are not provided, they will default to ``SNR>=3``. Note that the value ``sii`` will be applied to both ``[SII 6718]`` and ``[SII 6732]``. return_figure (bool): If ``True``, it also returns the matplotlib figure_ of the BPT diagram plot, which can be used to modify the style of the plot. use_oi (bool): If ``True``, uses the OI diagnostic diagram for spaxel classification. Returns: bpt_return: ``bpt_kewley06`` returns a dictionary of dictionaries of classification masks. The classification masks (not to be confused with bitmasks) are boolean arrays with the same shape as the Maps or Cube (without the spectral dimension) that can be used to select spaxels belonging to a certain excitation process (e.g., star forming). The returned dictionary has the following keys: ``'sf'`` (star forming), ``'comp'`` (composite), ``'agn'``, ``'seyfert'``, ``'liner'``, ``'invalid'`` (spaxels that are masked out at the DAP level), and ``'ambiguous'`` (good spaxels that do not fall in any classification or fall in more than one). Each key provides access to a new dictionary with keys ``'nii'`` (for the constraints in the diagram NII/Halpha vs OIII/Hbeta), ``'sii'`` (SII/Halpha vs OIII/Hbeta), ``'oi'`` (OI/Halpha vs OIII/Hbeta; only if ``use_oi=True``), and ``'global'``, which applies all the previous constraints at once. The ``'ambiguous'`` mask only contains the ``'global'`` subclassification, while the ``'comp'`` dictionary only contains ``'nii'``. ``'nii'`` is not available for ``'seyfert'`` and ``'liner'``. All the global masks are unique (a spaxel can only belong to one of them) with the exception of ``'agn'``, which intersects with ``'seyfert'`` and ``'liner'``. Additionally, if ``return_figure=True``, ``bpt_kewley06`` will also return the matplotlib figure for the generated plot, and a list of axes for each one of the subplots. Example: >>> maps_8485_1901 = Maps(plateifu='8485-1901') >>> bpt_masks, fig, axes = bpt_kewley06(maps_8485_1901) Gets the global mask for star forming spaxels >>> sf = bpt_masks['sf']['global'] Gets the seyfert mask based only on the SII/Halpha vs OIII/Hbeta diagnostics >>> seyfert_sii = bpt_masks['seyfert']['sii'] """ if 'snr' in kwargs: warnings.warn('snr is deprecated. Use snr_min instead. ' 'snr will be removed in a future version of marvin', MarvinDeprecationWarning) snr_min = kwargs.pop('snr') if len(kwargs.keys()) > 0: raise MarvinError('unknown keyword {0}'.format(list(kwargs.keys())[0])) # Gets the necessary emission line maps oiii = get_masked(maps, 'oiii_5008', snr=get_snr(snr_min, 'oiii')) nii = get_masked(maps, 'nii_6585', snr=get_snr(snr_min, 'nii')) ha = get_masked(maps, 'ha_6564', snr=get_snr(snr_min, 'ha')) hb = get_masked(maps, 'hb_4862', snr=get_snr(snr_min, 'hb')) oi = get_masked(maps, 'oi_6302', snr=get_snr(snr_min, 'oi')) sii_6718 = get_masked(maps, 'sii_6718', snr=get_snr(snr_min, 'sii')) sii_6732 = get_masked(maps, 'sii_6732', snr=get_snr(snr_min, 'sii')) sii = sii_6718 + sii_6732 # Calculate masked logarithms log_oiii_hb = np.ma.log10(oiii / hb) log_nii_ha = np.ma.log10(nii / ha) log_sii_ha = np.ma.log10(sii / ha) log_oi_ha = np.ma.log10(oi / ha) # Calculates masks for each emission mechanism according to the paper boundaries. # The log_nii_ha < 0.05, log_sii_ha < 0.32, etc are necessary because the classification lines # diverge and we only want the region before the asymptota. sf_mask_nii = ((log_oiii_hb < kewley_sf_nii(log_nii_ha)) & (log_nii_ha < 0.05)).filled(False) sf_mask_sii = ((log_oiii_hb < kewley_sf_sii(log_sii_ha)) & (log_sii_ha < 0.32)).filled(False) sf_mask_oi = ((log_oiii_hb < kewley_sf_oi(log_oi_ha)) & (log_oi_ha < -0.59)).filled(False) sf_mask = sf_mask_nii & sf_mask_sii & sf_mask_oi if use_oi else sf_mask_nii & sf_mask_sii comp_mask = ((log_oiii_hb > kewley_sf_nii(log_nii_ha)) & (log_nii_ha < 0.05)).filled(False) & \ ((log_oiii_hb < kewley_comp_nii(log_nii_ha)) & (log_nii_ha < 0.465)).filled(False) comp_mask &= (sf_mask_sii & sf_mask_oi) if use_oi else sf_mask_sii agn_mask_nii = ((log_oiii_hb > kewley_comp_nii(log_nii_ha)) | (log_nii_ha > 0.465)).filled(False) agn_mask_sii = ((log_oiii_hb > kewley_sf_sii(log_sii_ha)) | (log_sii_ha > 0.32)).filled(False) agn_mask_oi = ((log_oiii_hb > kewley_sf_oi(log_oi_ha)) | (log_oi_ha > -0.59)).filled(False) agn_mask = agn_mask_nii & agn_mask_sii & agn_mask_oi if use_oi else agn_mask_nii & agn_mask_sii seyfert_mask_sii = agn_mask & (kewley_agn_sii(log_sii_ha) < log_oiii_hb).filled(False) seyfert_mask_oi = agn_mask & (kewley_agn_oi(log_oi_ha) < log_oiii_hb).filled(False) seyfert_mask = seyfert_mask_sii & seyfert_mask_oi if use_oi else seyfert_mask_sii liner_mask_sii = agn_mask & (kewley_agn_sii(log_sii_ha) > log_oiii_hb).filled(False) liner_mask_oi = agn_mask & (kewley_agn_oi(log_oi_ha) > log_oiii_hb).filled(False) liner_mask = liner_mask_sii & liner_mask_oi if use_oi else liner_mask_sii # The invalid mask is the combination of spaxels that are invalid in all of the emission maps invalid_mask_nii = ha.mask | oiii.mask | nii.mask | hb.mask invalid_mask_sii = ha.mask | oiii.mask | sii.mask | hb.mask invalid_mask_oi = ha.mask | oiii.mask | oi.mask | hb.mask invalid_mask = ha.mask | oiii.mask | nii.mask | hb.mask | sii.mask if use_oi: invalid_mask |= oi.mask # The ambiguous mask are spaxels that are not invalid but don't fall into any of the # emission mechanism classifications. ambiguous_mask = ~(sf_mask | comp_mask | seyfert_mask | liner_mask) & ~invalid_mask sf_classification = {'global': sf_mask, 'nii': sf_mask_nii, 'sii': sf_mask_sii} comp_classification = {'global': comp_mask, 'nii': comp_mask} agn_classification = {'global': agn_mask, 'nii': agn_mask_nii, 'sii': agn_mask_sii} seyfert_classification = {'global': seyfert_mask, 'sii': seyfert_mask_sii} liner_classification = {'global': liner_mask, 'sii': liner_mask_sii} invalid_classification = {'global': invalid_mask, 'nii': invalid_mask_nii, 'sii': invalid_mask_sii} ambiguous_classification = {'global': ambiguous_mask} if use_oi: sf_classification['oi'] = sf_mask_oi agn_classification['oi'] = agn_mask_oi seyfert_classification['oi'] = seyfert_mask_oi liner_classification['oi'] = liner_mask_oi invalid_classification['oi'] = invalid_mask_oi bpt_return_classification = {'sf': sf_classification, 'comp': comp_classification, 'agn': agn_classification, 'seyfert': seyfert_classification, 'liner': liner_classification, 'invalid': invalid_classification, 'ambiguous': ambiguous_classification} if not return_figure: return bpt_return_classification # Does all the plotting with plt.style.context('seaborn-darkgrid'): fig, grid_bpt, gal_bpt = _get_kewley06_axes(use_oi=use_oi) sf_kwargs = {'marker': 's', 's': 12, 'color': 'c', 'zorder': 50, 'alpha': 0.7, 'lw': 0.0, 'label': 'Star-forming'} sf_handler = grid_bpt[0].scatter(log_nii_ha[sf_mask], log_oiii_hb[sf_mask], **sf_kwargs) grid_bpt[1].scatter(log_sii_ha[sf_mask], log_oiii_hb[sf_mask], **sf_kwargs) comp_kwargs = {'marker': 's', 's': 12, 'color': 'g', 'zorder': 45, 'alpha': 0.7, 'lw': 0.0, 'label': 'Composite'} comp_handler = grid_bpt[0].scatter(log_nii_ha[comp_mask], log_oiii_hb[comp_mask], **comp_kwargs) grid_bpt[1].scatter(log_sii_ha[comp_mask], log_oiii_hb[comp_mask], **comp_kwargs) seyfert_kwargs = {'marker': 's', 's': 12, 'color': 'r', 'zorder': 40, 'alpha': 0.7, 'lw': 0.0, 'label': 'Seyfert'} seyfert_handler = grid_bpt[0].scatter(log_nii_ha[seyfert_mask], log_oiii_hb[seyfert_mask], **seyfert_kwargs) grid_bpt[1].scatter(log_sii_ha[seyfert_mask], log_oiii_hb[seyfert_mask], **seyfert_kwargs) liner_kwargs = {'marker': 's', 's': 12, 'color': 'm', 'zorder': 35, 'alpha': 0.7, 'lw': 0.0, 'label': 'LINER'} liner_handler = grid_bpt[0].scatter(log_nii_ha[liner_mask], log_oiii_hb[liner_mask], **liner_kwargs) grid_bpt[1].scatter(log_sii_ha[liner_mask], log_oiii_hb[liner_mask], **liner_kwargs) amb_kwargs = {'marker': 's', 's': 12, 'color': '0.6', 'zorder': 30, 'alpha': 0.7, 'lw': 0.0, 'label': 'Ambiguous '} amb_handler = grid_bpt[0].scatter(log_nii_ha[ambiguous_mask], log_oiii_hb[ambiguous_mask], **amb_kwargs) grid_bpt[1].scatter(log_sii_ha[ambiguous_mask], log_oiii_hb[ambiguous_mask], **amb_kwargs) if use_oi: grid_bpt[2].scatter(log_oi_ha[sf_mask], log_oiii_hb[sf_mask], **sf_kwargs) grid_bpt[2].scatter(log_oi_ha[comp_mask], log_oiii_hb[comp_mask], **comp_kwargs) grid_bpt[2].scatter(log_oi_ha[seyfert_mask], log_oiii_hb[seyfert_mask], **seyfert_kwargs) grid_bpt[2].scatter(log_oi_ha[liner_mask], log_oiii_hb[liner_mask], **liner_kwargs) grid_bpt[2].scatter(log_oi_ha[ambiguous_mask], log_oiii_hb[ambiguous_mask], **amb_kwargs) # Creates the legend grid_bpt[0].legend([sf_handler, comp_handler, seyfert_handler, liner_handler, amb_handler], ['Star-forming', 'Composite', 'Seyfert', 'LINER', 'Ambiguous'], ncol=2, loc='upper left', frameon=True, labelspacing=0.1, columnspacing=0.1, handletextpad=0.1, fontsize=9) # Creates a RGB image of the galaxy, and sets the colours of the spaxels to match the # classification masks gal_rgb = np.zeros((ha.shape[0], ha.shape[1], 3), dtype=np.uint8) for ii in [1, 2]: # Cyan gal_rgb[:, :, ii][sf_mask] = 255 gal_rgb[:, :, 1][comp_mask] = 128 # Green gal_rgb[:, :, 0][seyfert_mask] = 255 # Red # Magenta gal_rgb[:, :, 0][liner_mask] = 255 gal_rgb[:, :, 2][liner_mask] = 255 for ii in [0, 1, 2]: gal_rgb[:, :, ii][invalid_mask] = 255 # White gal_rgb[:, :, ii][ambiguous_mask] = 169 # Grey # Shows the image. gal_bpt.imshow(gal_rgb, origin='lower', aspect='auto', interpolation='nearest') gal_bpt.set_xlim(0, ha.shape[1] - 1) gal_bpt.set_ylim(0, ha.shape[0] - 1) gal_bpt.set_xlabel('x [spaxels]') gal_bpt.set_ylabel('y [spaxels]') axes = grid_bpt.axes_all + [gal_bpt] # Adds custom method to create figure for ax in axes: setattr(ax.__class__, 'bind_to_figure', _bind_to_figure) return (bpt_return_classification, fig, axes)
def _bind_to_figure(self, fig=None): """Copies axes to a new figure. Uses ``marvin.utils.plot.utils.bind_to_figure`` with a number of tweaks. """ new_figure = bind_to_figure(self, fig=fig) if new_figure.axes[0].get_ylabel() == '': new_figure.axes[0].set_ylabel('log([OIII]/H$\\beta$)') return new_figure