Merge branches 'emiss_line' and 'master' of https://github.com/lineto…

…ols/linetools into emiss_line

docs/xspectrum1d.rst

@@ -69,7 +69,6 @@ Here are a series of example calls to generate the class::
     >>> sp = XSpectrum1D.from_file('PH957_f.fits')      # From a FITS file
     >>> sp = XSpectrum1D.from_file('q0002m422.txt.gz')  # From an ASCII table
     >>> sp = xspec1.copy()                              # From an XSpectrum1D object
-    >>> sp = XSpectrum1D.from_list(xspec1, xspec2)      # From a list of XSpectrum1D objects
     >>> sp = XSpectrum1D.from_tuple((wa, fl, sig), verbose=False)
 
 There are a number of methods to write a file, e.g.
@@ -88,6 +87,13 @@ data format of XSpectrum1D.  Here are some examples::
     >>> sp.write('QSO.ascii')                   # Same
 
 
+One can collate a list of XSpectrum1D objects into one with collate::
+
+    >>> sp1 = XSpectrum1D.from_file('PH957_f.fits')
+    >>> sp2 = XSpectrum1D.from_file('q0002m422.txt.gz')
+    >>> sp = linetools.spectra.utils.collate([sp1,sp2])
+
+
 Plotting
 --------
 
@@ -109,6 +115,11 @@ correlated::
 
     >>> newspec = sp.rebin(new_wv, do_sig=True) # doctest: +SKIP
 
+If the XSpectrum1D object containts multiple spectra, you can rebin
+all of them to the new wavelength array as well::
+
+    >>> newspec = sp.rebin(new_wv, do_sig=True, all=True) # doctest: +SKIP
+
 
 Continuum fitting
 -----------------

linetools/spectra/tests/test_xspec_io.py

@@ -20,7 +20,8 @@ def spec2():
 
 @pytest.fixture
 def specm(spec,spec2):
-    return XSpectrum1D.from_list([spec,spec2])
+    specm = ltsu.collate([spec,spec2])
+    return specm
 
 def data_path(filename):
     data_dir = os.path.join(os.path.dirname(__file__), 'files')

linetools/spectra/tests/test_xspec_utils.py

@@ -20,7 +20,8 @@ def spec2():
 
 @pytest.fixture
 def specm(spec,spec2):
-    return XSpectrum1D.from_list([spec,spec2])
+    specm = ltsu.collate([spec,spec2])
+    return specm
 
 def data_path(filename):
     data_dir = os.path.join(os.path.dirname(__file__), 'files')
@@ -40,7 +41,7 @@ def test_rebin_to_rest(spec,spec2):
     rest_spec = ltsu.rebin_to_rest(coll_spec, zarr, 100*u.km/u.s, debug=False)
     # Test
     assert rest_spec.totpix == 3716
-    np.testing.assert_allclose(rest_spec.wvmin.value, 986.021738877745, rtol=1e-5)
+    np.testing.assert_allclose(rest_spec.wvmin.value, 986.3506403, rtol=1e-5)
 
 
 def test_smash_spectra(spec,spec2):
@@ -54,7 +55,7 @@ def test_smash_spectra(spec,spec2):
     stack = ltsu.smash_spectra(rest_spec, method='average')
     # Test
     assert stack.totpix == 3716
-    np.testing.assert_allclose(stack.flux[0].value, -1.19753563, rtol=1e-5)
+    np.testing.assert_allclose(stack.flux[1].value, -3.3213510, rtol=1e-5)
 
 
 def test_airtovac_andback(spec):
@@ -72,7 +73,7 @@ def test_airtovac_andback(spec):
     assert spec.meta['airvac'] == 'air'
 
 
-def test_rebin(spec):
+def test_rebin(spec, specm):
     # Rebin
     new_wv = np.arange(3000., 9000., 5) * u.AA
     newspec = spec.rebin(new_wv)
@@ -81,18 +82,15 @@ def test_rebin(spec):
     assert newspec.flux.unit is u.dimensionless_unscaled
     # With sigma
     newspec = spec.rebin(new_wv, do_sig=True)
