Skip to content

Commit

Permalink
distinguish between resampling fluxes and flux densities
Browse files Browse the repository at this point in the history
  • Loading branch information
@ndrory
ndrory committed Oct 9, 2024
1 parent ead8fdb commit da0f284
Show file tree
Hide file tree
Showing 3 changed files with 55 additions and 96 deletions.
60 changes: 6 additions & 54 deletions examples/lco_com/testrebin.ipynb
View file

Large diffs are not rendered by default.

50 changes: 28 additions & 22 deletions python/lvmdrp/core/resample.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -335,7 +335,7 @@ def template_logwl_resample(spectrum, template, wblue=None, wred=None,
wave_array = numpy.power(10., log_wave_array) * spectrum.spectral_axis.unit


def resample_flux(xout, x, flux, ivar=None, extrapolate=False):
def resample_flux_density(xout, x, flux, ivar=None, extrapolate=False):
"""Returns a flux conserving resampling of an input flux density.
The total integrated flux is conserved.
Expand All @@ -360,26 +360,6 @@ def resample_flux(xout, x, flux, ivar=None, extrapolate=False):
This interpolation conserves flux such that, on average,
output_flux_density = input_flux_density
The input flux density outside of the range defined by the edges of the first
and last bins is considered null. The bin size of bin 'i' is given by (x[i+1]-x[i-1])/2
except for the first and last bin where it is (x[1]-x[0]) and (x[-1]-x[-2])
so flux density is zero for x<x[0]-(x[1]-x[0])/2 and x>x[-1]-(x[-1]-x[-2])/2
The input is interpreted as the nodes positions and node values of
a piece-wise linear function::
y(x) = sum_i y_i * f_i(x)
with::
f_i(x) = (x_{i-1}<x<=x_{i})*(x-x_{i-1})/(x_{i}-x_{i-1})
+ (x_{i}<x<=x_{i+1})*(x-x_{i+1})/(x_{i}-x_{i+1})
the output value is the average flux density in a bin::
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})
"""
if ivar is None:
return _unweighted_resample(xout, x, flux, extrapolate=extrapolate)
Expand All @@ -397,7 +377,33 @@ def resample_flux(xout, x, flux, ivar=None, extrapolate=False):

return outflux, outivar

def _unweighted_resample(output_x,input_x,input_flux_density, extrapolate=False) :
def resample_flux(xout, x, flux, extrapolate=False):
"""Returns a flux conserving resampling of an input flux.
The total integrated flux is conserved.
Args:
- xout: output SORTED vector, not necessarily linearly spaced
- x: input SORTED vector, not necessarily linearly spaced
- flux: input flux sampled at x
both x and xout must represent the same quantity with the same unit
Options:
- extrapolate: extrapolate using edge values of input array, default is False,
in which case values outside of input array are set to zero.
Returns:
returns outflux
This interpolation conserves flux such that, on average,
output_flux_density = input_flux_density
"""
flux_dens = flux/numpy.gradient(x)
f = _unweighted_resample(xout, x, flux_dens, extrapolate=extrapolate)
return f*numpy.gradient(xout)


def _unweighted_resample(output_x, input_x, input_flux_density, extrapolate=False) :
"""Returns a flux conserving resampling of an input flux density.
The total integrated flux is conserved.
Expand Down
41 changes: 21 additions & 20 deletions python/lvmdrp/core/rss.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +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.core.resample import resample_flux_density, rebin_spectra

from lvmdrp import __version__ as drpver

Expand Down Expand Up @@ -390,6 +390,7 @@ def from_channels(cls, rss_b, rss_r, rss_z, use_weights=True):
log.info(f"new wavelength sampling: min = {sampling.min():.2f}, max = {sampling.max():.2f}")

# define interpolators
# TODO: why are we rebinning again?
log.info("interpolating RSS data in new wavelength array")
for rss in rsss:
f = rebin_spectra(new_wave, rss._wave, rss._data, fill=numpy.nan)
Expand Down Expand Up @@ -1851,35 +1852,35 @@ def rectify_wave(self, wave=None, wave_range=None, wave_disp=None):

# fit and evaluate interpolators
for ifiber in range(rss._fibers):
f = resample_flux(wave, rss._wave[ifiber], rss._data[ifiber])
f = resample_flux_density(wave, rss._wave[ifiber], rss._data[ifiber])
new_rss._data[ifiber] = f.astype("float32")
if rss._error is not None:
f = resample_flux(wave, rss._wave[ifiber], rss._error[ifiber])
f = resample_flux_density(wave, rss._wave[ifiber], rss._error[ifiber])
new_rss._error[ifiber] = f.astype("float32")
f = resample_flux(wave, rss._wave[ifiber], rss._mask[ifiber])
f = resample_flux_density(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 = resample_flux(wave, rss._wave[ifiber], rss._lsf[ifiber])
f = resample_flux_density(wave, rss._wave[ifiber], rss._lsf[ifiber])
new_rss._lsf[ifiber] = f.astype("float32")
if rss._sky is not None:
f = resample_flux(wave, rss._wave[ifiber], rss._sky[ifiber])
f = resample_flux_density(wave, rss._wave[ifiber], rss._sky[ifiber])
new_rss._sky[ifiber] = f.astype("float32")
if rss._sky_error is not None:
f = resample_flux(wave, rss._wave[ifiber], rss._sky_error[ifiber])
f = resample_flux_density(wave, rss._wave[ifiber], rss._sky_error[ifiber])
new_rss._sky_error[ifiber] = f.astype("float32")
if rss._sky_east is not None:
f = resample_flux(wave, rss._wave[ifiber], rss._sky_east[ifiber])
f = resample_flux_density(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 = resample_flux(wave, rss._wave[ifiber], rss._sky_east_error[ifiber])
f = resample_flux_density(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 = resample_flux(wave, rss._wave[ifiber], rss._sky_west[ifiber])
f = resample_flux_density(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 = resample_flux(wave, rss._wave[ifiber], rss._sky_west_error[ifiber])
f = resample_flux_density(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:
Expand Down Expand Up @@ -1970,29 +1971,29 @@ 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 = resample_flux(wave[ifiber], rss._wave, rss._data[ifiber])
f = resample_flux_density(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])
f = resample_flux_density(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)
f = resample_flux_density(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 = resample_flux(wave[ifiber], rss._wave, rss._sky[ifiber])
f = resample_flux_density(wave[ifiber], rss._wave, rss._sky[ifiber])
new_rss._sky[ifiber] = f.astype("float32")
if rss._sky_error is not None:
f = resample_flux(wave[ifiber], rss._wave, rss._sky_error[ifiber])
f = resample_flux_density(wave[ifiber], rss._wave, rss._sky_error[ifiber])
new_rss._sky_error[ifiber] = f.astype("float32")
if rss._sky_east is not None:
f = resample_flux(wave[ifiber], rss._wave, rss._sky_east[ifiber])
f = resample_flux_density(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 = resample_flux(wave[ifiber], rss._wave, rss._sky_east_error[ifiber])
f = resample_flux_density(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 = resample_flux(wave[ifiber], rss._wave, rss._sky_west[ifiber])
f = resample_flux_density(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 = resample_flux(wave[ifiber], rss._wave, rss._sky_west_error[ifiber])
f = resample_flux_density(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"):
Expand Down

0 comments on commit da0f284

Please sign in to comment.