Repository: spacetelescope/synphot_refactor

Issue:
TST: Enable doctest
Examples in documentation should be tested. Care needs to be given such that floating point comparison and remote data access is marked appropriately. Also might have to skip some (e.g., some printed properties that will change from run to run).

Discussion:

Patch:
diff --git a/azure-pipelines.yml b/azure-pipelines.yml
index 149fb08..491518c 100644
--- a/azure-pipelines.yml
+++ b/azure-pipelines.yml
@@ -76,6 +76,25 @@ stages:
         codecov_token: $(CODECOV_TOKEN)
       displayName: 'Calculate coverage'
 
+  - job: 'LinkCheck'
+    pool:
+      vmImage: 'Ubuntu-16.04'
+    strategy:
+      maxParallel: 4
+
+    steps:
+    - task: UsePythonVersion@0
+
+    - script: |
+        python -m pip install --upgrade pip setuptools
+        pip install sphinx-astropy
+        python setup.py install
+
+    - bash: |
+        cd docs
+        make linkcheck
+      displayName: 'Run docs link check'
+
   - template: azure-templates.yml
     parameters:
       name: 'Older_Versions'
diff --git a/docs/conf.py b/docs/conf.py
index c888de6..8b9e306 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -149,3 +149,9 @@ intersphinx_mapping.update({
         'astropy': ('http://docs.astropy.org/en/latest/', None)})
 #intersphinx_mapping.update({
 #        'stsynphot': ('http://stsynphot.readthedocs.org/en/latest/', None)})
+
+# -- Options for linkcheck output -------------------------------------------
+linkcheck_retry = 5
+linkcheck_ignore = ['https://hsthelp.stsci.edu']
+linkcheck_timeout = 180
+linkcheck_anchors = False
diff --git a/docs/index.rst b/docs/index.rst
index 1f45347..006db38 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -1,5 +1,3 @@
-.. doctest-skip-all
-
 .. _astropy_synphot:
 
 ******************************
@@ -80,7 +78,7 @@ via HTTP, create a local directory where you plan to store the data files
 (e.g., ``/my/local/dir/cdbs``) and run the following:
 
     >>> from synphot.utils import download_data
-    >>> file_list = download_data('/my/local/dir/cdbs')
+    >>> file_list = download_data('/my/local/dir/cdbs')  # doctest: +SKIP
 
 Then copy `synphot.cfg <https://github.com/spacetelescope/synphot_refactor/blob/master/synphot/synphot.cfg>`_
 to your ``$HOME/.astropy/config/`` directory, and replace every instance of
@@ -103,6 +101,11 @@ can modify this line in your ``$HOME/.astropy/config/synphot.cfg`` file::
 
     johnson_v_file = /my/other/dir/my_johnson_v.fits
 
+.. testsetup::
+
+    >>> from synphot.config import conf
+    >>> conf.johnson_v_file = '/my/local/dir/cdbs/comp/nonhst/johnson_v_004_syn.fits'
+
 Alternately, you can also take advantage of :ref:`astropy:astropy_config`
 to manage **synphot** data files. For example, you can also temporarily use a
 different Johnson *V* throughput file::
@@ -116,6 +119,10 @@ different Johnson *V* throughput file::
     >>> print(conf.johnson_v_file)
     /my/local/dir/cdbs/comp/nonhst/johnson_v_004_syn.fits
 
+.. testsetup::
+
+    >>> conf.reset('johnson_v_file')
+
 
 .. _synphot_getting_started:
 
@@ -133,8 +140,41 @@ in Astropy. Despite this, some models like
 work, so they can be used directly. When in doubt, see if a model is in
 `synphot.models` first before using Astropy's.
 
+::
+
+    >>> from astropy import units as u
+    >>> from synphot import units, SourceSpectrum
+    >>> from synphot.models import (BlackBodyNorm1D, GaussianAbsorption1D,
+    ...                             GaussianFlux1D)
+    >>> from synphot.spectrum import BaseUnitlessSpectrum
+
+Create a Gaussian absorption profile that absorps 80% of the flux at
+4000 Angstrom with a sigma of 20 Angstrom::
+
+    >>> g_abs = BaseUnitlessSpectrum(GaussianAbsorption1D, amplitude=0.8,
+    ...                              mean=4000, stddev=20)
+
+Create a Gaussian emission line with the given total flux
+centered at 3000 Angstrom with FWHM of 100 Angstrom::
+
+    >>> g_em = SourceSpectrum(GaussianFlux1D,
+    ...                       total_flux=3.5e-13*u.erg/(u.cm**2 * u.s),
+    ...                       mean=3000, fwhm=100)
+
+Create a blackbody source spectrum with a temperature of 6000 K::
+
+    >>> bb = SourceSpectrum(BlackBodyNorm1D, temperature=6000)
+
+Combine the above components to create a source spectrum that is twice
+the original blackbody flux with the Gaussian emission and absorption lines::
+
+    >>> sp = g_abs * (g_em + 2 * bb)
+
+Plot the spectrum, zooming in on the line features::
+
+    >>> sp.plot(left=1, right=7000)  # doctest: +SKIP
+
 .. plot::
-    :include-source:
 
     from astropy import units as u
     from synphot import units, SourceSpectrum
@@ -161,37 +201,37 @@ work, so they can be used directly. When in doubt, see if a model is in
 
 Sample the spectrum at 0.3 micron::
 
-    >>> sp(0.3 * u.micron)
-    <Quantity 0.0012685 PHOTLAM>
+    >>> sp(0.3 * u.micron)  # doctest: +FLOAT_CMP
+    <Quantity 0.00129536 PHOTLAM>
 
 Or sample the same thing but in a different flux unit::
 
-    >>> sp(0.3 * u.micron, flux_unit=units.FLAM)
-    <Quantity 8.57716634e-15 FLAM>
+    >>> sp(0.3 * u.micron, flux_unit=units.FLAM)  # doctest: +FLOAT_CMP
+    <Quantity 8.57718043e-15 FLAM>
 
 Sample the spectrum at its native wavelength set::
 
-    >>> sp(sp.waveset)
+    >>> sp(sp.waveset)  # doctest: +FLOAT_CMP +ELLIPSIS
     <Quantity [0.00000000e+00, 0.00000000e+00, 0.00000000e+00, ...,
                4.47483287e-05, 4.34264666e-05, 4.21423394e-05] PHOTLAM>
 
 Models that built the spectrum::
 
-    >>> print(sp)
+    >>> print(sp)  # doctest: +ELLIPSIS
     SourceSpectrum at z=0.0
-    Model: CompoundModel2
+    Model: CompoundModel...
     Inputs: ('x',)
     Outputs: ('y',)
     Model set size: 1
     Expression: ([0] + ([1] | [2])) * [3]
     Components:
