Skip to content

Commit

Permalink
try to use new resampling in drp
Browse files Browse the repository at this point in the history
  • Loading branch information
@ndrory
ndrory committed Oct 8, 2024
1 parent 915a537 commit 634c0b5
Show file tree
Hide file tree
Showing 4 changed files with 168 additions and 109 deletions.
81 changes: 73 additions & 8 deletions examples/lco_com/testrebin.ipynb
View file

Large diffs are not rendered by default.

80 changes: 36 additions & 44 deletions python/lvmdrp/core/resample.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,7 @@
from scipy.sparse import csc_array
from scipy.signal import correlate
from scipy.signal.windows import tukey
from scipy import interpolate

def _normalize_for_template_matching(s1, s2):
"""
Expand Down Expand Up @@ -333,29 +334,6 @@ def template_logwl_resample(spectrum, template, wblue=None, wred=None,
# Build the corresponding wavelength array
wave_array = numpy.power(10., log_wave_array) * spectrum.spectral_axis.unit

# Resample spectrum and template into wavelength array so built
resampled_spectrum = resampler(spectrum, wave_array)
resampled_template = resampler(template, wave_array)

# Resampler leaves Nans on flux bins that aren't touched by it.
# We replace with zeros. This has the net effect of zero-padding
# the spectrum and/or template so they exactly match each other,
# wavelengthwise.
clean_spectrum_flux = numpy.nan_to_num(resampled_spectrum.flux.value) * resampled_spectrum.flux.unit
clean_template_flux = numpy.nan_to_num(resampled_template.flux.value) * resampled_template.flux.unit

clean_spectrum = Spectrum1D(spectral_axis=resampled_spectrum.spectral_axis,
flux=clean_spectrum_flux,
uncertainty=resampled_spectrum.uncertainty,
velocity_convention='optical',
rest_value=spectrum.rest_value)
clean_template = Spectrum1D(spectral_axis=resampled_template.spectral_axis,
flux=clean_template_flux,
uncertainty=resampled_template.uncertainty,
velocity_convention='optical',
rest_value=template.rest_value)

return clean_spectrum, clean_template

def resample_flux(xout, x, flux, ivar=None, extrapolate=False):
"""Returns a flux conserving resampling of an input flux density.
Expand Down Expand Up @@ -459,21 +437,36 @@ def _unweighted_resample(output_x,input_x,input_flux_density, extrapolate=False)
flux_out(j) = int_{x>(x_{j-1}+x_j)/2}^{x<(x_j+x_{j+1})/2} y(x) dx / 0.5*(x_{j+1}+x_{j-1})
"""
def interpolate_mask(x, y, mask, kind="linear", fill_value=0):
if not numpy.any(mask):
return y
if numpy.all(mask):
return y
known_x, known_v = x[~mask], y[~mask]
missing_x = x[mask]
missing_idx = numpy.where(mask)

f = interpolate.interp1d(known_x, known_v, kind=kind, fill_value=fill_value, bounds_error=False)
yy = y.copy()
yy[missing_idx] = f(missing_x)

return yy

# shorter names
ix=input_x
iy=input_flux_density
ox=output_x
ix = input_x
ox = output_x

iy = interpolate_mask(input_x, input_flux_density, ~numpy.isfinite(input_flux_density))

# boundary of output bins
bins=numpy.zeros(ox.size+1)
bins[1:-1]=(ox[:-1]+ox[1:])/2.
bins[0]=1.5*ox[0]-0.5*ox[1] # = ox[0]-(ox[1]-ox[0])/2
bins[-1]=1.5*ox[-1]-0.5*ox[-2] # = ox[-1]+(ox[-1]-ox[-2])/2
bins = numpy.zeros(ox.size+1)
bins[1:-1] = (ox[:-1]+ox[1:])/2.
bins[0] = 1.5*ox[0]-0.5*ox[1] # = ox[0]-(ox[1]-ox[0])/2
bins[-1] = 1.5*ox[-1]-0.5*ox[-2] # = ox[-1]+(ox[-1]-ox[-2])/2

# make a temporary node array including input nodes and output bin bounds
# first the boundaries of output bins
tx=bins.copy()
tx = bins.copy()

# if we do not extrapolate,
# because the input is a considered a piece-wise linear function, i.e. the sum of triangles f_i(x),
Expand All @@ -488,32 +481,31 @@ def _unweighted_resample(output_x,input_x,input_flux_density, extrapolate=False)
iy = numpy.append(iy, 0.)

# this sets values left and right of input range to first and/or last input values
# first and last values are=0 if we are not extrapolating
ty=numpy.interp(tx,ix,iy)
# first and last values are = 0 if we are not extrapolating
ty = numpy.interp(tx,ix,iy)

# add input nodes which are inside the node array
k=numpy.where((ix>=tx[0])&(ix<=tx[-1]))[0]
k = numpy.where((ix >= tx[0])&(ix <= tx[-1]))[0]
if k.size :
tx=numpy.append(tx,ix[k])
ty=numpy.append(ty,iy[k])
tx = numpy.append(tx,ix[k])
ty = numpy.append(ty,iy[k])

# sort this node array
p = tx.argsort()
tx=tx[p]
ty=ty[p]
p = tx.argsort()
tx = tx[p]
ty = ty[p]

# now we do a simple integration in each bin of the piece-wise
# linear function of the temporary nodes

# integral of individual trapezes
trapeze_integrals=(ty[1:]+ty[:-1])*(tx[1:]-tx[:-1])/2.
trapeze_integrals = (ty[1:]+ty[:-1])*(tx[1:]-tx[:-1])/2.

# output flux
# for each bin, we sum the trapeze_integrals that belong to that bin
# and divide by the bin size

trapeze_centers=(tx[1:]+tx[:-1])/2.
binsize = bins[1:]-bins[:-1]
trapeze_centers = (tx[1:]+tx[:-1])/2.
binsize = bins[1:]-bins[:-1]

if numpy.any(binsize<=0) :
raise ValueError("Zero or negative bin size")
Expand Down
114 changes: 58 additions & 56 deletions python/lvmdrp/core/rss.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,6 +23,7 @@
from lvmdrp.core.spectrum1d import Spectrum1D, find_continuum, wave_little_interpol
from lvmdrp.core import dataproducts as dp
from lvmdrp.core.fit_profile import polyfit2d, polyval2d
from lvmdrp.core.resample import resample_flux, rebin_spectra

from lvmdrp import __version__ as drpver

Expand Down Expand Up @@ -391,32 +392,32 @@ def from_channels(cls, rss_b, rss_r, rss_z, use_weights=True):
# define interpolators
log.info("interpolating RSS data in new wavelength array")
for rss in rsss:
f = interpolate.interp1d(rss._wave, rss._data, axis=1, bounds_error=False, fill_value=numpy.nan)
fluxes.append(f(new_wave).astype("float32"))
f = interpolate.interp1d(rss._wave, rss._error, axis=1, bounds_error=False, fill_value=numpy.nan)
errors.append(f(new_wave).astype("float32"))
f = interpolate.interp1d(rss._wave, rss._mask, axis=1, kind="nearest", bounds_error=False, fill_value=0)
f = rebin_spectra(new_wave, rss._wave, rss._data, fill=numpy.nan)
fluxes.append(f.astype("float32"))
f = rebin_spectra(new_wave, rss._wave, rss._error, fill=numpy.nan)
errors.append(f.astype("float32"))
f = rebin_spectra(new_wave, rss._wave, rss._mask, axis=1, kind="nearest", bounds_error=False, fill_value=0)
masks.append(f(new_wave).astype("uint8"))
f = interpolate.interp1d(rss._wave, rss._lsf, axis=1, bounds_error=False, fill_value=numpy.nan)
lsfs.append(f(new_wave).astype("float32"))
f = rebin_spectra(new_wave, rss._wave, rss._lsf, fill=numpy.nan)
lsfs.append(f.astype("float32"))
if rss._sky is not None:
f = interpolate.interp1d(rss._wave, rss._sky, axis=1, bounds_error=False, fill_value=numpy.nan)
skies.append(f(new_wave).astype("float32"))
f = rebin_spectra(new_wave, rss._wave, rss._sky, fill=numpy.nan)
skies.append(f.astype("float32"))
if rss._sky_error is not None:
f = interpolate.interp1d(rss._wave, rss._sky_error, axis=1, bounds_error=False, fill_value=numpy.nan)
sky_errors.append(f(new_wave).astype("float32"))
f = rebin_spectra(new_wave, rss._wave, rss._sky_error, fill=numpy.nan)
sky_errors.append(f.astype("float32"))
if rss._sky_east is not None:
f = interpolate.interp1d(rss._wave, rss._sky_east, axis=1, bounds_error=False, fill_value=numpy.nan)
skies_e.append(f(new_wave).astype("float32"))
f = rebin_spectra(new_wave, rss._wave, rss._sky_east, fill=numpy.nan)
skies_e.append(f.astype("float32"))
if rss._sky_east_error is not None:
f = interpolate.interp1d(rss._wave, rss._sky_east_error, axis=1, bounds_error=False, fill_value=numpy.nan)
sky_e_errors.append(f(new_wave).astype("float32"))
f = rebin_spectra(new_wave, rss._wave, rss._sky_east_error, fill=numpy.nan)
sky_e_errors.append(f.astype("float32"))
if rss._sky_west is not None:
f = interpolate.interp1d(rss._wave, rss._sky_west, axis=1, bounds_error=False, fill_value=numpy.nan)
skies_w.append(f(new_wave).astype("float32"))
f = rebin_spectra(new_wave, rss._wave, rss._sky_west, fill=numpy.nan)
skies_w.append(f.astype("float32"))
if rss._sky_west_error is not None:
f = interpolate.interp1d(rss._wave, rss._sky_west_error, axis=1, bounds_error=False, fill_value=numpy.nan)
sky_w_errors.append(f(new_wave).astype("float32"))
f = rebin_spectra(new_wave, rss._wave, rss._sky_west_error, fill=numpy.nan)
sky_w_errors.append(f.astype("float32"))
fluxes = numpy.asarray(fluxes)
errors = numpy.asarray(errors)
masks = numpy.asarray(masks)
Expand Down Expand Up @@ -1861,38 +1862,39 @@ def rectify_wave(self, wave=None, wave_range=None, wave_disp=None, method="linea
else:
raise ValueError(f"Invalid interpolation {method = }")

# TODO: when do we have fluxes vs flux densities? what do we want?
# fit and evaluate interpolators
for ifiber in range(rss._fibers):
f = interpolate.interp1d(rss._wave[ifiber], rss._data[ifiber], kind=method, bounds_error=False, fill_value=numpy.nan)
new_rss._data[ifiber] = f(wave).astype("float32")
f = resample_flux(wave, rss._wave[ifiber], rss._data[ifiber])
new_rss._data[ifiber] = f.astype("float32")
if rss._error is not None:
f = interpolate.interp1d(rss._wave[ifiber], rss._error[ifiber], kind=method, bounds_error=False, fill_value=numpy.nan)
new_rss._error[ifiber] = f(wave).astype("float32")
f = interpolate.interp1d(rss._wave[ifiber], rss._mask[ifiber], kind="nearest", bounds_error=False, fill_value=1)
new_rss._mask[ifiber] = f(wave).astype("bool")
f = resample_flux(wave, rss._wave[ifiber], rss._error[ifiber])
new_rss._error[ifiber] = f.astype("float32")
f = resample_flux(wave, rss._wave[ifiber], rss._mask[ifiber])
new_rss._mask[ifiber] = f.astype("bool")
new_rss._mask[ifiber] |= numpy.isnan(new_rss._data[ifiber])|(new_rss._data[ifiber]==0)
new_rss._mask[ifiber] |= ~numpy.isfinite(new_rss._error[ifiber])
if rss._lsf is not None:
f = interpolate.interp1d(rss._wave[ifiber], rss._lsf[ifiber], kind=method, bounds_error=False, fill_value=numpy.nan)
new_rss._lsf[ifiber] = f(wave).astype("float32")
f = resample_flux(wave, rss._wave[ifiber], rss._lsf[ifiber])
new_rss._lsf[ifiber] = f.astype("float32")
if rss._sky is not None:
f = interpolate.interp1d(rss._wave[ifiber], rss._sky[ifiber], kind=method, bounds_error=False, fill_value=numpy.nan)
new_rss._sky[ifiber] = f(wave).astype("float32")
f = resample_flux(wave, rss._wave[ifiber], rss._sky[ifiber])
new_rss._sky[ifiber] = f.astype("float32")
if rss._sky_error is not None:
f = interpolate.interp1d(rss._wave[ifiber], rss._sky_error[ifiber], kind=method, bounds_error=False, fill_value=numpy.nan)
new_rss._sky_error[ifiber] = f(wave).astype("float32")
f = resample_flux(wave, rss._wave[ifiber], rss._sky_error[ifiber])
new_rss._sky_error[ifiber] = f.astype("float32")
if rss._sky_east is not None:
f = interpolate.interp1d(rss._wave[ifiber], rss._sky_east[ifiber], kind=method, bounds_error=False, fill_value=numpy.nan)
new_rss._sky_east[ifiber] = f(wave).astype("float32")
f = resample_flux(wave, rss._wave[ifiber], rss._sky_east[ifiber])
new_rss._sky_east[ifiber] = f.astype("float32")
if rss._sky_east_error is not None:
f = interpolate.interp1d(rss._wave[ifiber], rss._sky_east_error[ifiber], kind=method, bounds_error=False, fill_value=numpy.nan)
new_rss._sky_east_error[ifiber] = f(wave).astype("float32")
f = resample_flux(wave, rss._wave[ifiber], rss._sky_east_error[ifiber])
new_rss._sky_east_error[ifiber] = f.astype("float32")
if rss._sky_west is not None:
f = interpolate.interp1d(rss._wave[ifiber], rss._sky_west[ifiber], kind=method, bounds_error=False, fill_value=numpy.nan)
new_rss._sky_west[ifiber] = f(wave).astype("float32")
f = resample_flux(wave, rss._wave[ifiber], rss._sky_west[ifiber])
new_rss._sky_west[ifiber] = f.astype("float32")
if rss._sky_west_error is not None:
f = interpolate.interp1d(rss._wave[ifiber], rss._sky_west_error[ifiber], kind=method, bounds_error=False, fill_value=numpy.nan)
new_rss._sky_west_error[ifiber] = f(wave).astype("float32")
f = resample_flux(wave, rss._wave[ifiber], rss._sky_west_error[ifiber])
new_rss._sky_west_error[ifiber] = f.astype("float32")
# add supersky information if available
if rss._supersky is not None:
new_rss.set_supersky(rss._supersky)
Expand Down Expand Up @@ -2000,30 +2002,30 @@ def to_native_wave(self, method="linear", interp_density=True, return_density=Fa

# interpolate data, error, mask and sky arrays from rectified grid to original grid
for ifiber in range(rss._fibers):
f = interpolate.interp1d(rss._wave, rss._data[ifiber], kind=method, bounds_error=False, fill_value=numpy.nan)
new_rss._data[ifiber] = f(wave[ifiber]).astype("float32")
f = interpolate.interp1d(rss._wave, rss._error[ifiber], kind=method, bounds_error=False, fill_value=numpy.nan)
new_rss._error[ifiber] = f(wave[ifiber]).astype("float32")
f = interpolate.interp1d(rss._wave, rss._mask[ifiber], kind="nearest", bounds_error=False, fill_value=1)
f = resample_flux(wave[ifiber], rss._wave, rss._data[ifiber])
new_rss._data[ifiber] = f.astype("float32")
f = resample_flux(wave[ifiber][ifiber], rss._wave, rss._error[ifiber])
new_rss._error[ifiber] = f.astype("float32")
f = resample_flux(wave[ifiber][ifiber], rss._wave, rss._mask[ifiber], kind="nearest", bounds_error=False, fill_value=1)
new_rss._mask[ifiber] = f(wave[ifiber]).astype("bool")
if rss._sky is not None:
f = interpolate.interp1d(rss._wave, rss._sky[ifiber], kind=method, bounds_error=False, fill_value=numpy.nan)
new_rss._sky[ifiber] = f(wave[ifiber]).astype("float32")
f = resample_flux(wave[ifiber], rss._wave, rss._sky[ifiber])
new_rss._sky[ifiber] = f.astype("float32")
if rss._sky_error is not None:
f = interpolate.interp1d(rss._wave, rss._sky_error[ifiber], kind=method, bounds_error=False, fill_value=numpy.nan)
new_rss._sky_error[ifiber] = f(wave[ifiber]).astype("float32")
f = resample_flux(wave[ifiber], rss._wave, rss._sky_error[ifiber])
new_rss._sky_error[ifiber] = f.astype("float32")
if rss._sky_east is not None:
f = interpolate.interp1d(rss._wave, rss._sky_east[ifiber], kind=method, bounds_error=False, fill_value=numpy.nan)
new_rss._sky_east[ifiber] = f(wave[ifiber]).astype("float32")
f = resample_flux(wave[ifiber], rss._wave, rss._sky_east[ifiber])
new_rss._sky_east[ifiber] = f.astype("float32")
if rss._sky_east_error is not None:
f = interpolate.interp1d(rss._wave, rss._sky_east_error[ifiber], kind=method, bounds_error=False, fill_value=numpy.nan)
new_rss._sky_east_error[ifiber] = f(wave[ifiber]).astype("float32")
f = resample_flux(wave[ifiber], rss._wave, rss._sky_east_error[ifiber])
new_rss._sky_east_error[ifiber] = f.astype("float32")
if rss._sky_west is not None:
f = interpolate.interp1d(rss._wave, rss._sky_west[ifiber], kind=method, bounds_error=False, fill_value=numpy.nan)
new_rss._sky_west[ifiber] = f(wave[ifiber]).astype("float32")
f = resample_flux(wave[ifiber], rss._wave, rss._sky_west[ifiber])
new_rss._sky_west[ifiber] = f.astype("float32")
if rss._sky_west_error is not None:
f = interpolate.interp1d(rss._wave, rss._sky_west_error[ifiber], kind=method, bounds_error=False, fill_value=numpy.nan)
new_rss._sky_west_error[ifiber] = f(wave[ifiber]).astype("float32")
f = resample_flux(wave[ifiber], rss._wave, rss._sky_west_error[ifiber])
new_rss._sky_west_error[ifiber] = f.astype("float32")

if not return_density and unit.endswith("/angstrom"):
dlambda = numpy.gradient(wave, axis=1)
Expand Down
2 changes: 1 addition & 1 deletion python/lvmdrp/functions/rssMethod.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1120,7 +1120,7 @@ def resample_wavelength(in_rss: str, out_rss: str, method: str = "linear",
log.info(f"using wavelength range {wave_range = } angstrom and {wave_disp = } angstrom pixel size")

# resample the wavelength solution
log.info(f"resampling the wavelength solution using {method = } interpolation")
log.info(f"resampling the spectra ...")
new_rss = rss.rectify_wave(wave_range=wave_range, wave_disp=wave_disp, method=method, return_density=convert_to_density)

# write output RSS
Expand Down

0 comments on commit 634c0b5

Please sign in to comment.