spec1 = Spectrum1D(spectral_axis=np.arange(1, 50) * u.micron, flux=np.random.randn(49)*u.Jy)
spec2 = Spectrum1D(spectral_axis=np.arange(51, 100) * u.micron, flux=(np.random.randn(49)+1)*u.Jy)

new_spectral_axis = np.concatenate([spec1.spectral_axis.value, spec2.spectral_axis.to_value(spec1.spectral_axis.unit)]) * spec1.spectral_axis.unit

resampler = LinearInterpolatedResampler(extrapolation_treatment='zero_fill')
new_spec1 = resampler(spec1, new_spectral_axis)
new_spec2 = resampler(spec2, new_spectral_axis)

final_spec = new_spec1 + new_spec2

# Yielding a spliced spectrum (the solid line below) composed of the splice of
# two other spectra (dashed lines)::

f, ax = plt.subplots()  # doctest: +IGNORE_OUTPUT
ax.step(final_spec.spectral_axis, final_spec.flux, where='mid', c='k', lw=2) # doctest: +IGNORE_OUTPUT
ax.step(spec1.spectral_axis, spec1.flux, ls='--', where='mid', lw=1) # doctest: +IGNORE_OUTPUT
ax.step(spec2.spectral_axis, spec2.flux, ls='--', where='mid', lw=1) # doctest: +IGNORE_OUTPUT