Initialization


In [1]:
%pylab inline


Populating the interactive namespace from numpy and matplotlib

In [2]:
import os
import matplotlib.gridspec as gridspec
from scipy.interpolate import interp1d
from astropy.io import fits

Configuration


In [3]:
# top level dir containg boss data, this notebook assumes a directory layout organized as follows: 
#   '<boss_root>/deltas/delta-<thing_id>.fits'
#   '<boss_root>/v5_7_0/spectra/lite/<plate>/spec-<plate>-<mjd>-<fiber>.fits'
boss_root = '/Users/daniel/data/boss'

BOSS Data Objects


In [4]:
class BOSSSpectrum:
    """
    Represents a BOSS spectrum
    """
    def __init__(self, loglam, flux, ivar):
        self.loglam = loglam
        self.flux = flux
        self.ivar = ivar

    def combine_pixels(self, num_combine, offset=None, ivar_weighted=True):
        """
        Combines groups of consecutive pixels.

        Args:
            num_combine (int): number of consecutive pixels to combine
            offset (None, int): pixel index to start groupings at (useful for aligning wavelength grids)
            ivar_weighted (bool): indicate whether or not to combine flux using an ivar weighted average

        Returns:
            A new BOSSSpectrum object with combined pixels of self.
        """
        # Find offset relative to fiducial wavelength grid
        if offset is None:
            fiducial_offset = self.get_fiducial_pixel_index_offset(self.loglam[0])
            offset = fiducial_offset % num_combine
        
        # Trim pixel arrays to multiple of num_combine
        num_pixels = len(self.loglam)
        num_leftover = num_pixels % num_combine
        if offset > 0 or num_leftover > 0:
            end_pixel = num_pixels - ((num_pixels - offset) % num_combine)
            trimmed = self.select_pixels(slice(offset, end_pixel))
        else:
            trimmed = self
        
        # Perform combination as a "reduction" of reshaped array dimension
        loglam = trimmed.loglam.reshape((-1, num_combine))
        flux = trimmed.flux.reshape((-1, num_combine))
        ivar = trimmed.ivar.reshape((-1, num_combine))
                
        loglam = np.average(loglam, axis=1)
        if ivar_weighted:
            flux = np.average(flux, weights=ivar, axis=1)
        else:
            flux = np.average(flux, axis=1)
        ivar = np.sum(ivar, axis=1)
                    
        return BOSSSpectrum(loglam, flux, ivar)
    
    def busca_combine(self, num_combine, z, forest_min=1040, forest_max=1200,
                      observed_min=3600, observed_max=5500, loglam_spacing=1e-4):
        
        loglam_min_rest = np.log10(forest_min)
        loglam_max_rest = np.log10(forest_max)
        
        loglam_min = np.log10(observed_min)
        loglam_max = np.log10(observed_max)
        
        dloglam = num_combine*loglam_spacing
        
        # calculate start and end wavelengths
        loglam_start = loglam_min_rest + np.log10(1+z)
        if loglam_start < loglam_min:
            loglam_start = loglam_min
            
        loglam_stop = loglam_max_rest + np.log10(1+z)
        if loglam_stop > loglam_max:
            loglam_stop = loglam_max
            
        # shift to coadd grid
        combined_bin_start = np.around((loglam_start-self.loglam[0])/dloglam).astype(int)
        loglam_start = self.loglam[0] + combined_bin_start*dloglam
        
        combined_bin_stop = np.around((loglam_stop-self.loglam[0])/dloglam).astype(int)
        loglam_stop = self.loglam[0] + combined_bin_stop*dloglam

        # calculate number of bins on combined grid
        num_combined_bins = np.around((loglam_stop-loglam_start)/dloglam).astype(int)+1
        
        #loglam = loglam_start + np.arange(num_combined_bins)*num_combined_bins
                
        # calculate new combined bin indices of uncombinned pixels
        bin_indices = np.around((self.loglam-loglam_start)/dloglam).astype(int);
        
        # trim spectrum to pixels that will go into combined bins
        bin_min = np.argmax(bin_indices >= 0)
        bin_max = np.argmax(bin_indices >= num_combined_bins)      
        trimmed = self.select_pixels(slice(bin_min, bin_max))
        
        # combine pixels
        return trimmed.combine_pixels(num_combine, offset=0)

    def select_pixels(self, selection):
        """
        Returns a copy of self after applying the selection provided.
        """
        return BOSSSpectrum(self.loglam[selection], self.flux[selection], self.ivar[selection])
    
    def get_forest_spectrum(self, z, forest_min=1040.0, forest_max=1200.0):
        """
        Returns a copy of self after selecting pixels in the lyman alpha forest
        """
        
        forest_min *= 1+z
        forest_max *= 1+z
        
        forest_pixels = np.where( \
            (self.loglam >= np.log10(forest_min)) & \
            (self.loglam <= np.log10(forest_max)) )[0]
        
        return self.select_pixels(forest_pixels)
    
    @staticmethod
    def get_fiducial_pixel_index_offset(coeff0, coeff1=1e-4, fid_lambda=3500.26):
        """
        Returns the pixel index offset from the start of the BOSS co-add fiducial wavelength grid.
        If more than 1% away from a fiducial pixel index, returns the exact offset (in pixels).

        Args:
            coeff0 (float): central wavelength (log10) of first pixel
            coeff1 (float, optional): log10 spacing between pixel centers

        Returns:
            pixel index offset from the start of the fiducial wavelength grid.
        """
        delta = (np.log10(fid_lambda)-coeff0)/coeff1
        # round to nearest pixel index
        offset = np.floor(delta+0.5).astype(int)
        # if more than %1 away from pixel index, return exact offset (in pixels)
        if np.any(np.fabs(delta-offset) > 0.01):
            return -delta
        # otherwise return nearest pixel index
        return -offset

