1- using Spectra: Spectra, AbstractSpectrum, SpectrumResampler, spectrum, wave, flux
1+ using Spectra: Spectra, AbstractSpectrum, SpectrumResampler, spectrum, spectral_axis, flux_axis
22using DataInterpolations: LinearInterpolation, ExtrapolationType
33using Unitful: @u_str, uconvert
44using UnitfulAstro
@@ -7,28 +7,16 @@ using Measurements: ±
77# TODO : See if it makes sense to have an exported API for resample/resample! even though we are using external interpolators.
88
99function resample(spec, new_wave)
10- interp = LinearInterpolation(flux (spec), wave (spec); extrapolation = ExtrapolationType. Constant)
10+ interp = LinearInterpolation(flux_axis (spec), spectral_axis (spec); extrapolation = ExtrapolationType. Constant)
1111 resampler = SpectrumResampler(spec, interp)
1212 return resampler(new_wave)
1313end
1414
1515function resample(spec1:: AbstractSpectrum , spec2:: AbstractSpectrum )
16- new_wave = wave (spec2)
16+ new_wave = spectral_axis (spec2)
1717 return resample(spec1, new_wave)
1818end
1919
20- # function resample!(spec::AbstractSpectrum, new_wave)
21- # spec_new = resample(spec, new_wave)
22- # spec.flux = flux(spec_new)
23- # spec.wave = wave(spec_new)
24- # return spec
25- # end
26-
27- # function resample!(spec1::AbstractSpectrum, spec2::AbstractSpectrum)
28- # new_wave = wave(spec2)
29- # resample!(spec1, new_wave)
30- # end
31-
3220function mock_spectrum(n:: Int = Int(1e3 ); use_units:: Bool = false )
3321 wave = range(1e4 , 3e4 , length = n)
3422 sigma = 0.1 .* sin.(wave)
4331
4432@testset " Resampler" begin
4533 spec = mock_spectrum()
46- interp = LinearInterpolation(flux(spec), wave(spec); extrapolation = ExtrapolationType. Constant)
34+ s, f = spectral_axis(spec), flux_axis(spec)
35+ interp = LinearInterpolation(f, s; extrapolation = ExtrapolationType. Constant)
4736 resampler = SpectrumResampler(spec, interp)
4837 expected = """
4938 SpectrumResampler(Float64, Measurements.Measurement{Float64})
50- spec: Spectrum(Float64, Measurements.Measurement{Float64})
51- name: Test Spectrum
52- interpolator: DataInterpolations.LinearInterpolation{Vector{Measurements.Measurement{Float64}}, Vector{Float64}, Vector{Measurements.Measurement{Float64}}, Vector{Measurements.Measurement{Float64}}, Measurements.Measurement{Float64}}"""
39+ spec: Spectra.SingleSpectrum{Float64, Measurements.Measurement{Float64}}
40+ interpolator: DataInterpolations.LinearInterpolation{Vector{Measurements.Measurement{Float64}}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, Vector{Measurements.Measurement{Float64}}, Vector{Measurements.Measurement{Float64}}, Measurements.Measurement{Float64}}"""
5341
5442 @test sprint(show, resampler) == expected
55- @test wave (resampler) == wave(spec)
56- @test flux (resampler) == flux(spec)
43+ @test spectral_axis (resampler) == s
44+ @test flux_axis (resampler) == f
5745end
5846
5947@testset " Resampling" begin
6048 spec = mock_spectrum()
61- new_wave = range(minimum(spec. wave), maximum(spec. wave); length = Integer(length(spec. wave) ÷ 2.4 ))
49+ s, f = spectral_axis(spec), flux_axis(spec)
50+ new_wave = range(minimum(s), maximum(s); length = Integer(length(s) ÷ 2.4 ))
6251 res_spec = resample(spec, new_wave)
6352 expected = """
64- Spectrum(Float64, Measurements.Measurement{Float64})
65- name: Test Spectrum"""
53+ SingleSpectrum(Float64, Measurements.Measurement{Float64})
54+ spectral axis (416,): 10000.0 .. 30000.0
55+ flux axis (416,): 67.0 ± 0.031 .. 0.827 ± 0.08
56+ meta: Dict{Symbol, Any}(:name => "Test Spectrum")"""
6657
6758 @test sprint(show, res_spec) == expected
68- @test wave(res_spec) == new_wave
69- @test length(flux(res_spec)) == length(new_wave)
70-
71- # resample!(spec, new_wave)
72- # @test wave(res_spec) == spec.wave
73- # @test flux(res_spec) == spec.flux
59+ @test spectral_axis(res_spec) == new_wave
60+ @test length(flux_axis(res_spec)) == length(new_wave)
7461
7562 # Unitful
76- spec = mock_spectrum(use_units= true )
77- new_wave = range(minimum(spec. wave), maximum(spec. wave), length= Integer(length(spec. wave) ÷ 2.4 ))
78- @assert unit(eltype(new_wave)) == unit(eltype(spec. wave))
63+ spec = mock_spectrum(; use_units = true )
64+ s, f = spectral_axis(spec), flux_axis(spec)
65+ new_wave = range(minimum(s), maximum(s); length = Integer(length(s) ÷ 2.4 ))
66+ @assert unit(eltype(new_wave)) == unit(eltype(s))
7967 res_spec = resample(spec, new_wave)
8068
81- @test wave(res_spec) == new_wave
82- @test length(flux(res_spec)) == length(new_wave)
83-
84- # resample!(spec, new_wave)
85- # @test wave(res_spec) == spec.wave
86- # @test flux(res_spec) == spec.flux
69+ @test spectral_axis(res_spec) == new_wave
70+ @test length(flux_axis(res_spec)) == length(new_wave)
8771
8872 # Test resampling to another Spectrum
8973 spec1 = mock_spectrum(Integer(1e4 ))
9074 spec2 = mock_spectrum(Integer(1e3 ))
9175 res_spec = resample(spec1, spec2)
92- @test wave(res_spec) == wave(spec2)
93-
94- # resample!(spec1, spec2)
95- # @test wave(spec1) == spec2.wave
96- # @test flux(spec1) == flux(res_spec)
76+ @test spectral_axis(res_spec) == spectral_axis(spec2)
9777
9878 # Unitful
99- spec1 = mock_spectrum(Integer(1e4 ), use_units= true )
100- spec2 = mock_spectrum(Integer(1e3 ), use_units= true )
79+ spec1 = mock_spectrum(Integer(1e4 ); use_units = true )
80+ spec2 = mock_spectrum(Integer(1e3 ); use_units = true )
10181 res_spec = resample(spec1, spec2)
10282
103- @test wave (res_spec) == spec2. wave
104- @test unit(eltype(flux (spec1))) == unit(eltype(flux (spec2))) == unit(eltype(flux (res_spec)))
83+ @test spectral_axis (res_spec) == spectral_axis( spec2)
84+ @test unit(eltype(flux_axis (spec1))) == unit(eltype(flux_axis (spec2))) == unit(eltype(flux_axis (res_spec)))
10585
10686 # Test when spectra2 has different units
10787 # spec1 = mock_spectrum(Integer(1e4), use_units=true)
10888 # spec2 = mock_spectrum(Integer(1e3), use_units=true)
10989 # spec1.wave = uconvert.(u"cm", spec1.wave)
11090 # resample!(spec1, spec2)
111- # @test ustrip.(wave (spec1)) ≈ ustrip.(unit(eltype(wave (spec1))), wave (spec2))
91+ # @test ustrip.(spectral_axis (spec1)) ≈ ustrip.(unit(eltype(spectral_axis (spec1))), spectral_axis (spec2))
11292
11393 # address bug when doing the same thing but affecting spec2
11494 # spec1 = mock_spectrum(Integer(1e4), use_units=true)
11595 # spec2 = mock_spectrum(Integer(1e3), use_units=true)
11696 # spec2.wave = uconvert.(u"cm", spec2.wave)
11797 # resample!(spec1, spec2)
118- # @test ustrip.(wave (spec1)) ≈ ustrip.(unit(eltype(wave (spec1))), wave (spec2))
98+ # @test ustrip.(spectral_axis (spec1)) ≈ ustrip.(unit(eltype(spectral_axis (spec1))), spectral_axis (spec2))
11999end
0 commit comments