docs: 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