-    """
-    i1 = np.argmin(np.abs(spec.wavelength-5000.*u.AA))
-    s2n_1 = spec.flux[i1] / spec.sig[i1]
-    i2 = np.argmin(np.abs(newspec.wavelength-5000.*u.AA))
-    s2n_2 = newspec.flux[i2] / newspec.sig[i2]
-    """
     imn = np.argmin(np.abs(newspec.wavelength-8055*u.AA))
     np.testing.assert_allclose(newspec.sig[imn].value, 0.0169634, rtol=1e-5)
     # With NANs
     spec.data['flux'][spec.select][100:110] = np.nan
     newspec = spec.rebin(new_wv)
     np.testing.assert_allclose(newspec.flux[1000], 0.9999280967617779)
+    # All
+    spec2 = specm.rebin(new_wv, all=True)
+    np.testing.assert_allclose(spec2.wvmax.value, 8995.0)
 
 
 def test_addnoise(spec):

linetools/spectra/utils.py

@@ -5,6 +5,7 @@
 # Import libraries
 import numpy as np
 import json
+import warnings
 import pdb
 
 from astropy import units as u
@@ -103,12 +104,18 @@ def splice_two(spec1, spec2, wvmx=None, scale=1., chk_units=True):
 
 
 def collate(spectra):
-    """ Pack list of XSpectrum1D files into one object
-    Meta data will be scrubbed..
+    """ Generate a single XSpectrum1D instance containing an array of
+    spectra from a list of individual XSpectrum1D spectra.
+    Each spectrum is padded with extra pixels so that the
+    wavelength ranges of all spectra are covered.
+    Padded pixels are masked.
+
+    Also note that masked pixels in the original data are ignored!
 
     Parameters
     ----------
-    spectra : list of XSpectrum1D
+    spectra : list
+      of XSpectrum1D
 
     Returns
     -------
@@ -118,39 +125,178 @@ def collate(spectra):
     from linetools.spectra.xspectrum1d import XSpectrum1D
     # Init
     maxpix = 0
-    flg_co = False
+    flg_co, flg_sig = False, False
     nspec = len(spectra)
     units = None
     for spec in spectra:
         maxpix = max(maxpix, spec.npix)
         if spec.co_is_set:
             flg_co = True
+        if spec.sig_is_set:
+            flg_sig = True
         # Check for identical units
         if units is None:
             units = spec.units
         else:
             assert spec.units == units
     # Generate data arrays
-    wave = np.zeros((nspec, maxpix))
-    flux = np.zeros_like(wave)
-    sig = np.zeros_like(wave)
+    wave = np.zeros((nspec, maxpix), dtype='float64')
+    flux = np.zeros_like(wave, dtype='float32')
+    if flg_sig:
+        sig = np.zeros_like(wave, dtype='float32')
+    else:
+        sig = None
     if flg_co:
         co = np.zeros_like(wave)
     else:
         co = None
     # Load
+    meta = dict(headers=[])
     for ii, spec in enumerate(spectra):
         wave[ii,:spec.npix] = spec.wavelength.value
         flux[ii,:spec.npix] = spec.flux.value
-        sig[ii,:spec.npix] = spec.sig.value
+        if flg_sig:
+            sig[ii,:spec.npix] = spec.sig.value
         if flg_co:
             co[ii,:spec.npix] = spec.co.value
+        # Meta
+        meta['headers'].append(spec.header)
     # Finish
-    new_spec = XSpectrum1D(wave, flux, sig=sig, co=co, units=units.copy())
-    #new_spec.meta = spec.meta.copy()
+    new_spec = XSpectrum1D(wave, flux, sig=sig, co=co, units=units.copy(),
+                           masking='edges', meta=meta)
     # Return
     return new_spec
 