-        [0]: <GaussianFlux1D(amplitude=...>
-
-        [1]: <BlackBodyNorm1D(temperature=6000.0)>
-
-        [2]: <Scale(factor=2.0)>
-
-        [3]: <GaussianAbsorption1D(amplitude=0.8, mean=4000.0, stddev=20.0)>
+        [0]: <GaussianFlux1D(...>
+    <BLANKLINE>
+        [1]: <BlackBodyNorm1D(...)>
+    <BLANKLINE>
+        [2]: <Scale(...)>
+    <BLANKLINE>
+        [3]: <GaussianAbsorption1D(...)>
     Parameters:
         ...
 
@@ -210,8 +250,8 @@ Normalize the source spectrum to 1 Jy in a given box bandpass and
 integrate it::
 
     >>> sp_rn = sp.normalize(1 * u.Jy, band=bp)
-    >>> sp_rn.integrate()
-    <Quantity 12901.13466265 ph / (cm2 s)>
+    >>> sp_rn.integrate()  # doctest: +FLOAT_CMP
+    <Quantity 12883.09903646 ph / (cm2 s)>
 
 Create an observation by passing the redshifted and normalized source spectrum
 through the box bandpass::
@@ -223,10 +263,10 @@ Calculate the count rate of the observation above for an 2-meter telescope:
 
     >>> import numpy as np
     >>> area = np.pi * (1 * u.m) ** 2
-    >>> area
+    >>> area  # doctest: +FLOAT_CMP
     <Quantity 3.14159265 m2>
-    >>> obs.countrate(area=area)
-    <Quantity 24219651.1854991 ct / s>
+    >>> obs.countrate(area=area)  # doctest: +FLOAT_CMP
+    <Quantity 24219650.71772057 ct / s>
 
 
 .. _synphot_using:
@@ -251,9 +291,7 @@ Using **synphot**
 A Brief History
 ===============
 
-A brief history: First, there was
-`STSDAS SYNPHOT <http://stsdas.stsci.edu/cgi-bin/gethelp.cgi?synphot.sys>`_.
-Then, there was
+A brief history: First, there was STSDAS SYNPHOT (IRAF). Then, there was
 ASTROLIB PYSYNPHOT (:ref:`Lim et al. 2015 <synphot-ref-lim2015>`), which
 aimed at replacing STSDAS SYNPHOT using Python. In order to take advantage of
 :ref:`astropy:astropy-modeling` and :ref:`astropy:astropy-units` and to
diff --git a/docs/synphot/bandpass.rst b/docs/synphot/bandpass.rst
index b7d676d..ccc419b 100644
--- a/docs/synphot/bandpass.rst
+++ b/docs/synphot/bandpass.rst
@@ -1,5 +1,3 @@
-.. doctest-skip-all
-
 .. _bandpass-main:
 
 Bandpass
@@ -33,8 +31,8 @@ To **evaluate** its transmission at a given wavelength, use its
     >>> from synphot import SpectralElement
     >>> from synphot.models import Box1D
     >>> bp = SpectralElement(Box1D, amplitude=1, x_0=5000, width=100)
-    >>> bp(500 * u.nm)
-    <Quantity 1.0>
+    >>> bp(500 * u.nm)  # doctest: +FLOAT_CMP
+    <Quantity 1.>
 
 Bandpass also has access to **photometric parameter** calculations akin to
 IRAF SYNPHOT ``bandpar`` task. (Also see :ref:`synphot_formulae` and respective
@@ -42,34 +40,34 @@ API documentations.) Some of these need the information of telescope collecting
 area, which will be set to HST value in the examples below::
 
     >>> area = 45238.93416 * (u.cm * u.cm)  # HST
-    >>> bp.avgwave()
-    <Quantity 4999.995000001091 Angstrom>
-    >>> bp.tlambda()
-    <Quantity 1.0>
-    >>> bp.tpeak()
-    <Quantity 1.0>
-    >>> bp.wpeak()  # For box, this is the first occurence of max throughput
-    <Quantity 4950.0 Angstrom>
-    >>> bp.efficiency()
-    <Quantity 0.020000686709115298>
-    >>> bp.equivwidth()
-    <Quantity 100.00000000218279 Angstrom>
-    >>> bp.rectwidth()
-    <Quantity 100.00000000218279 Angstrom>
-    >>> bp.rmswidth()
-    <Quantity 28.867513315773838 Angstrom>
-    >>> bp.photbw()
-    <Quantity 28.867032163097033 Angstrom>
-    >>> bp.fwhm()
-    <Quantity 67.97666597821402 Angstrom>
-    >>> bp.pivot()
-    <Quantity 4999.911663668414 Angstrom>
-    >>> bp.barlam()
-    <Quantity 4999.744993503255 Angstrom>
-    >>> bp.unit_response(area)
-    <Quantity 8.78202683842001e-19 FLAM>
-    >>> bp.emflx(area)
-    <Quantity 8.782026838611703e-17 FLAM>
+    >>> bp.avgwave()  # doctest: +FLOAT_CMP
+    <Quantity 4999.995 Angstrom>
+    >>> bp.tlambda()  # doctest: +FLOAT_CMP
+    <Quantity 1.>
+    >>> bp.tpeak()  # doctest: +FLOAT_CMP
+    <Quantity 1.>
+    >>> bp.wpeak()  # For box, this is the first occurence of max throughput # doctest: +FLOAT_CMP
+    <Quantity 4950. Angstrom>
+    >>> bp.efficiency()  # doctest: +FLOAT_CMP
+    <Quantity 0.02000069>
+    >>> bp.equivwidth()  # doctest: +FLOAT_CMP
+    <Quantity 100. Angstrom>
+    >>> bp.rectwidth()  # doctest: +FLOAT_CMP
+    <Quantity 100. Angstrom>
+    >>> bp.rmswidth()  # doctest: +FLOAT_CMP
+    <Quantity 28.86751332 Angstrom>
+    >>> bp.photbw()  # doctest: +FLOAT_CMP
+    <Quantity 28.86703216 Angstrom>
+    >>> bp.fwhm()  # doctest: +FLOAT_CMP
+    <Quantity 67.97666598 Angstrom>
+    >>> bp.pivot()  # doctest: +FLOAT_CMP
+    <Quantity 4999.91166367 Angstrom>
+    >>> bp.barlam()  # doctest: +FLOAT_CMP
+    <Quantity 4999.7449935 Angstrom>
+    >>> bp.unit_response(area)  # doctest: +FLOAT_CMP
+    <Quantity 8.78202761e-19 FLAM>
+    >>> bp.emflx(area)  # doctest: +FLOAT_CMP
+    <Quantity 8.78202761e-17 FLAM>
 
 To check the **overlap**, i.e., whether the wavelength range of another
 spectrum is defined everywhere within the a bandpass, you can use
@@ -82,24 +80,24 @@ any potential inaccurate representation of the result::
     >>> bp = SpectralElement(
     ...     Empirical1D, points=[2999.9, 3000, 6000, 6000.1],
     ...     lookup_table=[0, 1, 1, 0])
-    # Source spectrum is fully defined within bandpass waveset
+    >>> # Source spectrum is fully defined within bandpass waveset
     >>> sp_full = SourceSpectrum(
     ...     Empirical1D, points=[999.9, 1000, 9000, 9000.1],
     ...     lookup_table=[0, 1, 1, 0])
     >>> bp.check_overlap(sp_full)
     'full'
-    # 99% of spectrum's flux is in the overlap (not a concern)
+    >>> # 99% of spectrum's flux is in the overlap (not a concern)
     >>> sp_most = SourceSpectrum(
     ...     Empirical1D, points=[3005, 3005.1, 6000.1, 6000.2],
     ...     lookup_table=[0, 1, 1, 0])
     >>> bp.check_overlap(sp_most)
     'partial_most'
-    # Source spectrum needs significant extrapolation (guessing)
+    >>> # Source spectrum needs significant extrapolation (guessing)
     >>> sp_notmost = SourceSpectrum(
     ...     Empirical1D, points=[3999.9, 4500.1], lookup_table=[1, 1])
     >>> bp.check_overlap(sp_notmost)
     'partial_notmost'
-    # No overlap at all
+    >>> # No overlap at all
     >>> sp_none = SourceSpectrum(
     ...     Empirical1D, points=[99.9, 100, 2999.8, 2999.9],
     ...     lookup_table=[0, 1, 1, 0])
@@ -225,14 +223,14 @@ appropriate filter names:
 
 The example below loads and plots bandpasses from Johnson *BI*::
 
-    >>> import matplotlib.pyplot as plt
+    >>> import matplotlib.pyplot as plt  # doctest: +SKIP
     >>> from synphot import SpectralElement
-    >>> b = SpectralElement.from_filter('johnson_b')
-    >>> i = SpectralElement.from_filter('johnson_i')
-    >>> plt.plot(b.waveset, b(b.waveset), 'b', i.waveset, i(i.waveset), 'r')
-    >>> plt.ylim(0, 1.1)
+    >>> b = SpectralElement.from_filter('johnson_b')  # doctest: +REMOTE_DATA
+    >>> i = SpectralElement.from_filter('johnson_i')  # doctest: +REMOTE_DATA
+    >>> plt.plot(b.waveset, b(b.waveset), 'b', i.waveset, i(i.waveset), 'r')  # doctest: +SKIP
+    >>> plt.ylim(0, 1.1)  # doctest: +SKIP
     >>> # Label comes from DESCRIP keyword from FITS header
-    >>> plt.legend([b.meta['header']['descrip'], i.meta['header']['descrip']])
+    >>> plt.legend([b.meta['header']['descrip'], i.meta['header']['descrip']])  # doctest: +SKIP
 
 .. image:: images/johnson_bi.png
     :width: 600px
diff --git a/docs/synphot/c_ext.rst b/docs/synphot/c_ext.rst
index 16cb9e0..a48e8e4 100644
--- a/docs/synphot/c_ext.rst
+++ b/docs/synphot/c_ext.rst
@@ -1,5 +1,3 @@
-.. doctest-skip-all
-
 .. _synphot-c-ext:
 
 calcbinflux
diff --git a/docs/synphot/formulae.rst b/docs/synphot/formulae.rst
index 731e64e..1200484 100644
--- a/docs/synphot/formulae.rst
+++ b/docs/synphot/formulae.rst
@@ -1,5 +1,3 @@
-.. doctest-skip-all
-
 .. _synphot_formulae:
 
 Photometric Properties
@@ -58,9 +56,9 @@ equivalent to IRAF SYNPHOT ``bandpar`` results for ``avglam``,
 
 Example::
 
-    >>> bp.avgwave()
+    >>> bp.avgwave()  # doctest: +FLOAT_CMP
     <Quantity 5367.903565874441 Angstrom>
-    >>> bp.tlambda()
+    >>> bp.tlambda()  # doctest: +FLOAT_CMP
     <Quantity 0.22807748889452203>
 
 
@@ -77,9 +75,9 @@ if peak value is not unique).
 
 Example::
 
-    >>> bp.tpeak()
+    >>> bp.tpeak()  # doctest: +FLOAT_CMP
     <Quantity 0.24144500494003296>
-    >>> bp.wpeak()
+    >>> bp.wpeak()  # doctest: +FLOAT_CMP
     <Quantity 5059.8212890625 Angstrom>
 
 
@@ -98,7 +96,7 @@ result for ``qtlam``.
 
 Example::
 
-    >>> bp.efficiency()
+    >>> bp.efficiency()  # doctest: +FLOAT_CMP
     <Quantity 0.05090165033079963>
 
 
@@ -118,9 +116,9 @@ IRAF SYNPHOT ``bandpar`` result for ``equvw``.
 
 Example::
 
-    >>> bp.equivwidth()
+    >>> bp.equivwidth()  # doctest: +FLOAT_CMP
     <Quantity 272.0108162945954 Angstrom>
-    >>> bp.integrate()
+    >>> bp.integrate()  # doctest: +FLOAT_CMP
     <Quantity 272.0108162945954 Angstrom>
 
 
@@ -139,7 +137,7 @@ the rectangular width. It is equivalent to IRAF SYNPHOT ``bandpar`` result for
 
 Example::
 
-    >>> bp.rectwidth()
+    >>> bp.rectwidth()  # doctest: +FLOAT_CMP
     <Quantity 1126.5953352903448 Angstrom>
 
 
@@ -159,7 +157,7 @@ the bandpass RMS width as defined in
 
 Example::
 
-    >>> bp.rmswidth()
+    >>> bp.rmswidth()  # doctest: +FLOAT_CMP
     <Quantity 359.56457676412236 Angstrom>
 
 
@@ -179,7 +177,7 @@ as :ref:`synphot-formula-rmswidth`.
 
 Example::
 
-    >>> bp.photbw()
+    >>> bp.photbw()  # doctest: +FLOAT_CMP
     <Quantity 357.17951791474843 Angstrom>
 
 
@@ -198,7 +196,7 @@ is :ref:`synphot-formula-bandw`.
 
 Example::
 
-    >>> bp.fwhm()
+    >>> bp.fwhm()  # doctest: +FLOAT_CMP
     <Quantity 841.0934884601406 Angstrom>
 
 
@@ -220,7 +218,7 @@ It is equivalent to IRAF SYNPHOT ``bandpar`` results for ``barlam``.
 
 Example::
 
-    >>> bp.barlam()
+    >>> bp.barlam()  # doctest: +FLOAT_CMP
     <Quantity 5331.8648495386 Angstrom>
 
 
@@ -240,7 +238,7 @@ It is equivalent to IRAF SYNPHOT ``bandpar`` result for ``uresp``.
 
 Example::
 
-    >>> bp.unit_response(area)
+    >>> bp.unit_response(area)  # doctest: +FLOAT_CMP
     <Quantity 3.007277127274156e-19 FLAM>
 
 
@@ -262,7 +260,7 @@ In the formula below, ``uresp``, ``equvw``, and :math:`\lambda_{0}` are
 
 Example::
 
-    >>> bp.emflx(area)
+    >>> bp.emflx(area)  # doctest: +FLOAT_CMP
     <Quantity 3.586552579909415e-16 FLAM>
 
 .. _synphot-formula-effstim:
@@ -283,15 +281,15 @@ The default binning behavior is to be consistent with ASTROLIB PYSYNPHOT.
 
 Example::
 
-    >>> obs.effstim()
+    >>> obs.effstim()  # doctest: +FLOAT_CMP
     <Quantity 0.00053744 PHOTLAM>
-    >>> obs.effstim('flam')
+    >>> obs.effstim('flam')  # doctest: +FLOAT_CMP
     <Quantity 1.99333435e-15 FLAM>
-    >>> obs.effstim('count', area=area)  # Not binned
+    >>> obs.effstim('count', area=area)  # Not binned  # doctest: +FLOAT_CMP
     <Quantity 6628.36886854 ct / s>
-    >>> obs.countrate(area=area, binned=False)
+    >>> obs.countrate(area=area, binned=False)  # doctest: +FLOAT_CMP
     <Quantity 6628.36886854 ct / s>
-    >>> obs.countrate(area=area)  # Binned
+    >>> obs.countrate(area=area)  # Binned  # doctest: +FLOAT_CMP
     <Quantity 6628.36888121 ct / s>
 
 
@@ -317,11 +315,11 @@ The default binning behavior is to be consistent with ASTROLIB PYSYNPHOT.
 
 Example::
 
-    >>> obs.effective_wavelength()  # Binned
-    <Quantity 5401.267857308841 Angstrom>
-    >>> obs.effective_wavelength(mode='efflphot')  # Deprecated
+    >>> obs.effective_wavelength()  # Binned  # doctest: +FLOAT_CMP
+    <Quantity 5401.26785017 Angstrom>
+    >>> obs.effective_wavelength(mode='efflphot')  # Deprecated  # doctest: +SKIP
     WARNING: AstropyDeprecationWarning: Usage of EFFLPHOT is deprecated. [...]
-    <Quantity 5424.929868234263 Angstrom>
+    <Quantity 5424.92986116 Angstrom>
 
 
 .. _synphot-formula-pivwv:
@@ -341,7 +339,7 @@ The formula shown applies to a bandpass. For a source, replace
 
 Example::
 
-    >>> bp.pivot()
+    >>> bp.pivot()  # doctest: +FLOAT_CMP
     <Quantity 5355.863596422958 Angstrom>
-    >>> obs.pivot()  # Not binned
+    >>> obs.pivot()  # Not binned  # doctest: +FLOAT_CMP
     <Quantity 5389.368734064575 Angstrom>
diff --git a/docs/synphot/from_pysyn_iraf.rst b/docs/synphot/from_pysyn_iraf.rst
index a0d6141..6081976 100644
--- a/docs/synphot/from_pysyn_iraf.rst
+++ b/docs/synphot/from_pysyn_iraf.rst
@@ -1,5 +1,3 @@
-.. doctest-skip-all
-
 .. _synphot-switcher:
 
 Switching from Legacy Software
diff --git a/docs/synphot/observation.rst b/docs/synphot/observation.rst
index 675904a..53cf947 100644
--- a/docs/synphot/observation.rst
+++ b/docs/synphot/observation.rst
@@ -1,5 +1,3 @@
-.. doctest-skip-all
-
 .. _synphot_observation:
 
 Observation
@@ -22,11 +20,27 @@ cases, the bandpass has an extra ``binset`` attribute (wavelength values for
 bin centers) that must be passed into an observation's constructor, as shown
 below::
 
-    >>> import stsynphot as STS
-    >>> from synphot import Observation
-    >>> sp = STS.Vega
-    >>> bp = STS.band('acs,wfc1,f555w')
-    >>> obs = Observation(sp, bp, binset=bp.binset)
+    >>> import stsynphot as stsyn  # doctest: +SKIP
+    >>> from astropy import units as u
+    >>> from synphot import Observation, SourceSpectrum
+    >>> from synphot.models import GaussianFlux1D
+    >>> sp = SourceSpectrum(GaussianFlux1D, amplitude=1*u.Jy,
+    ...                     mean=6000*u.AA, fwhm=100*u.AA)
+    >>> bp = stsyn.band('acs,hrc,f555w')  # doctest: +SKIP
+    >>> obs = Observation(sp, bp, binset=bp.binset)  # doctest: +SKIP
+
+.. testsetup::
+
+    >>> import os
+    >>> import numpy as np
+    >>> from astropy.utils.data import get_pkg_data_filename
+    >>> from synphot import SpectralElement
+    >>> filename = get_pkg_data_filename(
+    ...     os.path.join('data', 'hst_acs_hrc_f555w.fits'),
+    ...     package='synphot.tests')
+    >>> bp = SpectralElement.from_file(filename)
+    >>> binset = np.arange(1000, 11001, dtype=np.float)
+    >>> obs = Observation(sp, bp, binset=binset)
 
 It has these main general components:
 
@@ -52,18 +66,18 @@ for flux conversion). To get binned flux values, use
 :meth:`~synphot.observation.Observation.sample_binned`, where you must provide
 the exact bin center(s)::
 
