Continuum Fitting with Python using specutils
At the end I settled on using a median filter from scipy (medfilt) and with a high enough kernel it removed the spikes completely.
The problem on my specutils implementation was that I needed to exclude the peaks in the spectrum window by doing this, althought I found this too hard to implement on 100-120 data measurement I had:
...
from specutils import SpectralRegion

spectrum = Spectrum1D(flux=y*u.Jy, spectral_axis=x*u.um)
g1_fit = fit_generic_continuum(spectrum, exclude_regions=[SpectralRegion(2 * u.um, 6 * u.um), SpectralRegion(58 * u.um, 65 * u.um)])
y_fit = g1_fit(x*u.um)