"""Manage pre-computed PDR models"""
import collections
import itertools
from copy import deepcopy
import astropy.units as u
import numpy as np
from astropy.table import Column, Table, unique, vstack
from .measurement import Measurement
from .utils import _OBS_UNIT_, get_table, model_dir
[docs]
class ModelSet:
"""Class for computed PDR Model Sets. :class:`ModelSet` provides interface to a directory containing the model FITS files and the ability to query details about.
Parameters
----------
name : str
Identifying name, e.g., 'wk2006'.
z : float
Metallicity in solar units.
medium : str
Medium type, e.g. 'constant density', 'clumpy', 'non-clumpy'.
losangle : float
The line-of-sight inclination angle (aka viewing angle), losangle=0 is
face-on, losangle=90 is edge-on. Valid values are 0, 30, 45, 60, 75, 90.
mass : float
Maximum clump mass (for KosmaTau models). Default: None (appropriate for
Wolfire/Kaufman models).
modelsetinfo : str or :class:`~astropy.table.Table`
For adding user specified ModelSet, this parameter specifies the file with
information about the ModelSet. It can be a pathname to an external tabular
file in an astropy-recognized format or an astropy Table object. If an
external tabular file, its format should be given by the `format` keyword.
The columns are:
.. code-block:: python
['PDRcode','name', 'version','path','filename','medium','z','inc', 'mass','description']
`z`, 'inc', and `mass` should be numbers, the rest should be strings.
The `name`, `version`, `medium`, `z`, 'inc', and `mass` columns contain
the values available in the input ModelSet as described above. `PDRcode`
is the originator of the code e.g., "KOSMA-tau", `version` is the code
version, `path` is the full-qualified pathname to the FITS files,
`filename` is the tabular file which contains descriptions of individual
FITS files. This must also be in the format specified by the `format`
keyword. `description` is a description of the input ModelSet.
format : str
If `modelsetinfo` is a file, give the format of file. Must be
`an astropy recognized Table format.
<https://docs.astropy.org/en/stable/io/unified.html#built-in-readers-writers>`_
e.g., ascii, ipac, votable. Default is `IPAC format
<https://docs.astropy.org/en/stable/api/astropy.io.ascii.Ipac.html#astropy.io.ascii.Ipac>`_
Raises
------
ValueError
If model set not recognized/found.
"""
def __init__(
self,
name,
z,
medium="constant density",
losangle=0,
mass=None,
modelsetinfo=None,
format="ipac",
debug=False,
):
if modelsetinfo is None:
# get the package default
self._all_models = get_table("all_models.tab")
elif isinstance(modelsetinfo, str):
self._all_models = Table.read(modelsetinfo, format=format)
else: # must be an Astropy Table
self._all_models = deepcopy(modelsetinfo)
self._all_models.add_index("name")
possible = dict()
if name not in self._all_models["name"]:
raise ValueError(
f"Unrecognized PDR model code {name:s}. Choices are: {set(list(self._all_models['name']))}"
)
if np.all(self._all_models.loc[name]["mass"].mask):
matching_rows = np.where(
(self._all_models["z"] == z)
& (self._all_models["medium"] == medium)
& (self._all_models["losangle"] == losangle)
)
possible["mass"] = None
else: # Kosma-tau
matching_rows = np.where(
(self._all_models["z"] == z)
& (self._all_models["medium"] == medium)
& (self._all_models["mass"] == mass)
)
possible["mass"] = self._all_models.loc[name]["mass"]
for key in ["z", "medium", "losangle"]:
possible[key] = self._all_models.loc[name][key]
# possible[] can be a Python native or a Column depending on how many rows matched.
# If only one row matches it will be a native, otherwise a Column.
# Check the type so we can create a sorted numpy array in either case.
for j in possible:
if possible[j] is None:
continue
if isinstance(possible[j], Column):
# convert Column to np.array
possible[j] = sorted(set(np.array(possible[j])))
else:
# convert native to np.array
possible[j] = sorted(set(np.array([possible[j]])))
if mass is None and possible["mass"] is not None:
raise ValueError(f"mass value is required for model {name:s}. Allowed values are {possible['mass']}")
if matching_rows[0].size == 0:
msg = (
f"Requested ModelSet not found in {name:s}. Check your input values. Allowed z are {possible['z']}. "
f" Allowed medium are {possible['medium']}. Allowed losangle are {possible['losangle']}."
)
if possible["mass"] is not None:
msg = msg + f" Allowed mass are {possible['mass']}."
raise ValueError(msg)
self._tabrow = self._all_models[matching_rows].loc[name]
tpath = self._tabrow["path"]
self._table = get_table(path=tpath, filename=self._tabrow["filename"], format=format)
self._table.add_index("ratio")
self._set_identifiers()
self._set_ratios()
self._default_unit = dict()
self._default_unit["ratio"] = u.dimensionless_unscaled
self._default_unit["intensity"] = _OBS_UNIT_
self._default_unit["emissivity"] = "erg / (cm3 ion s)"
self._user_added_models = dict()
self._av_table = get_table("av.tab", format="ascii")
self._av_table.add_index("LOSANGLE")
@property
def description(self):
"""The description of this model
Returns
-------
str
"""
s = f", Z={self.z:2.1f}, losangle={int(self.losangle):d}"
if self.avlos is not None:
s += f", avlos={self.avlos:.2f}"
if self.avperp is not None:
s += f", avperp={self.avperp:.2f}"
return self._tabrow["description"] + s
@property
def code(self):
"""The PDR code that created this ModelSet, e.g. 'KOSMA-tau'
Returns
-------
str
"""
return self._tabrow["PDRcode"]
@property
def name(self):
"""The name of this ModelSet
Returns
-------
str
"""
return self._tabrow["name"]
@property
def version(self):
"""The version of this ModelSet
Returns
-------
str
"""
return self._tabrow["version"]
@property
def medium(self):
"""The medium type of this model, e.g. 'constant density', 'clumpy', 'non-clumpy'
Returns
-------
str
"""
return self._tabrow["medium"]
@property
def z(self):
"""The metallicity of this ModelSet
Returns
-------
float
"""
return self.metallicity
@property
def metallicity(self):
"""The metallicity of this ModelSet
Returns
-------
float
"""
return self._tabrow["z"]
@property
def losangle(self):
"""The inclination (viewing) angle with respect to the line of sight, in degrees. A value of 0 means face-on.
Returns
-------
float
"""
return self._tabrow["losangle"]
@property
def viewangle(self):
"""The viewing angle (inclination) with respect to the line of sight, in degrees. A value of 0 means face-on.
Returns
-------
float
"""
return self.losangle
@property
def avperp(self):
"""The PDR thickness for inclined viewing angle models, expressed in visual magnitudes.
Returns ``None`` for face-on geometry (losangle=0) and for non-wk2020 models.
Returns
-------
float or None
"""
if not self.is_wk2020:
return None
avperp = float(self._av_table.loc[int(self.losangle)]["AVPERP"])
if avperp == 0.0:
return None
return avperp
@property
def avlos(self):
"""The visual extinction along the line of sight, expressed in magnitudes.
Returns ``None`` for non-wk2020 models.
Returns
-------
float or None
"""
if not self.is_wk2020:
return None
return float(self._av_table.loc[int(self.losangle)]["AVLOS"])
@property
def table(self):
"""The table containing details of the models in this ModelSet.
Returns
-------
:class:`astropy.table.Table`
"""
return self._table
@property
def identifiers(self):
"""Table of lines and continuum that are included in ratio models of this ModelSet. Only lines and continuum that are part of ratios are included in this list. For a separate list of line and continuum intensity models see :meth:`~pdrtpy.modelset.ModelSet.supported_intensities`.
Returns
-------
:class:`astropy.table.Table`
"""
return self._identifiers
@property
def supported_intensities(self):
"""Table of lines and continuum that are covered by this ModelSet and have models separate from
the any ratio model they might be in.
Returns
-------
:class:`astropy.table.Table`
"""
return self._supported_lines
@property
def supported_ratios(self):
"""The emission ratios that are covered by this ModelSet
Returns
-------
:class:`astropy.table.Table`
"""
return self._supported_ratios
@property
def user_added_models(self):
"""Show which models have been added to this ModelSet by the user
Returns
-------
list
"""
return list(self._user_added_models.keys())
[docs]
def ratiocount(self, m):
"""The number of valid ratios in this ModelSet, given a list of observation (:class:`~pdrtpy.measurement.Measurement`) identifiers.
Parameters
----------
m : list
List of string :class:`~pdrtpy.measurement.Measurement` IDs, e.g. ["CII_158","OI_145","FIR"].
Returns
-------
int
The number of model ratios found for the given list of measurement IDs.
"""
return len(self._get_ratio_elements(m))
[docs]
def find_pairs(self, m):
"""Find the valid model ratios labels in this ModelSet for a given list of measurement IDs
Parameters
----------
m : list
List of string `Measurement` IDs, e.g. ["CII_158","OI_145","FIR"].
Returns
-------
iterator
An iterator of model ratio labels for the given list of measurement IDs.
"""
if not isinstance(m, collections.abc.Iterable) or isinstance(m, (str, bytes)):
raise Exception("m must be an array of strings")
for q in itertools.product(m, m):
if q[0] == "FIR" and (q[1] == "OI_145" or q[1] == "OI_63") and "CII_158" in m:
s = q[1] + "+CII_158/" + q[0]
else:
s = q[0] + "/" + q[1]
if s in self.table["ratio"]:
yield s
[docs]
def find_files(self, m, ext="fits"):
"""Find the valid model ratios files in this ModelSet for a given list of measurement IDs. See :meth:`~pdrtpy.measurement.Measurement.id`
Parameters
----------
m : list
List of string :class:`~pdrtpy.measurement.Measurement` IDs, e.g. ["CII_158","OI_145","FIR"].
ext : str
File extension. Default: "fits".
Returns
-------
iterator
An iterator of model ratio files for the given list of Measurement IDs.
"""
if not isinstance(m, collections.abc.Iterable) or isinstance(m, (str, bytes)):
raise Exception("m must be an array of strings")
for q in itertools.product(m, m):
# must deal with OI+CII/FIR models. Note we must check for FIR first, since
# if you check q has OI,CII and m has FIR order you'll miss OI/CII.
if q[0] == "FIR" and (q[1] == "OI_145" or q[1] == "OI_63") and "CII_158" in m:
s = q[1] + "+CII_158/" + q[0]
else:
s = q[0] + "/" + q[1]
if s in self.table["ratio"]:
fullpath = self._tabrow["path"] + self.table.loc[s]["filename"] + "." + ext
tup = (s, fullpath)
yield tup
[docs]
def model_ratios(self, m):
"""Return the model ratios that match the input Measurement ID list. You must provide at least 2 Measurements IDs
Parameters
----------
m : list
List of string :class:`~pdrtpy.measurement.Measurement` IDs, e.g., ["CII_158","OI_145","FIR"].
Returns
-------
list
List of string identifiers of ratio IDs, e.g., ['OI_145/CII_158', 'OI_145+CII_158/FIR'].
"""
ratios = list()
if len(m) < 2:
raise Exception("m most contain at least two strings")
for p in self.find_files(m):
ratios.append(p[0])
return ratios
[docs]
def model_intensities(self, m):
"""Return the model intensities in this ModelSet that match the input Measurement ID list.
This method will return the intersection of the input list and the list of supported lines.
Parameters
----------
m : list
List of string :class:`~pdrtpy.measurement.Measurement` IDs, e.g., ["CII_158","OI_145","CS_21"].
Returns
-------
list
List of string identifiers of intensity IDs, e.g., ['CII_158','OI_145'].
"""
# get intersection of input list and supported lines
return list(set(m) & set(self._supported_lines["intensity label"]))
[docs]
def get_model(self, identifier, unit=None, ext="fits"):
"""Get a specific model by its identifier
Parameters
----------
identifier : str
A :class:`~pdrtpy.measurement.Measurement` ID. It can be an intensity or a ratio, e.g., "CII_158","CI_609/FIR".
Returns
-------
:class:`~pdrtpy.measurement.Measurement`
The model matching the identifier.
Raises
------
KeyError
If identifier not found in this ModelSet.
"""
if identifier in self._user_added_models:
return self._user_added_models[identifier]
if identifier not in self.table["ratio"]:
raise KeyError(f"{identifier} is not in this ModelSet")
_filename = self.table.loc[identifier]["filename"] + "." + ext
d = model_dir()
_thefile = d + self._tabrow["path"] + _filename
_title = self._table.loc[identifier]["title"]
# @TODO Fix this: see issues 66 & 67
if unit is None or unit == "":
if identifier == "TS":
unit = "K"
modeltype = "intensity" # ??'temperature'
# make a guess at the unit
elif "/" in identifier:
unit = self._default_unit["ratio"]
modeltype = "ratio"
else:
unit = self._default_unit["intensity"]
modeltype = "intensity" # this is wrong for emissivity modeltypes
elif unit == u.dimensionless_unscaled:
modeltype = "ratio"
else:
modeltype = "intensity"
_model = Measurement.read(_thefile, title=_title, unit=unit, identifier=identifier)
# if _model.unit=="":
# _model.unit = u.Unit("adu")#self._default_unit["ratio"]
_wcs = _model.wcs
if "MODELTYP" not in _model.header:
_model.header["MODELTYP"] = modeltype
_model.modeltype = modeltype
# @todo this is messy. clean up by doing if wcs.. first?
if self.is_wk2006 or self.name == "smc":
# fix WK2006 model headerslisthd
if _wcs.wcs.cunit[0] == "":
_model.header["CUNIT1"] = "cm^-3"
_wcs.wcs.cunit[0] = u.Unit("cm^-3")
else:
_model.header["CUNIT1"] = str(_wcs.wcs.cunit[0])
if _wcs.wcs.cunit[1] == "":
_model.header["CUNIT2"] = "Habing"
# Raises UnitScaleError:
# "The FITS unit format is not able to represent scales that are not powers of 10. Multiply your data by 1.600000e-03."
# This causes all sorts of downstream problems. Workaround in LineRatioFit.read_models().
# _model.wcs.wcs.cunit[1] = habing_unit
elif self.code == "KOSMA-tau":
# fix KosmaTau model headers
if _wcs.wcs.cunit[0] == "":
_model.header["CUNIT1"] = "cm^-3"
_wcs.wcs.cunit[0] = u.Unit("cm^-3")
else:
_model.header["CUNIT1"] = str(_wcs.wcs.cunit[0])
if _wcs.wcs.cunit[1] == "":
_model.header["CUNIT2"] = "Draine"
else:
_model.header["CUNIT2"] = str(_wcs.wcs.cunit[1])
else:
# copy wcs cunit to header. used later.
_model.header["CUNIT1"] = str(_wcs.wcs.cunit[0])
_model.header["CUNIT2"] = str(_wcs.wcs.cunit[1])
return _model
[docs]
def get_models(self, identifiers, model_type="ratio", ext="fits"):
"""Get the models from this ModelSet that match the input list of identifiers
Parameters
----------
identifiers : list
List of string :class:`~pdrtpy.measurement.Measurement` IDs, e.g., ["CII_158","OI_145","CS_21"].
model_type : str
Indicates which type of model is requested: one of 'ratio' or 'intensity'.
Returns
-------
list
The matching models as a list of :class:`~pdrtpy.measurement.Measurement`.
Raises
------
KeyError
If identifiers not found in this ModelSet.
"""
# if identifier not in self._identifiers["ID"]:
# raise Exception("There is no model in ModelSet %s with the identifier %s"%(self.name,identifier))
if model_type != "intensity" and model_type != "ratio" and model_type != "both":
raise ValueError("Unrecognized model_type: must be one of 'intensity', 'ratio', or 'both'")
models = dict()
a = list()
self._table.remove_indices("ratio")
self._table.add_index("ratio")
if model_type == "intensity" or model_type == "both":
a.extend(self.model_intensities(identifiers))
if model_type == "ratio" or model_type == "both":
a.extend(self.model_ratios(identifiers))
if model_type == "intensity" or model_type == "ratio":
_unit = self._default_unit[model_type]
else:
_unit = None
for k in a:
if k == "TS":
# kluge. we need to support model_type = "temperature"
models[k] = self.get_model(k, unit="K", ext=ext)
else:
models[k] = self.get_model(k, unit=_unit, ext=ext)
return models
[docs]
def add_model(self, identifier, model, title, overwrite=False):
r"""Add your own model to this ModelSet.
Parameters
----------
identifier : str
A :class:`~pdrtpy.measurement.Measurement` ID. It can be an intensity or a ratio, e.g., "CII_158","CI_609/FIR".
model : str or :class:`~pdrtpy.measurement.Measurement`
The model to add. If a string, this must be the fully-qualified path of a FITS file.
If a :class:`~pdrtpy.measurement.Measurement` it must have the same CTYPEs and CUNITs as the models in the ModelSet.
title : str
A formatted string (e.g., LaTeX) describing this observation that can be used for plotting.
Python r-strings are accepted, e.g., r'$^{13}$CO(3-2)' would give :math:`^{13}{\rm CO(3-2)}`.
overwrite : bool
Whether to overwrite the model if the identifier already exists in the ModelSet or has been
previously added. Default: False.
"""
if identifier not in self.table["ratio"] and identifier not in self._user_added_models:
self._really_add_model(identifier, model, title)
elif identifier in self._user_added_models and not overwrite:
raise Exception(
f"{identifier} was previously added to this ModelSet. If you wish to overwrite it, use overwrite=True"
)
elif identifier in self.table["ratio"] and not overwrite:
raise Exception(
f"{identifier} is already in the {self.name} ModelSet. If you wish to overwrite it, use overwrite=True"
)
else:
# print(f"Overwriting {identifier}.")
self._really_add_model(identifier, model, title)
def _really_add_model(self, identifier, model, title):
print(f"Adding user model {identifier}")
if type(model) is str:
m = Measurement.read(model, identifier=identifier)
else:
m = model
m._title = title
self._user_added_models[identifier] = m
if "/" in identifier: # it's a ratio
if identifier in self._supported_ratios["ratio label"]:
# ack, there should be a way just to replace title but I can't get Table.loc to work.
index = np.where(self._supported_ratios["ratio label"] == identifier)[0][0]
self._supported_ratios.remove_row(index)
self._supported_ratios.add_row([title, identifier])
numerator, denominator = identifier.split("/")
fakefilename = "user-" + numerator.replace("_", "") + "_" + denominator.replace("_", "")
else:
if identifier in self._supported_lines["intensity label"]:
index = np.where(self._supported_lines["intensity label"] == identifier)[0][0]
self._supported_lines.remove_row(index)
self._supported_lines.add_row([title, identifier])
numerator = identifier
denominator = 1
fakefilename = "user-" + numerator.replace("_", "")
# numerator denominator ratio filename z title
self.table.add_row([numerator, denominator, identifier, fakefilename, self.z, title, str(self.losangle)])
def _find_ratio_elements(self, m):
# TODO handle case of OI+CII/FIR so it is not special cased in lineratiofit.py
"""Find the valid model numerator,denominator pairs in this ModelSet for a given list of measurement IDs. See :meth:`~pdrtpy.measurement.Measurement.id`
Parameters
----------
m : list
List of string :class:`~pdrtpy.measurement.Measurement` IDs, e.g. ["CII_158","OI_145","FIR"].
Returns
-------
iterator
A dictionary iterator where key is numerator and value is denominator.
"""
if not isinstance(m, collections.abc.Iterable) or isinstance(m, (str, bytes)):
raise Exception("m must be an array of strings")
for q in itertools.product(m, m):
s = q[0] + "/" + q[1]
z = dict()
if s in self.table["ratio"]:
z = {"numerator": self.table.loc[s]["numerator"], "denominator": self.table.loc[s]["denominator"]}
yield z
def _get_ratio_elements(self, m):
"""Get the valid model numerator,denominator pairs in this ModelSet for a given list of measurement IDs. See :meth:`~pdrtpy.measurement.Measurement.id`
Parameters
----------
m : list
List of string :class:`~pdrtpy.measurement.Measurement` IDs, e.g. ["CII_158","OI_145","FIR"].
Returns
-------
dict
A dictionary where key is numerator and value is denominator.
"""
if not isinstance(m, collections.abc.Iterable) or isinstance(m, (str, bytes)):
raise Exception("m must be an array of strings")
k = list()
for q in itertools.product(m, m):
s = q[0] + "/" + q[1]
if s in self.table["ratio"]:
z = {"numerator": self.table.loc[s]["numerator"], "denominator": self.table.loc[s]["denominator"]}
k.append(z)
self._get_oi_cii_fir(m, k)
return k
def _get_oi_cii_fir(self, m, k):
"""Utility method for determining ratio elements, to handle special
case of ([O I] 63 micron + [C II] 158 micron)/I_FIR.
Parameters
----------
m : list
List of string :class:`~pdrtpy.measurement.Measurement` IDs, e.g. ["CII_158","OI_145","FIR"].
k : list
List of existing ratios to append to.
"""
if "CII_158" in m and "FIR" in m:
if "OI_63" in m:
num = "OI_63+CII_158"
den = "FIR"
# lab="OI_63+CII_158/FIR"
z = {"numerator": num, "denominator": den}
k.append(z)
if "OI_145" in m:
num = "OI_145+CII_158"
den = "FIR"
# lab="OI_145+CII_158/FIR"
z = {"numerator": num, "denominator": den}
k.append(z)
def _set_ratios(self):
"""make a useful table of ratios covered by this model"""
self._supported_ratios = Table([self.table["title"], self.table["denominator"], self.table["ratio"]], copy=True)
matching_rows = np.where(self._supported_ratios["denominator"] == "1")[0]
self._supported_lines = Table(self._supported_ratios[matching_rows], copy=True)
self._supported_lines.remove_column("denominator")
self._supported_ratios.remove_rows(matching_rows)
self._supported_ratios["title"].unit = None
self._supported_ratios["ratio"].unit = None
self._supported_ratios.remove_column("denominator")
self._supported_ratios.rename_column("ratio", "ratio label")
self._supported_lines.rename_column("ratio", "intensity label")
def _set_identifiers(self):
"""make a useful table of identifiers of lines covered by ratios in this ModelSet"""
# remove the single line intensity models from the list.
matching_rows = np.where(self._table["denominator"] != "1")[0]
n = deepcopy(self._table["numerator"][matching_rows])
n.name = "ID"
d = deepcopy(self._table["denominator"][matching_rows])
d.name = "ID"
t1 = Table([self._table["title"][matching_rows], n], copy=True)
# discard the summed fluxes as user would input them individually
for id in ["OI_145+CII_158", "OI_63+CII_158"]:
a = np.where(t1["ID"] == id)[0]
for z in a:
t1.remove_row(z)
# now remove denominator from title (everything from / onwards)
for i in range(len(t1["title"])):
if "/" in t1["title"][i]:
t1["title"][i] = t1["title"][i][0 : t1["title"][i].index("/")]
t2 = Table([self._table["title"][matching_rows], d], copy=True)
# remove numerator from title (everything before and including /)
for i in range(len(t2["title"])):
if "/" in t2["title"][i]:
t2["title"][i] = t2["title"][i][t2["title"][i].index("/") + 1 :]
t = vstack([t1, t2])
t = unique(t, keys=["ID"], keep="first", silent=True)
t["title"].unit = None
t["ID"].unit = None
t.rename_column("title", "canonical name")
self._identifiers = t
@property
def is_wk2006(self):
"""Indicate this is a wk2006 model, to deal with quirks of that modelset.
Returns
-------
bool
True if it is.
"""
return self.name == "wk2006"
@property
def is_wk2020(self):
"""Indicate this is a wk2020 model, to deal with quirks of that modelset.
Returns
-------
bool
True if it is.
"""
return self.name == "wk2020"
# ============= Static Methods =============
[docs]
@staticmethod
def list():
"""Print the names and descriptions of available ModelSets (not just this one)"""
ModelSet.all_sets().pprint_all(align="<")
[docs]
@staticmethod
def all_sets(debug=False):
"""Return a table of the names and descriptions of available ModelSets (not just this one)
Returns
-------
:class:`~astropy.table.Table`
"""
t = get_table("all_models.tab")
t.remove_column("path")
t.remove_column("filename")
return t