-    >>> obs(6000.5)  # Native (not binned) sampling; Angstrom
-    <Quantity 160.95438045640773 PHOTLAM>
-    >>> obs.sample_binned([6000, 6001])  # Binned flux in given centers
-    <Quantity [ 161.19216632, 160.65695961] PHOTLAM>
+    >>> obs(6000.5)  # Native (not binned) sampling; Angstrom  # doctest: +FLOAT_CMP
+    <Quantity 0.03080553 PHOTLAM>
+    >>> obs.sample_binned([6000, 6001])  # Binned flux in given centers  # doctest: +FLOAT_CMP
+    <Quantity [0.03084873, 0.0307453 ] PHOTLAM>
 
 To calculate **bin properties** such as covered wavelength or pixel ranges,
 you can use its :meth:`~synphot.observation.Observation.binned_waverange` and
 :meth:`~synphot.observation.Observation.binned_pixelrange` as follows::
 
     >>> # Wavelength range covered by 10 pixels centered at 5500 Angstrom
-    >>> obs.binned_waverange(5500, 10)
-    <Quantity [ 5495.5, 5505.5] Angstrom>
+    >>> obs.binned_waverange(5500, 10)  # doctest: +FLOAT_CMP
+    <Quantity [5495.5, 5505.5] Angstrom>
     >>> # Pixel range covered by above wavelength range
     >>> obs.binned_pixelrange([5495.5, 5505.5])
     10
