Repository: spacetelescope/synphot_refactor

Issue:
Use separate package for extinction curve
As brought up by @karllark over at spacetelescope/pysynphot#54
<hr/>
There are external packages that provide extinction curves for python. One of these packages could be used instead of the numerical tables currently used.
Pros: offload the updating of the extinction curves to the "community"
Cons: Less control of the provided curves (but could be mitigated by automated testing against existing numerical results?)

Possible packages I know of:
https://github.com/kbarbary/extinction
https://github.com/karllark/dust_extinction
https://github.com/karllark/dust_attenuation
https://github.com/eblur/dust

In the interest of full disclose: I'm responsible for the 2nd package that has the aspiration of being an astropy affiliated package. The goals of my package are to provide dust extinction curves to the community and for use in dust extinction research (mine and others). Feedback is always welcomed.

...

I agree that only the dust_extinction package uses the astropy.modeling framework. I'm finding using this framework very powerful and useful. I've also been using the astropy.units framework. Very nice as well.

Discussion:

Patch:
diff --git a/CHANGES.rst b/CHANGES.rst
index e7e583a..0e11d85 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -10,6 +10,7 @@
   compatible with Python 3.6 or later and ``astropy`` 3.x or later. [#243]
 - Added support for ``RickerWavelet1D`` model that is the renamed version
   of ``MexicanHat1D`` model to be consistent with ``astropy`` 4.0. [#250]
+- Added support for extinction curve from ``dust-extinction``. [#251]
 - Added option for ``synphot.utils.download_data()`` to download to the cache
   instead of a specific location. Please note that new option is not fully
   compatible with customization using ``synphot.cfg``. [#211]
diff --git a/docs/conf.py b/docs/conf.py
index a74a7f9..1dd3449 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -128,7 +128,8 @@ man_pages = [('index', project.lower(), project + u' Documentation',
 intersphinx_mapping.update({
     'astropy': ('https://docs.astropy.org/en/latest/', None),
     'stsynphot': ('https://stsynphot.readthedocs.io/en/latest/', None),
-    'specutils': ('https://specutils.readthedocs.io/en/latest/', None)})
+    'specutils': ('https://specutils.readthedocs.io/en/latest/', None),
+    'dust-extinction': ('https://dust-extinction.readthedocs.io/en/latest/', None)})
 
 # -- Options for linkcheck output ---------------------------------------------
 linkcheck_retry = 5
diff --git a/docs/synphot/spectrum.rst b/docs/synphot/spectrum.rst
index 0fe28e7..36a12eb 100644
--- a/docs/synphot/spectrum.rst
+++ b/docs/synphot/spectrum.rst
@@ -104,7 +104,9 @@ de-reddens the spectrum), and then multiply it to the source::
 .. |xg_ext_ref| replace:: :ref:`Calzetti et al. (2000) <synphot-ref-extinction-calzetti2000>`
 
 For extinction due to Lyman-alpha forest, see :ref:`tutorial_lyman_alpha`
-tutorial.
+tutorial. To use a model from the
+:ref:`dust-extinction <dust-extinction:dust-extinction-mainpage>` package, see
+:ref:`tutorial_dust_extinction`.
 
 You can **redshift** a source spectrum in several ways (shown in example
 below), either by setting its ``z`` attribute or passing in a ``z`` keyword
diff --git a/docs/synphot/tutorials.rst b/docs/synphot/tutorials.rst
index bcbdce1..7adee12 100644
--- a/docs/synphot/tutorials.rst
+++ b/docs/synphot/tutorials.rst
@@ -98,6 +98,35 @@ continuum-normalized spectrum.
     plt.axhline(1, ls='--', color='k')
 
 
+.. _tutorial_dust_extinction:
+
+Using dust-extinction model
+---------------------------
+
+In this tutorial, you will learn how to apply an extinction curve using a
+model from the
+:ref:`dust-extinction <dust-extinction:dust-extinction-mainpage>` package:
+
+.. doctest-requires:: dust-extinction
+
+    >>> import numpy as np
+    >>> from astropy import units as u
+    >>> from dust_extinction.parameter_averages import CCM89
+    >>> from synphot import SourceSpectrum, ReddeningLaw
+    >>> from synphot.models import BlackBodyNorm1D
+    >>> ccm89_model = CCM89(Rv=3.1)
+    >>> wav = np.arange(0.1, 3, 0.001) * u.micron
+    >>> ebv = 0.1  # E(B-V)
+    >>> redlaw = ReddeningLaw(ccm89_model)
+    >>> extcurve = redlaw.extinction_curve(ebv, wavelengths=wav)
+    >>> bb = SourceSpectrum(BlackBodyNorm1D, temperature=5000 * u.K)
+    >>> bb.integrate()  # doctest: +FLOAT_CMP
+    <Quantity 9.31004974 ph / (cm2 s)>
+    >>> sp = bb * extcurve
+    >>> sp.integrate()  # doctest: +FLOAT_CMP
+    <Quantity 8.27106961 ph / (cm2 s)>
+
+
 .. _tutorial_fit_ew:
 
 Fitting, Equivalent Width
diff --git a/setup.cfg b/setup.cfg
index 8e36b72..068051d 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -10,6 +10,7 @@ xfail_strict = true
 filterwarnings =
     ignore:can't resolve package from __spec__
     ignore:numpy.ufunc size changed:RuntimeWarning
+    ignore:Class CCM89 defines class attributes
 
 [metadata]
 name = synphot
@@ -46,6 +47,7 @@ python_requires = >=3.6
 [options.extras_require]
 all =
     specutils>=0.7
+    dust-extinction
 test =
     pytest-astropy
 docs =
diff --git a/synphot/compat.py b/synphot/compat.py
index bcd3c94..4c75f8b 100644
--- a/synphot/compat.py
+++ b/synphot/compat.py
@@ -11,6 +11,14 @@ except ImportError:
 else:
     HAS_SPECUTILS = True
 
-__all__ = ['ASTROPY_LT_4_0', 'HAS_SPECUTILS']
+try:
+    import dust_extinction  # noqa
+except ImportError:
+    HAS_DUST_EXTINCTION = False
+else:
+    HAS_DUST_EXTINCTION = True
+
+
+__all__ = ['ASTROPY_LT_4_0', 'HAS_SPECUTILS', 'HAS_DUST_EXTINCTION']
 
 ASTROPY_LT_4_0 = not minversion(astropy, '4.0')
diff --git a/synphot/models.py b/synphot/models.py
index 0ec273a..f124dd9 100644
--- a/synphot/models.py
+++ b/synphot/models.py
@@ -388,9 +388,15 @@ class Empirical1D(Tabular1D):
         # Assume NaN at both ends need to be extrapolated based on
         # nearest end point.
         if self.fill_value is np.nan:
-            # Cannot use sampleset() due to ExtinctionModel1D
+            # Cannot use sampleset() due to ExtinctionModel1D.
             x = np.squeeze(self.points)
 
+            # np.squeeze may throw unit away.
+            if (isinstance(self.points, tuple) and
+                    isinstance(self.points[0], u.Quantity) and
+                    not isinstance(x, u.Quantity)):
+                x = x * self.points[0].unit
+
             if np.isscalar(y):  # pragma: no cover
                 if inputs < x[0]:
                     y = self.lookup_table[0]
diff --git a/synphot/reddening.py b/synphot/reddening.py
index e67e403..029cb80 100644
--- a/synphot/reddening.py
+++ b/synphot/reddening.py
@@ -10,6 +10,7 @@ from astropy import units as u
 
 # LOCAL
 from . import exceptions, specio, units
+from .compat import HAS_DUST_EXTINCTION
 from .config import Conf
 from .models import Empirical1D
 from .spectrum import BaseUnitlessSpectrum
@@ -74,11 +75,16 @@ class ReddeningLaw(BaseUnitlessSpectrum):
         elif not isinstance(ebv, numbers.Real):
             raise exceptions.SynphotError('E(B-V)={0} is invalid.'.format(ebv))
 
-        x = self._validate_wavelengths(wavelengths).value
-        y = 10 ** (-0.4 * self(x).value * ebv)
-        header = {
-            'E(B-V)': ebv,
-            'ReddeningLaw': self.meta.get('expr', 'unknown')}
+        x = self._validate_wavelengths(wavelengths)
+        header = {'E(B-V)': ebv}
+
+        # Duck-typing dust-extinction package API.
+        if HAS_DUST_EXTINCTION and hasattr(self.model, 'extinguish'):
+            y = self.model.extinguish(x, Ebv=ebv)
+            header['ReddeningLaw'] = '{!r}'.format(self.model)
+        else:
+            y = 10 ** (-0.4 * self(x).value * ebv)
+            header['ReddeningLaw'] = self.meta.get('expr', 'unknown')
 
         return ExtinctionCurve(ExtinctionModel1D, points=x, lookup_table=y,
                                meta={'header': header})


Test patch:
diff --git a/synphot/tests/test_reddening.py b/synphot/tests/test_reddening.py
index 7fca735..ec70d44 100644
--- a/synphot/tests/test_reddening.py
+++ b/synphot/tests/test_reddening.py
@@ -13,12 +13,14 @@ import pytest
 # ASTROPY
 from astropy import units as u
 from astropy.io import fits
+from astropy.tests.helper import assert_quantity_allclose
 from astropy.utils import minversion
 from astropy.utils.data import get_pkg_data_filename
 from astropy.utils.exceptions import AstropyUserWarning
 
 # LOCAL
 from .. import exceptions, units
+from ..compat import HAS_DUST_EXTINCTION  # noqa
 from ..models import ConstFlux1D, Empirical1D
 from ..observation import Observation
 from ..reddening import ReddeningLaw, etau_madau
@@ -65,6 +67,22 @@ class TestExtinction(object):
         np.testing.assert_array_equal(w, self.redlaw.waveset)
         np.testing.assert_allclose(y, 10 ** (-0.4 * self.redlaw(w) * 0.3))
 
+    @pytest.mark.skipif('not HAS_DUST_EXTINCTION')
+    def test_extcurve_call_2(self):
+        from dust_extinction.parameter_averages import CCM89
+
+        Av = 2
+        wav = np.arange(0.1, 3, 0.001) * u.micron
+        ext = CCM89(Rv=3.1)
+        ans = ext.extinguish(wav, Av)
+
+        ebv = Av / ext.Rv
+        redlaw = ReddeningLaw(ext)
+        extcurve = redlaw.extinction_curve(ebv, wavelengths=wav)
+
+        assert_quantity_allclose(extcurve(wav), ans)
+        assert extcurve.meta['header']['ReddeningLaw'] == '<CCM89(Rv=3.1)>'
+
     def test_mul_spec(self):
         """Apply extinction curve in inverse micron to flat spectrum in
         Angstrom.