Merge pull request #284 from drvdputt/units-standard

Standardize units (extended)

pahfit/features/features.py

@@ -21,16 +21,33 @@
 import numpy as np
 from astropy.table import vstack, Table, TableAttribute
 from astropy.io.misc.yaml import yaml
-import astropy.units as u
 from importlib import resources
 from pahfit.errors import PAHFITFeatureError
 from pahfit.features.features_format import BoundedMaskedColumn, BoundedParTableFormatter
+import pahfit.units
+
+# Feature kinds and associated parameters
+KIND_PARAMS = {'starlight': {'temperature', 'tau'},
+               'dust_continuum': {'temperature', 'tau'},
+               'line': {'wavelength', 'power'},  # 'fwhm', Instrument Pack detail!
+               'dust_feature': {'wavelength', 'fwhm', 'power'},
+               'attenuation': {'model', 'tau', 'geometry'},
+               'absorption': {'wavelength', 'fwhm', 'tau', 'geometry'}}
+
+# Parameter default units: flux density/intensity/power (other units determined on fit)
+# Note: power is actually amplitude for now. Change this to
+# intensity_power when the power features are implemented.
+
+PARAM_UNITS = {'temperature': pahfit.units.temperature,
+               'wavelength': pahfit.units.wavelength,
+               'fwhm': pahfit.units.wavelength,
+               'power': pahfit.units.intensity}
 
 
 class UniqueKeyLoader(yaml.SafeLoader):
     def construct_mapping(self, node, deep=False):
         mapping = set()
-        for key_node, value_node in node.value:
+        for key_node, _ in node.value:
             key = self.construct_object(key_node, deep=deep)
             if key in mapping:
                 raise PAHFITFeatureError(f"Duplicate {key!r} key found in YAML.")
@@ -102,26 +119,24 @@ def value_bounds(val, bounds):
 
 
 class Features(Table):