@@ -73,14 +87,14 @@ and :ref:`synphot-formula-effwave`, which can be calculated in a way
 that is consistent with ASTROLIB PYSYNPHOT::
 
     >>> # Effective stimulus in FLAM
-    >>> obs.effstim(flux_unit='flam')
-    <Quantity 3.78664384e-09 FLAM>
+    >>> obs.effstim(flux_unit='flam')  # doctest: +FLOAT_CMP
+    <Quantity 4.12569567e-14 FLAM>
     >>> # Effective wavelength for binned sampling in FLAM
-    >>> obs.effective_wavelength()
-    <Quantity 5332.62717945 Angstrom>
+    >>> obs.effective_wavelength()  # doctest: +FLOAT_CMP
+    <Quantity 5991.99794956 Angstrom>
     >>> # Repeat for "native" sampling
-    >>> obs.effective_wavelength(binned=False)
-    <Quantity 5332.62724394 Angstrom>
+    >>> obs.effective_wavelength(binned=False)  # doctest: +FLOAT_CMP
+    <Quantity 5991.99798818 Angstrom>
 
 :meth:`~synphot.observation.Observation.countrate` is probably the most often
 used method for an observation.
@@ -90,8 +104,8 @@ some binning. By default, it uses ``binset``, which should be defined such that
 one wavelength bin corresponds to one detector pixel::
 
     >>> area = 45238.93416  # HST, in cm^2
-    >>> obs.countrate(area)
-    <Quantity 1.9249653e+10 ct / s>
+    >>> obs.countrate(area)  # doctest: +FLOAT_CMP
+    <Quantity 137190.19332899 ct / s>
 
 An observation can be converted to a **regular source spectrum** containing
 only the wavelength set and sampled flux (binned by default) by using its
