Merge branch 'master' into loaders-fix-tabular

.travis.yml

@@ -90,29 +90,24 @@ matrix:
       env: SETUP_CMD='test' EVENT_TYPE='cron'
 
     # Do a coverage test.
-    - os: linux
-      stage: Initial tests
+    - stage: Initial tests
       env: SETUP_CMD='test --coverage'
 
     # Check for sphinx doc build warnings - we do this first because it
     # may run for a long time
-    - os: linux
-      env: SETUP_CMD='build_docs -w'
+    - env: SETUP_CMD='build_docs -w'
         CONDA_DEPENDENCIES=$CONDA_DEPENDENCIES_DOC
 
     # check on somewhat older versions LTS version
-    - os: linux
-      env: PYTHON_VERSION=3.6 ASTROPY_VERSION=lts NUMPY_VERSION=1.17
+    - env: PYTHON_VERSION=3.6 ASTROPY_VERSION=lts NUMPY_VERSION=1.17
 
     # Add a job that runs from cron only and tests against astropy dev and
     # numpy dev to give a change for early discovery of issues and feedback
     # for both developer teams.
-    - os: linux
-      stage: Cron tests
+    - stage: Cron tests
       env: ASTROPY_VERSION=development NUMPY_VERSION=development
         EVENT_TYPE='cron'
-    - os: linux
-      env: ASTROPY_VERSION=development
+    - env: ASTROPY_VERSION=development
         EVENT_TYPE='pull_request push cron'
 
     # Try on Windows.
@@ -123,26 +118,24 @@ matrix:
     # Note we are not doing any numpy versions older than 1.16 since those are not compatible with Astropy 4.0
 
     # Do a PEP8 test with flake8
-    # - os: linux
-    #   stage: Initial tests
+    # - stage: Initial tests
     #   env: MAIN_CMD='flake8 specutils --count --show-source --statistics $FLAKE8_OPT' SETUP_CMD=''
     # Do a PEP8 test with pycodestyle
     - stage: Initial tests
-      os: linux
       env: MAIN_CMD='pycodestyle specutils --count'
            SETUP_CMD=''
 
     # Try development version of Astropy and GWCS
-    - os: linux
-      env: ASTROPY_VERSION=dev EVENT_TYPE='pull_request push cron'
+    - env: ASTROPY_VERSION=dev EVENT_TYPE='pull_request push cron'
         PIP_DEPENDENCIES="$GWCS_DEV $ASDF_DEV scipy matplotlib"
 
+    - env: PIP_DEPENDENCIES="$PIP_DEPENDENCIES git+https://github.com/spacetelescope/jwst@stable"
+
   allow_failures:
   # Do a PEP8 test with flake8
   # (do allow to fail unless your code completely compliant)
-   - os: linux
-     stage: Initial tests
-     env: MAIN_CMD='flake8 specutils --count --show-source --statistics $FLAKE8_OPT' SETUP_CMD=''
+   # - stage: Initial tests
+   #   env: MAIN_CMD='flake8 specutils --count --show-source --statistics $FLAKE8_OPT' SETUP_CMD=''
 
 install:
 

CHANGES.rst

@@ -18,15 +18,47 @@ Bug Fixes
 1.0
 ---
 
-API Changes
-^^^^^^^^^^^
-
 New Features
 ^^^^^^^^^^^^
 
