Apply the filter to the given [W, F], or [W, F, E] spectrum

        Parameters
        ----------
        spectrum: array-like
            The wavelength [um] and flux of the spectrum
            to apply the filter to
        plot: bool
            Plot the original and filtered spectrum

        Returns
        -------
        np.ndarray
            The filtered spectrum and error

def apply(self, spectrum, plot=False):
        """
        Apply the filter to the given [W, F], or [W, F, E] spectrum

        Parameters
        ----------
        spectrum: array-like
            The wavelength [um] and flux of the spectrum
            to apply the filter to
        plot: bool
            Plot the original and filtered spectrum

        Returns
        -------
        np.ndarray
            The filtered spectrum and error

        """
        # Convert to filter units if possible
        f_units = 1.
        if hasattr(spectrum[0], 'unit'):
            spectrum[0] = spectrum[0].to(self.wave_units)
        if hasattr(spectrum[1], 'unit'):
            spectrum[1] = spectrum[1].to(self.flux_units)
            f_units = self.flux_units
        if len(spectrum) >= 3 and hasattr(spectrum[2], 'unit'):
            spectrum[2] = spectrum[2].to(self.flux_units)

        # Make into iterable arrays
        wav, flx, *err = [np.asarray(i) for i in spectrum]

        # Check for error array
        if len(err) == 0:
            err = np.ones_like(flx)*np.nan
            unc = False
        else:
            err = err[0]
            unc = True

        # Make flux 2D
        if len(flx.shape) == 1:
            flx = np.expand_dims(flx, axis=0)
            err = np.expand_dims(err, axis=0)

        # Make throughput 3D
        rsr = np.copy(self.rsr)

        # Make empty filtered arrays
        filtered_flx = np.zeros((rsr.shape[0], flx.shape[0], rsr.shape[2]))
        filtered_err = np.zeros_like(filtered_flx)

        # Rebin the input spectra to the filter wavelength array
        # and apply the RSR curve to the spectrum
        for i, bn in enumerate(rsr):
            for j, (f, e) in enumerate(zip(flx, err)):
                filtered_flx[i][j] = np.interp(bn[0], wav, f, left=np.nan, right=np.nan)*bn[1]
                filtered_err[i][j] = np.interp(bn[0], wav, e, left=np.nan, right=np.nan)*bn[1]

        # Propagate the filter systematic uncertainties
        if unc:
            filtered_err += filtered_flx*self.systematics

        if plot:

            # Make the figure
            COLORS = color_gen('Category10')
            xlab = 'Wavelength [{}]'.format(self.wave_units)
            ylab = 'Flux Density [{}]'.format(self.flux_units)
            fig = figure(title=self.filterID, x_axis_label=xlab, y_axis_label=ylab)

            # Plot the unfiltered spectrum
            fig.line(wav, flx[0], legend='Input spectrum', color='black')
 
            # Plot the uncertainties
            if unc:
                band_x = np.append(wav, wav[::-1])
                band_y = np.append(flx-err, (flx+err)[::-1])
                fig.patch(band_x, band_y, color='black', fill_alpha=0.1, line_alpha=0)

            # Plot each spectrum bin
            for wav, bn, bne in zip(self.wave, filtered_flx, filtered_err):
                color = next(COLORS)
                fig.line(wav, bn[0], color=color)

                # Plot the uncertainties
                if unc:
                    band_x = np.append(wav, wav[::-1])
                    band_y = np.append(bn[0]-bne[0], (bn[0]+bne[0])[::-1])
                    fig.patch(band_x, band_y, color=color, fill_alpha=0.1, line_alpha=0)

            show(fig)

        return filtered_flx.squeeze()*f_units, filtered_err.squeeze()*f_units