diff --git a/docs/synphot/overview.rst b/docs/synphot/overview.rst
index b852a20..a555964 100644
--- a/docs/synphot/overview.rst
+++ b/docs/synphot/overview.rst
@@ -1,5 +1,3 @@
-.. doctest-skip-all
-
 .. _synphot_overview:
 
 Overview
diff --git a/docs/synphot/spectrum.rst b/docs/synphot/spectrum.rst
index 8fa8a5a..281ef09 100644
--- a/docs/synphot/spectrum.rst
+++ b/docs/synphot/spectrum.rst
@@ -1,5 +1,3 @@
-.. doctest-skip-all
-
 .. _source-spectrum-main:
 
 Source Spectrum
@@ -43,12 +41,12 @@ for flux conversion)::
     >>> from synphot.models import ConstFlux1D
     >>> sp = SourceSpectrum(ConstFlux1D, amplitude=1)  # PHOTLAM
     >>> wave = [1000, 10000]  # Angstrom
-    >>> sp(wave)
-    <Quantity [ 1., 1.] PHOTLAM>
-    >>> sp(wave, flux_unit=units.FNU)
-    <Quantity [  6.62606957e-24,  6.62606957e-23] FNU>
+    >>> sp(wave)  # doctest: +FLOAT_CMP
+    <Quantity [1., 1.] PHOTLAM>
+    >>> sp(wave, flux_unit=units.FNU)  # doctest: +FLOAT_CMP
+    <Quantity [6.62607015e-24, 6.62607015e-23] FNU>
     >>> area = 45238.93416 * units.AREA  # HST
-    >>> sp(wave, flux_unit=units.OBMAG, area=area)
+    >>> sp(wave, flux_unit=units.OBMAG, area=area)  # doctest: +FLOAT_CMP
     <Quantity [-21.52438718,-21.52438718] OBMAG>
 
 .. _synphot_reddening:
@@ -61,19 +59,19 @@ Astropy model) and then
 extinction curve with a given :math:`E(B-V)` value (negative value effectively
 de-reddens the spectrum), and then multiply it to the source::
 
-    >>> import matplotlib.pyplot as plt
+    >>> import matplotlib.pyplot as plt  # doctest: +SKIP
     >>> from synphot import SourceSpectrum, ReddeningLaw
     >>> from synphot.models import BlackBodyNorm1D
     >>> em = SourceSpectrum(BlackBodyNorm1D, temperature=5000)
     >>> ext = ReddeningLaw.from_extinction_model(
-    ...     'lmcavg').extinction_curve(0.1)
-    >>> sp = em * ext
+    ...     'lmcavg').extinction_curve(0.1)  # doctest: +REMOTE_DATA
+    >>> sp = em * ext  # doctest: +REMOTE_DATA
     >>> wave = em.waveset
-    >>> plt.plot(wave, em(wave), 'b', wave, sp(wave), 'r')
-    >>> plt.xlim(1000, 30000)
-    >>> plt.xlabel('Wavelength (Angstrom)')
-    >>> plt.ylabel('Flux (PHOTLAM)')
-    >>> plt.legend(['E(B-V)=0', 'E(B-V)=0.1'], loc='upper right')
+    >>> plt.plot(wave, em(wave), 'b', wave, sp(wave), 'r')  # doctest: +SKIP
+    >>> plt.xlim(1000, 30000)  # doctest: +SKIP
+    >>> plt.xlabel('Wavelength (Angstrom)')  # doctest: +SKIP
+    >>> plt.ylabel('Flux (PHOTLAM)')  # doctest: +SKIP
+    >>> plt.legend(['E(B-V)=0', 'E(B-V)=0.1'], loc='upper right')  # doctest: +SKIP
 
 .. image:: images/bb5000_lmcavg.png
    :width: 600px
@@ -167,20 +165,20 @@ bandpass using its :meth:`~synphot.spectrum.BaseSourceSpectrum.normalize`
 method. The resultant spectrum is basically the source multiplied with a factor
 necessary to achieve the desired normalization::
 
-    >>> import matplotlib.pyplot as plt
+    >>> import matplotlib.pyplot as plt  # doctest: +SKIP
     >>> from synphot import SourceSpectrum, SpectralElement, units
     >>> from synphot.models import BlackBodyNorm1D
     >>> sp = SourceSpectrum(BlackBodyNorm1D, temperature=5000)
-    >>> bp = SpectralElement.from_filter('johnson_v')
-    >>> vega = SourceSpectrum.from_vega()  # For unit conversion
-    >>> sp_norm = sp.normalize(17 * units.VEGAMAG, bp, vegaspec=vega)
+    >>> bp = SpectralElement.from_filter('johnson_v')  # doctest: +REMOTE_DATA
+    >>> vega = SourceSpectrum.from_vega()  # For unit conversion  # doctest: +REMOTE_DATA
+    >>> sp_norm = sp.normalize(17 * units.VEGAMAG, bp, vegaspec=vega)  # doctest: +REMOTE_DATA
     >>> wave = sp.waveset
-    >>> plt.plot(wave, sp(wave), 'b', wave, sp_norm(wave), 'r')
-    >>> plt.xlim(1000, 30000)
-    >>> plt.xlabel('Wavelength (Angstrom)')
-    >>> plt.ylabel('Flux (PHOTLAM)')
-    >>> plt.title(sp.meta['expr'])
-    >>> plt.legend(['Original', 'Normalized'], loc='upper right')
+    >>> plt.plot(wave, sp(wave), 'b', wave, sp_norm(wave), 'r')  # doctest: +SKIP
+    >>> plt.xlim(1000, 30000)  # doctest: +SKIP
+    >>> plt.xlabel('Wavelength (Angstrom)')  # doctest: +SKIP
+    >>> plt.ylabel('Flux (PHOTLAM)')  # doctest: +SKIP
+    >>> plt.title(sp.meta['expr'])  # doctest: +SKIP
+    >>> plt.legend(['Original', 'Normalized'], loc='upper right')  # doctest: +SKIP
 
 .. image:: images/bb5000_renorm.png
    :width: 600px
@@ -192,13 +190,14 @@ method. It uses trapezoid integration (but could be expanded to perform
 analytical calculations instead in the future when that is supported by
 Astropy). By default, integration is done in internal units::
 
+    >>> from astropy import units as u
     >>> from synphot import SourceSpectrum, units
     >>> from synphot.models import GaussianFlux1D
     >>> sp = SourceSpectrum(GaussianFlux1D, mean=6000, fwhm=10,
     ...                     total_flux=1*u.erg/(u.cm**2 * u.s))
-    >>> sp.integrate()
-    <Quantity 3.02046763e+11 ph / (cm2 s)
-    >>> sp.integrate(flux_unit=units.FLAM)
+    >>> sp.integrate()  # doctest: +FLOAT_CMP
+    <Quantity 3.02046763e+11 ph / (cm2 s)>
+    >>> sp.integrate(flux_unit=units.FLAM)  # doctest: +FLOAT_CMP
     <Quantity 0.99999972 erg / (cm2 s)>
 
 
@@ -469,8 +468,8 @@ It is loaded from ``synphot.conf.vega_file`` using
 The example below loads and plots the built-in Vega spectrum::
 
     >>> from synphot import SourceSpectrum
-    >>> sp = SourceSpectrum.from_vega()
-    >>> sp.plot(right=12000, flux_unit='flam', title=sp.meta['expr'])
+    >>> sp = SourceSpectrum.from_vega()  # doctest: +REMOTE_DATA
+    >>> sp.plot(right=12000, flux_unit='flam', title=sp.meta['expr'])  # doctest: +SKIP
 
 .. image:: images/vega_spec.png
    :width: 600px
diff --git a/docs/synphot/tutorials.rst b/docs/synphot/tutorials.rst
index f0c4dcc..b35c4f5 100644
--- a/docs/synphot/tutorials.rst
+++ b/docs/synphot/tutorials.rst
@@ -1,5 +1,3 @@
-.. doctest-skip-all
-
 .. _synphot-tutorials:
 
 Tutorials
@@ -15,7 +13,7 @@ Emission Line
 -------------
 
 This tutorial is adapted from
-`Exposure Time Calculator User's Guide on a similar topic <http://etc.stsci.edu/etcstatic/users_guide/1_ref_8.5_emlines.html#using-synphot-to-modify-emission-lines>`_.
+`Exposure Time Calculator User's Guide on a similar topic <http://etc.stsci.edu/etcstatic/users_guide/1_ref_8.5_emlines.html#using-pysynphot-to-modify-emission-lines>`_.
 In this tutorial, you will learn how to manipulate and superimpose
 an emission line to a continuum spectrum.
 
@@ -28,6 +26,7 @@ Create a continuum spectrum of a 5500 K blackbody with z=0.6::
 Create a Gaussian emission line with 8E-14 FLAM total flux,
 FWHM of 100 Angstrom, and centered at 7000 Angstrom::
 
+    >>> from astropy import units as u
     >>> from synphot import units
     >>> from synphot.models import GaussianFlux1D
     >>> em = SourceSpectrum(
@@ -43,12 +42,12 @@ to the composite spectrum::
 
     >>> from synphot import ReddeningLaw
     >>> ext = ReddeningLaw.from_extinction_model(
-    ...     'lmcavg').extinction_curve(1.3)
-    >>> my_spec = sp * ext
+    ...     'lmcavg').extinction_curve(1.3)  # doctest: +REMOTE_DATA
+    >>> my_spec = sp * ext  # doctest: +REMOTE_DATA
 
 Plot the result::
 
-    >>> my_spec.plot(right=45000)
+    >>> my_spec.plot(right=45000)  # doctest: +SKIP
 
 .. image:: images/tutorial_em_line.png
    :width: 600px
@@ -114,17 +113,17 @@ Read in the real data. If your own data has a different format, you need to
 adjust the example accordingly::
 
     >>> from astropy.io import fits
-    >>> with fits.open('/path/to/combined_13330_G130M_v40_bin4.fits') as pf:
-    ...     dat = pf[1].data
-    ...     wave = dat.field('WAVELENGTH').flatten()  # Angstrom
-    ...     flux = dat.field('FLUX').flatten()  # FLAM
+    >>> with fits.open('/path/to/combined_13330_G130M_v40_bin4.fits') as pf:  # doctest: +SKIP
+    ...     dat = pf[1].data  # doctest: +SKIP
+    ...     wave = dat.field('WAVELENGTH').flatten()  # Angstrom  # doctest: +SKIP
+    ...     flux = dat.field('FLUX').flatten()  # FLAM  # doctest: +SKIP
 
 For a good fit, only use data around the feature of interest. In this example,
 the feature is between 1202 and 1211 Angstrom::
 
-    >>> mask = (wave >= 1202) & (wave <= 1211)
-    >>> x = wave[mask]
-    >>> y = flux[mask]
+    >>> mask = (wave >= 1202) & (wave <= 1211)  # doctest: +SKIP
+    >>> x = wave[mask]  # doctest: +SKIP
+    >>> y = flux[mask]  # doctest: +SKIP
 
 Create a composite model with some initial parameters close to the desired
 result (usually sufficient to guess from looking at a plot of the data) and
@@ -136,17 +135,17 @@ iterations are required for a good fit)::
     >>> gs = models.Gaussian1D(amplitude=3.5E-14, mean=1206, stddev=1)
     >>> init_model = bg - gs
     >>> fitter = fitting.LevMarLSQFitter()
-    >>> fit_model = fitter(init_model, x, y)
-    >>> y_fit = fit_model(x)
-    >>> print(fit_model)
-    Model: CompoundModel1
+    >>> fit_model = fitter(init_model, x, y)  # doctest: +SKIP
+    >>> y_fit = fit_model(x)  # doctest: +SKIP
+    >>> print(fit_model)  # doctest: +SKIP
+    Model: CompoundModel...
     Inputs: ('x',)
     Outputs: ('y',)
     Model set size: 1
     Expression: [0] - [1]
     Components:
         [0]: <Const1D(amplitude=3.5e-14)>
-
+    <BLANKLINE>
         [1]: <Gaussian1D(amplitude=3.5e-14, mean=1206.0, stddev=1.0)>
     Parameters:
            amplitude_0       amplitude_1        mean_1       stddev_1
@@ -155,15 +154,15 @@ iterations are required for a good fit)::
 
 Plot the fitted model on top of input data::
 
-    >>> import matplotlib.pyplot as plt
-    >>> from matplotlib import ticker
-    >>> fig, ax = plt.subplots()
-    >>> ax.plot(x, y, 'b', x, y_fit, 'r')
+    >>> import matplotlib.pyplot as plt  # doctest: +SKIP
+    >>> from matplotlib import ticker  # doctest: +SKIP
+    >>> fig, ax = plt.subplots()  # doctest: +SKIP
+    >>> ax.plot(x, y, 'b', x, y_fit, 'r')  # doctest: +SKIP
     >>> ax.get_xaxis().set_major_formatter(
-    ...     ticker.FuncFormatter(ticker.FormatStrFormatter('%.0f')))
-    >>> ax.set_xlabel('Wavelength (Angstrom)')
-    >>> ax.set_ylabel('Flux (FLAM)')
-    >>> ax.legend(['Data', 'Fit'], loc='lower right')
+    ...     ticker.FuncFormatter(ticker.FormatStrFormatter('%.0f')))  # doctest: +SKIP
+    >>> ax.set_xlabel('Wavelength (Angstrom)')  # doctest: +SKIP
+    >>> ax.set_ylabel('Flux (FLAM)')  # doctest: +SKIP
+    >>> ax.legend(['Data', 'Fit'], loc='lower right')  # doctest: +SKIP
 
 .. image:: images/tutorial_fit_ab.png
    :width: 600px
@@ -173,9 +172,9 @@ Calculate equivalent width using the fitted model::
 
     >>> import math
     >>> area = (math.sqrt(2 * math.pi) * fit_model.amplitude_1 *
-    ...     fit_model.stddev_1)  # Area under curve
-    >>> height = fit_model.amplitude_0  # Continuum level
-    >>> print('EW = {:.4f} Angstrom'.format(area / height))
+    ...     fit_model.stddev_1)  # Area under curve  # doctest: +SKIP
+    >>> height = fit_model.amplitude_0  # Continuum level  # doctest: +SKIP
+    >>> print('EW = {:.4f} Angstrom'.format(area / height))  # doctest: +SKIP
     EW = 0.5937 Angstrom
 
 
diff --git a/docs/synphot/units.rst b/docs/synphot/units.rst
index e21c989..1ca6a61 100644
--- a/docs/synphot/units.rst
+++ b/docs/synphot/units.rst
@@ -1,5 +1,3 @@
-.. doctest-skip-all
-
 .. _synphot_units:
 
 Units and Constants
@@ -216,15 +214,15 @@ for HST::
 
     >>> from synphot import units
     >>> area = 45238.93416 * units.AREA
-    >>> sp([499, 500, 501, 502] * u.nm, flux_unit='count', area=area)
-    <Quantity [ 219.25781807, 220.04014043, 220.81868137, 221.59342806] ct>
+    >>> sp([499, 500, 501, 502] * u.nm, flux_unit='count', area=area)  # doctest: +FLOAT_CMP
+    <Quantity [219.3785428 , 220.16129659, 220.94026689, 221.71544085] ct>
 
 Sample the source at 5E+15 Hz and obtain flux in FLAM::
 
-    >>> sp([5E+15] * u.Hz, flux_unit=units.FLAM)
-    <Quantity [  3.52271822e-29] FLAM>
+    >>> sp([5E+15] * u.Hz, flux_unit=units.FLAM)  # doctest: +FLOAT_CMP
+    <Quantity [3.52467344e-29] FLAM>
 
 Sample the source in internal units (Angstrom and PHOTLAM)::
 
-    >>> sp(6000)
-    <Quantity 0.0006152610208167509 PHOTLAM>
+    >>> sp(6000)  # doctest: +FLOAT_CMP
+    <Quantity 0.0006156 PHOTLAM>
diff --git a/setup.cfg b/setup.cfg
index c3dc8d5..60c7aec 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -8,9 +8,12 @@ upload-dir = docs/_build/html
 show-response = 1
 
 [tool:pytest]
-minversion = 3.0
+minversion = 3.1
 norecursedirs = build docs/_build synphot/src
-addopts = -p no:warnings
+doctest_plus = enabled
+text_file_format = rst
+addopts = --doctest-rst
+xfail_strict = true
 
 [ah_bootstrap]
 auto_use = True
diff --git a/synphot/models.py b/synphot/models.py
index 51d78ac..23249c4 100644
--- a/synphot/models.py
+++ b/synphot/models.py
@@ -276,7 +276,7 @@ class Empirical1D(Tabular1D):
 
     .. note::
 
-        This model requires `SciPy <http://www.scipy.org>`_ 0.14
+        This model requires `SciPy <https://www.scipy.org>`_ 0.14
         or later to be installed.
 
     Parameters


Test patch:
diff --git a/synphot/tests/test_models.py b/synphot/tests/test_models.py
index 8f25f32..15286aa 100644
--- a/synphot/tests/test_models.py
+++ b/synphot/tests/test_models.py
@@ -12,6 +12,7 @@
 
 # STDLIB
 import os
+import warnings
 
 # THIRD-PARTY
 import numpy as np
@@ -23,6 +24,7 @@ from astropy import units as u
 from astropy.modeling.models import Const1D
 from astropy.utils import minversion
 from astropy.utils.data import get_pkg_data_filename
+from astropy.utils.exceptions import AstropyUserWarning
 
 # LOCAL
 from .. import specio, units
@@ -180,8 +182,12 @@ class TestEmpirical1D(object):
         [(True, [-1.1, 0, 1.1]),
          (False, [0, 0, 1.1])])
     def test_neg_array(self, keep_neg, ans):
-        m2 = Empirical1D(points=[1, 2, 3], lookup_table=[-1.1, 0, 1.1],
-                         keep_neg=keep_neg)
+        with warnings.catch_warnings():
+            warnings.filterwarnings(
+                'ignore', message=r'.*contained negative flux or throughput.*',
+                category=AstropyUserWarning)
+            m2 = Empirical1D(points=[1, 2, 3], lookup_table=[-1.1, 0, 1.1],
+                             keep_neg=keep_neg)
         np.testing.assert_array_equal(m2([1, 2, 3]), ans)
         if not keep_neg:
             assert 'NegativeFlux' in m2.meta['warnings']
@@ -191,8 +197,12 @@ class TestEmpirical1D(object):
         [(True, -1),
          (False, 0)])
     def test_neg_scalar(self, keep_neg, ans):
-        m2 = Empirical1D(points=[1, 2, 3], lookup_table=[-1, 0, 1],
-                         keep_neg=keep_neg)
+        with warnings.catch_warnings():
+            warnings.filterwarnings(
+                'ignore', message=r'.*contained negative flux or throughput.*',
+                category=AstropyUserWarning)
+            m2 = Empirical1D(points=[1, 2, 3], lookup_table=[-1, 0, 1],
+                             keep_neg=keep_neg)
         np.testing.assert_array_equal(m2(1), ans)
         if not keep_neg:
             assert 'NegativeFlux' in m2.meta['warnings']
diff --git a/synphot/tests/test_observation.py b/synphot/tests/test_observation.py
index 000f7d7..dcf9254 100644
--- a/synphot/tests/test_observation.py
+++ b/synphot/tests/test_observation.py
@@ -3,6 +3,7 @@
 
 # STDLIB
 import os
+import warnings
 
 # THIRD-PARTY
 import numpy as np
@@ -14,7 +15,8 @@ from astropy.modeling.models import Const1D
 from astropy.tests.helper import catch_warnings
 from astropy.utils import minversion
 from astropy.utils.data import get_pkg_data_filename
-from astropy.utils.exceptions import AstropyDeprecationWarning
+from astropy.utils.exceptions import (AstropyDeprecationWarning,
+                                      AstropyUserWarning)
 
 # LOCAL
 from .test_units import _area
@@ -144,7 +146,10 @@ class TestObservation(object):
         sp = SourceSpectrum(
             GaussianFlux1D, mean=5000, total_flux=(1 * tf_unit), fwhm=10)
         bp = SpectralElement(Const1D, amplitude=1)
-        obs2 = Observation(sp, bp, force='extrap')
+        with pytest.warns(AstropyUserWarning,
+                          match=r'Source spectrum will be evaluated outside '
+                          r'pre-defined waveset'):
+            obs2 = Observation(sp, bp, force='extrap')
         np.testing.assert_array_equal(obs2.binset, sp.waveset)
 
     def test_undefined_binset(self):
@@ -177,7 +182,11 @@ class TestInitWithForce(object):
         [('taper', 0),
          ('extrap', 0.75)])
     def test_force(self, force_type, ans):
-        obs = Observation(self.sp, self.bp, force=force_type)
+        with warnings.catch_warnings():
+            warnings.filterwarnings(
+                'ignore', message=r'.*Source spectrum.*',
+                category=AstropyUserWarning)
+            obs = Observation(self.sp, self.bp, force=force_type)
         assert obs(4005).value == ans
         assert force_type in obs.warnings['PartialOverlap']
 
@@ -387,8 +396,11 @@ class TestCountRate(object):
          ([999, 1016], 172.75)])
     def test_force(self, w, ans):
         """Force calculation for partial overlap."""