+- Implement ``SpectralCoord`` object. [#524] 
+
+- Implement cross-correlation for finding redshift/radial velocity. [#544]
+
+- Improve FITS file identification in default_loaders. [#545]
+
+- Support ``len()`` for ``SpectrumCollection`` objects. [#575]
+
+- Improved 1D JWST loader and allow parsing into an ``SpectrumCollection`` object. [#579]
+
+- Implemented 2D and 3D data loaders for JWST products. [#595] 
+
+- Include documentation on how to use dust_extinction in specutils. [#594]
+
+- Include example of spectrum shifting in docs. [#600]
+
+- Add new default excise_regions exciser function and improve subregion handling. [#609]
+
+- Implement use of ``SpectralCoord`` in ``Spectrum1D`` objects. [#610]
+
 Bug Fixes
 ^^^^^^^^^
 
+- Fix stacking and unit treatment in ``SpectrumCollection.from_spectra``. [#578]
+
+- Fix spectral axis unit retrieval. [#581]
+
+- Fix bug in subspectrum fitting. [#586] 
+
+- Fix uncertainty to weight conversion to match astropy assumptions. [#594]
+
+- Fix warnings and log messages from ASDF serialization tests. [#597]
+
+Other Changes and Additions
+^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+- Remove spectral_resolution stub from Spectrum1D. [#606]
+
 
 0.7
 ---

docs/analysis.rst

@@ -20,7 +20,7 @@ sub-sections below - a gaussian-profile line with flux of 5 GHz Jy.  See
     >>> from astropy.modeling import models
     >>> from specutils import Spectrum1D, SpectralRegion
     >>> np.random.seed(42)
-    >>> spectral_axis = np.linspace(0., 10., 200) * u.GHz
+    >>> spectral_axis = np.linspace(11., 1., 200) * u.GHz
     >>> spectral_model = models.Gaussian1D(amplitude=5*(2*np.pi*0.8**2)**-0.5*u.Jy, mean=5*u.GHz, stddev=0.8*u.GHz)
     >>> flux = spectral_model(spectral_axis)
     >>> flux += np.random.normal(0., 0.05, spectral_axis.shape) * u.Jy
@@ -43,29 +43,29 @@ spectrum:
     >>> from specutils.analysis import snr
     >>> snr(noisy_gaussian)  # doctest:+FLOAT_CMP
     <Quantity 2.47730726>
-    >>> snr(noisy_gaussian, SpectralRegion(4*u.GHz, 6*u.GHz))  # doctest:+FLOAT_CMP
-    <Quantity 9.84136331>
+    >>> snr(noisy_gaussian, SpectralRegion(6*u.GHz, 4*u.GHz))  # doctest:+FLOAT_CMP
+    <Quantity 9.8300873>
 
 
 A second method to calculate SNR does not require the uncertainty defined
-on the `~specutils.Spectrum1D` object. This computes the signal to noise 
-ratio DER_SNR following the definition set forth by the Spectral 
+on the `~specutils.Spectrum1D` object. This computes the signal to noise
+ratio DER_SNR following the definition set forth by the Spectral
 Container Working Group of ST-ECF, MAST and CADC. This is based on the
 code at http://www.stecf.org/software/ASTROsoft/DER_SNR/.
 
 .. code-block:: python
 
     >>> from specutils.analysis import snr_derived
     >>> snr_derived(noisy_gaussian)  # doctest:+FLOAT_CMP
-    <Quantity 1.0448650809155666>
-    >>> snr_derived(noisy_gaussian, SpectralRegion(4*u.GHz, 6*u.GHz))  # doctest:+FLOAT_CMP
-    <Quantity 40.53528208761309>
+    <Quantity 1.13359867>
+    >>> snr_derived(noisy_gaussian, SpectralRegion(6*u.GHz, 4*u.GHz))  # doctest:+FLOAT_CMP
+    <Quantity 42.10020601>
 
-The conditions on the data for this implementation for it to be an unbiased estimator of the SNR 
+The conditions on the data for this implementation for it to be an unbiased estimator of the SNR
 are strict.  In particular:
 
   * the noise is uncorrelated in wavelength bins spaced two pixels apart
-  * for large wavelength regions, the signal over the scale of 5 or more pixels can be approximated by a straight line 
+  * for large wavelength regions, the signal over the scale of 5 or more pixels can be approximated by a straight line
 
 
 Line Flux Estimates
@@ -91,9 +91,9 @@ of a spectrum.  Both are demonstrated below:
 
     >>> from specutils.analysis import line_flux
     >>> line_flux(noisy_gaussian).to(u.erg * u.cm**-2 * u.s**-1)  # doctest:+FLOAT_CMP
-    <Quantity 4.97826405e-14 erg / (cm2 s)>
-    >>> line_flux(noisy_gaussian, SpectralRegion(3*u.GHz, 7*u.GHz))  # doctest:+FLOAT_CMP
-    <Quantity 4.92933252 GHz Jy>
+    <Quantity 4.97826284e-14 erg / (cm2 s)>
+    >>> line_flux(noisy_gaussian, SpectralRegion(7*u.GHz, 3*u.GHz))  # doctest:+FLOAT_CMP
+    <Quantity 4.93254052 GHz Jy>
 
 For the equivalent width, note the need to add a continuum level:
 
@@ -102,9 +102,9 @@ For the equivalent width, note the need to add a continuum level:
     >>> from specutils.analysis import equivalent_width
     >>> noisy_gaussian_with_continuum = noisy_gaussian + 1*u.Jy
     >>> equivalent_width(noisy_gaussian_with_continuum)  # doctest:+FLOAT_CMP
-    <Quantity -4.97826405 GHz>
-    >>> equivalent_width(noisy_gaussian_with_continuum, regions=SpectralRegion(3*u.GHz, 7*u.GHz))  # doctest:+FLOAT_CMP
-    <Quantity -4.92933252 GHz>
+    <Quantity -4.97826284 GHz>
+    >>> equivalent_width(noisy_gaussian_with_continuum, regions=SpectralRegion(7*u.GHz, 3*u.GHz))  # doctest:+FLOAT_CMP
+    <Quantity -4.93254052 GHz>
 
 
 Centroid
@@ -116,8 +116,8 @@ estimate the center of a spectral feature:
 .. code-block:: python
 
     >>> from specutils.analysis import centroid
-    >>> centroid(noisy_gaussian, SpectralRegion(3*u.GHz, 7*u.GHz))  # doctest:+FLOAT_CMP
-    <Quantity 4.99881315 GHz>
+    >>> centroid(noisy_gaussian, SpectralRegion(7*u.GHz, 3*u.GHz))  # doctest:+FLOAT_CMP
+    <Quantity 4.99909151 GHz>
 
 While this example is "pre-subtracted", this function only performs well if the
 contiuum has already been subtracted, as for the other functions above and
@@ -155,16 +155,16 @@ Each of the width analysis functions are applied to this spectrum below:
 
    >>> from specutils.analysis import gaussian_sigma_width, gaussian_fwhm, fwhm, fwzi
    >>> gaussian_sigma_width(noisy_gaussian) # doctest: +FLOAT_CMP
-   <Quantity 0.76925064 GHz>
+   <Quantity 0.74075431 GHz>
    >>> gaussian_fwhm(noisy_gaussian) # doctest: +FLOAT_CMP
-   <Quantity 1.81144683 GHz>
+   <Quantity 1.74434311 GHz>
    >>> fwhm(noisy_gaussian) # doctest: +FLOAT_CMP
-   <Quantity 1.91686888 GHz>
+   <Quantity 1.86047666 GHz>
    >>> fwzi(noisy_gaussian) # doctest: +FLOAT_CMP
-   <Quantity 89.99999455 GHz>
+   <Quantity 94.99997484 GHz>
 
 
-Template Comparison
+Template comparison
 -------------------
 
 The ~`specutils.analysis.template_comparison.template_match` function takes an
@@ -225,6 +225,63 @@ An example of how to do template matching with an unknown redshift is:
    >>> tm_result = template_comparison.template_match(observed_spectrum=observed_spectrum, spectral_templates=spectral_template, resample_method=resample_method, redshift=rs_values) # doctest:+FLOAT_CMP
 
 
+Dust extinction
+---------------
+
+Dust extinction can be applied to Spectrum1D instances via their internal arrays, using
+the ``dust_extinction`` package (http://dust-extinction.readthedocs.io/en/latest)
+
+Below is an example of how to apply extinction.
+
+.. code-block:: python
+
+    from astropy.modeling.blackbody import blackbody_lambda
+    from dust_extinction.parameter_averages import F99
+
+    wave = np.logspace(np.log10(1000), np.log10(3e4), num=10) * u.AA
+    flux = blackbody_lambda(wave, 10000 * u.K)
+    spec = Spectrum1D(spectral_axis=wave, flux=flux)
+
+    # define the model
+    ext = F99(Rv=3.1)
+
+    # extinguish (redden) the spectrum
+    flux_ext = spec.flux * ext.extinguish(spec.spectral_axis, Ebv=0.5)
+    spec_ext = Spectrum1D(spectral_axis=wave, flux=flux_ext)
+
+
+Template Cross-correlation
+--------------------------
+
+The cross-correlation function between an observed spectrum and a template spectrum that both share a common spectral
+axis can be calculated with the function `~template_correlate` in the `~specutils.analysis` module.
+
+An example of how to get the cross correlation follows. Note that the observed spectrum must have a rest wavelength
+value set.
+
+.. code-block:: python
+
+    >>> from specutils.analysis import correlation
+    >>> size = 200
+    >>> spec_axis = np.linspace(4500., 6500., num=size) * u.AA
+    >>> f1 = np.random.randn(size)*0.5 * u.Jy
+    >>> f2 = np.random.randn(size)*0.5 * u.Jy
+    >>> rest_value = 6000. * u.AA
+    >>> mean1 = 5035. * u.AA
+    >>> mean2 = 5015. * u.AA
+    >>> g1 = models.Gaussian1D(amplitude=30 * u.Jy, mean=mean1, stddev=10. * u.AA)
+    >>> g2 = models.Gaussian1D(amplitude=30 * u.Jy, mean=mean2, stddev=10. * u.AA)
+    >>> flux1 = f1 + g1(spec_axis)
+    >>> flux2 = f2 + g2(spec_axis)
+    >>> uncertainty = StdDevUncertainty(0.2*np.ones(size)*u.Jy)
+    >>> ospec = Spectrum1D(spectral_axis=spec_axis, flux=flux1, uncertainty=uncertainty, velocity_convention='optical', rest_value=rest_value)
+    >>> tspec = Spectrum1D(spectral_axis=spec_axis, flux=flux2, uncertainty=uncertainty)
+    >>> corr, lag = correlation.template_correlate(ospec, tspec)
+
+The lag values are reported in km/s units. The correlation values are computed after the template spectrum is
+normalized in order to have the same total flux as the observed spectrum.
+
+
 Reference/API
 -------------
 .. automodapi:: specutils.analysis

docs/manipulation.rst

@@ -251,7 +251,7 @@ known uncertainty:
 
     >>> from astropy.modeling import models
     >>> np.random.seed(42)
-    >>> spectral_axis = np.linspace(0., 10., 200) * u.GHz
+    >>> spectral_axis = np.linspace(10, 1, 200) * u.GHz
     >>> spectral_model = models.Gaussian1D(amplitude=3*u.Jy, mean=5*u.GHz, stddev=0.8*u.GHz)
     >>> flux = spectral_model(spectral_axis)
     >>> flux += np.random.normal(0., 0.2, spectral_axis.shape) * u.Jy
@@ -264,10 +264,10 @@ the line:
 
     >>> from specutils import SpectralRegion
     >>> from specutils.manipulation import noise_region_uncertainty
-    >>> noise_region = SpectralRegion([(0, 3), (7, 10)]*u.GHz)
+    >>> noise_region = SpectralRegion([(10, 7), (3, 0)] * u.GHz)
     >>> spec_w_unc = noise_region_uncertainty(noisy_gaussian, noise_region)
     >>> spec_w_unc.uncertainty # doctest: +ELLIPSIS
-    StdDevUncertainty([0.18461457, ..., 0.18461457])
+    StdDevUncertainty([0.18823157, ..., 0.18823157])
 
 Or similarly, expressed in pixels:
 
@@ -276,7 +276,7 @@ Or similarly, expressed in pixels:
     >>> noise_region = SpectralRegion([(0, 25), (175, 200)]*u.pix)
     >>> spec_w_unc = noise_region_uncertainty(noisy_gaussian, noise_region)
     >>> spec_w_unc.uncertainty # doctest: +ELLIPSIS
-    StdDevUncertainty([0.18714535, ..., 0.18714535])
+    StdDevUncertainty([0.18739524, ..., 0.18739524])
 
 S/N Threshold Mask
 ------------------
@@ -315,6 +315,49 @@ and True elsewhere. It is this way to be consistent with ``astropy.nddata``.
 .. note:: The mask attribute is the only attribute modified by ``snr_threshold()``. To
              retrieve the masked flux data use ``spectrum.masked.flux_masked``.
 
+Shifting
+--------
+
+In addition to resampling, you may sometimes wish to simply shift the
+``spectral_axis`` of a spectrum (a la the ``specshift`` iraf task).
+There is no explicit function for this because it is a basic transform of
+the ``spectral_axis``. Therefore one can use a construct like this:
+
+.. code-block:: python
+
+    >>> from specutils import Spectrum1D
+    >>> np.random.seed(42)
+    >>> wavelengths = np.arange(0, 10) * u.um
+    >>> flux = 100 * np.abs(np.random.randn(10)) * u.Jy
+    >>> spectrum = Spectrum1D(spectral_axis=wavelengths, flux=flux)
+    >>> spectrum  #doctest:+ELLIPSIS
+    <Spectrum1D(flux=<Quantity [ 49.6714153 ,  13.82643012,  64.76885381, 152.30298564,
+                23.41533747,  23.41369569, 157.92128155,  76.74347292,
+                46.94743859,  54.25600436] Jy>,
+     spectral_axis=<SpectralCoord [0., 1., 2., 3., 4., 5., 6., 7., 8., 9.] um,
+       radial_velocity=0.0 km / s,
+       redshift=0.0,
+       doppler_rest=0.0 Angstrom,
+       doppler_convention=None,
+       observer=None,
+       target=None>)>
+
+    >>> shift = 12300 * u.AA
+    >>> new_spec = Spectrum1D(spectral_axis=spectrum.spectral_axis + shift, flux=spectrum.flux)
+    >>> new_spec #doctest:+ELLIPSIS
+    <Spectrum1D(flux=<Quantity [ 49.6714153 ,  13.82643012,  64.76885381, 152.30298564,
+                23.41533747,  23.41369569, 157.92128155,  76.74347292,
+                46.94743859,  54.25600436] Jy>, spectral_axis=<SpectralCoord [ 1.23,  2.23,  3.23,  4.23,  5.23,  6.23,  7.23,  8.23,
+                     9.23, 10.23] um,
+       radial_velocity=0.0 km / s,
+       redshift=0.0,
+       doppler_rest=0.0 Angstrom,
+       doppler_convention=None,
+       observer=None,
+       target=None>)>
+
+
+
 Reference/API
 -------------