In [11]:
class DeltaLOS():
    """
    Represents a delta file (delta field along a line of sight)
    """
    def __init__(self, thing_id):
        self.thing_id = thing_id

        hdulist = fits.open(os.path.join(boss_root,'deltas','delta-%s.fits' % thing_id)) 
        data = hdulist[1].data
        self.loglam = data['loglam']
        self.delta = data['delta']
        self.weight = data['weight']
        self.r_comov = data['r_comov']
        self.cont = data['cont']
        self.msha = data['msha']
        self.mabs = data['mabs']
        self.ivar = data['ivar']

        self.flux = (1+self.delta)*self.cont*self.msha*self.mabs

        header = hdulist[1].header
        self.z = header['z']
        self.plate = header['plate']
        self.mjd = header['mjd']
        self.fiber = header['fiberid']    
        self.ra = header['ra']
        self.dec = header['dec']

In [18]:
class LiteSpec():
    """
    Represents a lite spec file (single target's observation, co-add only)
    """
    def __init__(self, plate, mjd, fiber):
        self.plate = plate
        self.mjd = mjd
        self.fiber = fiber
        
        filename = 'spec-%s-%s-%04d.fits' % (plate, mjd, fiber)
        hdulist = fits.open(os.path.join(boss_root,'v5_7_0','spectra','lite',str(plate),filename))
        data = hdulist[1].data
        self.loglam = data['loglam']
        self.flux = data['flux']
        self.ivar = data['ivar']
        self.ra = hdulist[2].data['ra']
        self.dec = hdulist[2].data['dec']

Compare Binning