+def rebin(spec, new_wv, do_sig=False, all=False, **kwargs):
+    """ Rebin a single spectrum in an XSpectrum1D object to a new wavelength array
+
+    Uses simple linear interpolation.  The default (and only)
+    option conserves counts (and flambda).
+
+    WARNING: Do not trust either edge pixel of the new array.
+      In fact the sig is set to 0 for each of these
+    Also be aware that neighboring pixels are likely to be
+    correlated in a manner that is not described by the error
+    array.
+
+    Parameters
+    ----------
+    new_wv : Quantity array
+      New wavelength array
+    do_sig : bool, optional
+      Rebin error too (if it exists).
+      S/N is only crudely conserved.
+      Rejected pixels are propagated.
+    all : bool, optional
+      Rebin all spectra in the XSpectrum1D object?
+
+    Returns
+    -------
+    XSpectrum1D of the rebinned spectrum
+    """
+    from linetools.spectra.xspectrum1d import XSpectrum1D
+    from scipy.interpolate import interp1d
+    # Save flux info to avoid unit issues
+    funit = spec.flux.unit
+    flux = spec.flux.value
+
+    # Deal with nan
+    badf = np.isnan(flux)
+    if np.sum(badf) > 0:
+        warnings.warn("Ignoring NAN in flux")
+    gdf = ~badf
+    flux = flux[gdf]
+
+    # Endpoints of original pixels
+    npix = len(spec.wavelength)
+    wvh = (spec.wavelength + np.roll(spec.wavelength, -1)) / 2.
+    wvh[npix - 1] = spec.wavelength[npix - 1] + \
+                    (spec.wavelength[npix - 1] - spec.wavelength[npix - 2]) / 2.
+    dwv = wvh - np.roll(wvh, 1)
+    dwv[0] = 2 * (wvh[0] - spec.wavelength[0])
+    med_dwv = np.median(dwv.value)
+
+    wvh = wvh[gdf]
+    dwv = dwv[gdf]
+
+    # Error
+    if do_sig:
+        var = spec.sig.value**2
+        var = var[gdf]
+    else:
+        var = np.ones_like(flux)
+
+    # Cumulative Sum
+    cumsum = np.cumsum(flux * dwv)
+    cumvar = np.cumsum(var * dwv)
+
+    # Interpolate (loses the units)
+    fcum = interp1d(wvh, cumsum, fill_value=0., bounds_error=False)
+    fvar = interp1d(wvh, cumvar, fill_value=0., bounds_error=False)
+
+    # Endpoints of new pixels
+    nnew = len(new_wv)
+    nwvh = (new_wv + np.roll(new_wv, -1)) / 2.
+    nwvh[nnew - 1] = new_wv[nnew - 1] + \
+                     (new_wv[nnew - 1] - new_wv[nnew - 2]) / 2.
+    # Pad starting point
+    bwv = np.zeros(nnew + 1) * new_wv.unit
+    bwv[0] = new_wv[0] - (new_wv[1] - new_wv[0]) / 2.
+    bwv[1:] = nwvh
+
+    # Evaluate and put unit back
+    newcum = fcum(bwv) * dwv.unit
+    newvar = fvar(bwv) * dwv.unit
+
+    # Endpoint
+    if (bwv[-1] > wvh[-1]):
+        newcum[-1] = cumsum[-1]
+        newvar[-1] = cumvar[-1]
+
+    # Rebinned flux, var
+    new_fx = (np.roll(newcum, -1) - newcum)[:-1]
+    new_var = (np.roll(newvar, -1) - newvar)[:-1]
+
+    # Normalize (preserve counts and flambda)
+    new_dwv = bwv - np.roll(bwv, 1)
+    #import pdb
+    # pdb.set_trace()
+    new_fx = new_fx / new_dwv[1:]
+    # Preserve S/N (crudely)
+    med_newdwv = np.median(new_dwv.value)
+    new_var = new_var / (med_newdwv/med_dwv) / new_dwv[1:]
+
+    # Return new spectrum
+    if do_sig:
+        # Create new_sig
+        new_sig = np.zeros_like(new_var)
+        gd = new_var > 0.
+        new_sig[gd] = np.sqrt(new_var[gd].value)
+        # Deal with bad pixels
+        bad = np.where(var <= 0.)[0]
+        for ibad in bad:
+            bad_new = np.where(np.abs(new_wv-spec.wavelength[ibad]) <
+                               (new_dwv[1:]+dwv[ibad])/2)[0]
+            new_sig[bad_new] = 0.
+        # Zero out edge pixels -- not to be trusted
+        new_sig[0] = 0.
+        new_sig[-1] = 0.
+    else:
+        new_sig = None
+
+    # update continuum
+    if spec.co_is_set:
+        x, y = spec._get_contpoints()
+        new_co = spec._interp_continuum(x, y, new_wv)
+    else:
+        new_co = None
+
+    # Finish
+    newspec = XSpectrum1D.from_tuple((new_wv, new_fx*funit,
+                                      new_sig, new_co),
+                                     meta=spec.meta.copy(), **kwargs)
+    # Return
+    return newspec
 
 def rebin_to_rest(spec, zarr, dv, debug=False):
     """ Shuffle an XSpectrum1D dataset to an array of