-    """A class for holding PAHFIT features and their associated
-    parameter information.  Note that each parameter has an associated
-    `kind', and that each kind has an associated set of allowable
-    parameters (see _kind_params, below).
+    """A class for holding a table of PAHFIT features and associated
+    parameter information.
+
+    Note that each parameter has an associated `kind', and that each
+    kind has an associated set of allowable parameters (see
+    `KIND_PARAMS`).
+
+    See Also
+    --------
+    `~astropy.table.Table`: The parent table class.
     """
 
     TableFormatter = BoundedParTableFormatter
     MaskedColumn = BoundedMaskedColumn
 
     param_covar = TableAttribute(default=[])
-    _kind_params = {'starlight': {'temperature', 'tau'},
-                    'dust_continuum': {'temperature', 'tau'},
-                    'line': {'wavelength', 'power'},  # 'fwhm', Instrument Pack detail!
-                    'dust_feature': {'wavelength', 'fwhm', 'power'},
-                    'attenuation': {'model', 'tau', 'geometry'},
-                    'absorption': {'wavelength', 'fwhm', 'tau', 'geometry'}}
-
-    _units = {'temperature': u.K, 'wavelength': u.um, 'fwhm': u.um}
+    _param_attrs = set(('value', 'bounds'))  # params can have these attributes
     _group_attrs = set(('bounds', 'features', 'kind'))  # group-level attributes
-    _param_attrs = set(('value', 'bounds'))  # Each parameter can have these attributes
     _no_bounds = set(('name', 'group', 'geometry', 'model'))  # String attributes (no bounds)
 
     @classmethod
@@ -178,7 +193,7 @@ def _read_scipack(cls, file):
                 raise PAHFITFeatureError(f"No kind found for {name}\n\t{file}")
 
             try:
-                valid_params = cls._kind_params[kind]
+                valid_params = KIND_PARAMS[kind]
             except KeyError:
                 raise PAHFITFeatureError(f"Unknown kind {kind} for {name}\n\t{file}")
             unknown_params = [x for x in keys
@@ -263,7 +278,7 @@ def _add_feature(cls, kind: str, t: dict, name: str, *,
         t[kind][name]['group'] = group
         t[kind][name]['kind'] = kind
         for (param, val) in pars.items():
-            if param not in cls._kind_params[kind]:
+            if param not in KIND_PARAMS[kind]:
                 continue
             if isinstance(val, dict):  # A param attribute dictionary
                 unknown_attrs = [x for x in val.keys() if x not in cls._param_attrs]
@@ -306,28 +321,38 @@ def _construct_table(cls, inp: dict):
         """
         tables = []
         for (kind, features) in inp.items():
-            kind_params = cls._kind_params[kind]  # All params for this kind
+            if kind == "_ratios":
+                continue
+            kp = KIND_PARAMS[kind]  # All params for this kind
             rows = []
             for (name, params) in features.items():
-                for missing in kind_params - params.keys():
+                for missing in kp - params.keys():
                     if missing in cls._no_bounds:
                         params[missing] = 0.0
                     else:
                         params[missing] = value_bounds(0.0, bounds=(0.0, None))
                 rows.append(dict(name=name, **params))
             table_columns = rows[0].keys()
             t = cls(rows, names=table_columns)
-            for p in cls._kind_params[kind]:
+            for p in KIND_PARAMS[kind]:
                 if p not in cls._no_bounds:
                     t[p].info.format = "0.4g"  # Nice format (customized by Formatter)
             tables.append(t)
         tables = vstack(tables)
         for cn, col in tables.columns.items():
-            if cn in cls._units:
-                col.unit = cls._units[cn]
-        tables.add_index('name')
+            if cn in PARAM_UNITS:
+                col.unit = PARAM_UNITS[cn]
+        cls._index_table(tables)
+
+        if '_ratios' in inp:
+            tables.meta['_ratios'] = inp['_ratios']
         return tables
 
+    @staticmethod
+    def _index_table(tbl):
+        for indx in ('name', 'group'):
+            tbl.add_index(indx)
+
     def mask_feature(self, name, mask_value=True):
         """Mask all the parameters of a feature.
 
@@ -344,7 +369,7 @@ def mask_feature(self, name, mask_value=True):
 
         """
         row = self.loc[name]
-        relevant_params = self._kind_params[row['kind']]
+        relevant_params = KIND_PARAMS[row['kind']]
         for col_name in relevant_params:
             if col_name in self._no_bounds:
                 # these are all strings, so can't mask

pahfit/helpers.py

@@ -1,11 +1,14 @@
 import os
 from importlib import resources
 
-from specutils import Spectrum1D
-
 from pahfit.component_models import BlackBody1D, S07_attenuation
+from pahfit import units
+
+from specutils import Spectrum1D
 from astropy.modeling.physical_models import Drude1D
 from astropy.modeling.functional_models import Gaussian1D
+from astropy import units as u
+from astropy.nddata import StdDevUncertainty
 
 __all__ = ["read_spectrum", "calculate_compounds"]
 
@@ -81,8 +84,18 @@ def read_spectrum(specfile, format=None):
     else:
         tformat = format
 
-    return Spectrum1D.read(specfile, format=tformat)
+    s = Spectrum1D.read(specfile, format=tformat)
+
+    # Convert to intensity units by assuming an arbitrary solid angle
+    # for now. To be removed when dual unit support (intensity and flux
+    # density) is supported.
+    if s.flux.unit.is_equivalent(units.flux_density):
+        solid_angle = (3 * u.arcsec)**2
+        alt_flux = (s.flux / solid_angle).to(units.intensity)
+        alt_unc_array = (s.uncertainty.array * s.flux.unit / solid_angle).to(units.intensity)
+        s = Spectrum1D(alt_flux, s.spectral_axis, uncertainty=StdDevUncertainty(alt_unc_array))
 
+    return s
 
 def calculate_compounds(obsdata, pmodel):
     """

pahfit/model.py

@@ -12,6 +12,7 @@
 from pahfit import instrument
 from pahfit.errors import PAHFITModelError
 from pahfit.component_models import BlackBody1D
+from pahfit import units
 
 
 class Model:
@@ -288,9 +289,12 @@ def amp_guess(row, fwhm):
 
     @staticmethod
     def _convert_spec_data(spec, z):
-        """
-        Turn astropy quantities stored in Spectrum1D into fittable
-        numbers.
+        """Convert Spectrum1D Quantities to fittable numbers.
+
+        The unit of the input spectrum has to be a multiple of MJy / sr,
+        the internal intensity unit. The output of this function
+        consists of simple unitless arrays (the numbers in these arrays
+        are assumed to be consistent with the internal units).
 
         Also corrects for redshift.
 
@@ -300,10 +304,13 @@ def _convert_spec_data(spec, z):
 
         xz, yz, uncz: wavelength in micron, flux, uncertainty
             corrected for redshift
+
         """
+        if not spec.flux.unit.is_equivalent(units.intensity):
+            raise PAHFITModelError("For now, PAHFIT only supports intensity units, i.e. convertible to MJy / sr.")
+        y = spec.flux.to(units.intensity).value
         x = spec.spectral_axis.to(u.micron).value
-        y = spec.flux.value
-        unc = spec.uncertainty.array
+        unc = (spec.uncertainty.array * spec.flux.unit).to(units.intensity).value
 
         # transform observed wavelength to "physical" wavelength
         xz = x / (1 + z)  # wavelength shorter

pahfit/units.py

@@ -0,0 +1,17 @@
+import astropy.units as u
+from astropy.units import CompositeUnit
+
+
+# Working/default unitParameter default units: flux density/intensity/power
+# These are PAHFITs default science packs parameter and output units
+temperature = u.K
+wavelength = u.um
+flux_density = u.mJy
+flux_power = CompositeUnit(1e-22, (u.W, u.m), (1, -2))
+solid_angle = u.sr
+intensity = u.MJy / u.sr
+intensity_power = CompositeUnit(1e-10, (u.W, u.m, u.sr), (1, -2, -1))
+
+# Note: integrated power units of 1e-22 W/m^2 (from flux) corresponds
+# to the unit 1e-10 W/m^2/sr (from intensity) if it occurs uniformly
+# over a solid angle 0.21" on a side (about a small JWST IFU pixel)