In [19]:
def compare_binning(thing_id, num_combine=3, offset=None):
    # Load delta los data and raw BOSS spectrum
    delta_los = DeltaLOS(thing_id)
    print 'RA, Dec: %f, %f' % (delta_los.ra, delta_los.dec)
    lite_spec = LiteSpec(delta_los.plate, delta_los.mjd, delta_los.fiber)
    print 'RA, Dec: %f, %f' % (lite_spec.ra, lite_spec.dec)
    boss_spectrum = BOSSSpectrum(lite_spec.loglam, lite_spec.flux, lite_spec.ivar)
    print 'Plate, mjd, fiber:', delta_los.plate, delta_los.mjd, delta_los.fiber
    
    # combine pixels from raw spectrum into "analysis" pixels
    analysis_spectrum = boss_spectrum.combine_pixels(num_combine, offset=offset, ivar_weighted=True)
    forest_spectrum = analysis_spectrum.get_forest_spectrum(delta_los.z)
    
    busca_spectrum = boss_spectrum.busca_combine(num_combine, delta_los.z)
    test_spectrum = busca_spectrum#.get_forest_spectrum(delta_los.z)

    # plot results 
    
    # plot forest/spec limits 
    z = delta_los.z
    forest_min = 1040*(1+z)
    forest_max = 1200*(1+z)
    observed_min = 10**boss_spectrum.loglam[0]
    observed_max = 10**boss_spectrum.loglam[-1]
    
    fig = plt.figure(figsize=(14, 6))
    fig.suptitle('THING_ID = %s, Z_VI = %.2f' % (thing_id, z))
    gs = gridspec.GridSpec(2, 2, height_ratios=[2,1], width_ratios=[1,1], hspace=0.1, wspace=0.05)

    ax1 = plt.subplot(gs[0])

    # plot the raw boss spectrum
    plt.plot(10**boss_spectrum.loglam, boss_spectrum.flux, 
        color='gray', label='Raw Spectrum', lw=0, marker='+')
    
    # plot busca forest
    s = test_spectrum.ivar > 0
    plt.scatter(10**test_spectrum.loglam[s], test_spectrum.flux[s],
        color='purple', marker='o', facecolors='none', label='Busca Pixels')
    
    # plot the delta los spectrum
    s = delta_los.ivar > 0
    plt.errorbar(10**delta_los.loglam[s], delta_los.flux[s], yerr=1/np.sqrt(delta_los.ivar)[s],
        color='blue', marker='x', fmt='o', alpha=.7, label='FPG Pixels')
    # plot the combined and trimmed forest
    s = forest_spectrum.ivar > 0
    plt.errorbar(10**forest_spectrum.loglam[s], forest_spectrum.flux[s], yerr=1/np.sqrt(forest_spectrum.ivar[s]),
        color='red', marker='x', fmt='o', alpha=.7, label='IVAR Weighted Pixels')

        
    plt.axvline(forest_min, color='black', linestyle='--')
    plt.axvline(forest_max, color='black', linestyle='--')
    plt.axvline(observed_min, color='black', linestyle='-.')
    plt.grid()

    plt.ylabel('Observed Flux')

    # hide x-axis ticklabels
    plt.setp(ax1.get_xticklabels(), visible=False)
    
    ax2 = plt.subplot(gs[1], sharey=ax1)
    
    # plot the raw boss spectrum
    plt.plot(10**boss_spectrum.loglam, boss_spectrum.flux, 
        color='gray', label='Raw Spectrum', lw=0, marker='+')
    
    # plot busca forest
    s = test_spectrum.ivar > 0
    plt.scatter(10**test_spectrum.loglam[s], test_spectrum.flux[s],
        color='purple', marker='o', facecolors='none', label='Busca Pixels')
    
    # plot the delta los spectrum
    s = delta_los.ivar > 0
    plt.errorbar(10**delta_los.loglam[s], delta_los.flux[s], yerr=1/np.sqrt(delta_los.ivar)[s],
        color='blue', marker='x', fmt='o', alpha=.7, label='FPG Pixels')
    # plot the combined and trimmed forest
    s = forest_spectrum.ivar > 0
    plt.errorbar(10**forest_spectrum.loglam[s], forest_spectrum.flux[s], yerr=1/np.sqrt(forest_spectrum.ivar[s]),
        color='red', marker='x', fmt='o', alpha=.7, label='IVAR Weighted Pixels')

    plt.axvline(forest_max, color='black', linestyle='--')
    plt.grid()
    
    plt.ylim(np.min(forest_spectrum.flux-1/np.sqrt(forest_spectrum.ivar))*0.9,
             np.max(forest_spectrum.flux+1/np.sqrt(forest_spectrum.ivar))*1.1)

    plt.setp(ax2.get_xticklabels(), visible=False)
    plt.setp(ax2.get_yticklabels(), visible=False)

    # switch to lower subplot, share x axis
    ax3 = plt.subplot(gs[2], sharex=ax1)

    # create interpolation of delta field fluxes so we can evaulate using grid of our combined pixels
    delta_los_interp = interp1d(delta_los.loglam, delta_los.flux, bounds_error=False)
    residual_flux = forest_spectrum.flux - delta_los_interp(forest_spectrum.loglam)

    # plot residual
    plt.step(10**forest_spectrum.loglam, residual_flux, 
        color='black', where='mid', lw=.5)
    plt.step(10**test_spectrum.loglam, test_spectrum.flux - delta_los_interp(test_spectrum.loglam),
        color='purple', where='mid', lw=.5)

    # plot forest/spec limits 
    plt.axvline(forest_min, color='black', linestyle='--')
    plt.axvline(forest_max, color='black', linestyle='--')
    plt.axvline(observed_min, color='black', linestyle='-.')
    plt.grid()

    plot_min = max(forest_min, observed_min)-10
    plot_max = min(plot_min+80, forest_max+20)

    plt.xlim(plot_min, plot_max)

    plt.xlabel('Observed Wavelength ($\AA$)')
    plt.ylabel('Residual Flux')
    
    ax4 = plt.subplot(gs[3], sharex=ax2, sharey=ax3)
    
    # plot residual
    plt.step(10**forest_spectrum.loglam, residual_flux, 
        color='black', where='mid', lw=.5)
    plt.step(10**test_spectrum.loglam, test_spectrum.flux - delta_los_interp(test_spectrum.loglam),
        color='purple', where='mid', lw=.5)

    # plot forest/spec limits 
    plt.axvline(forest_min, color='black', linestyle='--')
    plt.axvline(forest_max, color='black', linestyle='--')
    plt.axvline(observed_min, color='black', linestyle='-.')
    plt.grid()

    plt.xlim(forest_max+10 - (plot_max - plot_min), forest_max+10)

    plt.xlabel('Observed Wavelength ($\AA$)')
    
    plt.setp(ax4.get_yticklabels(), visible=False)

In [20]:
compare_binning(66692618, num_combine=3)


RA, Dec: 0.000100, -0.023126
RA, Dec: 0.005746, -1.325009
Plate, mjd, fiber: 4354 55810 646

In [17]:
compare_binning(401440280, num_combine=3)


RA, Dec: 6.283173, 0.605993
Plate, mjd, fiber: 7145 56567 240

In [10]:
compare_binning(401440280, num_combine=3, offset=2)


Plate, mjd, fiber: 7145 56567 240

In [10]: