#!/usr/bin/env python3
# -*- coding: utf-8 -*-
from __future__ import annotations
from typing import Literal, Optional
from numpy.typing import ArrayLike
from copy import copy
from astropy.constants import c
from astropy import units as u
import math
import matplotlib.pyplot as plt
import numpy as np
from matplotlib import gridspec, rcParams
import matplotlib as mpl
from matplotlib.ticker import StrMethodFormatter
from q3dfit.contfit import readcf
from q3dfit.exceptions import InitializationError
from q3dfit.q3din import q3din
from q3dfit.spectConvol import spectConvol
from q3dfit.q3dout import q3dout
from q3dfit.q3dutil import lmlabel
[docs]
def plotcont(q3do: q3dout,
savefig: bool=False,
outfile: Optional[str]=None,
argssavefig: dict={'bbox_inches': 'tight', 'dpi': 300},
questfit: bool=False,
q3di: Optional[q3din]=None,
compspec: Optional[np.ndarray]=None,
complabs: Optional[list]=None,
compcols: Optional[list]=None,
title: Optional[str]=None,
xran: Optional[list]=None,
yran: Optional[list]=None,
zeroymin: bool=False,
xstyle: Literal['log', 'lin']='lin',
ystyle: Literal['log', 'lin']='lin',
figsize: tuple=(10, 5),
waveunit_in: Literal['micron','Angstrom']='micron',
waveunit_out: Optional[str]=None,
fluxunit_out: Literal['flambda', 'lambdaflambda', 'fnu'] = 'flambda',
mode: Literal['dark', 'light'] = 'dark'
):
'''
Created on Tue Jun 1 13:32:37 2021
@author: annamurphree
Plots continuum fit of optical data (fit by fitqsohost or ppxf)
or IR data (fit by questfit).
Parameters
----------
q3do
:py:class:`~q3dfit.q3dout.q3dout` object containing results of fit.
savefig
Optional. If True, saves the plot to a file. Defaults to False.
outfile
Optional. Full path and name of output plot. Defaults to None, which
means no output file is created.
argssavefig
Optional. Dictionary of arguments to pass to
:py:meth:`~matplotlib.pyplt.savefig()`. Defaults to
{'bbox_inches': 'tight', 'dpi': 300}.
questfit
Optional. If True, indicates that the fit is an IR spectrum fit by
questfit. Defaults to False.
q3di
Optional. :py:class:`~q3dfit.q3din.q3din` object containing fit initialization.
Only needed for IR spectra fit by questfit. Defaults to None.
compspec
Optional. Array of component spectra to overplot on continuum data.
If None, no components are overplotted.
complabs
Optional. List of labels for component spectra. Defaults to None, which
means labels are automatically generated as 'Component 1', 'Component 2', etc.
compcols
Optional. List of colors for component spectra. Defaults to None, which
means colors are automatically generated as 'C0', 'C1', etc.
title
Optional. Title for the plot. Defaults to None.
xran
Optional. Range of x-axis (wavelength) to plot. If None, applies the
`fitrange` attribute of :py:class:`~q3dfit.q3dout.q3dout`.
yran
Optional. Range of y-axis (flux) to plot. If None, the y-axis limits
are determined automatically based on the data.
zeroymin
Optional. If True, sets the minimum y-axis value to 0. Defaults to False.
xstyle
Optional. Style of x-axis scale, either 'log' or 'lin'. Defaults to 'lin'.
ystyle
Optional. Style of y-axis scale, either 'log' or 'lin'. Defaults to 'lin'.
figsize
Optional. Size of the figure in inches, specified as a tuple (width, height).
Defaults to (10, 5).
mode
Optional. Plot style, either 'dark' or 'light'. Defaults to 'dark'.
This affects the background color and text color of the plot.
waveunit_in
Optional. Input wavelength unit, either 'micron' or 'Angstrom'.
Defaults to 'micron'.
waveunit_out
Optional. Output wavelength unit, either 'micron' or 'Angstrom'.
If None, defaults to waveunit_in.
fluxunit_out
Optional. Output flux unit, either 'flambda', 'lambdaflambda', or 'fnu'.
Defaults to 'flambda'. The input flux unit is always 'flambda'. Choosing
'lambdaflambda' multiplies the flux by wavelength (in the units of
waveunit_out), and 'fnu' multiplies the flux by wavelength^2/c,
where c is the speed of light in the units of waveunit_out.
'''
rcParamsOrig = rcParams.copy()
# dark mode just for fun:
if mode == 'dark':
pltstyle = 'dark_background'
dcolor = 'w'
else:
pltstyle = 'seaborn-v0_8-ticks'
dcolor = 'k'
wave = q3do.wave.copy()
specstars = copy(q3do.cont_dat)
modstars = copy(q3do.cont_fit)
if waveunit_out is None:
# if no output wavelength unit is specified, use the input wavelength unit
waveunit_out = waveunit_in
wavein = wave.copy()*getattr(u,waveunit_in)
waveout = wavein.copy().to(waveunit_out)
# for optical spectra fit by fitqsohost or ppxf:
if not questfit:
if compspec is not None:
ccompspec = copy(compspec)
if len(ccompspec) > 1:
ncomp = len(ccompspec)
else:
ncomp = 1
if compcols is None:
compcols = ['C' + str(i) for i in range(ncomp)]
if complabs is None:
complabs = ['Component ' + str(i + 1) for i in range(ncomp)]
else:
ncomp = 0
if xran is None:
xran = copy(q3do.fitrange)
if waveunit_in == 'Angstrom' and waveunit_out == 'micron':
# convert angstrom to microns
xran = list(np.divide(xran, 10**4))
elif waveunit_in == 'micron' and waveunit_out == 'Angstrom':
# convert microns to angstroms
xran = list(np.multiply(xran, 10**4))
if fluxunit_out == 'lambdaflambda':
# multiply the flux by wavelength
specstars = list(np.multiply(specstars, wave))
modstars = list(np.multiply(modstars, wave))
if ncomp > 0:
for i in range(0, ncomp):
ccompspec[i] = list(np.multiply(ccompspec[i], wave))
ytit = '$\lambda$F$_\lambda$'
elif fluxunit_out == 'fnu':
# multiply the flux by wavelength^2/c
specstars = \
list(np.multiply(specstars,
np.divide(np.multiply(wavein.value, waveout.value),
c.to(waveunit_out+'/s').value)))
modstars = \
list(np.multiply(modstars,
np.divide(np.multiply(wavein.value, waveout.value),
c.to(waveunit_out+'/s').value)))
if ncomp > 0:
for i in range(0, ncomp):
ccompspec[i] = \
list(np.multiply(ccompspec[i],
np.divide(np.multiply(wavein.value, waveout.value),
c.to(waveunit_out+'/s').value)))
ytit = 'F$_\u03BD$'
else:
ytit = 'F$_\lambda$'
# plot on a log scale:
if xstyle == 'log' or ystyle == 'log':
plt.style.use(pltstyle)
# CB: Otherwise the background becomes black and the axes ticks
# unreadable when saving the figure
if mode == 'light':
rcParams['savefig.facecolor'] = 'white'
fig = plt.figure(figsize=figsize)
# fig = plt.figure()
plt.axis('off') # so the subplots don't share a y-axis
fig.add_subplot(1, 1, 1)
ydat = specstars
ymod = modstars
# plotting
plt.xlim(xran)
if yran is not None:
plt.ylim(yran)
fig.axes[0].axis('off') # so the subplots don't share a y-axis
fig.axes[1].axis('off') # so the subplots don't share a y-axis
gs = fig.add_gridspec(4, 1)
ax1 = fig.add_subplot(gs[:3, :])
# ax1.legend(ncol=2)
if xstyle == 'log':
ax1.set_xscale('log')
# ax1.set_xticklabels([])
if ystyle == 'log':
ax1.set_yscale('log')
ax1.set_ylabel(ytit, fontsize=20)
#if title == 'QSO':
# ax1.set_ylim(10e-7)
# actually plotting
plt.plot(waveout.value, ydat, dcolor, linewidth=1)
plt.plot(waveout.value, ymod, 'r', linewidth=2, label='Total')
if ncomp > 0:
for i in range(0, ncomp):
plt.plot(waveout.value, ccompspec[i], compcols[i], linewidth=2,
label=complabs[i])
# tick formatting
yticks_used = ax1.get_yticks()
ylim_used = ax1.get_ylim()
yticks_used = np.append(np.append(ylim_used[0], yticks_used),
ylim_used[1])
ax1.set_yticks(yticks_used)
ax1.set_ylim(ylim_used)
ax1.minorticks_on()
#ax1.tick_params(which='major', length=20, pad=10, labelsize=10)
#ax1.tick_params(which='minor', length=7, labelsize=8)
l = ax1.legend(loc='upper right', fontsize=16)
for text in l.get_texts():
text.set_color(dcolor)
ax2 = fig.add_subplot(gs[-1, :], sharex=ax1)
ax2.plot(waveout.value, np.divide(specstars, modstars), color=dcolor)
ax2.axhline(1, color='grey', linestyle='--', alpha=0.7, zorder=0)
ax2.set_ylabel('Data/Model', fontsize=19)
# ax2.tick_params(which='major', length=20, pad=20, labelsize=9)
# ax2.tick_params(which='minor', length=7, labelsize=8)
if waveunit_out == 'micron':
ax2.set_xlabel('Wavelength ($\mu$m)', fontsize=20)
elif waveunit_out == 'Angstrom':
ax2.set_xlabel('Wavelength ($\AA$)', fontsize=20)
gs.update(wspace=0.0, hspace=0.05)
plt.gcf().subplots_adjust(bottom=0.1)
if title is not None:
plt.suptitle(title, fontsize=20)
if savefig and outfile is not None:
plt.savefig(outfile[0], **argssavefig)
elif xstyle == 'lin' or ystyle == 'lin':
dxran = xran[1] - xran[0]
xran1 = [xran[0], xran[0] + np.around(dxran/3.0, 3)]
xran2 = [xran[0] + np.around(dxran/3.0, 3),
xran[0] + 2.0 * np.around(dxran/3.0, 3)]
xran3 = [xran[0] + 2.0 * np.around(dxran/3.0, 3),
xran[1]]
i1 = [None]
i2 = [None]
i3 = [None]
i1.pop(0)
i2.pop(0)
i3.pop(0)
ydat = specstars
ymod = modstars
for i in range(0, len(waveout.value)):
if waveout.value[i] > xran1[0] and waveout.value[i] < xran1[1]:
i1.append(i)
if waveout.value[i] > xran2[0] and waveout.value[i] < xran2[1]:
i2.append(i)
if waveout.value[i] > xran3[0] and waveout.value[i] < xran3[1]:
i3.append(i)
maxthresh = 0.2
ntop = 20
nbottom = 20
if len(waveout.value) < 100:
ntop = 10
nbottom = 10
++ntop
--nbottom
if waveunit_out == 'micron':
xtit = 'Observed Wavelength ($\mu$m)'
elif waveunit_out == 'Angstrom':
xtit = 'Observed Wavelength ($\AA$)'
plt.style.use(pltstyle)
fig = plt.figure(figsize=figsize)
plt.axis('off') # so the subplots don't share a y-axis
maximum = 0
minimum = 0
''
idict = {1: i1, 2: i2, 3: i3}
xrans = {1: xran1, 2: xran2, 3: xran3}
for group in range(1, 4):
if len(idict[group]) > 0:
fig.add_subplot(3, 1, group)
# finding min/max values at indices from idict
dat_et_mod = np.concatenate((ydat[idict[group]],
ymod[idict[group]]))
maximum = np.nanmax(dat_et_mod)
minimum = np.nanmin(dat_et_mod)
# set min and max in yran
yranpan = None
if yran is None:
yranpan = [minimum, maximum]
# finding yran[1] aka max
ydi = np.zeros(len(idict[group]))
ydi = np.array(ydat)[idict[group]]
ymodi = np.zeros(len(idict[group]))
ymodi = np.array(ymod)[idict[group]]
y = np.array(ydi - ymodi)
ny = len(y)
iysort = np.argsort(y)
ysort = np.array(y)[iysort]
ymodisort = ymodi[iysort]
if ysort[ny - ntop] < ysort[ny - 1] * maxthresh:
yranpan[1] = np.nanmax(ysort[0:ny - ntop] +
ymodisort[0:ny - ntop])
else:
yranpan = copy(yran)
if zeroymin:
yranpan[0] = 0.
# plotting
plt.xlim(xrans[group][0], xrans[group][1])
plt.ylim(yranpan)
plt.ylabel(ytit, fontsize=15)
if group == 3:
plt.xlabel(xtit, fontsize=15, labelpad=10)
if ystyle == 'log':
plt.yscale('log')
# tick formatting
plt.minorticks_on()
plt.tick_params(which='major', length=10, pad=5)
plt.tick_params(which='minor', length=5)
#if waveunit_out == 'micron':
# xticks = np.arange(np.around(xrans[group][0],1)-0.025,
# np.around(xrans[group][1],1), 0.025)[:-1]
# plt.xticks(xticks, fontsize=10)
#elif waveunit_out == 'Angstrom':
# xticks = np.arange(math.floor(xrans[group][0]/100.0)*100,
# (math.floor(xrans[group][1]/100)*100)+100, 100)
# plt.xticks(xticks, fontsize=10)
if np.nanmin(ydat) > 1e-10:
# this will fail if fluxes are very low (<~1e-10)
plt.yticks(np.arange(yranpan[0], yranpan[1],
np.around((yranpan[1] - yranpan[0])/5.,
decimals=2)), fontsize=10)
else:
plt.yticks()
# actually plotting
plt.plot(waveout.value, ydat, dcolor, linewidth=1)
if ncomp > 0:
for i in range(0, ncomp):
plt.plot(waveout.value, ccompspec[i], compcols[i],
linewidth=2, label=complabs[i])
plt.plot(waveout.value, ymod, 'r', linewidth=2, label=title)
if group == 1:
plt.legend(loc='upper right')
# more formatting
plt.subplots_adjust(hspace=0.25)
#plt.tight_layout(pad=5)
#plt.gcf().subplots_adjust(bottom=0.1)
if title is not None:
plt.suptitle(title, fontsize=20)
if savefig and outfile is not None:
if len(outfile[0])>1:
plt.savefig(outfile[0], **argssavefig)
else:
plt.savefig(outfile, **argssavefig)
plt.show()
# for IR spectra fit with questfit:
else:
comp_best_fit = q3do.ct_coeff['comp_best_fit']
if xstyle == 'log' or ystyle == 'log':
fig = plt.figure(figsize=figsize)
gs = fig.add_gridspec(4,1)
ax1 = fig.add_subplot(gs[:3, :])
MIRgdlambda = wave #[q3do.ct_indx]
MIRgdflux = q3do.spec #[q3do.ct_indx]
MIRcontinuum = modstars #[q3do.ct_indx]
if waveunit_in =='micron' and waveunit_out == 'Angstrom':
# convert microns to angstroms
MIRgdlambda = list(np.multiply(MIRgdlambda, 10**4))
elif waveunit_in =='Angstrom' and waveunit_out == 'micron':
# convert angstroms to microns
MIRgdlambda = list(np.divide(MIRgdlambda, 10**4))
if fluxunit_out == 'lambdaflambda':
# multiply the flux by wavelength
MIRgdflux = list(np.multiply(MIRgdflux, MIRgdlambda))
MIRcontinuum = list(np.multiply(MIRcontinuum, MIRgdlambda))
if len(comp_best_fit.keys()) > 0:
for i in range(0, len(comp_best_fit.keys())):
comp_best_fit[list(comp_best_fit.keys())[i]] = \
np.multiply(comp_best_fit[list(comp_best_fit.keys())[i]],
MIRgdlambda)
ytit = '$\lambda$F$_\lambda$'
elif fluxunit_out == 'fnu':
# multiply the flux by wavelength^2/c
MIRgdflux = list(np.multiply(MIRgdflux, MIRgdlambda))
MIRcontinuum = list(np.multiply(MIRgdflux, MIRgdlambda))
ytit = 'F$_\u03BD$'
else:
ytit = 'F$_\lambda$'
plt.style.use(pltstyle)
ax1.plot(MIRgdlambda, MIRgdflux, label='Data',color=dcolor)
ax1.plot(MIRgdlambda, MIRcontinuum, label='Model', color='r')
if 'global_ext_model' in q3di.argscontfit:
for i in np.arange(0,len(comp_best_fit.keys())-2,1):
ax1.plot(MIRgdlambda,
np.multiply(comp_best_fit[list(comp_best_fit.keys())[i]],
np.multiply(comp_best_fit[list(comp_best_fit.keys())[-2]],
comp_best_fit[list(comp_best_fit.keys())[-1]])),
label=list(comp_best_fit.keys())[i],
linestyle='--',alpha=0.5)
else:
for i in np.arange(0,len(comp_best_fit.keys()),3):
ax1.plot(MIRgdlambda,
np.multiply(comp_best_fit[list(comp_best_fit.keys())[i]],
np.multiply(comp_best_fit[list(comp_best_fit.keys())[i+1]],
comp_best_fit[list(comp_best_fit.keys())[i+2]])),
label=list(comp_best_fit.keys())[i],
linestyle='--',alpha=0.5)
for comp_i in comp_best_fit.keys():
if 'ext' not in comp_i and 'abs' not in comp_i:
spec_out = comp_best_fit[comp_i]
if comp_i+'_ext' in comp_best_fit.keys():
spec_out *= comp_best_fit[comp_i+'_ext']
if comp_i+'_abs' in comp_best_fit.keys():
spec_out *= comp_best_fit[comp_i+'_abs']
plt.plot(MIRgdlambda, spec_out, label=comp_i,linestyle='--',alpha=0.5)
#ax1.legend(ncol=2)
ax1.legend(loc='upper right',bbox_to_anchor=(1.15, 1),prop={'size': 10})
if xstyle == 'log':
ax1.set_xscale('log')
if ystyle == 'log':
ax1.set_yscale('log')
ax1.set_ylim(1e-4)
ax1.set_ylabel(ytit, fontsize=12)
ax2 = fig.add_subplot(gs[-1, :], sharex=ax1)
ax2.plot(MIRgdlambda,np.divide(MIRgdflux,MIRcontinuum),color=dcolor)
ax2.axhline(1, color='grey', linestyle='--', alpha=0.7, zorder=0)
ax2.set_ylabel('Data/Model', fontsize=12)
if waveunit_out == 'Angstrom':
ax2.set_xlabel('Wavelength ($\AA$)', fontsize=12)
elif waveunit_out == 'micron':
ax2.set_xlabel('Wavelength ($\mu$m)', fontsize=12)
gs.update(wspace=0.0, hspace=0.05)
plt.suptitle('Total', fontsize=30)
elif xstyle == 'lin' or ystyle == 'lin':
if xran is None:
xran = q3do.fitrange
MIRgdlambda = wave #[q3do.ct_indx]
MIRgdflux = q3do.spec #[q3do.ct_indx]
MIRcontinuum = modstars #[q3do.ct_indx]
xtit = ''
if waveunit_in == 'microns' and waveunit_out == 'Angstrom':
# convert wave list from microns to angstroms
MIRgdlambda = list(np.multiply(MIRgdlambda, 10**4))
xtit = 'Observed Wavelength ($\AA$)'
elif waveunit_in == 'Angstrom' and waveunit_out == 'micron':
# convert wave list from angstroms to microns
MIRgdlambda = list(np.divide(MIRgdlambda, 10**4))
xtit = 'Observed Wavelength ($\mu$m)'
if fluxunit_out == 'lambdaflambda':
# multiply the flux by wavelength
MIRgdflux = list(np.multiply(MIRgdflux, MIRgdlambda))
MIRcontinuum = list(np.multiply(MIRcontinuum, MIRgdlambda))
if len(comp_best_fit.keys()) > 0:
for i in range(0, len(comp_best_fit.keys())):
comp_best_fit[list(comp_best_fit.keys())[i]] = \
list(np.multiply(comp_best_fit[list(comp_best_fit.keys())[i]],
MIRgdlambda))
ytit = '$\lambda$F$_\lambda$'
elif fluxunit_out == 'fnu':
# multiply the flux by wavelength
MIRgdflux = list(np.multiply(MIRgdflux, MIRgdlambda))
MIRcontinuum = list(np.multiply(MIRgdflux, MIRgdlambda))
ytit = 'F$_\u03BD$'
else:
ytit = 'F$_\lambda$'
wave = MIRgdlambda
ydat = MIRgdflux
ymod = MIRcontinuum
dxran = xran[1] - xran[0]
xran1 = [xran[0], xran[0] + np.around(dxran/3.0,3)]
xran2 = [xran[0] + np.around(dxran/3.0,3), xran[0] + 2.0 * np.around(dxran/3.0,3)]
xran3 = [xran[0] + 2.0 * np.around(dxran/3.0,3), xran[1]]
i1 = [None]
i2 = [None]
i3 = [None]
i1.pop(0)
i2.pop(0)
i3.pop(0)
for i in range(0, len(wave)):
if wave[i] > xran1[0] and wave[i] < xran1[1]:
i1.append(i)
if wave[i] > xran2[0] and wave[i] < xran2[1]:
i2.append(i)
if wave[i] > xran3[0] and wave[i] < xran3[1]:
i3.append(i)
maxthresh = 0.2
ntop = 20
nbottom = 20
if len(wave) < 100:
ntop = 10
nbottom = 10
++ntop
--nbottom
plt.style.use(pltstyle)
fig = plt.figure(figsize=figsize)
#fig = plt.figure()
plt.axis('off') # so the subplots don't share a y-axis
maximum = 0
minimum = 0
idict = {1:i1, 2:i2, 3:i3}
xrans = {1:xran1, 2:xran2, 3:xran3}
for group in range(1,4):
if len(idict[group]) > 0:
fig.add_subplot(3, 1, group)
ax = plt.subplot(3, 1, group)
# shrink current axis by 10% to fit legend on side
box = ax.get_position()
ax.set_position([box.x0, box.y0, box.width * 0.8, box.height])
# finding max value between ydat and ymod at indices from i1
for i in idict[group]:
bigboy = np.nanmax([ydat[i], ymod[i]])
if bigboy > maximum:
maximum = bigboy
# finding min
for i in idict[group]:
smallboy = np.nanmin([ydat[i], ymod[i]])
if smallboy < minimum:
minimum = smallboy
# set min and max in yran
yranpan = None
if yran is None:
yranpan = [minimum, maximum]
# finding yran[1] aka max
ydi = np.zeros(len(idict[group]))
ydi = np.array(ydat)[idict[group]]
ymodi = np.zeros(len(idict[group]))
ymodi = np.array(ymod)[idict[group]]
y = np.array(ydi - ymodi)
ny = len(y)
iysort = np.argsort(y)
ysort = np.array(y)[iysort]
ymodisort = ymodi[iysort]
if ysort[ny - ntop] < ysort[ny - 1] * maxthresh:
yranpan[1] = np.nanmax(ysort[0:ny - ntop] +
ymodisort[0:ny - ntop])
else:
yranpan = copy(yran)
if zeroymin:
yranpan[0] = 0.
# plotting
plt.xlim(xrans[group][0], xrans[group][1])
plt.ylim(yranpan)
plt.ylabel(ytit, fontsize=15)
if group == 3:
plt.xlabel(xtit, fontsize=15, labelpad=10)
if ystyle == 'log':
plt.yscale('log')
# tick formatting
plt.minorticks_on()
plt.tick_params(which='major', length=10, pad=5)
plt.tick_params(which='minor', length=5)
if waveunit_out == 'micron':
xticks = np.arange(np.around(xrans[group][0]), np.around(xrans[group][1]), 1)
plt.xticks(xticks, fontsize=10)
elif waveunit_out == 'Angstrom':
xticks = np.arange(math.floor(xrans[group][0]/1000.0)*1000,
(math.floor(xrans[group][1]/1000.0)*1000)+1000, 10000)
plt.xticks(xticks, fontsize=10)
if fluxunit_out != 'fnu':
# this will fail if fluxes are very low (<~1e-10)
plt.yticks(np.arange(yranpan[0], yranpan[1],
np.around((yranpan[1] - yranpan[0])/5.,
decimals=2)),
fontsize=10)
else:
plt.yticks()
# actually plotting
plt.plot(MIRgdlambda, MIRgdflux, label='Data',
color=dcolor)
plt.plot(MIRgdlambda, MIRcontinuum, label='Model',
color='red')
if 'global_ext_model' in q3di.argscontfit:
for i in np.arange(0,len(comp_best_fit.keys())-2,1):
plt.plot(MIRgdlambda,
np.multiply(comp_best_fit[list(comp_best_fit.keys())[i]],
np.multiply(comp_best_fit[list(comp_best_fit.keys())[-2]],
comp_best_fit[list(comp_best_fit.keys())[-1]])),
label=list(comp_best_fit.keys())[i],linestyle='--',alpha=0.5)
else:
for comp_i in comp_best_fit.keys():
if 'ext' not in comp_i and 'abs' not in comp_i:
spec_out = comp_best_fit[comp_i]
if comp_i+'_ext' in comp_best_fit.keys():
spec_out *= comp_best_fit[comp_i+'_ext']
if comp_i+'_abs' in comp_best_fit.keys():
spec_out *= comp_best_fit[comp_i+'_abs']
plt.plot(MIRgdlambda, spec_out, label=comp_i,linestyle='--',alpha=0.5)
if group == 1:
ax.legend(loc='upper right',bbox_to_anchor=(1.22, 1),prop={'size': 10})
# more formatting
plt.subplots_adjust(hspace=0.25)
plt.tight_layout(pad=15)
plt.gcf().subplots_adjust(bottom=0.1)
plt.gcf().subplots_adjust(right=0.85)
if title is not None:
plt.suptitle(title, fontsize=20)
if savefig and outfile is not None:
if len(outfile[0])>1:
plt.savefig(outfile[0], **argssavefig)
else:
plt.savefig(outfile, **argssavefig)
plt.show()
rcParams.update(rcParamsOrig)
[docs]
def plotline(q3do: q3dout,
nx: int=1,
ny: int=1,
figsize: tuple=(16,13),
line: Optional[ArrayLike]=None,
center_obs: Optional[ArrayLike]=None,
center_rest: Optional[ArrayLike]=None,
waveunit_in: Literal['micron','Angstrom']='micron',
waveunit_out: Optional[str]=None,
specConv: Optional[spectConvol]=None,
yran: Optional[list]=None,
zeroymin: bool=False,
size: float|ArrayLike=300.,
savefig: bool=False,
outfile: Optional[str]=None,
argssavefig: dict={'bbox_inches': 'tight', 'dpi': 300},
mode: Literal['dark', 'light'] = 'dark'):
"""
Plot emission line fit and output to JPG
Parameters
----------
q3do
:py:class:`~q3dfit.q3dout.q3dout` object containing the output of the
fit.
nx
Number of columns in the plot grid. Defaults to 1.
ny
Number of rows in the plot grid. Defaults to 1.
figsize
Size of the figure in inches, specified as a tuple (width, height).
line
Optional. List of lines in which to center the subplots. If None,
the center of the plot window is determined from center_obs or
center_rest.
center_obs
Optional. List of wavelengths in the observed frame to center the
subplots. If None, the center of the plot window is determined from
line or center_rest.
center_rest
Optional. List of wavelengths in the rest frame to center the
subplots. If None, the center of the plot window is determined from
line or center_obs.
waveunit_in
Optional. Input wavelength unit, either 'micron' or 'Angstrom'.
Defaults to 'micron'.
waveunit_out
Optional. Output wavelength unit, either 'micron' or 'Angstrom'.
If None, defaults to waveunit_in.
specConv
Optional. :py:class:`~q3dfit.spectConvol.spectConvol` object for
spectral convolution. If None, no convolution is applied.
yran
Optional. Range of y-axis (flux) to plot. If None, the y-axis limits
are determined automatically based on the data.
zeroymin
Optional. If True, sets the minimum y-axis value to 0. Defaults to False.
size
Optional. Size of the plot window in Angstroms. If a single float is
savefig
Optional. If True, saves the plot to a file. Defaults to False.
outfile
Optional. Full path and name of output plot. Defaults to None, which
means no output file is created.
argssavefig
Optional. Dictionary of arguments to pass to
:py:meth:`~matplotlib.pyplt.savefig()`. Defaults to
{'bbox_inches': 'tight', 'dpi': 300}.
mode
Optional. Plotting mode, either 'dark' or 'light'. Defaults to 'dark'
Changes the plot background and text colors for better visibility.
"""
rcParamsOrig = rcParams.copy()
#setting mode parameters
if mode == 'dark':
pltstyle = 'dark_background'
dcolor = 'w'
elif mode == 'light':
pltstyle = 'seaborn-v0_8-ticks'
dcolor = 'k'
ncomp = q3do.maxncomp
colors = ['Magenta', 'Green', 'Orange', 'Teal']
if waveunit_out is None:
# if no output wavelength unit is specified, use the input wavelength unit
waveunit_out = waveunit_in
wave = q3do.wave.copy()
spectot = q3do.spec
specstars = q3do.cont_dat
modstars = q3do.cont_fit
modlines = q3do.line_fit
modtot = modstars + modlines
if waveunit_in == 'Angstrom' and waveunit_out == 'micron':
# convert angstrom to microns
wave = list(np.divide(wave, 10**4))
elif waveunit_in == 'micron' and waveunit_out == 'Angstrom':
# convert microns to angstroms
wave = list(np.multiply(wave, 10**4))
# To-do: Allow output wavelengths in Angstrom
#'waveunit_out' = 'micron'
# if 'waveunit_out' in pltpar:
# if pltpar['waveunit_out = 'Angstrom':
# waveunit_out = 'Angstrom'
# To-do: Get masking code from pltcont
# lines
linelist = q3do.linelist['lines']
linelabel = q3do.linelist['name']
linetext = q3do.linelist['linelab']
# Sort in wavelength order
isortlam = np.argsort(linelist)
linelist = linelist[isortlam]
linelabel = linelabel[isortlam]
linetext = linetext[isortlam]
#
# Plotting parameters
#
# Look for line list, then determine center of plot window from fitted
# wavelength
if line is not None:
sub_linlab = line
linwav = np.empty(len(sub_linlab), dtype='float32')
for i in range(0, len(sub_linlab)):
# Get wavelength from zeroth component
if sub_linlab[i] != '':
lmline = lmlabel(sub_linlab[i])
# if ncomp > 0
if f'{lmline.lmlabel}_0_cwv' in q3do.param.keys():
linwav[i] = q3do.param[f'{lmline.lmlabel}_0_cwv']
# otherwise
else:
idx = np.where(q3do.linelist['name'] == sub_linlab[i])
if len(idx) > 0:
linwav[i] = q3do.linelist['lines'][idx] * \
(1. + q3do.zstar)
else:
raise InitializationError(f'Line {sub_linlab[i]} not fit.')
else:
linwav[i] = 0.
# If linelist not present, get cwavelength enter of plot window from list
# first option: wavelength center specified in observed (plotted) frame
elif center_obs is not None:
linwav = np.array(center_obs)
# second option: wavelength center specified in rest frame, then converted
# to observed (plotted) frame
elif center_rest is not None:
linwav = np.array(center_rest) * q3do.zstar
else:
raise InitializationError('LINE, CENTER_OBS, or CENTER_REST ' +
'list not given in ARGSPLTLIN dictionary')
nlin = len(linwav)
# Size of plot in wavelength, in observed frame
# case of single size for all panels
if isinstance(size, float):
size = np.full(nlin, size) # default size currently 300 A ... fix for
# case of array of sizes
else:
size = np.array(size)
# other units!
off = np.array([-1.*size/2., size/2.])
off = off.transpose()
plt.style.use(pltstyle)
fig = plt.figure(figsize=figsize)
for i in range(0, nlin):
outer = gridspec.GridSpec(ny, nx, wspace=0.2, hspace=0.2)
inner = \
gridspec.GridSpecFromSubplotSpec(2, 1,
subplot_spec=outer[i],
wspace=0.1, hspace=0,
height_ratios=[4, 2],
width_ratios=None)
# create xran and ind
linwavtmp = linwav[i]
offtmp = off[i, :]
xran = linwavtmp + offtmp
ind = np.array([0])
for h in range(0, len(wave)):
if wave[h] > xran[0] and wave[h] < xran[1]:
ind = np.append(ind, h)
ind = np.delete(ind, [0])
ct = len(ind)
if ct > 0:
# create subplots
ax0 = plt.Subplot(fig, inner[0])
ax1 = plt.Subplot(fig, inner[1])
fig.add_subplot(ax0)
fig.add_subplot(ax1)
# create x-ticks
xticks = np.linspace(xran[0],xran[1],num=5,endpoint=False)
xticks = np.delete(xticks, [0])
if waveunit_out == 'Angstrom':
plt.gca().xaxis.set_major_formatter(StrMethodFormatter('{x:.0f}'))
elif waveunit_out == 'micron':
plt.gca().xaxis.set_major_formatter(StrMethodFormatter('{x:.2f}'))
# xticks = xticks * 1e4
# create minor x-ticks
xmticks = np.linspace(xran[0],xran[1],num=25,endpoint=False)
xmticks = np.delete(xmticks, [0])
#if waveunit_out == 'Angstrom':
# xmticks = xticks * 1e4
# set ticks
ax0.set_xticks(xticks)
ax1.set_xticks(xticks)
ax0.set_xticks(xmticks, minor=True)
ax1.set_xticks(xmticks, minor=True)
ax0.tick_params('x', which='major', direction='in', length=7,
width=2, color=dcolor)
ax0.tick_params('x', which='minor', direction='in', length=5,
width=1, color=dcolor)
ax1.tick_params('x', which='major', direction='in', length=7,
width=2, color=dcolor)
ax1.tick_params('x', which='minor', direction='in', length=5,
width=1, color=dcolor)
# create yran
ydat = spectot
ymod = modtot
ydattmp = np.zeros((ct), dtype=float)
ymodtmp = np.zeros((ct), dtype=float)
for j in range(0, len(ind)):
ydattmp[j] = ydat[(ind[j])]
ymodtmp[j] = ymod[(ind[j])]
ydatmin = min(ydattmp)
ymodmin = min(ymodtmp)
if ydatmin <= ymodmin:
yranmin = ydatmin
else:
yranmin = ymodmin
ydatmax = max(ydattmp)
ymodmax = max(ymodtmp)
if ydatmax >= ymodmax:
yranmax = ydatmax
else:
yranmax = ymodmax
if yran is None:
yranpan = [yranmin, yranmax]
if zeroymin:
yranpan[0]=0.
icol = (float(i))/(float(nx))
if icol % 1 == 0:
ytit = 'Fit'
else:
ytit = ''
ax0.set(ylabel=ytit)
ax0.set_xlim([xran[0], xran[1]])
ax0.set_ylim([yranpan[0], yranpan[1]])
# plots on ax0
ax0.plot(wave, ydat, color=dcolor, linewidth=1)
if waveunit_out == 'micron':
xtit = 'Wavelength ($\mu$m)'
elif waveunit_out == 'Angstrom':
xtit = 'Wavelength ($\AA$)'
ytit = ''
ax0.plot(wave, ymod, color='Red', linewidth=2)
# Plot all lines visible in plot range
for j in range(0, ncomp):
ylaboff = 0.07
for k, line in enumerate(linelabel):
lmline = lmlabel(line)
if f'{lmline.lmlabel}_{j}_cwv' in q3do.param.keys():
refwav = q3do.param[f'{lmline.lmlabel}_{j}_cwv']
else:
irefwav = np.where(q3do.linelist['name'] == line)
refwav = q3do.linelist['lines'][irefwav] * \
(1. + q3do.zstar)
if refwav >= xran[0] and refwav <= xran[1]:
if f'{lmline.lmlabel}_{j}_cwv' in \
q3do.param.keys():
flux = q3do.cmplin(line, j)
#import pdb; pdb.set_trace()
if specConv is not None:
conv = specConv.spect_convolver(wave, flux,
wavecen=refwav)
else:
conv = flux
ax0.plot(wave, yranpan[0] + conv, color=colors[j],
linewidth=2, linestyle='dashed')
ax0.annotate(linetext[k], (0.05, 1. - ylaboff),
xycoords='axes fraction',
va='center', fontsize=10)
ylaboff += 0.07
# if nmasked > 0:
# for r in range(0,nmasked):
# ax0.plot([masklam[r,0], masklam[r,1]], [yran[0], yran[0]],linewidth=8, color='Cyan')
# set new value for yran
ydat = specstars
ymod = modstars
ydattmp = np.zeros((len(ind)), dtype=float)
ymodtmp = np.zeros((len(ind)), dtype=float)
for j in range(0, len(ind)):
ydattmp[j] = ydat[(ind[j])]
ymodtmp[j] = ymod[(ind[j])]
ydatmin = min(ydattmp)
ymodmin = min(ymodtmp)
if ydatmin <= ymodmin:
yranmin = ydatmin
else:
yranmin = ymodmin
ydatmax = max(ydattmp)
ymodmax = max(ymodtmp)
if ydatmax >= ymodmax:
yranmax = ydatmax
else:
yranmax = ymodmax
yranpan = [yranmin, yranmax]
if icol % 1 == 0:
ytit = 'Residual'
else:
ytit = ''
ax1.set(ylabel=ytit)
# plots on ax1
ax1.set_xlim([xran[0], xran[1]])
ax1.set_ylim([yranpan[0], yranpan[1]])
ax1.plot(wave, ydat, linewidth=1)
ax1.plot(wave, ymod, color='Red')
# title
if waveunit_out == 'micron':
xtit = 'Wavelength ($\mu$m)'
elif waveunit_out == 'Angstrom':
xtit = 'Wavelength ($\AA$)'
fig.suptitle(xtit, fontsize=20)
if savefig and outfile is not None:
if len(outfile[0])>1:
fig.savefig(outfile[0], **argssavefig)
else:
fig.savefig(outfile, **argssavefig)
plt.show()
rcParams.update(rcParamsOrig)
[docs]
def plotcontcomponents(q3do: q3dout,
plottype: Literal['line', 'stackplot'] = 'line',
totals_plot: bool=True,
totals_plot_components: Optional[ArrayLike] = None,
totals_plot_labels: Optional[ArrayLike] = None,
totals_plot_colors: Optional[ArrayLike] = None,
components_plot: Optional[bool] = True,
min_weight: Optional[float] = 0,
sortvar: Optional[Literal['medians', 'ages', 'zs', 'weights']] = 'medians',
sortorder: Optional[bool] = True,
linalpha: Optional[float] = 0.8,
linwidth: Optional[float] = 1,
mode: Optional[Literal['dark', 'light']] = 'dark',
figsize: tuple = (12,12),
savefig: bool=False,
outfile: Optional[str]='none',
argssavefig: dict={'bbox_inches': 'tight', 'dpi': 300}):
"""
Plot stellar components of fit and output to JPG
Parameters
----------
q3do
:py:class:`~q3dfit.q3dout.q3dout` object containing the output of the
fit.
plottype
Select type of plot for components, currently only "line" and "stackplot", defaults to "line"
totals_plot
Optional. If True plots the component sum on a seperate plot, defaults to True
totals_plot_components
Optional. Additional data to plot on the totals plot.
totals_plot_labels
Optional. Labels for additional data to plot on the totals plot.
totals_plot_colors
Optional. Colors for additional data to plot on the totals plot.
components_plot
Optional. If True plots the stellar components, defaults to True
compcmap
Optional. Matplotlib colormap to use for stellar components plot, defaults to 'cividis'
min_weight
Optional. Minimum weight for a component to be plotted. Defaults to 0, which means all components are plotted.
sortvar
Optional. Variable by which to sort the stellar components, either 'medians', 'ages', 'zs', or 'weights'. Defaults to 'medians'.
sortorder
Optional. If True plots stellar components in decending order, if False plots in ascending order. Defaults to True
linalpha
Optional. Sets the alpha value of plotted lines.
linwidth
Optional. Sets the line width of plotted lines.
mode
set plotting mode, "light" or "dark", defaults to "dark"
figsize
Size of the figure in inches, specified as a tuple (width, height).
savefig
Optional. If True, saves the plot to a file. Defaults to False.
outfile
Optional. Full path and name of output plot. Defaults to None, which
means no output file is created.
argssavefig
Optional. Dictionary of arguments to pass to
:py:meth:`~matplotlib.pyplt.savefig()`. Defaults to
{'bbox_inches': 'tight', 'dpi': 300}.
mode
Optional. Plotting mode, either 'dark' or 'light'. Defaults to 'dark'
Changes the plot background and text colors for better visibility.
"""
rcParamsorig = rcParams.copy()
# Dark Mode!
wave = q3do.wave
if mode == 'dark':
pltstyle = 'dark_background'
dcolor = 'w'
elif mode == 'light':
pltstyle = 'seaborn-v0_8-ticks'
dcolor = 'k'
else:
raise ValueError("Invalid mode. Choose 'dark' or 'light'.")
# Configuring plot layout based on enabled subplots
plt.style.use(pltstyle)
if totals_plot and components_plot:
fig, ax = plt.subplots(2, 1, sharex=True, figsize=figsize)
elif components_plot and not totals_plot:
fig, ax = plt.subplots(1, 1, sharex=True, figsize=figsize)
ax = [ax]
elif totals_plot and not components_plot:
fig, ax = plt.subplots(1, 1, sharex=True, figsize=figsize)
ax = [None, ax]
else:
raise ValueError("At least one of totals_plot or components_plot must be True.")
# Initializing variables
ages = q3do.component_templates['age']
zs = q3do.component_templates['zs']
templates_list = q3do.component_templates['convolved']
weights = q3do.ct_coeff['stelweights']
indecies = q3do.component_templates['index']
if 'flux_fraction' in q3do.ct_coeff:
flux_percentage = [f' | Flux: {q3do.ct_coeff['flux_fraction'][i]*100:.2f}%' for i in indecies]
else:
flux_percentage = [''] * len(indecies)
if 'mass_fraction' in q3do.ct_coeff:
mass_percentage = [f' | Mass: {q3do.ct_coeff['mass_fraction'][i]*100:.2f}%' for i in indecies]
else:
mass_percentage = [''] * len(indecies)
# Create a list of relavent templates
templates_list = [templates_list[:, i] for i in range(len(indecies))]
labels_list = [f'Age: {(ages[i]/1e9):.2f} Gyr | Z: {zs[i]:.3f}{flux_percentage[i]}{mass_percentage[i]}' for i in range(len(indecies))]
# Plotting the compoonents plot
if components_plot:
# Sort templates by their mean value
if sortvar != None:
if sortvar == 'medians':
templates_medians = [np.mean(template) for template in templates_list]
sorted_indices = np.argsort(templates_medians)
if sortvar == 'ages':
sorted_indices = np.argsort(ages)
if sortvar == 'zs':
sorted_indices = np.argsort(zs)
if sortvar == 'weights':
sorted_indices = np.argsort(weights[indecies])
if sortvar == 'maxima':
sorted_indices = np.argsort([np.max(template) for template in templates_list])
if sortorder:
sorted_indices = sorted_indices[::-1] # Sort in descending order
templates_list = [templates_list[i] for i in sorted_indices]
labels_list = [labels_list[i] for i in sorted_indices]
# Generating colors for each component
#n_lines = len(templates_list)
#cmap = plt.get_cmap(compcmap)
#color = cmap(np.linspace(0, 1, n_lines))
# Plot each component in the upper panel
if plottype == 'line':
for i in range(len(templates_list)):
if weights[i] >= min_weight:
ax[0].plot(wave, templates_list[i], label=labels_list[i], alpha=linalpha, lw=linwidth)#, color=color[i])
if plottype == 'stackplot':
ax[0].stackplot(wave, templates_list, labels=labels_list, alpha=linalpha, lw=linwidth)# colors=color)
ax[0].legend()
ax[0].set_title('Stellar fit components')
ax[0].set_xlabel('Wavelength (Angstrom)')
ax[0].set_ylabel('Flux')
# Plotting the totals plot
# Finding the total fit.
if totals_plot:
tempSum = q3do.component_templates['convolved'].sum(1) + q3do.polymod
# Establish comppnents for totals plot
if totals_plot_components == None:
totals_plot_components = [q3do.cont_dat, q3do.cont_fit]
if q3do.add_poly_degree >= 0:
totals_plot_components += [q3do.polymod, q3do.stelmod]
if totals_plot_labels == None:
totals_plot_labels = ['Continuum Data', 'Continuum fit']
if q3do.add_poly_degree >= 0:
totals_plot_labels += [f'Ord. {q3do.add_poly_degree} Legendre poly', 'Sum of Convolved Templates']
if totals_plot_colors == None:
totals_plot_colors = [dcolor, 'red']
if q3do.add_poly_degree >= 0:
totals_plot_colors += ['magenta', 'cyan']
totals_plot_components.insert(-1, tempSum)
totals_plot_labels.insert(-1, 'Total Component Sum (Templates + Polynomial)')
totals_plot_colors.insert(-1, 'blue')
# Plot each component for totals plot
for i in range(len(totals_plot_components)):
ax[1].plot(q3do.wave, totals_plot_components[i], label=totals_plot_labels[i], lw=1, zorder=0+i, color=totals_plot_colors[i])
ax[1].legend()
ax[1].set_title('Total stellar fit')
ax[1].set_xlabel('Wavelength (Angstrom)')
ax[1].set_ylabel('Flux')
fig.suptitle('Stellar Population Component Decomposition', fontsize=16)
if savefig and outfile is not None:
if len(outfile[0])>1:
fig.savefig(outfile[0], **argssavefig)
else:
fig.savefig(outfile, **argssavefig)
plt.show()
rcParams.update(rcParamsorig)
[docs]
def plotpopheatmap(q3do: q3dout,
startempfile: str = None,
stelweights: ArrayLike = None,
mode: Optional[Literal['dark', 'light']] = 'dark',
savefig: Optional[bool] = False,
outfile: Optional[str] = None,
argssavefig: Optional[dict] = {'bbox_inches': 'tight', 'dpi': 300}):
'''
Plot a heatmap of the stellar population weights as a function of age and metallicity.
Parameters
----------
q3do
:py:class:`~q3dfit.q3dout.q3dout` object containing the output of the fit.
startempfile
Optional. Path to a .npz file containing the stellar population templates.
Used to load the ages and metallicities of the templates. If None, uses data saved in the q3do object.
stelweights
Optional. Array of stellar population weights. If None, uses the weights from the q3do object.
mode
Optional. Plotting mode, either 'dark' or 'light'. Defaults to 'dark'.
Changes the plot background and text colors for better visibility.
savefig
Boolean indicating whether to save the figure.
outfile
String specifying the path and filename for the saved figure.
argssavefig
Dictionary of arguments to pass to the savefig function.
'''
rcParamsorig = rcParams.copy()
if mode == 'dark':
pltstyle = 'dark_background'
dcolor = 'w'
elif mode == 'light':
pltstyle = 'seaborn-v0_8-ticks'
dcolor = 'k'
else:
raise ValueError("Invalid mode. Choose 'dark' or 'light'.")
plt.style.use(pltstyle)
if type(startempfile) != type(None):
templates = np.load(startempfile, allow_pickle=True)[()]
ages = np.asarray(templates['ages'])
zs = np.asarray(templates['zs'])
# Getting unique ages and metallicities from templates
unique_log_ages = np.unique(np.log10(ages))
unique_zs = np.unique(zs)
else:
# Getting unique ages and metallicities from q3do
unique_log_ages = np.log10(q3do.component_templates['unique_ages'])
unique_zs = q3do.component_templates['unique_zs']
if type(stelweights) != type(None):
# getting weights from stelweights
weights = np.asarray(stelweights)
else:
# getting weights from q3do
weights = np.asarray(q3do.ct_coeff['stelweights'])
n_ages = len(unique_log_ages)
n_metal = len(unique_zs)
# Sanity check and reshape
if weights.size != n_metal * n_ages:
raise ValueError(f"weights length {weights.size} != n_metal*n_ages ({n_metal}*{n_ages}={n_metal*n_ages})")
weights_grid = weights.reshape(n_metal, n_ages)
#weights_grid = weights_grid[:, -20:]
#unique_log_ages = unique_log_ages[-20:]
#nages = len(unique_log_ages)
# Plotting heatmap
fig, ax = plt.subplots(figsize=(10, 5))
im = ax.imshow(weights_grid, origin='lower', aspect='auto', cmap='viridis')
# set x-ticks at log-age intervals
tick_interval = 0.2
half_ints = np.arange(unique_log_ages[0], unique_log_ages[-1], tick_interval)
xtick_positions = []
xtick_labels = []
for val in half_ints:
idx = np.where(np.isclose(unique_log_ages, val))[0]
if idx.size:
xtick_positions.append(int(idx[0]))
xtick_labels.append(f"{val:.1f}")
ax.set_xticks(xtick_positions)
ax.set_xticklabels(xtick_labels, rotation=45)
ax.set_yticks(np.arange(n_metal))
ax.set_yticklabels([f"{z:.3f}" for z in unique_zs])
ax.set_xlabel('Stellar Age (log10 Age [yr])')
ax.set_ylabel('Metallicity (Z)')
ax.set_title('Stellar Population Heatmap')
fig.colorbar(im, ax=ax, label='Weight')
fig.tight_layout()
if savefig and outfile is not None:
if len(outfile[0])>1:
fig.savefig(outfile[0], **argssavefig)
else:
fig.savefig(outfile, **argssavefig)
fig.show()
rcParams.update(rcParamsorig)
def adjust_ax(ax,
fig,
fs=20,
minor=False):
'''
CB: Function defined to adjust the sizes of xlabel, ylabel, and the
ticklabels (in an inelegant way for the latter).
Presently just a utility function to be used in plotquest.
Further documentation pending more testing and development.
Parameters
-----
ax: matplotlib axis object
ax object of the plot you want to adjust
fig: matplotlib fig object
fig object that contains the ax object
returns
-------
Nothing
'''
fig.canvas.draw()
xlabel = ax.get_xlabel()
ylabel = ax.get_ylabel()
ax.set_xlabel(xlabel, fontsize=fs)
ax.set_ylabel(ylabel, fontsize=fs)
ax.tick_params(labelsize=fs-3)
# -- Trying to prune xtickslabels if increasing the fontsize made them overlap
xticks_old = ax.get_xticks()
if minor:
xticks_old = ax.get_xticks(minor=True)
xfigsize = fig.get_size_inches()[0] # in inches
textstrlen = len(ax.get_xticklabels()[0]._text.replace('\\mathdefault', '')) # length of tick labels depends on nr of decimals specified
textwidth_inch = textstrlen * (fs-3)*0.7 / 72. # Assume width of number in text = 0.7* height. Matplotlib uses 72 Points per inch (ppi): https://stackoverflow.com/questions/47633546/relationship-between-dpi-and-figure-size
if (len(xticks_old)+1)*textwidth_inch > 0.9* xfigsize * ax.get_position().width:
xticks_new = np.array([])
for i in range(len(xticks_old)):
if i%2==1:
xticks_new = np.append(xticks_new, xticks_old[i])
if not minor:
ax.set_xticks(xticks_new, fontsize=fs-3)
else:
ax.set_xticks(xticks_new, fontsize=fs-3, minor=True)
ax.set_xticklabels(ax.get_xticks(), fontsize=fs-3)
ax.tick_params(axis='x', which='both', labelsize=fs-3)
fig.tight_layout()
def plotdecomp(q3do,
q3di,
savefig=True,
outfile=None,
templ_mask=[],
do_lines=False,
show=False,
mode='light',
ymin=-1,
ymax=-1,
try_adjust_ax=True):
'''
Calls plotquest to plot the quasar-host galaxy decomposition. Not sure what
the difference is between this and plotquest.
Further documentation pending more testing and development.
'''
wave = q3do.wave
specstars = q3do.cont_dat
modstars = q3do.cont_fit
MIRgdlambda = wave
MIRgdflux = q3do.spec
MIRcontinuum = modstars
if outfile is None:
outfile=q3do.filelab + '_decomp'
if do_lines:
plotquest(q3do.wave, q3do.spec, q3do.cont_fit, q3do.ct_coeff, q3di, zstar=q3do.zstar, savefig=savefig, outfile=outfile,
templ_mask=templ_mask, lines=q3do.linelist['lines'], linespec=q3do.line_fit, show=show, mode=mode, ymin=ymin, ymax=ymax,
try_adjust_ax=try_adjust_ax, row=q3do.row, col=q3do.col)
else:
plotquest(q3do.wave, q3do.spec, q3do.cont_fit, q3do.ct_coeff, q3di, zstar=q3do.zstar, savefig=savefig, outfile=outfile,
templ_mask=templ_mask, show=show, mode=mode, ymin=ymin, ymax=ymax, try_adjust_ax=try_adjust_ax, row=q3do.row, col=q3do.col)
[docs]
def plotquest(MIRgdlambda,
MIRgdflux,
MIRcontinuum,
ct_coeff,
q3di,
zstar=0.,
savefig=True,
outfile=None,
templ_mask=[],
lines=[],
linespec=[],
show=False,
mode='light',
ymin=-1,
ymax=-1,
try_adjust_ax=True,
row=-1,
col=-1):
'''
Plot the fit to the residual of the quasar-host galaxy decomposition, if
refit is done with questfit. This function is presently only called
internally by :py:func:`~q3dfit.contfit.fitqsohost`.
Further documentation pending more testing and development.
'''
rcParamsOrig = rcParams.copy()
# dark mode just for fun:
if mode == 'dark':
pltstyle = 'dark_background'
dcolor = 'w'
else:
pltstyle = 'seaborn-v0_8-ticks'
dcolor = 'k'
plt.style.use(pltstyle)
# CB: Otherwise the background becomes black and the axes ticks
# unreadable when saving the figure
if mode == 'light':
rcParams['savefig.facecolor'] = 'white'
comp_best_fit = ct_coeff['comp_best_fit']
plot_noext = False # Remove dust contribution and plot intrinstic components
if 'plot_decomp' in q3di.argscontfit:
config_file = readcf(q3di.argscontfit['config_file'])
global_extinction = False
for key in config_file:
try:
if 'global' in config_file[key][3]:
global_extinction = True
except:
continue
fig = plt.figure(figsize=(6, 9))
gs = fig.add_gridspec(6,1, top=0.95, bottom=0.08, left=0.2)
ax1 = fig.add_subplot(gs[:5, :])
ax1.plot(MIRgdlambda, MIRgdflux,color='black')
if len(lines)==0:
ax1.plot(MIRgdlambda, MIRcontinuum, color='r')
else:
ax1.plot(MIRgdlambda, MIRcontinuum + linespec, color='darkorange')
if len(templ_mask)>0:
MIRgdlambda_temp = MIRgdlambda[templ_mask]
else:
MIRgdlambda_temp = MIRgdlambda
if len(lines)>0:
for line_i in lines:
ax1.axvline(line_i * (1. + zstar), color='grey', linestyle='--', alpha=0.7, zorder=0)
#ax1.axvspan(line_i-max(q3di.siglim_gas), line_i+max(q3di.siglim_gas))
ax1.plot(MIRgdlambda, linespec, color='r', linestyle='-', alpha=0.7, linewidth=1.5)
colour_list = ['dodgerblue', 'mediumblue', 'salmon', 'palegreen', 'orange', 'purple', 'forestgreen', 'darkgoldenrod', 'mediumblue', 'magenta', 'plum', 'yellowgreen']
if global_extinction:
str_global_ext = list(comp_best_fit.keys())[-2]
str_global_ice = list(comp_best_fit.keys())[-1]
# global_ext is a multi-dimensional array
if len(comp_best_fit[str_global_ext].shape) > 1:
comp_best_fit[str_global_ext] = comp_best_fit[str_global_ext] [:,0,0]
# global_ice is a multi-dimensional array
if len(comp_best_fit[str_global_ice].shape) > 1:
comp_best_fit[str_global_ice] = comp_best_fit[str_global_ice] [:,0,0]
count = 0
for i, el in enumerate(comp_best_fit):
if (el != str_global_ext) and (el != str_global_ice):
if len(comp_best_fit[el].shape) > 1: # component is a multi-dimensional array
comp_best_fit[el] = comp_best_fit[el] [:,0,0]
if plot_noext:
if count>len(colour_list)-1:
ax1.plot(MIRgdlambda_temp, comp_best_fit[el]/comp_best_fit[str_global_ext]/comp_best_fit[str_global_ice], label=el,linestyle='--',alpha=0.5)
else:
ax1.plot(MIRgdlambda_temp, comp_best_fit[el]/comp_best_fit[str_global_ext]/comp_best_fit[str_global_ice], color=colour_list[count], label=el,linestyle='--',alpha=0.5)
else:
if count>len(colour_list)-1:
ax1.plot(MIRgdlambda_temp, comp_best_fit[el], label=el,linestyle='--',alpha=0.5)
else:
ax1.plot(MIRgdlambda_temp, comp_best_fit[el], color=colour_list[count], label=el,linestyle='--',alpha=0.5)
count += 1
else:
count = 0
for i, el in enumerate(comp_best_fit):
if len(comp_best_fit[el].shape) > 1:
comp_best_fit[el] = comp_best_fit[el] [:,0,0]
if not ('_ext' in el or '_abs' in el):
spec_i = comp_best_fit[el]
label_i = el
if not plot_noext:
if el+'_ext' in comp_best_fit.keys():
spec_i = spec_i*comp_best_fit[el+'_ext']
if el+'_abs' in comp_best_fit.keys():
spec_i = spec_i*comp_best_fit[el+'_abs']
if count>len(colour_list)-1:
ax1.plot(MIRgdlambda_temp, spec_i, label=label_i,linestyle='--',alpha=0.5)
else:
ax1.plot(MIRgdlambda_temp, spec_i, label=label_i, color=colour_list[i], linestyle='--',alpha=0.5)
count += 1
ax1.legend(ncol=2)
ax1.set_xscale('log')
ax1.set_yscale('log')
#ax1.set_ylim(1e-5,1e2)
ax1.set_ylabel('Flux')
if try_adjust_ax:
adjust_ax(ax1, fig, minor=True)
ax1.tick_params(axis='x', which='both', bottom=False, top=False, labelbottom=False) # turn off major & minor ticks on the x-axis
ax2 = fig.add_subplot(gs[5:6, :], sharex=ax1)
if len(lines)>=1:
ax1.set_ylim(min(MIRcontinuum)/1e3, 3*max(MIRcontinuum + linespec))
ax2.plot(MIRgdlambda,MIRgdflux/(MIRcontinuum + linespec),color='black')
else:
ax1.set_ylim(min(MIRcontinuum)/1e3, 3*max(max(MIRgdflux), max(MIRcontinuum)))
ax2.plot(MIRgdlambda,MIRgdflux/MIRcontinuum,color='black')
if ymin>0.:
ax1.set_ylim(bottom=ymin)
if ymax>0.:
ax1.set_ylim(top=ymax)
ax2.axhline(1, color='grey', linestyle='--', alpha=0.7, zorder=0)
ax2.set_ylabel('Data/Model')
ax2.set_xlabel('Wavelength [micron]')
from matplotlib.ticker import ScalarFormatter
ax2.xaxis.set_major_formatter(ScalarFormatter())
ax2.xaxis.set_minor_formatter(ScalarFormatter())
ax2.ticklabel_format(style='plain')
if row>-1 and col>-1:
ax1.set_title('Spaxel [{}, {}]'.format(col, row), fontsize=20)
gs.update(wspace=0.0, hspace=0.05)
adjust_ax(ax2, fig)
if savefig and outfile is not None:
if len(outfile[0])>1:
plt.savefig(outfile[0]+'.jpg')
else:
plt.savefig(outfile+'.jpg')
else:
fig.savefig(outfile + '.jpg')
if show:
plt.show()
rcParams.update(rcParamsOrig)