Merge pull request #101 from pllim/fix-norm-overlap

 Fixed normalization requiring extrapolation

synphot/observation.py

@@ -50,7 +50,9 @@ class Observation(BaseSourceSpectrum):
         and bandpass do not fully overlap:
 
             * `None` or 'none' - Source must encompass bandpass (default)
-            * 'extrap' - Extrapolate source spectrum
+            * 'extrap' - Extrapolate source spectrum (this changes the
+              underlying model of ``spec`` to always extrapolate,
+              if applicable)
             * 'taper' - Taper source spectrum
 
     Raises
@@ -105,9 +107,7 @@ def __init__(self, spec, band, binset=None, force='none'):
                 warn['PartialOverlap'] = msg
 
             elif force.startswith('extrap'):
-                if isinstance(spec.model, Empirical1D):
-                    spec.model.method = 'nearest'
-                    spec.model.fill_value = None
+                if spec.force_extrapolation():
                     msg = ('Source spectrum will be extrapolated (at constant '
                            'value for empirical model).')
                 else:

synphot/spectrum.py

@@ -627,6 +627,37 @@ def check_overlap(self, other, wavelengths=None, threshold=0.01):
 
         return result
 
+    def force_extrapolation(self):
+        """Force the underlying model to extrapolate.
+
+        An example where this is useful: You read in a source spectrum
+        from file (by default, the model returns zero flux if given
+        wavelength is out of range) and you wish to force the underlying
+        empirical model to extrapolate based on nearest point.
+
+        .. note::
+
+            This is only applicable to `~synphot.models.Empirical1D` model
+            and should still work even if the source spectrum has been
+            redshifted.
+
+        Returns
+        -------
+        is_forced : bool
+            `True` if the model is successfully forced to be extrapolated,
+            else `False`.
+
+        """
+        # We use _model here in case the spectrum is redshifted.
+        if isinstance(self._model, Empirical1D):
+            self._model.method = 'nearest'
+            self._model.fill_value = None
+            is_forced = True
+        else:
+            is_forced = False
+
+        return is_forced
+
     def taper(self, wavelengths=None):
         """Taper the spectrum or bandpass.
 
@@ -864,7 +895,9 @@ def normalize(self, renorm_val, band=None, wavelengths=None, force=False,
             By default (`False`), renormalization is only done
             when band wavelength limits are within ``self``
             or at least 99% of the flux is within the overlap.
-            Set to `True` to force renormalization for partial overlap.
+            Set to `True` to force renormalization for partial overlap
+            (this changes the underlying model of ``self`` to always
+            extrapolate, if applicable).
             Disjoint bandpass raises an exception regardless.
 
         area, vegaspec
@@ -903,29 +936,32 @@ def normalize(self, renorm_val, band=None, wavelengths=None, force=False,
                 raise exceptions.DisjointError(
                     'Spectrum and renormalization band are disjoint.')
 
-            elif stat == 'partial_most':
-                warn_str = (
-                    'Spectrum is not defined everywhere in renormalization' +
-                    'bandpass. At least 99% of the band throughput has' +
-                    'data. Spectrum will be extrapolated at constant value.')
-                warndict['PartialRenorm'] = warn_str
-                warnings.warn(warn_str, AstropyUserWarning)
-
-            elif stat == 'partial_notmost':
-                if force:
-                    warn_str = (
-                        'Spectrum is not defined everywhere in '
-                        'renormalization bandpass. Less than 99% of the ' +
-                        'band throughput has data. Spectrum will be ' +
-                        'extrapolated at constant value.')
-                    warndict['PartialRenorm'] = warn_str
-                    warnings.warn(warn_str, AstropyUserWarning)
+            elif 'partial' in stat:
+                if stat == 'partial_most':
+                    warn_str = 'At least'
+                elif stat == 'partial_notmost' and force:
+                    warn_str = 'Less than'
                 else:
                     raise exceptions.PartialOverlap(
                         'Spectrum and renormalization band do not fully '
                         'overlap. You may use force=True to force the '
                         'renormalization to proceed.')
 
+                warn_str = (
+                    'Spectrum is not defined everywhere in renormalization '
+                    'bandpass. {0} 99% of the band throughput has '
+                    'data. Spectrum will be').format(warn_str)
+
+                if self.force_extrapolation():
+                    warn_str = ('{0} extrapolated at constant '
+                                'value.').format(warn_str)
+                else:
+                    warn_str = ('{0} evaluated outside pre-defined '
+                                'waveset.').format(warn_str)
+
+                warnings.warn(warn_str, AstropyUserWarning)
+                warndict['PartialRenorm'] = warn_str
+
             elif stat != 'full':  # pragma: no cover
                 raise exceptions.SynphotError(
                     'Overlap result of {0} is unexpected.'.format(stat))
@@ -946,12 +982,14 @@ def normalize(self, renorm_val, band=None, wavelengths=None, force=False,
             stdflux = 1.0
         else:
             totalflux = sp.integrate(wavelengths=wavelengths)
+
             # VEGAMAG
             if renorm_unit_name == units.VEGAMAG.to_string():
                 if not isinstance(vegaspec, SourceSpectrum):
                     raise exceptions.SynphotError(
                         'Vega spectrum is missing.')
                 stdspec = vegaspec
+
             # Magnitude flux-density units
             elif renorm_val.unit in (u.STmag, u.ABmag):
                 stdspec = SourceSpectrum(

synphot/tests/test_spectrum.py

@@ -527,6 +527,33 @@ def test_exceptions(self):
             self.bp.check_overlap(1)
 
 
+class TestForceExtrap(object):
+    """Test forcing extrapolation on a source spectrum."""
+
+    @pytest.mark.skipif('not HAS_SCIPY')
+    @pytest.mark.parametrize('z', [0, 0.03])
+    def test_empirical(self, z):
+        sp = SourceSpectrum(Empirical1D, points=[1000, 2000, 3000, 4000],
+                            lookup_table=[0.5, 0.6, 10.6, 1.5])
+        sp.z = z
+        w = [900, 4300]
+        np.testing.assert_allclose(sp(w).value, 0)  # No extrapolation
+
+        is_forced = sp.force_extrapolation()  # Force extrapolation
+        assert is_forced
+        np.testing.assert_allclose(sp(w).value, [0.5, 1.5])
+
+    def test_analytical(self):
+        """Forcing is not possible."""
+        sp = SourceSpectrum(GaussianFlux1D, mean=5500, total_flux=1, fwhm=10)
+        w = [100, 10000]
+        np.testing.assert_allclose(sp(w).value, 0)
+
+        is_forced = sp.force_extrapolation()
+        assert not is_forced
+        np.testing.assert_allclose(sp(w).value, 0)
+
+
 @pytest.mark.skipif('not HAS_SCIPY')
 class TestNormalize(object):
     """Test source spectrum normalization."""