-        np.testing.assert_allclose(
-            self.obs.countrate(_area, waverange=w, force=True).value, ans)
+        with pytest.warns(AstropyUserWarning, match=r'Count rate calculated '
+                          r'only for wavelengths in the overlap between '
+                          r'observation and given range'):
+            np.testing.assert_allclose(
+                self.obs.countrate(_area, waverange=w, force=True).value, ans)
 
         # Must raise error without force
         with pytest.raises(exceptions.PartialOverlap):
@@ -412,7 +424,11 @@ class TestCountRateNegFlux(object):
         [(True, 1510.219531414891),
          (False, 1627.8250215634343)])
     def test_neg_handling(self, keep_neg, ans):
-        sp = SourceSpectrum.from_file(self.spfile, keep_neg=keep_neg)
+        with warnings.catch_warnings():
+            warnings.filterwarnings(
+                'ignore', message=r'.*contained negative flux or throughput.*',
+                category=AstropyUserWarning)
+            sp = SourceSpectrum.from_file(self.spfile, keep_neg=keep_neg)
         obs = Observation(sp, self.bp)
         c = obs.countrate(_area)
         np.testing.assert_allclose(c.value, ans, rtol=1e-4)
diff --git a/synphot/tests/test_reddening.py b/synphot/tests/test_reddening.py
index 6b277a5..7fca735 100644
--- a/synphot/tests/test_reddening.py
+++ b/synphot/tests/test_reddening.py
@@ -15,6 +15,7 @@ from astropy import units as u
 from astropy.io import fits
 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
@@ -85,8 +86,9 @@ class TestExtinction(object):
             os.path.join('data', 'qso_template_001.dat')))
         extcurve = self.redlaw.extinction_curve(1.0 * u.mag)
         spext = qso * extcurve
-        sp = spext.normalize(25 * u.STmag, bp)
-        obs = Observation(sp, bp, force='taper')
+        with pytest.warns(AstropyUserWarning):
+            sp = spext.normalize(25 * u.STmag, bp)
+            obs = Observation(sp, bp, force='taper')
         area = 45238.93416  # HST cm^2
         c = obs.countrate(area)
         ans = 1.104404103799421e-07  # From similar setup in Astrolib PySynphot
diff --git a/synphot/tests/test_specio.py b/synphot/tests/test_specio.py
index b1fd4d5..ed4f5cb 100644
--- a/synphot/tests/test_specio.py
+++ b/synphot/tests/test_specio.py
@@ -15,6 +15,7 @@ from astropy import units as u
 from astropy.io import fits
 from astropy.tests.helper import assert_quantity_allclose
 from astropy.utils.data import get_pkg_data_filename
+from astropy.utils.exceptions import AstropyUserWarning
 
 # LOCAL
 from .. import exceptions, specio, units
@@ -67,11 +68,12 @@ class TestReadWriteFITS(object):
         outfile = os.path.join(self.outdir, 'outspec1.fits')
 
         # Write it out
-        specio.write_fits_spec(
-            outfile, self.wave.value, self.flux.value, pri_header=self.prihdr,
-            ext_header=self.scihdr, trim_zero=False, pad_zero_ends=False,
-            precision='single', wave_unit=self.wave.unit,
-            flux_unit=self.flux.unit)
+        with pytest.warns(AstropyUserWarning, match=r'rows are thrown out'):
+            specio.write_fits_spec(
+                outfile, self.wave.value, self.flux.value,
+                pri_header=self.prihdr, ext_header=self.scihdr,
+                trim_zero=False, pad_zero_ends=False, precision='single',
+                wave_unit=self.wave.unit, flux_unit=self.flux.unit)
 
         # Read it back in and check values (flux_unit should be ignored)
         hdr, wave, flux = specio.read_spec(outfile, flux_unit='foo')
diff --git a/synphot/tests/test_spectrum.py b/synphot/tests/test_spectrum.py
index 5d7bf9f..d23c574 100644
--- a/synphot/tests/test_spectrum.py
+++ b/synphot/tests/test_spectrum.py
@@ -5,6 +5,7 @@
 import os
 import shutil
 import tempfile
+import warnings
 
 # THIRD-PARTY
 import numpy as np
@@ -20,6 +21,7 @@ from astropy.modeling.models import (
     PowerLaw1D, RedshiftScaleFactor)
 from astropy.utils import minversion
 from astropy.utils.data import get_pkg_data_filename
+from astropy.utils.exceptions import AstropyUserWarning
 
 # LOCAL
 from .test_units import _area, _wave, _flux_jy, _flux_photlam, _flux_vegamag
@@ -167,8 +169,10 @@ class TestEmpiricalSourceFromFile(object):
 
     def test_neg_flux(self):
         w = [1000, 5000, 9000]
-        sp = SourceSpectrum(
-            Empirical1D, points=w, lookup_table=[100, -45, 5e-17])
+        with pytest.warns(AstropyUserWarning,
+                          match=r'contained negative flux or throughput'):
+            sp = SourceSpectrum(
+                Empirical1D, points=w, lookup_table=[100, -45, 5e-17])
         np.testing.assert_array_equal(sp(w).value, [100, 0, 5e-17])
         assert 'NegativeFlux' in sp.warnings
 
@@ -646,7 +650,11 @@ class TestNormalize(object):
     def _compare_countrate(self, rn_sp, ans_countrate):
         # Observation is needed to compare with expected count rate
         # although it is tested in test_observation.py
-        obs = Observation(rn_sp, self.acs, force='extrap')
+        with warnings.catch_warnings():
+            warnings.filterwarnings(
+                'ignore', message=r'.*Source spectrum will be evaluated '
+                r'outside pre-defined waveset.*', category=AstropyUserWarning)
+            obs = Observation(rn_sp, self.acs, force='extrap')
         ct_rate = obs.countrate(_area)
 
         # 0.1% agreement with IRAF SYNPHOT COUNTRATE
@@ -706,7 +714,9 @@ class TestNormalize(object):
         """Test force=True for 'partial_notmost' overlap."""
         sp = SourceSpectrum(Empirical1D, points=[5000, 6000],
                             lookup_table=[1, 1])
-        rn_sp = sp.normalize(1e-23 * u.Jy, band=self.acs, force=True)
+        with pytest.warns(AstropyUserWarning,
+                          match=r'Spectrum is not defined everywhere'):
+            rn_sp = sp.normalize(1e-23 * u.Jy, band=self.acs, force=True)
         assert 'PartialRenorm' in rn_sp.warnings
         assert 'PartialRenorm' not in sp.warnings
 
@@ -717,7 +727,9 @@ class TestNormalize(object):
     def test_renorm_partial_most(self):
         """Test 'partial_most' overlap."""
         bp = SpectralElement(Box1D, amplitude=1, x_0=5600, width=870)
-        rn_sp = self.em.normalize(1e-23 * u.Jy, band=bp)
+        with pytest.warns(AstropyUserWarning,
+                          match=r'Spectrum is not defined everywhere'):
+            rn_sp = self.em.normalize(1e-23 * u.Jy, band=bp)
         assert 'PartialRenorm' in rn_sp.warnings
         assert 'PartialRenorm' not in self.em.warnings
         assert '99%' in rn_sp.warnings['PartialRenorm']