Source code for species.plot.plot_comparison

"""
Module with functions for plotting results from a spectral
analysis that compares data with a library of empirical
spectra or a grid of model spectra.
"""

import configparser
import os
import warnings

from typing import List, Optional, Tuple

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

from matplotlib.ticker import AutoMinorLocator
from scipy.interpolate import interp1d, RegularGridInterpolator
from typeguard import typechecked

from species.core import constants
from species.read.read_model import ReadModel
from species.read.read_object import ReadObject
from species.util.dust_util import ism_extinction
from species.util.plot_util import create_model_label, update_labels
from species.util.spec_util import smooth_spectrum


[docs] @typechecked def plot_statistic( tag: str, xlim: Optional[Tuple[float, float]] = None, ylim: Optional[Tuple[float, float]] = None, title: Optional[str] = None, offset: Optional[Tuple[float, float]] = None, figsize: Optional[Tuple[float, float]] = (4.0, 2.5), output: Optional[str] = None, ) -> mpl.figure.Figure: """ Function for plotting the goodness-of-fit statistic from a comparison with an empirical spectral library with :class:`~species.fit.compare_spectra.CompareSpectra.spectral_type` that enables a determination of the spectral type Parameters ---------- tag : str Database tag where the results from the empirical comparison with :class:`~species.fit.compare_spectra.CompareSpectra.spectral_type` are stored. xlim : tuple(float, float) Limits of the spectral type axis in numbers (i.e. 0=M0, 5=M5, 10=L0, etc.). ylim : tuple(float, float) Limits of the goodness-of-fit axis. title : str Plot title. offset : tuple(float, float) Offset for the label of the x- and y-axis. figsize : tuple(float, float) Figure size. output : str, None Output filename for the plot. The plot is shown in an interface window if the argument is set to ``None``. Returns ------- NoneType None """ if output is None: print("Plotting goodness-of-fit statistic...", end="") else: print(f"Plotting goodness-of-fit statistic: {output}...", end="") config_file = os.path.join(os.getcwd(), "species_config.ini") config = configparser.ConfigParser() config.read(config_file) db_path = config["species"]["database"] with h5py.File(db_path, "r") as hdf5_file: dset = hdf5_file[f"results/empirical/{tag}/names"] names = np.array(dset) sptypes = np.array(hdf5_file[f"results/empirical/{tag}/sptypes"]) g_fit = np.array(hdf5_file[f"results/empirical/{tag}/goodness_of_fit"]) plt.rcParams["font.family"] = "serif" plt.rcParams["mathtext.fontset"] = "dejavuserif" plt.rcParams["axes.axisbelow"] = False fig = plt.figure(figsize=figsize) gridsp = mpl.gridspec.GridSpec(1, 1) gridsp.update(wspace=0, hspace=0, left=0, right=1, bottom=0, top=1) ax = plt.subplot(gridsp[0, 0]) ax.tick_params( axis="both", which="major", colors="black", labelcolor="black", direction="in", width=1, length=5, labelsize=12, top=True, bottom=True, left=True, right=True, ) ax.tick_params( axis="both", which="minor", colors="black", labelcolor="black", direction="in", width=1, length=3, labelsize=12, top=True, bottom=True, left=True, right=True, ) ax.xaxis.set_minor_locator(AutoMinorLocator(5)) ax.yaxis.set_minor_locator(AutoMinorLocator(5)) ax.set_xlabel("Spectral type", fontsize=13) ax.set_ylabel(r"G$_\mathregular{k}$", fontsize=13) if offset is not None: ax.get_xaxis().set_label_coords(0.5, offset[0]) ax.get_yaxis().set_label_coords(offset[1], 0.5) else: ax.get_xaxis().set_label_coords(0.5, -0.1) ax.get_yaxis().set_label_coords(-0.1, 0.5) if title is not None: ax.set_title(title, y=1.02, fontsize=13) ax.set_xticks(np.linspace(0.0, 30.0, 7, endpoint=True)) ax.set_xticklabels(["M0", "M5", "L0", "L5", "T0", "T5", "Y0"]) if xlim is None: ax.set_xlim(0.0, 30.0) else: ax.set_xlim(xlim[0], xlim[1]) if ylim is not None: ax.set_ylim(ylim[0], ylim[1]) sptype_num = np.zeros(names.shape[0]) for i, item in enumerate(sptypes): for j in range(10): if not isinstance(item, str): item = item.decode("utf-8") if item == f"M{j}": sptype_num[i] = float(j) elif item == f"L{j}": sptype_num[i] = float(10 + j) elif item == f"T{j}": sptype_num[i] = float(20 + j) ax.plot( sptype_num, g_fit, "s", ms=3.0, mew=0.5, color="lightgray", markeredgecolor="darkgray", ) if output is None: plt.show() else: plt.savefig(output, bbox_inches="tight") print(" [DONE]") return fig
[docs] @typechecked def plot_empirical_spectra( tag: str, n_spectra: Optional[int] = None, flux_offset: Optional[float] = None, label_pos: Optional[Tuple[float, float]] = None, xlim: Optional[Tuple[float, float]] = None, ylim: Optional[Tuple[float, float]] = None, title: Optional[str] = None, offset: Optional[Tuple[float, float]] = None, figsize: Optional[Tuple[float, float]] = (4.0, 2.5), output: Optional[str] = None, ) -> mpl.figure.Figure: """ Function for plotting the results from the empirical spectrum comparison. Parameters ---------- tag : str Database tag where the results from the empirical comparison with :class:`~species.fit.compare_spectra.CompareSpectra.spectral_type` are stored. n_spectra : int, None The number of spectra with the lowest goodness-of-fit statistic that will be plotted in comparison with the data. All spectra are selected if the argument is set to ``None``. label_pos : tuple(float, float), None Position for the name labels. Should be provided as (x, y) for the lowest spectrum. The ``flux_offset`` will be applied to the remaining spectra. The labels are only plotted if the argument of both ``label_pos`` and ``flux_offset`` are not ``None``. flux_offset : float, None Offset to be applied such that the spectra do not overlap. No offset is applied if the argument is set to ``None``. xlim : tuple(float, float) Limits of the spectral type axis. ylim : tuple(float, float) Limits of the goodness-of-fit axis. title : str Plot title. offset : tuple(float, float) Offset for the label of the x- and y-axis. figsize : tuple(float, float) Figure size. output : str, None Output filename for the plot. The plot is shown in an interface window if the argument is set to ``None``. Returns ------- matplotlib.figure.Figure The ``Figure`` object that can be used for further customization of the plot. """ if output is None: print("Plotting empirical spectra comparison...", end="") else: print(f"Plotting empirical spectra comparison: {output}...", end="") if flux_offset is None: flux_offset = 0.0 config_file = os.path.join(os.getcwd(), "species_config.ini") config = configparser.ConfigParser() config.read(config_file) db_path = config["species"]["database"] with h5py.File(db_path, "r") as hdf5_file: dset = hdf5_file[f"results/empirical/{tag}/names"] object_name = dset.attrs["object_name"] spec_library = dset.attrs["spec_library"] n_spec_name = dset.attrs["n_spec_name"] spec_name = [] for i in range(n_spec_name): spec_name.append(dset.attrs[f"spec_name{i}"]) names = np.array(dset) flux_scaling = np.array(hdf5_file[f"results/empirical/{tag}/flux_scaling"]) av_ext = np.array(hdf5_file[f"results/empirical/{tag}/av_ext"]) rad_vel = np.array(hdf5_file[f"results/empirical/{tag}/rad_vel"]) rad_vel *= 1e3 # (m s-1) if n_spectra is None: n_spectra = names.size plt.rcParams["font.family"] = "serif" plt.rcParams["mathtext.fontset"] = "dejavuserif" plt.rcParams["axes.axisbelow"] = False fig = plt.figure(figsize=figsize) gridsp = mpl.gridspec.GridSpec(1, 1) gridsp.update(wspace=0, hspace=0, left=0, right=1, bottom=0, top=1) ax = plt.subplot(gridsp[0, 0]) ax.tick_params( axis="both", which="major", colors="black", labelcolor="black", direction="in", width=1, length=5, labelsize=12, top=True, bottom=True, left=True, right=True, ) ax.tick_params( axis="both", which="minor", colors="black", labelcolor="black", direction="in", width=1, length=3, labelsize=12, top=True, bottom=True, left=True, right=True, ) ax.xaxis.set_minor_locator(AutoMinorLocator(5)) ax.yaxis.set_minor_locator(AutoMinorLocator(5)) ax.set_xlabel("Wavelength (\N{GREEK SMALL LETTER MU}m)", fontsize=13) if flux_offset == 0.0: ax.set_ylabel( r"$\mathregular{F}_\lambda$" + " (W m$^{-2}$ \N{GREEK SMALL LETTER MU}m$^{-1}$)", fontsize=11, ) else: ax.set_ylabel( r"$\mathregular{F}_\lambda$ (W m$^{-2}$" + " \N{GREEK SMALL LETTER MU}m$^{-1}$) + offset", fontsize=11, ) if xlim is not None: ax.set_xlim(xlim[0], xlim[1]) if ylim is not None: ax.set_ylim(ylim[0], ylim[1]) if offset is not None: ax.get_xaxis().set_label_coords(0.5, offset[0]) ax.get_yaxis().set_label_coords(offset[1], 0.5) else: ax.get_xaxis().set_label_coords(0.5, -0.1) ax.get_yaxis().set_label_coords(-0.1, 0.5) if title is not None: ax.set_title(title, y=1.02, fontsize=13) read_obj = ReadObject(object_name) obj_spec = [] obj_res = [] for item in spec_name: obj_spec.append(read_obj.get_spectrum()[item][0]) obj_res.append(read_obj.get_spectrum()[item][3]) if flux_offset == 0.0: for spec_item in obj_spec: ax.plot(spec_item[:, 0], spec_item[:, 1], "-", lw=0.5, color="black") for i in range(n_spectra): if isinstance(names[i], str): name_item = names[i] else: name_item = names[i].decode("utf-8") with h5py.File(db_path, "r") as hdf5_file: dset = hdf5_file[f"spectra/{spec_library}/{name_item}"] sptype = dset.attrs["sptype"] spectrum = np.asarray(dset) if flux_offset != 0.0: for spec_item in obj_spec: ax.plot( spec_item[:, 0], (n_spectra - i - 1) * flux_offset + spec_item[:, 1], "-", lw=0.5, color="black", ) for j, spec_item in enumerate(obj_spec): ism_ext = ism_extinction(av_ext[i], 3.1, spectrum[:, 0]) ext_scaling = 10.0 ** (-0.4 * ism_ext) wavel_shifted = ( spectrum[:, 0] + spectrum[:, 0] * rad_vel[i] / constants.LIGHT ) flux_smooth = smooth_spectrum( wavel_shifted, spectrum[:, 1] * ext_scaling, spec_res=obj_res[j], force_smooth=True, ) interp_spec = interp1d( spectrum[:, 0], flux_smooth, fill_value="extrapolate" ) indices = np.where( (spec_item[:, 0] > np.amin(spectrum[:, 0])) & (spec_item[:, 0] < np.amax(spectrum[:, 0])) )[0] flux_resample = interp_spec(spec_item[indices, 0]) ax.plot( spec_item[indices, 0], (n_spectra - i - 1) * flux_offset + flux_scaling[i][j] * flux_resample, color="tomato", lw=0.5, ) if label_pos is not None and flux_offset != 0.0: label_text = name_item + ", " + sptype if av_ext[i] != 0.0: label_text += r", A$_\mathregular{V}$ = " + f"{av_ext[i]:.1f}" ax.text( label_pos[0], label_pos[1] + (n_spectra - i - 1) * flux_offset, label_text, fontsize=8.0, ha="left", ) if output is None: plt.show() else: plt.savefig(output, bbox_inches="tight") print(" [DONE]") return fig
[docs] @typechecked def plot_grid_statistic( tag: str, upsample: bool = False, xlim: Optional[Tuple[float, float]] = None, ylim: Optional[Tuple[float, float]] = None, title: Optional[str] = None, offset: Optional[Tuple[float, float]] = None, figsize: Optional[Tuple[float, float]] = (4.0, 2.5), output: Optional[str] = None, extra_param: Optional[str] = None, nlevels_main: int = 20, nlevels_extra: int = 10, ) -> mpl.figure.Figure: """ Function for plotting the results from the comparison with a grid of empirical or model spectra Parameters ---------- tag : str Database tag where the results from the comparison with :class:`~species.fit.compare_spectra.CompareSpectra` are stored. upsample : bool Upsample the goodness-of-fit grid to a higher resolution for a smoother appearance. xlim : tuple(float, float), None Limits of the x-axis (spectral type or effective temperature). ylim : tuple(float, float), None Limits of the y-axis. title : str, None Title that is shown above the plot. offset : tuple(float, float), None Offset for the label for the x- and y-axis. figsize : tuple(float, float), None Figure size. output : str, None Output filename for the plot. The plot is shown in an interface window if the argument is set to ``None``. extra_param : str, None Extra parameter to be overplotted with contours. The argument can be set to any of the atmospheric parameters that were used for the comparison, for example, 'teff', 'logg', 'feh', or 'radius'. Optionally, the argument can be set to 'ism_ext' in case the ``av_points`` parameter of the :func:`~species.fit.compare_spectra.CompareSpectra.compare_model` method was used. Extra contours are not plotted if the argument is set to ``None``. nlevels_main : int Number of contour levels for the main plot. nlevels_extra : int Number of contour levels for the optional extra parameter. Returns ------- matplotlib.figure.Figure The ``Figure`` object that can be used for further customization of the plot. """ if output is None: print("Plotting goodness-of-fit of model grid...", end="") else: print(f"Plotting goodness-of-fit of model grid: {output}...", end="") config_file = os.path.join(os.getcwd(), "species_config.ini") config = configparser.ConfigParser() config.read(config_file) db_path = config["species"]["database"] with h5py.File(db_path, "r") as hdf5_file: dset = hdf5_file[f"results/comparison/{tag}/goodness_of_fit"] goodness_fit = np.array(dset) n_param = dset.attrs["n_param"] parallax = dset.attrs["parallax"] flux_scaling = np.array(hdf5_file[f"results/comparison/{tag}/flux_scaling"]) radius = ( np.sqrt(flux_scaling * (constants.PARSEC / (1e-3 * parallax)) ** 2) / constants.R_JUP ) # if 'extra_scaling' in hdf5_file[f'results/comparison/{tag}']: # extra_scaling = np.array(hdf5_file[f'results/comparison/{tag}/extra_scaling']) # else: # extra_scaling = None model_param = [] coord_points = [] for i in range(n_param): model_param.append(dset.attrs[f"parameter{i}"]) coord_points.append( np.array(hdf5_file[f"results/comparison/{tag}/coord_points{i}"]) ) object_name = dset.attrs["object_name"] read_obj = ReadObject(object_name) n_wavel = 0 for item in read_obj.get_spectrum().values(): n_wavel += item[0].shape[0] # Set the coordinate for the x-axis to Teff coord_x = coord_points[0] coord_y = None param_y = None for i, item in enumerate(coord_points[1:]): if len(item) > 1: # Set the coordinate for the y-axis to the 1st axis after # Teff that has a length larger than 1 coord_y = item param_y = model_param[i + 1] break plt.rcParams["font.family"] = "serif" plt.rcParams["mathtext.fontset"] = "dejavuserif" plt.rcParams["axes.axisbelow"] = False fig = plt.figure(figsize=figsize) if coord_y is None: # Create a line plot if there is not a parameter for the y-axis gridsp = mpl.gridspec.GridSpec(1, 1) gridsp.update(wspace=0, hspace=0, left=0, right=1, bottom=0, top=1) ax = plt.subplot(gridsp[0, 0]) else: # Create a contour plot if there is a second parameter to show gridsp = mpl.gridspec.GridSpec(1, 2, width_ratios=[4.0, 0.25]) gridsp.update(wspace=0.07, hspace=0, left=0, right=1, bottom=0, top=1) ax = plt.subplot(gridsp[0, 0]) ax_cb = plt.subplot(gridsp[0, 1]) ax.tick_params( axis="both", which="major", colors="black", labelcolor="black", direction="in", width=1, length=5, labelsize=11, top=True, bottom=True, left=True, right=True, pad=5, ) ax.tick_params( axis="both", which="minor", colors="black", labelcolor="black", direction="in", width=1, length=3, labelsize=11, top=True, bottom=True, left=True, right=True, pad=5, ) ax.xaxis.set_minor_locator(AutoMinorLocator(5)) ax.yaxis.set_minor_locator(AutoMinorLocator(5)) ax.set_xlabel(r"$T_\mathregular{eff}$ (K)", fontsize=13.0) if param_y is None: ax.set_ylabel(r"$\Delta\log\,G_k$", fontsize=13.0) elif param_y == "ism_ext": ax.set_ylabel(r"$\mathregular{A}_\mathregular{V}$", fontsize=13.0) elif param_y == "logg": ax.set_ylabel(r"$\mathregular{log}\,g$", fontsize=13.0) elif param_y == "feh": ax.set_ylabel("[Fe/H]", fontsize=13.0) elif param_y == "c_o_ratio": ax.set_ylabel("C/O", fontsize=13.0) elif param_y == "fsed": ax.set_ylabel(r"$f_\mathregular{sed}$", fontsize=13.0) if xlim is not None: ax.set_xlim(xlim[0], xlim[1]) if ylim is not None: ax.set_ylim(ylim[0], ylim[1]) if offset is not None: ax.get_xaxis().set_label_coords(0.5, offset[0]) ax.get_yaxis().set_label_coords(offset[1], 0.5) else: ax.get_xaxis().set_label_coords(0.5, -0.11) ax.get_yaxis().set_label_coords(-0.1, 0.5) if title is not None: ax.set_title(title, y=1.02, fontsize=14.0) # Sum/collapse over parameters with a single value for i, item in enumerate(coord_points): if len(item) == 1: goodness_fit = np.sum(goodness_fit, axis=i) # Indices of the best-fit model best_index = np.unravel_index(np.nanargmin(goodness_fit), goodness_fit.shape) if len(model_param) > 2 or extra_param == "radius": if extra_param is None: # Select all axes beyond the 1st and 2nd axis ax_list = [] for i in range(len(coord_points) - 2): ax_list.append(2 + i) # Select minimum G_k along the axes of ax_list # This creates a 3D array of goodness_fit goodness_fit = np.nanmin(goodness_fit, axis=tuple(ax_list)) else: if extra_param == "radius": goodness_full = goodness_fit.copy() # Select all axes beyond the 2nd axis ax_list = [] for i in range(len(coord_points) - 2): ax_list.append(2 + i) # Select minimum G_k along the axes of ax_list # This creates a 2D array from goodness_fit goodness_fit = np.nanmin(goodness_fit, axis=tuple(ax_list)) extra_map = np.zeros(goodness_fit.shape[:2]) for i in range(goodness_fit.shape[0]): for j in range(goodness_fit.shape[1]): # Get the indices in the goodness_full array # for the values from the collapsed (2D) # goodness_fit array min_idx = np.argwhere(goodness_fit[i, j] == goodness_full) if len(min_idx) == 0: extra_map[i, j] = np.nan continue if len(min_idx) > 1: warnings.warn( f"Found {len(min_idx)} positions in the " "goodness-of-fit grid with the value " f"{goodness_fit[i, j]:.2f}. Using the " "first position in the list but this " "warning is not expected to have occurred." ) extra_map[i, j] = radius[tuple(min_idx[0])] else: extra_idx = model_param.index(extra_param) if extra_idx != 2: goodness_fit = np.swapaxes(goodness_fit, extra_idx, 2) coord_points_new = [ coord_points[0], coord_points[1], coord_points[extra_idx], ] for coord_idx in range(len(coord_points)): if coord_idx in [0, 1, extra_idx]: continue # Add remaining coordinate points but skip # 1st, 2nd, and extra axis that have already # been added before the for loop coord_points_new.append(coord_points[coord_idx]) coord_points = coord_points_new.copy() # Set the extra axis to the 3rd axis after the swap extra_idx = 2 if len(model_param) > 3: # Select all axes beyond the axis of extra_param ax_list = [] for i in range(len(coord_points) - 3): ax_list.append(3 + i) # Select minimum G_k along the axes of ax_list # This creates a 3D array of goodness_fit goodness_fit = np.nanmin(goodness_fit, axis=tuple(ax_list)) # Indices with the minimum G_k for the tested values # of the 3rd axis, that is, the extra_param axis # This creates a 2D array with the shape of the # 1st and 2nd axis of goodness_fit indices = np.argmin(goodness_fit, axis=extra_idx) # Select minimum G_k along the 3rd axis # This creates a 2D array of goodness_fit goodness_fit = np.nanmin(goodness_fit, axis=extra_idx) # Create 2D array with the extra map that will # be plotted as contours over the main map extra_map = np.zeros(goodness_fit.shape) for i in range(extra_map.shape[0]): for j in range(extra_map.shape[1]): extra_map[i, j] = coord_points[extra_idx][indices[i, j]] # Transpose for plot so make Teff the x axis goodness_fit = np.transpose(goodness_fit) if extra_param is not None: extra_map = np.transpose(extra_map) if coord_y is not None: if upsample: fit_interp = RegularGridInterpolator((coord_y, coord_x), goodness_fit) x_new = np.linspace(coord_x[0], coord_x[-1], 50) y_new = np.linspace(coord_y[0], coord_y[-1], 50) x_grid, y_grid = np.meshgrid(x_new, y_new) goodness_fit = fit_interp((y_grid, x_grid)) else: x_grid, y_grid = np.meshgrid(coord_x, coord_y) goodness_fit = np.log10(goodness_fit) goodness_fit -= np.nanmin(goodness_fit) nan_points = np.sum(np.isnan(goodness_fit)) if nan_points > 0: warnings.warn( f"Found {nan_points} NaN values in the " "goodness-of-fit grid. These points will be set " "to zero in the contour map." ) goodness_fit = np.nan_to_num(goodness_fit) if coord_y is None: ax.plot( coord_x, goodness_fit[0,], ) else: c = ax.contourf(x_grid, y_grid, goodness_fit, levels=nlevels_main) cb = mpl.colorbar.Colorbar( ax=ax_cb, mappable=c, orientation="vertical", ticklocation="right", format="%.1f", ) cb.ax.minorticks_on() cb.ax.tick_params( which="major", width=0.8, length=5, labelsize=12, direction="in", color="black", ) cb.ax.tick_params( which="minor", width=0.8, length=3, labelsize=12, direction="in", color="black", ) cb.ax.set_ylabel( r"$\Delta\log\,G_k$", rotation=270, labelpad=22, fontsize=13.0, ) if extra_param is not None: if upsample: extra_interp = RegularGridInterpolator((coord_y, coord_x), extra_map) extra_map = extra_interp((y_grid, x_grid)) cs = ax.contour( x_grid, y_grid, extra_map, levels=nlevels_extra, colors="white", linewidths=0.7, ) else: cs = ax.contour( coord_x, coord_y, extra_map, levels=nlevels_extra, colors="white", linewidths=0.7, ) # manual = [(2350, 0.8), (2500, 0.8), (2600, 0.8), (2700, 0.8), # (2800, 0.8), (2950, 0.8), (3100, 0.8), (3300, 0.8)] ax.clabel(cs, cs.levels, inline=True, fontsize=8, fmt="%1.1f") # if extra_scaling is not None and len(coord_points[2]) > 1: # ratio = np.transpose(flux_scaling[:, 0, :])/np.transpose(extra_scaling[:, 0, :, 0]) # # cs = ax.contour(coord_x, coord_y, ratio, levels=nlevels_extra, colors='white', # linestyles='-', linewidths=0.7) # # ax.clabel(cs, cs.levels, inline=True, fontsize=8, fmt='%1.1f') ax.plot( coord_x[best_index[0]], coord_y[best_index[1]], marker="X", ms=10.0, color="#eb4242", mfc="#eb4242", mec="black", ) # best_param = (coord_x[best_index[0]], coord_y[best_index[1]]) # # par_key, par_unit, par_label = quantity_unit(model_param, object_type='planet') # # par_text = f'{par_label[0]} = {best_param[0]:.0f} {par_unit[0]}\n' \ # f'{par_label[1]} = {best_param[1]:.1f}' # # ax.annotate(par_text, (best_param[0]+50., best_param[1]), ha='left', va='center', # color='white', fontsize=12.) if extra_param is not None: extra_label = update_labels([extra_param])[0] ax.plot([], [], ls="-", lw=1.2, color="white", label=extra_label) ax.legend(loc="best", frameon=False, labelcolor="linecolor", fontsize=12.0) if output is None: plt.show() else: plt.savefig(output, bbox_inches="tight") print(" [DONE]") return fig
[docs] @typechecked def plot_model_spectra( tag: str, n_spectra: Optional[int] = None, flux_offset: Optional[float] = None, label_pos: Optional[Tuple[float, float]] = None, xlim: Optional[Tuple[float, float]] = None, ylim: Optional[Tuple[float, float]] = None, title: Optional[str] = None, offset: Optional[Tuple[float, float]] = None, figsize: Optional[Tuple[float, float]] = (4.0, 2.5), output: Optional[str] = None, leg_param: Optional[List[str]] = None, ) -> mpl.figure.Figure: """ Function for plotting the results from comparing a spectrum with a grid of model spectra. Parameters ---------- tag : str Database tag where the results from the model comparison with :func:`~species.fit.compare_spectra.CompareSpectra.compare_model` are stored. n_spectra : int, None The number of spectra with the lowest goodness-of-fit statistic that will be plotted in comparison with the data. All spectra are selected if the argument is set to ``None``. label_pos : tuple(float, float), None Position for the name labels. Should be provided as (x, y) for the lowest spectrum. The ``flux_offset`` will be applied to the remaining spectra. The labels are only plotted if the argument of both ``label_pos`` and ``flux_offset`` are not ``None``. flux_offset : float, None Offset to be applied such that the spectra do not overlap. No offset is applied if the argument is set to ``None``. xlim : tuple(float, float), None Limits of the wavelength axis. ylim : tuple(float, float), None Limits of the flux axis. title : str Plot title. offset : tuple(float, float) Offset for the label of the x- and y-axis. figsize : tuple(float, float) Figure size. output : str, None Output filename for the plot. The plot is shown in an interface window if the argument is set to ``None``. leg_param : list(str), None List with the parameters to include in the legend of the model spectra. Apart from atmospheric parameters (e.g. 'teff', 'logg', 'radius') also parameters such as 'mass' and 'luminosity' can be included. The default atmospheric parameters are included in the legend if the argument is set to ``None``. Returns ------- matplotlib.figure.Figure The ``Figure`` object that can be used for further customization of the plot. """ if output is None: print("Plotting model spectra comparison...", end="") else: print(f"Plotting model spectra comparison: {output}...", end="") if flux_offset is None: flux_offset = 0.0 config_file = os.path.join(os.getcwd(), "species_config.ini") config = configparser.ConfigParser() config.read(config_file) db_path = config["species"]["database"] with h5py.File(db_path, "r") as hdf5_file: dset = hdf5_file[f"results/comparison/{tag}/goodness_of_fit"] object_name = dset.attrs["object_name"] n_spec_name = dset.attrs["n_spec_name"] model_name = dset.attrs["model"] n_param = dset.attrs["n_param"] n_scale_spec = dset.attrs["n_scale_spec"] parallax = dset.attrs["parallax"] n_inc_phot = dset.attrs["n_inc_phot"] spec_name = [] for i in range(n_spec_name): spec_name.append(dset.attrs[f"spec_name{i}"]) inc_phot = [] for i in range(n_inc_phot): inc_phot.append(dset.attrs[f"inc_phot{i}"]) model_param = [] coord_points = [] for i in range(n_param): model_param.append(dset.attrs[f"parameter{i}"]) coord_points.append( np.array(hdf5_file[f"results/comparison/{tag}/coord_points{i}"]) ) scale_spec = [] for i in range(n_scale_spec): scale_spec.append(dset.attrs[f"scale_spec{i}"]) goodness_fit = np.array(dset) goodness_fit = np.log10(goodness_fit) goodness_fit -= np.nanmin(goodness_fit) sort_idx = np.unravel_index( np.argsort(goodness_fit, axis=None), goodness_fit.shape ) sort_idx = list(sort_idx) for i, item in enumerate(sort_idx): sort_idx[i] = item[:n_spectra] flux_scaling = np.array(hdf5_file[f"results/comparison/{tag}/flux_scaling"]) if n_scale_spec > 0: extra_scaling = np.array( hdf5_file[f"results/comparison/{tag}/extra_scaling"] ) else: extra_scaling = None radius = ( np.sqrt(flux_scaling * (constants.PARSEC / (1e-3 * parallax)) ** 2) / constants.R_JUP ) plt.rcParams["font.family"] = "serif" plt.rcParams["mathtext.fontset"] = "dejavuserif" plt.rcParams["axes.axisbelow"] = False fig = plt.figure(figsize=figsize) gridsp = mpl.gridspec.GridSpec(1, 1) gridsp.update(wspace=0, hspace=0, left=0, right=1, bottom=0, top=1) ax = plt.subplot(gridsp[0, 0]) ax.tick_params( axis="both", which="major", colors="black", labelcolor="black", direction="in", width=1, length=5, labelsize=12, top=True, bottom=True, left=True, right=True, ) ax.tick_params( axis="both", which="minor", colors="black", labelcolor="black", direction="in", width=1, length=3, labelsize=12, top=True, bottom=True, left=True, right=True, ) ax.xaxis.set_minor_locator(AutoMinorLocator(5)) ax.yaxis.set_minor_locator(AutoMinorLocator(5)) ax.set_xlabel("Wavelength (\N{GREEK SMALL LETTER MU}m)", fontsize=13) if flux_offset == 0.0: ax.set_ylabel( r"$\mathregular{F}_\lambda$" + " (W m$^{-2}$ \N{GREEK SMALL LETTER MU}m$^{-1}$)", fontsize=11, ) else: ax.set_ylabel( r"$\mathregular{F}_\lambda$ (W m$^{-2}$" + " \N{GREEK SMALL LETTER MU}m$^{-1}$) + offset", fontsize=11, ) if xlim is not None: ax.set_xlim(xlim[0], xlim[1]) if ylim is not None: ax.set_ylim(ylim[0], ylim[1]) if offset is not None: ax.get_xaxis().set_label_coords(0.5, offset[0]) ax.get_yaxis().set_label_coords(offset[1], 0.5) else: ax.get_xaxis().set_label_coords(0.5, -0.1) ax.get_yaxis().set_label_coords(-0.1, 0.5) if title is not None: ax.set_title(title, y=1.02, fontsize=13) read_obj = ReadObject(object_name) obj_spec = read_obj.get_spectrum() if flux_offset == 0.0: for spec_key, spec_item in obj_spec.items(): if spec_key not in spec_name: continue ax.plot(spec_item[0][:, 0], spec_item[0][:, 1], "-", lw=0.5, color="black") model_reader = ReadModel(model_name) for i in range(n_spectra): param_select = {"parallax": parallax} idx_select = [] for param_idx, param_item in enumerate(model_param): param_select[param_item] = coord_points[param_idx][sort_idx[param_idx][i]] idx_select.append(sort_idx[param_idx][i]) param_select["radius"] = radius[tuple(idx_select)] if flux_offset != 0.0: for spec_key, spec_item in obj_spec.items(): if spec_key not in spec_name: continue if spec_key in scale_spec: spec_idx = scale_spec.index(spec_key) scaling_idx = np.append(idx_select, spec_idx) data_scaling = extra_scaling[tuple(scaling_idx)] else: data_scaling = 1.0 ax.plot( spec_item[0][:, 0], (n_spectra - i - 1) * flux_offset + data_scaling * spec_item[0][:, 1], "-", lw=0.5, color="black", ) if n_inc_phot > 0 and n_spec_name == 0: model_box = model_reader.get_data( model_param=param_select, spec_res=100.0, ) ax.plot( model_box.wavelength, (n_spectra - i - 1) * flux_offset + model_box.flux, color="tomato", lw=0.5, ) for spec_key, spec_item in obj_spec.items(): if spec_key not in spec_name: continue model_box = model_reader.get_data( model_param=param_select, spec_res=spec_item[3], wavel_resample=spec_item[0][:, 0], ) ax.plot( model_box.wavelength, (n_spectra - i - 1) * flux_offset + model_box.flux, color="tomato", lw=0.5, ) if label_pos is not None and flux_offset != 0.0: if leg_param is None: leg_param = model_param.copy() leg_param.append("radius") label_text = create_model_label( model_param=param_select, model_name=model_name, inc_model_name=False, object_type="planet", leg_param=leg_param, ) label_text = ( rf"$\Delta\log\,G_k = " rf"{goodness_fit[tuple(idx_select)]:.2f}$: {label_text}" ) ax.text( label_pos[0], label_pos[1] + (n_spectra - i - 1) * flux_offset, label_text, fontsize=8.0, ha="left", ) if output is None: plt.show() else: plt.savefig(output, bbox_inches="tight") print(" [DONE]") return fig