Source code for pdrtpy.plot.lineratioplot

#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')