#todo:
# keywords for show_both need to be arrays. ugh.
#
# allow levels to be percent?
#
# Look into seaborn https://seaborn.pydata.org
# Also https://docs.bokeh.org/en
# especially for coloring and style
from copy import deepcopy
import warnings
import numpy as np
import scipy.stats as stats
from matplotlib.lines import Line2D
import astropy.wcs as wcs
import astropy.units as u
from astropy.units import UnitsWarning
from .plotbase import PlotBase
from .modelplot import ModelPlot
from .. import pdrutils as utils
[docs]class LineRatioPlot(PlotBase):
"""LineRatioPlot plots the results of :class:`~pdrtpy.tools.lineratiofit.LineRatioFit`. It can plot maps of fit results, observations with errors on top of models, chi-square and confidence intervals and more.
:Keyword Arguments:
The methods of this class can take a variety of optional keywords. See the general `Plot Keywords`_ documentation
"""
def __init__(self,tool):
"""Init method
:param tool: The line ratio fitting tool that is to be plotted.
:type tool: `~pdrtpy.tool.LineRatioFit`
"""
super().__init__(tool)
self._figure = None
self._axis = None
self._modelplot = ModelPlot(self._tool._modelset,self._figure,self._axis)
self._ratiocolor=[]
[docs] def modelintensity(self,id,**kwargs):
"""Plot one of the model intensities
:param id: the intensity identifier, such as `CO_32``.
:type id: str
:param \**kwargs: see class documentation above
:raises KeyError: if is id not in existing model intensities
"""
ms = self._tool.modelset
if id not in ms.supported_intensities["intensity label"]:
raise KeyError(f"{id} is not in the ModelSet of your LineRatioFit")
model = ms.get_models([id],model_type="intensity")
kwargs_opts = {'title': self._tool._modelset.table.loc[id]["title"], 'colorbar':True}
kwargs_opts.update(kwargs)
self._modelplot._plot_no_wcs(model[id],**kwargs_opts)
self._figure = self._modelplot.figure
self._axis = self._modelplot.axis
[docs] def modelratio(self,id,**kwargs):
"""Plot one of the model ratios
:param id: the ratio identifier, such as ``CII_158/CO_32``.
:type id: str
:param \**kwargs: see class documentation above
:raises KeyError: if is id not in existing model ratios
"""
if self._tool._modelratios[id].shape == (1,): #does this ever occur??
return self._tool._modelratios[id]
kwargs_opts = {'title': self._tool._modelset.table.loc[id]["title"], 'units': u.dimensionless_unscaled , 'colorbar':True}
kwargs_opts.update(kwargs)
self._modelplot._plot_no_wcs(self._tool._modelratios[id],**kwargs_opts)
self._figure = self._modelplot.figure
self._axis = self._modelplot.axis
[docs] def observedratio(self,id,**kwargs):
"""Plot one of the observed ratios
:param id: the ratio identifier, such as ``CII_158/CO_32``.
:type id: - str
:raises KeyError: if id is not in existing observed ratios
"""
if self._tool._observedratios[id].shape == (1,0) or self._tool.has_vectors:
return self._tool._observedratios[id]
kwargs_opts = {'title': self._tool._modelset.table.loc[id]["title"], 'units': u.dimensionless_unscaled , 'colorbar':False}
kwargs_opts.update(kwargs)
self._plot(data=self._tool._observedratios[id],**kwargs_opts)
[docs] def density(self,**kwargs):
'''Plot the hydrogen nucleus volume density map that was computed by :class:`~pdrtpy.tool.lineratiofit.LineRatioFit` tool. Default units: cm :math:`^{-3}`
'''
kwargs_opts = {'units': 'cm^-3',
'aspect': 'equal',
'image':True,
'contours': False,
'label': False,
'linewidths': 1.0,
'levels': None,
'norm': None,
'title': None}
kwargs_opts.update(kwargs)
# handle single pixel or multi-pixel non-map cases.
if self._tool._density.shape == (1,) or self._tool.has_vectors:
return utils.to(kwargs_opts['units'],self._tool._density)
tunit=u.Unit(kwargs_opts['units'])
kwargs_opts['title'] = r"n [{0:latex_inline}]".format(tunit)
self._plot(self._tool._density,**kwargs_opts)
[docs] def radiation_field(self,**kwargs):
'''Plot the radiation field map that was computed by :class:`~pdrtpy.tool.lineratiofit.LineRatioFit` tool. Default units: Habing.
'''
kwargs_opts = {'units': 'Habing',
'aspect': 'equal',
'image':True,
'contours': False,
'label': False,
'linewidths': 1.0,
'levels': None,
'norm': None,
'title': None}
kwargs_opts.update(kwargs)
# handle single pixel or multi-pixel non-map cases.
if self._tool.radiation_field.shape == (1,) or self._tool.has_vectors:
return utils.to(kwargs_opts['units'],self._tool.radiation_field)
if kwargs_opts['title'] is None:
rad_title = utils.get_rad(kwargs_opts['units'])
tunit=u.Unit(kwargs_opts['units'])
kwargs_opts['title'] = r"{0} [{1:latex_inline}]".format(rad_title,tunit)
self._plot(self._tool.radiation_field,**kwargs_opts)
#@TODO refactor this method with reduced_chisq()
[docs] def chisq(self, **kwargs):
'''Plot the :math:`\chi^2` map that was computed by the
:class:`~pdrtpy.tool.lineratiofit.LineRatioFit` tool.
'''
kwargs_opts = {'units': None,
'aspect': 'equal',
'image':True,
'contours': True,
'label': False,
'colors': ['white'],
'linewidths': 1.0,
'norm': 'simple',
'stretch': 'linear',
'xaxis_unit': None,
'yaxis_unit': None,
'legend': None,
'bbox_to_anchor':None,
'loc':'upper center',
'title': None}
kwargs_opts.update(kwargs)
# make a sensible choice about contours if image is not shown
if not kwargs_opts['image'] and kwargs_opts['colors'][0] == 'white':
kwargs_opts['colors'][0] = 'black'
if self._tool.has_vectors:
raise NotImplementedError("Plotting of chi-square is not yet implemented for vector Measurements.")
if self._tool.has_maps:
data = self._tool.chisq(min=True)
if 'title' not in kwargs:
kwargs_opts['title'] = r'$\chi^2$ (dof=%d)'%self._tool._dof
self._plot(data,**kwargs_opts)
else:
data = self._tool.chisq(min=False)
self._modelplot._plot_no_wcs(data,header=None,**kwargs_opts)
# Put a crosshair where the chisq minimum is.
if False:
# To do this we first get the array index of the minimum
# then use WCS to translate to world coordinates.
[row,col] = np.where(self._tool._chisq==self._tool._chisq_min.value)
mywcs = wcs.WCS(data.header)
# Suppress WCS warning about 1/cm3 not being FITS
warnings.simplefilter('ignore',category=UnitsWarning)
logn,logrf = mywcs.array_index_to_world(row,col)
warnings.resetwarnings()
n = (10**logn.value[0])*u.Unit(logn.unit.to_string())
rf = (10**logrf.value[0])*u.Unit(logrf.unit.to_string())
# since the minimum can now be between pixels, we use the value
# of radiation field and density directly. (Why didn't we do this before?)
if kwargs_opts['xaxis_unit'] is not None:
x = utils.to(self._tool._density,kwargs_opts['xaxis_unit']).value
else:
x = self._tool._density.value
if kwargs_opts['yaxis_unit'] is not None:
y = utils.to(kwargs_opts['yaxis_unit'],self._tool._radiation_field).value
else:
y = self._tool._radiation_field.value
#print("Plot x,y %s,%s"%(x,y))
if kwargs_opts['title'] is None:
kwargs_opts['title'] = r'$\chi^2$ (dof=%d)'%self._tool._dof
label = r'$\chi_{min}^2$ = %.2g @ (n,FUV) = (%.2g,%.2g)'%(self._tool._chisq_min.value,x,y)
self._modelplot._axis[0].scatter(x,y,c='r',marker='+',s=200,linewidth=2,label=label)
# handle legend locally
if kwargs_opts['legend']:
legend = self._modelplot._axis[0].legend(title=kwargs_opts['title'],
bbox_to_anchor=kwargs_opts['bbox_to_anchor'],
loc=kwargs_opts['loc'])
self._figure = self._modelplot.figure
self._axis = self._modelplot.axis
[docs] def reduced_chisq(self, **kwargs):
'''Plot the reduced :math:`\chi^2` map that was computed by the
:class:`~pdrtpy.tool.lineratiofit.LineRatioFit` tool.
'''
kwargs_opts = {'units': None,
'aspect': 'equal',
'image':True,
'contours': True,
'label': False,
'colors': ['white'],
'linewidths': 1.0,
'norm': 'simple',
'stretch': 'linear',
'xaxis_unit': None,
'yaxis_unit': None,
'legend': None,
'bbox_to_anchor':None,
'loc':'upper center',
'title': None
}
kwargs_opts.update(kwargs)
if self._tool.has_vectors:
raise NotImplementedError("Plotting of chi-square is not yet implemented for vector Measurements.")
if self._tool.has_maps:
if 'title' not in kwargs:
kwargs_opts['title'] = r'$\chi_\nu^2$ (dof=%d)'%self._tool._dof
data = self._tool.reduced_chisq(min=True)
self._plot(data,**kwargs_opts)
# doesn't make sense to point out minimum chisq on a spatial-spatial map,
# so no legend
else:
data = self._tool.reduced_chisq(min=False)
self._modelplot._plot_no_wcs(data,header=None,**kwargs_opts)
# Put a crosshair where the chisq minimum is.
if False:
# To do this we first get the array index of the minimum
# then use WCS to translate to world coordinates.
[row,col] = np.where(self._tool._reduced_chisq==self._tool._reduced_chisq_min.value)
mywcs = wcs.WCS(data.header)
# Suppress WCS warning about 1/cm3 not being FITS
warnings.simplefilter('ignore',category=UnitsWarning)
logn,logrf = mywcs.array_index_to_world(row,col)
warnings.resetwarnings()
# logn, logrf are Quantities of the log(density) and log(radiation field),
# respectively. The model default units are cm^-2 and erg/s/cm^-2.
# These must be converted to plot units based on user input
# xaxis_unit and yaxis_unit.
# Note: multiplying by n.unit causes the ValueError:
# "The unit '1/cm3' is unrecognized, so all arithmetic operations with it are invalid."
# Yet by all other measures this appears to be a valid unit.
# The workaround is to us to_string() method.
n = (10**logn.value[0])*u.Unit(logn.unit.to_string())
rf = (10**logrf.value[0])*u.Unit(logrf.unit.to_string())
# since the minimum can now be between pixels, we use the value
# of radiation field and density directly. (Why didn't we do this before?)
if kwargs_opts['xaxis_unit'] is not None:
x = utils.to(self._tool._density,kwargs_opts['xaxis_unit']).value
else:
x = self._tool._density.value
if kwargs_opts['yaxis_unit'] is not None:
y = utils.to(kwargs_opts['yaxis_unit'],self._tool._radiation_field).value
else:
y = self._tool._radiation_field.value
#print("Plot x,y %s,%s"%(x,y))
if kwargs_opts['title'] is None:
kwargs_opts['title'] = r'$\chi_\nu^2$ (dof=%d)'%self._tool._dof
label = r'$\chi_{\nu,min}^2$ = %.2g @ (n,FUV) = (%.2g,%.2g)'%(self._tool._reduced_chisq_min.value,x,y)
self._modelplot.axis[0].scatter(x,y,c='r',marker='+',s=200,linewidth=2,label=label)
# handle legend locally
if kwargs_opts['legend']:
legend = self._modelplot.axis[0].legend(title=kwargs_opts['title'],
bbox_to_anchor=kwargs_opts['bbox_to_anchor'],
loc=kwargs_opts['loc'])
self._figure = self._modelplot.figure
self._axis = self._modelplot.axis
[docs] def show_both(self,units = ['Habing','cm^-3'], **kwargs):
'''Plot both radiation field and volume density maps computed by the
:class:`~pdrtpy.tool.lineratiofit.LineRatioFit` tool in a 1x2 panel subplot. Default units: ['Habing','cm^-3']
'''
_index = [1,2]
_reset = [True,False]
kwargs_opts = {'image':True,
'aspect': 'equal',
'contours': False,
'label': False,
'levels': None,
'norm': None,
'title': None,
'nrows': 1,
'ncols': 2,
'index': _index[0],
'reset': _reset[0]
}
kwargs_opts.update(kwargs)
rf = self.radiation_field(units=units[0],**kwargs_opts)
kwargs_opts['index'] = _index[1]
kwargs_opts['reset'] = _reset[1]
d = self.density(units=units[1],**kwargs_opts)
# @todo don't return for plots. print for non-plots
return (rf,d)
[docs] def confidence_intervals(self,**kwargs):
'''Plot the confidence intervals from the :math:`\chi^2` map computed by the
:class:`~pdrtpy.tool.lineratiofit.LineRatioFit` tool. Default levels: [50., 68., 80., 95., 99.]
**Currently only works for single-pixel Measurements**
'''
if self._tool.has_maps or self._tool.has_vectors:
raise NotImplementedError("Plotting of confidence intervals is not yet implemented for maps or vectors.")
kwargs_opts = {'units': None,
'aspect': 'auto',
'image':False,
'contours': True,
'label': True,
'levels': [50., 68., 80., 95., 99.],
'colors': ['black'],
'linewidths': 1.0,
'norm': 'simple',
'stretch': 'linear',
'xaxis_unit': None,
'yaxis_unit': None,
'title': "Confidence Intervals"}
# ensure levels are in ascending order
kwargs_opts['levels'] = sorted(kwargs_opts['levels'])
kwargs_opts.update(kwargs)
confidence = deepcopy(self._tool._chisq)
confidence.data = 100*stats.distributions.chi2.cdf(confidence.data,self._tool._dof)
self._tool.confidence = confidence
self._modelplot._plot_no_wcs(data=confidence,header=None,**kwargs_opts)
self._figure = self._modelplot.figure
self._axis = self._modelplot.axis
[docs] def overlay_all_ratios(self,**kwargs):
'''Overlay all the measured ratios and their errors on the :math:`(n,F_{FUV})` space.
This only works for single-valued Measurements; an overlay for multi-pixel doesn't make sense.
'''
#NB: could have position and area though.
if self._tool.has_maps or self._tool.has_vectors:
raise NotImplementedError("Plotting of ratio overlays is not yet implemented for maps or vectors.")
kwargs_opts = {'units': None,
'image':False,
'contours': False,
'meas_color': self._CB_color_cycle,
'levels' : None,
'label': False,
'linewidths': 1.0,
'ncols': 1,
'norm': None,
'title': None,
'reset': True,
'legend': True,
'bbox_to_anchor':None,
'loc': 'upper center'}
kwargs_opts.update(kwargs)
# force this as ncols !=1 makes no sense.
kwargs_opts['ncols'] = 1
i = 0
_measurements = list()
_models = list()
# get_model will correctly raise exception if m.id not in ModelSet
meas_passed = False
if kwargs_opts.get('measurements',None) is not None:
# avoid modifying a passed parameter
_measurements = deepcopy(kwargs_opts['measurements'])
meas_passed = True
for m in _measurements:
if i > 0:
kwargs_opts['reset']=False
val = self._tool.modelset.get_model(m.id)
_models.append(val)
kwargs_opts['measurements'] = [utils.convert_if_necessary(m)]
self._modelplot._plot_no_wcs(_models[i],header=None,colorcounter=i,**kwargs_opts)
i = i+1
km = kwargs_opts.pop('measurements')
#for m in _measurements:
# print("A ",type(m))
for key,val in self._tool._modelratios.items():
#kwargs_opts['measurements'] = [self._tool._observedratios[key]]
#print(" K ",type(self._tool._observedratios[key]))
#_measurements.append(self._tool._observedratios[key])
if i > 0:
kwargs_opts['reset']=False
# pass the index of the contour color to use via the "secret" colorcounter keyword.
self._modelplot._plot_no_wcs(val,header=None,colorcounter=i,**kwargs_opts,measurements=[self._tool._observedratios[key]])
i = i+1
#for m in _measurements:
# print("B ",type(m))
if kwargs_opts['legend']:
lines = [Line2D([0], [0], color=c, linewidth=3, linestyle='-') for c in kwargs_opts['meas_color'][0:i]]
labels = list()
if meas_passed:
# THIS WILL EXCEPT IF m.id NOT IN MODELSET WHICH IS VERY POSSIBLE!
for m in _measurements:
try:
tt = self._tool.modelset.get_model(m.id).title
except Exception:
tt = m.id
labels.append(tt)
#labels = [m.id for m in _measurements]
title = "Observed Ratios and Intensities"
else:
title = "Observed Ratios"
labels.extend([self._tool._modelratios[k].title for k in self._tool._modelratios])
#print("LABELS ",labels)
self._plt.legend(lines, labels,loc=kwargs_opts['loc'],
bbox_to_anchor=kwargs_opts['bbox_to_anchor'],
title=title)
self._figure = self._modelplot.figure
self._axis = self._modelplot.axis
[docs] def ratios_on_models(self,**kwargs):
'''Overlay all the measured ratios and their errors on the individual models for those ratios. Plots are displayed in multi-column format, controlled the `ncols` keyword. Default: ncols=2
**Currently only works for single-pixel Measurements**
'''
if self._tool.has_maps or self._tool.has_vectors:
raise NotImplementedError("Plotting of ratio overlays is not yet implemented for maps or vectors.")
kwargs_opts = {'units': None,
'image':True,
'colorbar': True,
'contours': True,
'colors': ['white'],
'levels' : None,
'label': False,
'linewidths': 1.0,
'meas_color': ['#4daf4a'],
'ncols': 2,
'norm': 'simple',
'stretch': 'linear',
'index': 1,
'reset': True,
'legend': True,
'bbox_to_anchor':None,
'loc':'upper center'}
kwargs_opts.update(kwargs)
kwargs_opts["ncols"] = min(kwargs_opts["ncols"],self._tool.ratiocount)
kwargs_opts["nrows"] = int(round(self._tool.ratiocount/kwargs_opts["ncols"]+0.49,0))
# defend against meas_color not being a list
if type(kwargs_opts['meas_color']) == str:
#warnings.warn("meas_color should be a list")
kwargs_opts['meas_color']=[kwargs_opts['meas_color']]
for key,val in self._tool._modelratios.items():
axidx = kwargs_opts['index']-1
if kwargs_opts['index'] > 1:
kwargs_opts['reset'] = False
m = self._tool._model_files_used[key]
kwargs_opts['measurements'] = [self._tool._observedratios[key]]
self._modelplot._plot_no_wcs(val,header=None,**kwargs_opts)
self._axis = self._modelplot._axis
self._figure = self._modelplot._figure
kwargs_opts['index'] = kwargs_opts['index'] + 1
if kwargs_opts['legend']:
if 'title' not in kwargs: # then it was None, and we customize it
_title = self._tool._modelratios[key].title
else:
_title = kwargs['title']
lines = list()
labels = list()
if kwargs_opts['contours']:
lines.append(Line2D([0], [0], color=kwargs_opts['colors'][0], linewidth=3, linestyle='-'))
labels.append("model")
lines.append(Line2D([0], [0], color=kwargs_opts['meas_color'][0], linewidth=3, linestyle='-'))
labels.append("observed")
#maybe loc should be 'best' but then it bounces around
self._axis[axidx].legend(lines, labels,
bbox_to_anchor=kwargs_opts['bbox_to_anchor'],
loc=kwargs_opts['loc'],
title=_title)
# Turn off subplots greater than the number of
# available ratios
for i in range(self._tool.ratiocount,len(self._axis)):
self._axis[i].axis('off')